home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1996 September / macformat-041.iso / mac / Shareware City / Graphics / MacSPD / Sources / libvec.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-04-04  |  15.9 KB  |  648 lines  |  [TEXT/MMCC]
  1. /*
  2.  * libvec.c - a library of vector operations and a random number generator
  3.  *
  4.  * Author:  Eric Haines, 3D/Eye, Inc.
  5.  *
  6.  * Modified: 1 October 1992
  7.  *           Alexander R. Enzmann
  8.  *          Object generation routines split off to keep file size down
  9.  *
  10.  * Modified: 17 Jan 1993
  11.  *         Eduard [esp] Schwan
  12.  *           Removed unused local variables, added <stdlib.h> include for
  13.  *          exit() call
  14.  *
  15.  * Modified: 1 November 1994
  16.  *           Alexander R. Enzmann
  17.  *          Added routines to invert matrices, transform normals, and
  18.  *          determine rotate/scale/translate from a transform matrix.
  19.  *          Modified computation of perspective view matrix to be a
  20.  *          little cleaner.
  21.  *
  22.  */
  23.  
  24. #include <stdio.h>
  25. #include <stdlib.h> /* exit() */
  26. #include <math.h>
  27. #include "libvec.h"
  28.  
  29. /*
  30.  * Normalize the vector (X,Y,Z) so that X*X + Y*Y + Z*Z = 1.
  31.  *
  32.  * The normalization divisor is returned.  If the divisor is zero, no
  33.  * normalization occurs.
  34.  *
  35.  */
  36. double
  37. lib_normalize_vector(cvec)
  38.     COORD3 cvec;
  39. {
  40.     double divisor;
  41.  
  42.     divisor = sqrt( (double)DOT_PRODUCT(cvec, cvec) );
  43.     if (divisor > 0.0) {
  44.        cvec[X] /= divisor;
  45.        cvec[Y] /= divisor;
  46.        cvec[Z] /= divisor;
  47.        }
  48.     return divisor;
  49. }
  50.  
  51.  
  52. /*
  53.  * Set all matrix elements to zero.
  54.  */
  55. void
  56. lib_zero_matrix(mx)
  57.     MATRIX mx;
  58. {
  59.     int i, j;
  60.  
  61.     for (i=0;i<4;++i)
  62.        for (j=0;j<4;++j)
  63.       mx[i][j] = 0.0;
  64. }
  65.  
  66. /*
  67.  * Create identity matrix.
  68.  */
  69. void
  70. lib_create_identity_matrix(mx)
  71.     MATRIX mx;
  72. {
  73.     int i;
  74.  
  75.     lib_zero_matrix(mx);
  76.     for (i=0;i<4;++i)
  77.        mx[i][i] = 1.0;
  78. }
  79.  
  80. /*
  81.  * Copy one matrix to another
  82.  */
  83. void
  84. lib_copy_matrix(mres, mx)
  85.     MATRIX mres, mx;
  86. {
  87.     int i, j;
  88.  
  89.     for (i=0;i<4;i++)
  90.        for (j=0;j<4;j++)
  91.       mres[i][j] = mx[i][j];
  92. }
  93.  
  94. /*
  95.  * Create a rotation matrix along the given axis by the given angle in radians.
  96.  */
  97. void
  98. lib_create_rotate_matrix(mx, axis, angle)
  99.     MATRIX mx;
  100.     int axis;
  101.     double angle;
  102. {
  103.     double cosine, sine;
  104.  
  105.     lib_zero_matrix(mx);
  106.     cosine = cos((double)angle);
  107.     sine = sin((double)angle);
  108.     switch (axis) {
  109.       case X_AXIS:
  110.      mx[0][0] = 1.0;
  111.      mx[1][1] = cosine;
  112.      mx[2][2] = cosine;
  113.      mx[1][2] = sine;
  114.      mx[2][1] = -sine;
  115.        break ;
  116.      case Y_AXIS:
  117.      mx[1][1] = 1.0;
  118.      mx[0][0] = cosine;
  119.      mx[2][2] = cosine;
  120.      mx[2][0] = sine;
  121.      mx[0][2] = -sine;
  122.        break;
  123.       case Z_AXIS:
  124.      mx[2][2] = 1.0;
  125.      mx[0][0] = cosine;
  126.      mx[1][1] = cosine;
  127.      mx[0][1] = sine;
  128.      mx[1][0] = -sine;
  129.        break;
  130.     }
  131.     mx[3][3] = 1.0;
  132. }
  133.  
  134. /*
  135.  * Create a rotation matrix along the given axis by the given angle in radians.
  136.  * The axis is a set of direction cosines.
  137.  */
  138. void
  139. lib_create_axis_rotate_matrix(mx, axis, angle)
  140.     MATRIX mx;
  141.     COORD3 axis;
  142.     double angle;
  143. {
  144.     double  cosine, one_minus_cosine, sine;
  145.  
  146.     cosine = cos((double)angle);
  147.     sine = sin((double)angle);
  148.     one_minus_cosine = 1.0 - cosine;
  149.  
  150.     mx[0][0] = SQR(axis[X]) + (1.0 - SQR(axis[X])) * cosine;
  151.     mx[0][1] = axis[X] * axis[Y] * one_minus_cosine + axis[Z] * sine;
  152.     mx[0][2] = axis[X] * axis[Z] * one_minus_cosine - axis[Y] * sine;
  153.     mx[0][3] = 0.0;
  154.  
  155.     mx[1][0] = axis[X] * axis[Y] * one_minus_cosine - axis[Z] * sine;
  156.     mx[1][1] = SQR(axis[Y]) + (1.0 - SQR(axis[Y])) * cosine;
  157.     mx[1][2] = axis[Y] * axis[Z] * one_minus_cosine + axis[X] * sine;
  158.     mx[1][3] = 0.0;
  159.  
  160.     mx[2][0] = axis[X] * axis[Z] * one_minus_cosine + axis[Y] * sine;
  161.     mx[2][1] = axis[Y] * axis[Z] * one_minus_cosine - axis[X] * sine;
  162.     mx[2][2] = SQR(axis[Z]) + (1.0 - SQR(axis[Z])) * cosine;
  163.     mx[2][3] = 0.0;
  164.  
  165.     mx[3][0] = 0.0;
  166.     mx[3][1] = 0.0;
  167.     mx[3][2] = 0.0;
  168.     mx[3][3] = 1.0;
  169. }
  170.  
  171. /* Create translation matrix */
  172. void
  173. lib_create_translate_matrix(mx, vec)
  174.     MATRIX mx;
  175.     COORD3 vec;
  176. {
  177.     lib_create_identity_matrix(mx);
  178.     mx[3][0] = vec[X];
  179.     mx[3][1] = vec[Y];
  180.     mx[3][2] = vec[Z];
  181. }
  182.  
  183. /* Create scaling matrix */
  184. void
  185. lib_create_scale_matrix(mx, vec)
  186.     MATRIX mx;
  187.     COORD3 vec;
  188. {
  189.     lib_zero_matrix(mx);
  190.     mx[0][0] = vec[X];
  191.     mx[1][1] = vec[Y];
  192.     mx[2][2] = vec[Z];
  193.     mx[3][3] = 1.0;
  194. }
  195.  
  196. /* Given a point and a direction, find the transform that brings a point
  197.    in a canonical coordinate system into a coordinate system defined by
  198.    that position and direction.    Both the forward and inverse transform
  199.    matrices are calculated. */
  200. void
  201. lib_create_canonical_matrix(trans, itrans, origin, up)
  202.     MATRIX trans, itrans;
  203.     COORD3 origin;
  204.     COORD3 up;
  205. {
  206.     MATRIX trans1, trans2;
  207.     COORD3 tmpv;
  208.     double ang;
  209.  
  210.     /* Translate "origin" to <0, 0, 0> */
  211.     SET_COORD3(tmpv, -origin[X], -origin[Y], -origin[Z]);
  212.     lib_create_translate_matrix(trans1, tmpv);
  213.  
  214.     /* Determine the axis to rotate about */
  215.     if (fabs(up[Z]) == 1.0)
  216.        SET_COORD3(tmpv, 1.0, 0.0, 0.0)
  217.     else
  218.        SET_COORD3(tmpv, -up[Y], up[X], 0.0)
  219.     lib_normalize_vector(tmpv);
  220.     ang = acos(up[Z]);
  221.  
  222.     /* Calculate the forward matrix */
  223.     lib_create_axis_rotate_matrix(trans2, tmpv, -ang);
  224.     lib_matrix_multiply(trans, trans1, trans2);
  225.  
  226.     /* Calculate the inverse transform */
  227.     lib_create_axis_rotate_matrix(trans1, tmpv, ang);
  228.     lib_create_translate_matrix(trans2, origin);
  229.     lib_matrix_multiply(itrans, trans1, trans2);
  230. }
  231.  
  232. #ifdef _DEBUG
  233. static void
  234. show_matrix(mat)
  235.    MATRIX mat;
  236. {
  237.    int i, j;
  238.    for (i=0;i<4;i++) {
  239.       for (j=0;j<4;j++)
  240.      fprintf(stderr, "%7.4g ", mat[i][j]);
  241.       fprintf(stderr, "\n");
  242.       }
  243. }
  244. #endif
  245.  
  246. /* Find the transformation that takes the current eye position and
  247.    orientation and puts it at {0, 0, 0} with the up vector aligned
  248.    with {0, 1, 0} */
  249. void
  250. lib_create_view_matrix(trans, from, at, up, xres, yres, angle, aspect)
  251.     MATRIX trans;
  252.     COORD3 from, at, up;
  253.     int xres, yres;
  254.     double angle, aspect;
  255. {
  256.     double xscale, yscale, view_dist;
  257.     COORD3 right, up_dir;
  258.     COORD3 Va;
  259.     MATRIX Tv, Tt, Tx;
  260.  
  261.     /* Change the view angle into radians */
  262.     angle = PI * angle / 180.0;
  263.  
  264.     /* Set the hither clipping plane */
  265.     view_dist = 0.001;
  266.  
  267.     /* Translate point of interest to the origin */
  268.     SET_COORD3(Va, -from[X], -from[Y], -from[Z]);
  269.     lib_create_translate_matrix(Tv, Va);
  270.  
  271.     /* Move the eye by the same amount */
  272.     SUB3_COORD3(Va, at, from)
  273.  
  274.     /* Make sure this is a valid sort of setup */
  275.     if (Va[X] == 0.0 && Va[Y] == 0.0 && Va[Z] == 0.0) {
  276.       fprintf(stderr, "Degenerate perspective transformation\n");
  277.      exit(1);
  278.     }
  279.  
  280.     /* Get the up vector to be perpendicular to the view vector */
  281.     lib_normalize_vector(Va);
  282.     COPY_COORD3(up_dir, up);
  283.     CROSS(right, up_dir, Va);
  284.     lib_normalize_vector(right);
  285.     CROSS(up_dir, Va, right);
  286.     lib_normalize_vector(up_dir);
  287.  
  288.     /* Create the orthogonal view transformation */
  289.     Tt[0][0] = right[0];
  290.     Tt[1][0] = right[1];
  291.     Tt[2][0] = right[2];
  292.     Tt[3][0] = 0;
  293.  
  294.     Tt[0][1] = up_dir[0];
  295.     Tt[1][1] = up_dir[1];
  296.     Tt[2][1] = up_dir[2];
  297.     Tt[3][1] = 0;
  298.  
  299.     Tt[0][2] = Va[0];
  300.     Tt[1][2] = Va[1];
  301.     Tt[2][2] = Va[2];
  302.     Tt[3][2] = 0;
  303.  
  304.     Tt[0][3] = 0;
  305.     Tt[1][3] = 0;
  306.     Tt[2][3] = 0;
  307.     Tt[3][3] = 1;
  308.     lib_matrix_multiply(Tx, Tv, Tt);
  309.     lib_copy_matrix(Tv, Tx);
  310.  
  311.     /* Now add in the perspective transformation */
  312.     lib_create_identity_matrix(Tt);
  313.     Tt[2][2] =  1.0 / (1.0 + view_dist);
  314.     Tt[3][2] =  view_dist / (1.0 + view_dist);
  315.     Tt[2][3] =  1.0;
  316.     Tt[3][3] =  0.0;
  317.     lib_matrix_multiply(Tx, Tv, Tt);
  318.     lib_copy_matrix(Tv, Tx);
  319.  
  320.     /* Determine how much to scale things by */
  321.     yscale = (double)yres / (2.0 * tan(angle/2.0));
  322.     xscale = yscale * (double)xres / ((double)yres * aspect);
  323.     SET_COORD3(Va, -xscale, -yscale, 1.0);
  324.     lib_create_scale_matrix(Tt, Va);
  325.     lib_matrix_multiply(Tx, Tv, Tt);
  326.     lib_copy_matrix(Tv, Tx);
  327.  
  328.     /* Translate the points so that <0,0> is at the center of
  329.        the screen */
  330.     SET_COORD3(Va, xres/2, yres/2, 0.0);
  331.     lib_create_translate_matrix(Tt, Va);
  332.  
  333.     /* We now have the final transformation matrix */
  334.     lib_matrix_multiply(trans, Tv, Tt);
  335. }
  336.  
  337. /*
  338.  * Multiply a vector by a matrix.
  339.  */
  340. void
  341. lib_transform_vector(vres, vec, mx)
  342.     COORD3 vres, vec;
  343.     MATRIX mx;
  344. {
  345.     COORD3 vtemp;
  346.     vtemp[X] = vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0];
  347.     vtemp[Y] = vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1];
  348.     vtemp[Z] = vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2];
  349.     COPY_COORD3(vres, vtemp);
  350. }
  351.  
  352. /*
  353.  * Multiply a normal by a matrix (i.e. multiply by transpose).
  354.  */
  355. void
  356. lib_transform_normal(vres, vec, mx)
  357.     COORD3 vres, vec;
  358.     MATRIX mx;
  359. {
  360.     COORD3 vtemp;
  361.     vtemp[X] = vec[X]*mx[0][0] + vec[Y]*mx[0][1] + vec[Z]*mx[0][2];
  362.     vtemp[Y] = vec[X]*mx[1][0] + vec[Y]*mx[1][1] + vec[Z]*mx[1][2];
  363.     vtemp[Z] = vec[X]*mx[2][0] + vec[Y]*mx[2][1] + vec[Z]*mx[2][2];
  364.     COPY_COORD3(vres, vtemp);
  365. }
  366.  
  367. /*
  368.  * Multiply a point by a matrix.
  369.  */
  370. void
  371. lib_transform_point(vres, vec, mx)
  372.     COORD3 vres, vec;
  373.     MATRIX mx;
  374. {
  375.     COORD3 vtemp;
  376.     vtemp[X] = vec[X]*mx[0][0] + vec[Y]*mx[1][0] + vec[Z]*mx[2][0] + mx[3][0];
  377.     vtemp[Y] = vec[X]*mx[0][1] + vec[Y]*mx[1][1] + vec[Z]*mx[2][1] + mx[3][1];
  378.     vtemp[Z] = vec[X]*mx[0][2] + vec[Y]*mx[1][2] + vec[Z]*mx[2][2] + mx[3][2];
  379.     COPY_COORD3(vres, vtemp);
  380. }
  381.  
  382. /*
  383.  * Multiply a 4 element vector by a matrix.  Typically used for
  384.  * homogenous transformation from world space to screen space.
  385.  */
  386. void
  387. lib_transform_coord(vres, vec, mx)
  388.     COORD4 vres, vec;
  389.     MATRIX mx;
  390. {
  391.     COORD4 vtemp;
  392.     vtemp[X] = vec[X]*mx[0][0]+vec[Y]*mx[1][0]+vec[Z]*mx[2][0]+vec[W]*mx[3][0];
  393.     vtemp[Y] = vec[X]*mx[0][1]+vec[Y]*mx[1][1]+vec[Z]*mx[2][1]+vec[W]*mx[3][1];
  394.     vtemp[Z] = vec[X]*mx[0][2]+vec[Y]*mx[1][2]+vec[Z]*mx[2][2]+vec[W]*mx[3][2];
  395.     vtemp[W] = vec[X]*mx[0][3]+vec[Y]*mx[1][3]+vec[Z]*mx[2][3]+vec[W]*mx[3][3];
  396.     COPY_COORD4(vres, vtemp);
  397. }
  398.  
  399. /* Determinant of a 3x3 matrix */
  400. static double
  401. det3x3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
  402.    double a1, a2, a3;
  403.    double b1, b2, b3;
  404.    double c1, c2, c3;
  405. {
  406.     double ans;
  407.  
  408.     ans = a1 * (b2 * c3 - b3 * c2) -
  409.       b1 * (a2 * c3 - a3 * c2) +
  410.       c1 * (a2 * b3 - a3 * b2);
  411.     return ans;
  412. }
  413.  
  414. /* Adjoint of a 4x4 matrix - used in the computation of the inverse
  415.    of a 4x4 matrix */
  416. static void
  417. adjoint(out_mat, in_mat)
  418.    MATRIX out_mat, in_mat;
  419. {
  420.    double a1, a2, a3, a4, b1, b2, b3, b4;
  421.    double c1, c2, c3, c4, d1, d2, d3, d4;
  422.  
  423.    a1 = in_mat[0][0]; b1 = in_mat[0][1];
  424.    c1 = in_mat[0][2]; d1 = in_mat[0][3];
  425.    a2 = in_mat[1][0]; b2 = in_mat[1][1];
  426.    c2 = in_mat[1][2]; d2 = in_mat[1][3];
  427.    a3 = in_mat[2][0]; b3 = in_mat[2][1];
  428.    c3 = in_mat[2][2]; d3 = in_mat[2][3];
  429.    a4 = in_mat[3][0]; b4 = in_mat[3][1];
  430.    c4 = in_mat[3][2]; d4 = in_mat[3][3];
  431.  
  432.  
  433.     /* row column labeling reversed since we transpose rows & columns */
  434.    out_mat[0][0]  =   det3x3( b2, b3, b4, c2, c3, c4, d2, d3, d4);
  435.    out_mat[1][0]  = - det3x3( a2, a3, a4, c2, c3, c4, d2, d3, d4);
  436.    out_mat[2][0]  =   det3x3( a2, a3, a4, b2, b3, b4, d2, d3, d4);
  437.    out_mat[3][0]  = - det3x3( a2, a3, a4, b2, b3, b4, c2, c3, c4);
  438.  
  439.    out_mat[0][1]  = - det3x3( b1, b3, b4, c1, c3, c4, d1, d3, d4);
  440.    out_mat[1][1]  =   det3x3( a1, a3, a4, c1, c3, c4, d1, d3, d4);
  441.    out_mat[2][1]  = - det3x3( a1, a3, a4, b1, b3, b4, d1, d3, d4);
  442.    out_mat[3][1]  =   det3x3( a1, a3, a4, b1, b3, b4, c1, c3, c4);
  443.  
  444.    out_mat[0][2]  =   det3x3( b1, b2, b4, c1, c2, c4, d1, d2, d4);
  445.    out_mat[1][2]  = - det3x3( a1, a2, a4, c1, c2, c4, d1, d2, d4);
  446.    out_mat[2][2]  =   det3x3( a1, a2, a4, b1, b2, b4, d1, d2, d4);
  447.    out_mat[3][2]  = - det3x3( a1, a2, a4, b1, b2, b4, c1, c2, c4);
  448.  
  449.    out_mat[0][3]  = - det3x3( b1, b2, b3, c1, c2, c3, d1, d2, d3);
  450.    out_mat[1][3]  =   det3x3( a1, a2, a3, c1, c2, c3, d1, d2, d3);
  451.    out_mat[2][3]  = - det3x3( a1, a2, a3, b1, b2, b3, d1, d2, d3);
  452.    out_mat[3][3]  =   det3x3( a1, a2, a3, b1, b2, b3, c1, c2, c3);
  453. }
  454.  
  455. /* Determinant of a 4x4 matrix */
  456. double
  457. lib_matrix_det4x4(mat)
  458.    MATRIX mat;
  459. {
  460.    double ans;
  461.    double a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4;
  462.  
  463.    a1 = mat[0][0]; b1 = mat[0][1]; c1 = mat[0][2]; d1 = mat[0][3];
  464.    a2 = mat[1][0]; b2 = mat[1][1]; c2 = mat[1][2]; d2 = mat[1][3];
  465.    a3 = mat[2][0]; b3 = mat[2][1]; c3 = mat[2][2]; d3 = mat[2][3];
  466.    a4 = mat[3][0]; b4 = mat[3][1]; c4 = mat[3][2]; d4 = mat[3][3];
  467.  
  468.    ans = a1 * det3x3(b2, b3, b4, c2, c3, c4, d2, d3, d4) -
  469.      b1 * det3x3(a2, a3, a4, c2, c3, c4, d2, d3, d4) +
  470.      c1 * det3x3(a2, a3, a4, b2, b3, b4, d2, d3, d4) -
  471.      d1 * det3x3(a2, a3, a4, b2, b3, b4, c2, c3, c4);
  472.    return ans;
  473. }
  474.  
  475. /* Find the inverse of a 4x4 matrix */
  476. void
  477. lib_invert_matrix(out_mat, in_mat)
  478.    MATRIX out_mat, in_mat;
  479. {
  480.    int i, j;
  481.    double det;
  482.  
  483.    adjoint(out_mat, in_mat);
  484.    det = lib_matrix_det4x4(in_mat);
  485.    if (fabs(det) < EPSILON) {
  486.       lib_create_identity_matrix(out_mat);
  487.       return;
  488.       }
  489.    for (i=0;i<4;i++)
  490.       for (j=0;j<4;j++)
  491.      out_mat[i][j] /= det;
  492. }
  493.  
  494. /*
  495.  * Compute transpose of matrix.
  496.  */
  497. void
  498. lib_transpose_matrix( mxres, mx )
  499.     MATRIX mxres, mx;
  500. {
  501.     int i, j;
  502.  
  503.     for (i=0;i<4;i++)
  504.        for (j=0;j<4;j++)
  505.        mxres[j][i] = mx[i][j];
  506. }
  507.  
  508. /*
  509.  * Multiply two 4x4 matrices.  Note that mxres had better not be
  510.  * the same as either mx1 or mx2 or bad results will be returned.
  511.  */
  512. void
  513. lib_matrix_multiply(mxres, mx1, mx2)
  514.     MATRIX mxres, mx1, mx2;
  515. {
  516.     int i, j;
  517.  
  518.     for (i=0;i<4;i++)
  519.        for (j=0;j<4;j++)
  520.       mxres[i][j] = mx1[i][0]*mx2[0][j] + mx1[i][1]*mx2[1][j] +
  521.             mx1[i][2]*mx2[2][j] + mx1[i][3]*mx2[3][j];
  522. }
  523.  
  524. /* Performs a 3D clip of a line segment from start to end against the
  525.     box defined by bounds.  The actual values of start and end are modified. */
  526. int
  527. lib_clip_to_box(start, end,  bounds)
  528.     COORD3 start, end;
  529.     double bounds[2][3];
  530. {
  531.     int i;
  532.     double d, t, dir, pos;
  533.     double tmin, tmax;
  534.     COORD3 D;
  535.  
  536.     SUB3_COORD3(D, end, start);
  537.     d = lib_normalize_vector(D);
  538.     tmin = 0.0;
  539.     tmax = d;
  540.     for (i=0;i<3;i++) {
  541.     pos = start[i];
  542.     dir = D[i];
  543.      if (dir < -EPSILON) {
  544.        t = (bounds[0][i] - pos) / dir;
  545.      if (t < tmin)
  546.            return 0;
  547.        if (t <= tmax)
  548.           tmax = t;
  549.        t = (bounds[1][i] - pos) / dir;
  550.      if (t >= tmin) {
  551.         if (t > tmax)
  552.            return 0;
  553.        tmin = t;
  554.        }
  555.    } else if (dir > EPSILON) {
  556.      t = (bounds[1][i] - pos) / dir;
  557.      if (t < tmin)
  558.            return 0;
  559.        if (t <= tmax)
  560.           tmax = t;
  561.        t = (bounds[0][i] - pos) / dir;
  562.      if (t >= tmin) {
  563.         if (t > tmax)
  564.            return 0;
  565.        tmin = t;
  566.        }
  567.    }
  568.        else if (pos < bounds[0][i] || pos > bounds[1][i])
  569.       return 0;
  570.     }
  571.     if (tmax < d) {
  572.     end[X] = start[X] + tmax * D[X];
  573.     end[Y] = start[Y] + tmax * D[Y];
  574.     end[Z] = start[Z] + tmax * D[Z];
  575.     }
  576.     if (tmin > 0.0) {
  577.       start[X] = start[X] + tmin * D[X];
  578.       start[Y] = start[Y] + tmin * D[Y];
  579.       start[Z] = start[Z] + tmin * D[Z];
  580.     }
  581.     return 1;
  582. }
  583.  
  584. /*
  585.  * Rotate a vector pointing towards the major-axis faces (i.e. the major-axis
  586.  * component of the vector is defined as the largest value) 90 degrees to
  587.  * another cube face.  Mod_face is a face number.
  588.  *
  589.  * If the routine is called six times, with mod_face=0..5, the vector will be
  590.  * rotated to each face of a cube.  Rotations are:
  591.  *   mod_face = 0 mod 3, +Z axis rotate
  592.  *   mod_face = 1 mod 3, +X axis rotate
  593.  *   mod_face = 2 mod 3, -Y axis rotate
  594.  */
  595. void
  596. lib_rotate_cube_face(vec, major_axis, mod_face)
  597.     COORD3 vec;
  598.     int major_axis, mod_face;
  599. {
  600.     double swap;
  601.  
  602.     mod_face = (mod_face+major_axis) % 3;
  603.     if (mod_face == 0) {
  604.     swap   = vec[X];
  605.     vec[X] = -vec[Y];
  606.     vec[Y] = swap;
  607.     } else if (mod_face == 1) {
  608.     swap   = vec[Y];
  609.     vec[Y] = -vec[Z];
  610.     vec[Z] = swap;
  611.     } else {
  612.     swap   = vec[X];
  613.     vec[X] = -vec[Z];
  614.     vec[Z] = swap;
  615.     }
  616. }
  617.  
  618. /*
  619.  * Portable gaussian random number generator (from "Numerical Recipes", GASDEV)
  620.  * Returns a uniform random deviate between 0.0 and 1.0.  'iseed' must be
  621.  * less than M1 to avoid repetition, and less than (2**31-C1)/A1 [= 300718]
  622.  * to avoid overflow.
  623.  */
  624. #define M1  134456
  625. #define IA1   8121
  626. #define IC1  28411
  627. #define RM1 1.0/M1
  628.  
  629. double
  630. lib_gauss_rand(iseed)
  631.     long iseed;
  632. {
  633.     double fac, v1, v2, r;
  634.     long     ix1, ix2;
  635.  
  636.     ix2 = iseed;
  637.     do {
  638.     ix1 = (IC1+ix2*IA1) % M1;
  639.        ix2 = (IC1+ix1*IA1) % M1;
  640.        v1 = ix1 * 2.0 * RM1 - 1.0;
  641.      v2 = ix2 * 2.0 * RM1 - 1.0;
  642.      r = v1*v1 + v2*v2;
  643.     } while ( r >= 1.0 );
  644.  
  645.     fac = sqrt((double)(-2.0 * log((double)r) / r));
  646.     return v1 * fac;
  647. }
  648.