home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch7_7 / gg4d / veclib4d.c < prev   
Encoding:
C/C++ Source or Header  |  1995-04-04  |  19.0 KB  |  1,018 lines

  1. #include <stdio.h>
  2. #include <stdlib.h>    /* for malloc() definition */
  3. #include <math.h>
  4. #include "GGems.h"
  5. #include "GGems4d.h"
  6.  
  7. /* 4D VECTOR LIBRARY.
  8.  *
  9.  * Steve Hill, Computing Laboratory, University of Kent, UK
  10.  *
  11.  * Email: S.A.Hill@ukc.ac.uk
  12.  *
  13.  * Included with Graphics Gems V, Academic Press ed. Alan Paeth.
  14.  *
  15.  * Extends the standard Graphics Gems library to four dimensions
  16.  * and provides functions to create and manipulate homogeneous
  17.  * transformations for 4-vectors (including projections).
  18.  */
  19.  
  20. /* V4SquaredLength()
  21.  *
  22.  * Compute the square of the magnitude of a 4-vector
  23.  */
  24.  
  25. double
  26. V4SquaredLength(v)
  27. Vector4 *v;
  28. {
  29.     return v->x * v->x +
  30.            v->y * v->y +
  31.            v->z * v->z +
  32.            v->w * v->w;
  33. }
  34.  
  35. /* V4Length()
  36.  *
  37.  * Compute the magnitude of a 4-vector
  38.  * If square of magnitude is required use V4SquaredLength().
  39.  */
  40.  
  41. double
  42. V4Length(v)
  43. Vector4    *v;
  44. {
  45.     return sqrt(V4SquaredLength(v));
  46. }
  47.  
  48. /* V4Negate()
  49.  *
  50.  * Negate vector elementwise.
  51.  */
  52.  
  53. Vector4 *
  54. V4Negate(v)
  55. Vector4    *v;
  56. {
  57.     v->x = -v->x;
  58.     v->y = -v->y;
  59.     v->z = -v->z;
  60.     v->w = -v->w;
  61.  
  62.     return v;
  63. }
  64.  
  65. /* V4Normalize()
  66.  *
  67.  * Convert vector to unit magnitude.
  68.  * A zero-length vector is unchanged.
  69.  */
  70.  
  71. Vector4 *
  72. V4Normalize(v)
  73. Vector4    *v;
  74. {
  75.     double    len = V4Length(v);
  76.  
  77.     if (len != 0.0)
  78.     {
  79.         v->x /= len;
  80.         v->y /= len;
  81.         v->z /= len;
  82.         v->w /= len;
  83.     }
  84.  
  85.     return v;
  86. }
  87.  
  88. /* V4Scale()
  89.  *
  90.  * Scale vector to requested magnitude.
  91.  * Zero length vectors unchanged.
  92.  */
  93.  
  94. Vector4 *
  95. V4Scale(v, newlen)
  96. Vector4    *v;
  97. double    newlen;
  98. {
  99.     double    len = V4Length(v);
  100.  
  101.     if (len != 0.0)
  102.     {
  103.         v->x *= newlen/len;
  104.         v->y *= newlen/len;
  105.         v->z *= newlen/len;
  106.         v->w *= newlen/len;
  107.     }
  108.  
  109.     return v;
  110. }
  111.  
  112. /* V4Add()
  113.  *
  114.  * Add (elementwise) two 4-vectors
  115.  */
  116.  
  117. Vector4 *
  118. V4Add(a, b, result)
  119. Vector4    *a, *b, *result;
  120. {
  121.     result->x = a->x + b-> x;
  122.     result->y = a->y + b-> y;
  123.     result->z = a->z + b-> z;
  124.     result->w = a->w + b-> w;
  125.  
  126.     return result;
  127. }
  128.  
  129. /* V4Sub()
  130.  *
  131.  * Subtract (elementwise) vector b from vector a.
  132.  */
  133.  
  134. Vector4 *
  135. V4Sub(a, b, result)
  136. Vector4    *a, *b, *result;
  137. {
  138.     result->x = a->x - b-> x;
  139.     result->y = a->y - b-> y;
  140.     result->z = a->z - b-> z;
  141.     result->w = a->w - b-> w;
  142.  
  143.     return result;
  144. }
  145.  
  146. /* V4Dot()
  147.  *
  148.  * Returns the 4-space dot-product of two vectors.
  149.  */
  150.  
  151. double
  152. V4Dot(a, b)
  153. Vector4    *a, *b;
  154. {
  155.     return    a->x * b->x +
  156.         a->y * b->y +
  157.         a->z * b->z +
  158.         a->w * b->w;
  159. }
  160.  
  161. /* V4Lerp()
  162.  *
  163.  * Calculates vector that is a linear interpolation between
  164.  * lo and hi, accoreding to alpha.
  165.  * See also GraphicsGems.h (GGems.h)
  166.  */
  167.  
  168. Vector4 *
  169. V4Lerp(lo, hi, alpha, result)
  170. Vector4    *lo, *hi, *result;
  171. double    alpha;
  172. {
  173.     result->x = LERP(alpha, lo->x, hi->x);
  174.     result->y = LERP(alpha, lo->y, hi->y);
  175.     result->z = LERP(alpha, lo->z, hi->z);
  176.     result->w = LERP(alpha, lo->w, hi->w);
  177.  
  178.     return result;
  179. }
  180.  
  181. /* V4Combine()
  182.  *
  183.  * Calculate a linear combination of two vectors.
  184.  */
  185.  
  186. Vector4 *
  187. V4Combine(a, b, result, ascl, bscl)
  188. Vector4    *a, *b, *result;
  189. double    ascl, bscl;
  190. {
  191.     result->x = ascl * a->x + bscl * b->x;
  192.     result->y = ascl * a->y + bscl * b->y;
  193.     result->z = ascl * a->z + bscl * b->z;
  194.     result->w = ascl * a->w + bscl * b->w;
  195.  
  196.     return result;
  197. }
  198.  
  199. /* V4Mul()
  200.  *
  201.  * Multiply (elementwise) two vectors.
  202.  */
  203.  
  204. Vector4 *
  205. V4Mul(a, b, result)
  206. Vector4    *a, *b, *result;
  207. {
  208.     result->x = a->x * b->x;
  209.     result->y = a->y * b->y;
  210.     result->z = a->z * b->z;
  211.     result->w = a->w * b->w;
  212.  
  213.     return result;
  214. }
  215.  
  216. /* V4DistanceBetween2Points()
  217.  *
  218.  * Calculates the distance between two 4-space points.
  219.  */
  220.  
  221. double
  222. V4DistanceBetween2Points(a, b)
  223. Point4    *a, *b;
  224. {
  225.     double    dx = a->x - b->x,
  226.         dy = a->y - b->y,
  227.         dz = a->z - b->z,
  228.         dw = a->w - b->w;
  229.  
  230.     return sqrt(dx*dx + dy*dy + dz*dz + dw*dw);
  231. }
  232.  
  233. /* V4Cross()
  234.  *
  235.  * Calculates the cross product of three 4-space vectors.
  236.  */
  237.  
  238. Vector4 *
  239. V4Cross(a, b, c, result)
  240. Vector4    *a, *b, *c, *result;
  241. {
  242.     double    d1, d2, d3, d4, d5, d6;
  243.  
  244.     d1 = (b->z * c->w) - (b->w * c->z);
  245.     d2 = (b->y * c->w) - (b->w * c->y);
  246.     d3 = (b->y * c->z) - (b->z * c->y);
  247.     d4 = (b->x * c->w) - (b->w * c->x);
  248.     d5 = (b->x * c->z) - (b->z * c->x);
  249.     d6 = (b->x * c->y) - (b->y * c->x);
  250.  
  251.     result->x = - a->y * d1 + a->z * d2 - a->w * d3;
  252.     result->y =   a->x * d1 - a->z * d4 + a->w * d5;
  253.     result->z = - a->x * d2 + a->y * d4 - a->w * d6;
  254.     result->w =   a->x * d3 - a->y * d5 + a->z * d6;
  255.  
  256.     return result;
  257. }
  258.  
  259. /* V4New()
  260.  *
  261.  * Allocate and initialise a new 4-vector
  262.  */
  263.  
  264. Vector4 *
  265. V4New(x, y, z, w)
  266. double    x, y, z, w;
  267. {
  268.     Vector4    *v = NEWTYPE(Vector4);
  269.  
  270.     v->x = x;
  271.     v->y = y;
  272.     v->z = z;
  273.     v->w = w;
  274.  
  275.     return v;
  276. }
  277.  
  278. /* V4Duplicate()
  279.  *
  280.  * Create a copy of a 4-vector.
  281.  */
  282.  
  283. Vector4 *
  284. V4Duplicate(a)
  285. Vector4    *a;
  286. {
  287.     Vector4    *v = NEWTYPE(Vector4);
  288.  
  289.     v->x = a->x;
  290.     v->y = a->y;
  291.     v->z = a->z;
  292.     v->w = a->w;
  293.  
  294.     return v;
  295. }
  296.  
  297. /* V4MatPrint()
  298.  *
  299.  * Diagnostic function to print out a 5x5 matrix
  300.  *
  301.  */
  302.  
  303. void
  304. V4MatPrint(file, mat)
  305. FILE    *file;
  306. Matrix5    *mat;
  307. {
  308.     int    i, j;
  309.  
  310.     for (i = 0; i < 5; i += 1)
  311.     {
  312.         for (j = 0; j < 5; j += 1)
  313.             fprintf(file, "%lf ", mat->element[i][j]);
  314.         putc('\n', file);
  315.     }
  316. }
  317.  
  318. /* V4MatCopy()
  319.  *
  320.  * Copy a 5x5 matrix.
  321.  */
  322.  
  323. Matrix5 *
  324. V4MatCopy(from, to)
  325. Matrix5    *from, *to;
  326. {
  327.     int    i, j;
  328.  
  329.     for (i = 0; i < 5; i += 1)
  330.     for (j = 0; j < 5; j += 1)
  331.         to->element[i][j] = from->element[i][j];
  332.  
  333.     return to;
  334. }
  335.  
  336. /* V4MulPointByMatrix()
  337.  *
  338.  * Apply 4x4 transformation matrix to 4-space point.
  339.  *
  340.  * In common with the standard Graphics Gems library, points/vectors
  341.  * are considered to be row vectors, hence multiplication is
  342.  * post-multiplication.  The transformation T o S o R represents
  343.  * a translation followed by a scale followed by a rotation.
  344.  * Where o is matrix multiplication (see V4MatMul() and V4MatMul2() below).
  345.  */
  346.  
  347. Point4 *
  348. V4MulPointByMatrix(pin, m, pout)
  349. Point4    *pin, *pout;
  350. Matrix4    *m;
  351. {
  352.     pout->x = pin->x * m->element[0][0] +
  353.           pin->y * m->element[1][0] +
  354.           pin->z * m->element[2][0] +
  355.           pin->w * m->element[3][0];
  356.     pout->y = pin->x * m->element[0][1] +
  357.           pin->y * m->element[1][1] +
  358.           pin->z * m->element[2][1] +
  359.           pin->w * m->element[3][1];
  360.     pout->z = pin->x * m->element[0][2] +
  361.           pin->y * m->element[1][2] +
  362.           pin->z * m->element[2][2] +
  363.           pin->w * m->element[3][2];
  364.     pout->w = pin->x * m->element[0][3] +
  365.           pin->y * m->element[1][3] +
  366.           pin->z * m->element[2][3] +
  367.           pin->w * m->element[3][3];
  368.  
  369.     return pout;
  370. }
  371.  
  372. /* V4MulPointByProjMatrix()
  373.  *
  374.  * Apply 5x5 (projection) matrix to 4-space point
  375.  *
  376.  * For example, the projection matrix:
  377.  *
  378.  *   [  1  0  0  0    0    ]
  379.  *   [  0  1  0  0    0    ]
  380.  *   [  0  0  0  0   1/dz  ]
  381.  *   [  0  0  0  0   1/dw  ]
  382.  *   [  0  0  0  0    1    ]
  383.  *
  384.  * represents a combined perspective/parallel projection along
  385.  * the Z and W axes.  As dz or dw tend to infinity the projections
  386.  * along those axes becomes parallel.
  387.  *
  388.  * 5x5 matrices may also be used to perform translations.
  389.  *
  390.  *   [  1  0  0  0  0 ]
  391.  *   [  0  1  0  0  0 ]
  392.  *   [  0  0  1  0  0 ]
  393.  *   [  0  0  0  1  0 ]
  394.  *   [  x  y  z  w  1 ]
  395.  */
  396.  
  397. Point4 *
  398. V4MulPointByProjMatrix(pin, m, pout)
  399. Point4    *pin, *pout;
  400. Matrix5    *m;
  401. {
  402.     double    v;
  403.  
  404.     pout->x = pin->x * m->element[0][0] +
  405.           pin->y * m->element[1][0] +
  406.           pin->z * m->element[2][0] +
  407.           pin->w * m->element[3][0] +
  408.           m->element[4][0];
  409.     pout->y = pin->x * m->element[0][1] +
  410.           pin->y * m->element[1][1] +
  411.           pin->z * m->element[2][1] +
  412.           pin->w * m->element[3][1] +
  413.           m->element[4][1];
  414.     pout->z = pin->x * m->element[0][2] +
  415.           pin->y * m->element[1][2] +
  416.           pin->z * m->element[2][2] +
  417.           pin->w * m->element[3][2] +
  418.           m->element[4][2];
  419.     pout->w = pin->x * m->element[0][3] +
  420.           pin->y * m->element[1][3] +
  421.           pin->z * m->element[2][3] +
  422.           pin->w * m->element[3][3] +
  423.           m->element[4][3];
  424.  
  425.     v = pin->x * m->element[0][4] +
  426.         pin->y * m->element[1][4] +
  427.         pin->z * m->element[2][4] +
  428.         pin->w * m->element[3][4] +
  429.         m->element[4][4];
  430.  
  431.     if (v != 0.0)
  432.     {
  433.         pout->x /= v;
  434.         pout->y /= v;
  435.         pout->z /= v;
  436.         pout->w /= v;
  437.     }
  438.  
  439.     return pout;
  440. }
  441.  
  442. /* V4MatMul()
  443.  *
  444.  * Multiply two 5x5 matrices.  Result matrix must be distinct
  445.  * from argument matrices.
  446.  */
  447.  
  448. Matrix5 *
  449. V4MatMul(a, b, result)
  450. Matrix5    *a, *b, *result;
  451. {
  452.     int    i, j, k;
  453.  
  454.     for (i = 0; i < 5; i++)
  455.     for (j = 0; j < 5; j++)
  456.     {
  457.         double t = 0.0;
  458.  
  459.         for (k = 0; k < 5; k++)
  460.             t += a->element[i][k] * b->element[k][j];
  461.  
  462.         result->element[i][j] = t;
  463.     }
  464.  
  465.     return result;
  466. }
  467.  
  468. /* V4MatMul2
  469.  *
  470.  * Variation of V4MatMul() where result matrix may be the same
  471.  * as matrices a or b.
  472.  */
  473.  
  474. Matrix5 *
  475. V4MatMul2(a, b, result)
  476. Matrix5 *a, *b, *result;
  477. {
  478.     Matrix5    tmp;
  479.  
  480.     V4MatMul(a, b, &tmp);
  481.     V4MatCopy(&tmp, result);
  482.     return result;
  483. }
  484.  
  485. /* 4D TRANSFORMATIONS
  486.  *
  487.  * The following functions provide various methods
  488.  * to apply and manipulate transformations.
  489.  */
  490.  
  491. /* V4MatZero()
  492.  *
  493.  * Clear a 5x5 matrix.
  494.  */
  495.  
  496. Matrix5 *
  497. V4MatZero(mat)
  498. Matrix5    *mat;
  499. {
  500.     int    i, j;
  501.  
  502.     for (i = 0; i < 5; i += 1)
  503.     for (j = 0; j < 5; j += 1)
  504.         mat->element[i][j] = 0.0;
  505. }
  506.  
  507. /* STANDARD TRANSFORMATIONS
  508.  *
  509.  * The following functions initialise the standard affine transformations
  510.  * ie. Identity, Rotation (6 kinds), Translation and Scaling.
  511.  */
  512.  
  513. /* V4MatIdentity()
  514.  *
  515.  * Initialises to identity matrix.
  516.  */
  517.  
  518. Matrix5 *
  519. V4MatIdentity(mat)
  520. Matrix5    *mat;
  521. {
  522.     int    i, j;
  523.  
  524.     for (i = 0; i < 5; i += 1)
  525.     for (j = 0; j < 5; j += 1)
  526.         mat->element[i][j] = i == j ? 1.0 : 0.0;
  527. }
  528.  
  529. /* V4MatRotationXY() etc.
  530.  *
  531.  * Initialises a matrix to a rotation by theta about the indicated hyperplane.
  532.  */
  533.  
  534. Matrix5 *
  535. V4MatRotationXY(mat, theta)
  536. Matrix5    *mat;
  537. double    theta;
  538. {
  539.     double    s, c;
  540.  
  541.     s = sin(theta);
  542.     c = cos(theta);
  543.  
  544.     V4MatIdentity(mat);
  545.     mat->element[0][0] = c;
  546.     mat->element[0][1] = s;
  547.     mat->element[1][0] = -s;
  548.     mat->element[1][1] = c;
  549.  
  550.     return mat;
  551. }
  552.  
  553. Matrix5 *
  554. V4MatRotationXW(mat, theta)
  555. Matrix5    *mat;
  556. double    theta;
  557. {
  558.     double    s, c;
  559.  
  560.     s = sin(theta);
  561.     c = cos(theta);
  562.  
  563.     V4MatIdentity(mat);
  564.     mat->element[0][0] = c;
  565.     mat->element[0][3] = s;
  566.     mat->element[3][0] = -s;
  567.     mat->element[3][3] = c;
  568.  
  569.     return mat;
  570. }
  571.  
  572. Matrix5 *
  573. V4MatRotationYZ(mat, theta)
  574. Matrix5    *mat;
  575. double    theta;
  576. {
  577.     double    s, c;
  578.  
  579.     s = sin(theta);
  580.     c = cos(theta);
  581.  
  582.     V4MatIdentity(mat);
  583.     mat->element[1][1] = c;
  584.     mat->element[1][2] = s;
  585.     mat->element[2][1] = -s;
  586.     mat->element[2][2] = c;
  587.  
  588.     return mat;
  589. }
  590.  
  591. Matrix5 *
  592. V4MatRotationYW(mat, theta)
  593. Matrix5    *mat;
  594. double    theta;
  595. {
  596.     double    s, c;
  597.  
  598.     s = sin(theta);
  599.     c = cos(theta);
  600.  
  601.     V4MatIdentity(mat);
  602.     mat->element[1][1] = c;
  603.     mat->element[1][3] = -s;
  604.     mat->element[3][1] = s;
  605.     mat->element[3][3] = c;
  606.  
  607.     return mat;
  608. }
  609.  
  610. Matrix5 *
  611. V4MatRotationZX(mat, theta)
  612. Matrix5    *mat;
  613. double    theta;
  614. {
  615.     double    s, c;
  616.  
  617.     s = sin(theta);
  618.     c = cos(theta);
  619.  
  620.     V4MatIdentity(mat);
  621.     mat->element[0][0] = c;
  622.     mat->element[0][2] = -s;
  623.     mat->element[2][0] = s;
  624.     mat->element[2][2] = c;
  625.  
  626.     return mat;
  627. }
  628.  
  629. Matrix5 *
  630. V4MatRotationZW(mat, theta)
  631. Matrix5    *mat;
  632. double    theta;
  633. {
  634.     double    s, c;
  635.  
  636.     s = sin(theta);
  637.     c = cos(theta);
  638.  
  639.     V4MatIdentity(mat);
  640.     mat->element[2][2] = c;
  641.     mat->element[2][3] = -s;
  642.     mat->element[3][2] = s;
  643.     mat->element[3][3] = c;
  644.  
  645.     return mat;
  646. }
  647.  
  648. /* V4MatTranslation()
  649.  *
  650.  * Initialise a matrix to the translation v
  651.  */
  652.  
  653. Matrix5 *
  654. V4MatTranslation(mat, v)
  655. Matrix5    *mat;
  656. Vector4    *v;
  657. {
  658.     V4MatIdentity(mat);
  659.     mat->element[4][0] = v->x;
  660.     mat->element[4][1] = v->y;
  661.     mat->element[4][2] = v->z;
  662.     mat->element[4][3] = v->w;
  663.  
  664.     return mat;
  665. }
  666.  
  667. /* V4MatScaling()
  668.  *
  669.  * Initialise a matrix to the scaling v
  670.  */
  671.  
  672. Matrix5 *
  673. V4MatScaling(mat, v)
  674. Matrix5    *mat;
  675. Vector4    *v;
  676. {
  677.     V4MatIdentity(mat);
  678.     mat->element[0][0] = v->x;
  679.     mat->element[1][1] = v->y;
  680.     mat->element[2][2] = v->z;
  681.     mat->element[3][3] = v->w;
  682.  
  683.     return mat;
  684. }
  685.  
  686. /* MATRIX OPERATIONS
  687.  *
  688.  * The following functions are provided for a fast compact method
  689.  * to create compound transformations.  Matrix multiplication has
  690.  * been folded into the functions thus eliminating the need for
  691.  * many redundant calculations.  For sparse matrices containing
  692.  * mainly just ones and zeros this represents a considerable saving.
  693.  *
  694.  * Each function takes the form F(M, args) and has the effect
  695.  * M := M o F1(args) where F1 is the transformation corresponding
  696.  * to F.
  697.  */
  698.  
  699. Matrix5    *
  700. V4MatRotateXY(mat, theta)
  701. Matrix5    *mat;
  702. double    theta;
  703. {
  704.     double    a, b, f, g, k, l, p, q, u, v;
  705.     double    st, ct;
  706.  
  707.     st = sin(theta);
  708.     ct = cos(theta);
  709.  
  710.     a = mat->element[0][0]; b = mat->element[0][1];
  711.     f = mat->element[1][0]; g = mat->element[1][1];
  712.     k = mat->element[2][0]; l = mat->element[2][1];
  713.     p = mat->element[3][0]; q = mat->element[3][1];
  714.     u = mat->element[4][0]; v = mat->element[4][1];
  715.  
  716.     mat->element[0][0] = a * ct - b * st;
  717.     mat->element[0][1] = a * st + b * ct;
  718.     mat->element[1][0] = f * ct - g * st;
  719.     mat->element[1][1] = f * st + g * ct;
  720.     mat->element[2][0] = k * ct - l * st;
  721.     mat->element[2][1] = k * st + l * ct;
  722.     mat->element[3][0] = p * ct - q * st;
  723.     mat->element[3][1] = p * st + q * ct;
  724.     mat->element[4][0] = u * ct - v * st;
  725.     mat->element[4][1] = u * st + v * ct;
  726.  
  727.     return mat;
  728. }
  729.  
  730. Matrix5    *
  731. V4MatRotateXW(mat, theta)
  732. Matrix5    *mat;
  733. double    theta;
  734. {
  735.     double    a, d, f, i, k, n, p, s, u, x;
  736.     double    st, ct;
  737.  
  738.     st = sin(theta);
  739.     ct = cos(theta);
  740.  
  741.     a = mat->element[0][0]; d = mat->element[0][3];
  742.     f = mat->element[1][0]; i = mat->element[1][3];
  743.     k = mat->element[2][0]; n = mat->element[2][3];
  744.     p = mat->element[3][0]; s = mat->element[3][3];
  745.     u = mat->element[4][0]; x = mat->element[4][3];
  746.  
  747.     mat->element[0][0] = a * ct - d * st;
  748.     mat->element[0][3] = a * st + d * ct;
  749.     mat->element[1][0] = f * ct - i * st;
  750.     mat->element[1][3] = f * st + i * ct;
  751.     mat->element[2][0] = k * ct - n * st;
  752.     mat->element[2][3] = k * st + n * ct;
  753.     mat->element[3][0] = p * ct - s * st;
  754.     mat->element[3][3] = p * st + s * ct;
  755.     mat->element[4][0] = u * ct - x * st;
  756.     mat->element[4][3] = u * st + x * ct;
  757.  
  758.     return mat;
  759. }
  760.  
  761. Matrix5    *
  762. V4MatRotateYZ(mat, theta)
  763. Matrix5    *mat;
  764. double    theta;
  765. {
  766.     double    b, c, g, h, l, m, q, r, v, w;
  767.     double st, ct;
  768.  
  769.     st = sin(theta);
  770.     ct = cos(theta);
  771.  
  772.     b = mat->element[0][1]; c = mat->element[0][2];
  773.     g = mat->element[1][1]; h = mat->element[1][2];
  774.     l = mat->element[2][1]; m = mat->element[2][2];
  775.     q = mat->element[3][1]; r = mat->element[3][2];
  776.     v = mat->element[4][1]; w = mat->element[4][2];
  777.  
  778.     mat->element[0][1] = b * ct - c * st;
  779.     mat->element[0][2] = b * st + c * ct;
  780.     mat->element[1][1] = g * ct - h * st;
  781.     mat->element[1][2] = g * st + h * ct;
  782.     mat->element[2][1] = l * ct - m * st;
  783.     mat->element[2][2] = l * st + m * ct;
  784.     mat->element[3][1] = q * ct - r * st;
  785.     mat->element[3][2] = q * st + r * ct;
  786.     mat->element[4][1] = v * ct - w * st;
  787.     mat->element[4][2] = v * st + w * ct;
  788.  
  789.     return mat;
  790. }
  791.  
  792. Matrix5    *
  793. V4MatRotateYW(mat, theta)
  794. Matrix5    *mat;
  795. double    theta;
  796. {
  797.     double    b, d, g, i, l, n, q, s, v, x;
  798.     double    st, ct;
  799.  
  800.     st = sin(theta);
  801.     ct = cos(theta);
  802.  
  803.     b = mat->element[0][1]; d = mat->element[0][3];
  804.     g = mat->element[1][1]; i = mat->element[1][3];
  805.     l = mat->element[2][1]; n = mat->element[2][3];
  806.     q = mat->element[3][1]; s = mat->element[3][3];
  807.     v = mat->element[4][1]; x = mat->element[4][3];
  808.  
  809.     mat->element[0][1] =   b * ct + d * st;
  810.     mat->element[0][3] = - b * st + d * ct;
  811.     mat->element[1][1] =   g * ct + i * st;
  812.     mat->element[1][3] = - g * st + i * ct;
  813.     mat->element[2][1] =   l * ct + n * st;
  814.     mat->element[2][3] = - l * st + n * ct;
  815.     mat->element[3][1] =   q * ct + s * st;
  816.     mat->element[3][3] = - q * st + s * ct;
  817.     mat->element[4][1] =   v * ct + x * st;
  818.     mat->element[4][3] = - v * st + x * ct;
  819.  
  820.     return mat;
  821. }
  822.  
  823. Matrix5    *
  824. V4MatRotateZX(mat, theta)
  825. Matrix5    *mat;
  826. double    theta;
  827. {
  828.     double    a, c, f, h, k, m, p, r, u, w;
  829.     double    st, ct;
  830.  
  831.     st = sin(theta);
  832.     ct = cos(theta);
  833.  
  834.     a = mat->element[0][0]; c = mat->element[0][2];
  835.     f = mat->element[1][0]; h = mat->element[1][2];
  836.     k = mat->element[2][0]; m = mat->element[2][2];
  837.     p = mat->element[3][0]; r = mat->element[3][2];
  838.     u = mat->element[4][0]; w = mat->element[4][2];
  839.  
  840.     mat->element[0][0] =   a * ct + c * st;
  841.     mat->element[0][2] = - a * st + c * ct;
  842.     mat->element[1][0] =   f * ct + h * st;
  843.     mat->element[1][2] = - f * st + h * ct;
  844.     mat->element[2][0] =   k * ct + m * st;
  845.     mat->element[2][2] = - k * st + m * ct;
  846.     mat->element[3][0] =   p * ct + r * st;
  847.     mat->element[3][2] = - p * st + r * ct;
  848.     mat->element[4][0] =   u * ct + w * st;
  849.     mat->element[4][2] = - u * st + w * ct;
  850.  
  851.     return mat;
  852. }
  853.  
  854.  
  855. Matrix5    *
  856. V4MatRotateZW(mat, theta)
  857. Matrix5    *mat;
  858. double    theta;
  859. {
  860.     double    c, d, h, i, m, n, r, s, w, x;
  861.     double    st, ct;
  862.  
  863.     st = sin(theta);
  864.     ct = cos(theta);
  865.  
  866.     c = mat->element[0][2]; d = mat->element[0][3];
  867.     h = mat->element[1][2]; i = mat->element[1][3];
  868.     m = mat->element[2][2]; n = mat->element[2][3];
  869.     r = mat->element[3][2]; s = mat->element[3][3];
  870.     w = mat->element[4][2]; x = mat->element[4][3];
  871.  
  872.     mat->element[0][2] =   c * ct + d * st;
  873.     mat->element[0][3] = - c * st + d * ct;
  874.     mat->element[1][2] =   h * ct + i * st;
  875.     mat->element[1][3] = - h * st + i * ct;
  876.     mat->element[2][2] =   m * ct + n * st;
  877.     mat->element[2][3] = - m * st + n * ct;
  878.     mat->element[3][2] =   r * ct + s * st;
  879.     mat->element[3][3] = - r * st + s * ct;
  880.     mat->element[4][2] =   w * ct + x * st;
  881.     mat->element[4][3] = - w * st + x * ct;
  882.  
  883.     return mat;
  884. }
  885.  
  886. Matrix5    *
  887. V4MatTranslate(mat, v)
  888. Matrix5    *mat;
  889. Vector4    *v;
  890. {
  891.     double    e, j, o, t, y;
  892.  
  893.     e = mat->element[0][4];
  894.     j = mat->element[1][4];
  895.     o = mat->element[2][4];
  896.     t = mat->element[3][4];
  897.     y = mat->element[4][4];
  898.  
  899.     if (e != 0.0)
  900.     {
  901.         mat->element[0][0] += e * v->x;
  902.         mat->element[0][1] += e * v->y;
  903.         mat->element[0][2] += e * v->z;
  904.         mat->element[0][3] += e * v->w;
  905.     }
  906.     if (j != 0.0)
  907.     {
  908.         mat->element[1][0] += j * v->x;
  909.         mat->element[1][1] += j * v->y;
  910.         mat->element[1][2] += j * v->z;
  911.         mat->element[1][3] += j * v->w;
  912.     }
  913.     if (o != 0.0)
  914.     {
  915.         mat->element[2][0] += o * v->x;
  916.         mat->element[2][1] += o * v->y;
  917.         mat->element[2][2] += o * v->z;
  918.         mat->element[2][3] += o * v->w;
  919.     }
  920.     if (t != 0.0)
  921.     {
  922.         mat->element[3][0] += t * v->x;
  923.         mat->element[3][1] += t * v->y;
  924.         mat->element[3][2] += t * v->z;
  925.         mat->element[3][3] += t * v->w;
  926.     }
  927.     if (y != 0.0)
  928.     {
  929.         mat->element[4][0] += y * v->x;
  930.         mat->element[4][1] += y * v->y;
  931.         mat->element[4][2] += y * v->z;
  932.         mat->element[4][3] += y * v->w;
  933.     }
  934.  
  935.     return mat;
  936. }
  937.  
  938. Matrix5 *
  939. V4MatScale(mat, v)
  940. Matrix5    *mat;
  941. Vector4    *v;
  942. {
  943.     if (v->x != 1.0)
  944.     {
  945.         mat->element[0][0] *= v->x;
  946.         mat->element[1][0] *= v->x;
  947.         mat->element[2][0] *= v->x;
  948.         mat->element[3][0] *= v->x;
  949.         mat->element[4][0] *= v->x;
  950.     }
  951.     if (v->y != 1.0)
  952.     {
  953.         mat->element[0][1] *= v->y;
  954.         mat->element[1][1] *= v->y;
  955.         mat->element[2][1] *= v->y;
  956.         mat->element[3][1] *= v->y;
  957.         mat->element[4][1] *= v->y;
  958.     }
  959.     if (v->z != 1.0)
  960.     {
  961.         mat->element[0][2] *= v->z;
  962.         mat->element[1][2] *= v->z;
  963.         mat->element[2][2] *= v->z;
  964.         mat->element[3][2] *= v->z;
  965.         mat->element[4][2] *= v->z;
  966.     }
  967.     if (v->w != 1.0)
  968.     {
  969.         mat->element[0][3] *= v->w;
  970.         mat->element[1][3] *= v->w;
  971.         mat->element[2][3] *= v->w;
  972.         mat->element[3][3] *= v->w;
  973.         mat->element[4][3] *= v->w;
  974.     }
  975.  
  976.     return mat;
  977. }
  978.  
  979. /* Postscript
  980.  *
  981.  * The elements of the 5x5 matrix are labelled in the above functions
  982.  * as follows:
  983.  *
  984.  * a = mat->element[0][0];
  985.  * b = mat->element[0][1];
  986.  * c = mat->element[0][2];
  987.  * d = mat->element[0][3];
  988.  * e = mat->element[0][4];
  989.  * f = mat->element[1][0];
  990.  * g = mat->element[1][1];
  991.  * h = mat->element[1][2];
  992.  * i = mat->element[1][3];
  993.  * j = mat->element[1][4];
  994.  * k = mat->element[2][0];
  995.  * l = mat->element[2][1];
  996.  * m = mat->element[2][2];
  997.  * n = mat->element[2][3];
  998.  * o = mat->element[2][4];
  999.  * p = mat->element[3][0];
  1000.  * q = mat->element[3][1];
  1001.  * r = mat->element[3][2];
  1002.  * s = mat->element[3][3];
  1003.  * t = mat->element[3][4];
  1004.  * u = mat->element[4][0];
  1005.  * v = mat->element[4][1];
  1006.  * w = mat->element[4][2];
  1007.  * x = mat->element[4][3];
  1008.  * y = mat->element[4][4];
  1009.  *
  1010.  * ie.
  1011.  * [ a b c d e ]
  1012.  * [ f g h i j ]
  1013.  * [ k l m n o ]
  1014.  * [ p q r s t ]
  1015.  * [ u v w x y ]
  1016.  *
  1017.  */
  1018.