home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Professional / OS2PRO194.ISO / os2 / graphic / dkb / source / quartics.c < prev    next >
C/C++ Source or Header  |  1992-07-28  |  19KB  |  522 lines
  1. /*****************************************************************************
  2. *
  3. *                                   quartics.c
  4. *
  5. *   from DKBTrace (c) 1990  David Buck
  6. *
  7. *  This module implements the code for the quartic (4th order) shapes.
  8. *
  9. *  This file was written by Alexander Enzmann.  He wrote the code for DKB 2.11
  10. *  QUARTICS (4th order shapes) and generously provided us these enhancements.
  11. *
  12. * This software is freely distributable. The source and/or object code may be
  13. * copied or uploaded to communications services so long as this notice remains
  14. * at the top of each file.  If any changes are made to the program, you must
  15. * clearly indicate in the documentation and in the programs startup message
  16. * who it was who made the changes. The documentation should also describe what
  17. * those changes were. This software may not be included in whole or in
  18. * part into any commercial package without the express written consent of the
  19. * author.  It may, however, be included in other public domain or freely
  20. * distributed software so long as the proper credit for the software is given.
  21. *
  22. * This software is provided as is without any guarantees or warranty. Although
  23. * the author has attempted to find and correct any bugs in the software, he
  24. * is not responsible for any damage caused by the use of the software.  The
  25. * author is under no obligation to provide service, corrections, or upgrades
  26. * to this package.
  27. *
  28. * Despite all the legal stuff above, if you do find bugs, I would like to hear
  29. * about them.  Also, if you have any comments or questions, you may contact me
  30. * at the following address:
  31. *
  32. *     David Buck
  33. *     22C Sonnet Cres.
  34. *     Nepean Ontario
  35. *     Canada, K2H 8W7
  36. *
  37. *  I can also be reached on the following bulleton boards:
  38. *
  39. *     OMX              (613) 731-3419
  40. *     Mystic           (613) 596-4249  or  (613) 596-4772
  41. *
  42. *  Fidonet:   1:163/109.9
  43. *  Internet:  dbuck@ccs.carleton.ca
  44. *  The "You Can Call Me RAY" BBS    (708) 358-5611
  45. *
  46. *  IBM Port by Aaron A. Collins. Aaron may be reached on the following BBS'es:
  47. *
  48. *     The "You Can Call Me RAY" BBS (708) 358-5611
  49. *     The Information Exchange BBS  (708) 945-5575
  50. *
  51. *****************************************************************************/
  52.  
  53. #include "frame.h"
  54. #include "vector.h"
  55. #include "dkbproto.h"
  56.  
  57.  
  58. /* Basic form of a quartic equation
  59.    a00*x^4 + a01*x^3*y + a02*x^3*z + a03*x^3 + a04*x^2*y^2 + 
  60.    a05*x^2*y*z + a06*x^2*y + a07*x^2*z^2 + a08*x^2*z + a09*x^2 + 
  61.    a10*x*y^3 + a11*x*y^2*z + a12*x*y^2 + a13*x*y*z^2 + a14*x*y*z + 
  62.    a15*x*y + a16*x*z^3 + a17*x*z^2 + a18*x*z + a19*x + a20*y^4 + 
  63.    a21*y^3*z + a22*y^3 + a23*y^2*z^2 + a24*y^2*z + a25*y^2 + a26*y*z^3 + 
  64.    a27*y*z^2 + a28*y*z + a29*y + a30*z^4 + a31*z^3 + a32*z^2 + a33*z + a34
  65. */
  66.  
  67. int factorials [5] = {1, 1, 2, 6, 24};    /* 0!, 1!, 2!, 3!, 4! */
  68. int term_counts [5] = {1, 4, 10, 20, 35};
  69. int terms4 [35] [4] = {
  70.    {4, 0, 0, 0 }, {3, 1, 0, 0 }, {3, 0, 1, 0 }, {3, 0, 0, 1 }, {2, 2, 0, 0 },
  71.    {2, 1, 1, 0 }, {2, 1, 0, 1 }, {2, 0, 2, 0 }, {2, 0, 1, 1 }, {2, 0, 0, 2 },
  72.    {1, 3, 0, 0 }, {1, 2, 1, 0 }, {1, 2, 0, 1 }, {1, 1, 2, 0 }, {1, 1, 1, 1 },
  73.    {1, 1, 0, 2 }, {1, 0, 3, 0 }, {1, 0, 2, 1 }, {1, 0, 1, 2 }, {1, 0, 0, 3 },
  74.    {0, 4, 0, 0 }, {0, 3, 1, 0 }, {0, 3, 0, 1 }, {0, 2, 2, 0 }, {0, 2, 1, 1 },
  75.    {0, 2, 0, 2 }, {0, 1, 3, 0 }, {0, 1, 2, 1 }, {0, 1, 1, 2 }, {0, 1, 0, 3 },
  76.    {0, 0, 4, 0 }, {0, 0, 3, 1 }, {0, 0, 2, 2 }, {0, 0, 1, 3 }, {0, 0, 0, 4 },
  77.    };
  78.  
  79. int terms3 [20] [4] = {
  80.    {3, 0, 0, 0 }, {2, 1, 0, 0 }, {2, 0, 1, 0 }, {2, 0, 0, 1 }, {1, 2, 0, 0 },
  81.    {1, 1, 1, 0 }, {1, 1, 0, 1 }, {1, 0, 2, 0 }, {1, 0, 1, 1 }, {1, 0, 0, 2 },
  82.    {0, 3, 0, 0 }, {0, 2, 1, 0 }, {0, 2, 0, 1 }, {0, 1, 2, 0 }, {0, 1, 1, 1 },
  83.    {0, 1, 0, 2 }, {0, 0, 3, 0 }, {0, 0, 2, 1 }, {0, 0, 1, 2 }, {0, 0, 0, 3 },
  84.    };
  85.  
  86. int terms2 [10] [4] = {
  87.    {2, 0, 0, 0 }, {1, 1, 0, 0 }, {1, 0, 1, 0 }, {1, 0, 0, 1 }, {0, 2, 0, 0 },
  88.    {0, 1, 1, 0 }, {0, 1, 0, 1 }, {0, 0, 2, 0 }, {0, 0, 1, 1 }, {0, 0, 0, 2 },
  89.    };
  90.  
  91. int terms1 [4] [4] = {
  92.    {1, 0, 0, 0 }, {0, 1, 0, 0 }, {0, 0, 1, 0 }, {0, 0, 0, 1 },
  93.    };
  94.  
  95.  
  96. METHODS Quartic_Methods =
  97.    { Object_Intersect, All_Quartic_Intersections,
  98.      Inside_Quartic, Quartic_Normal, Copy_Quartic,
  99.      Translate_Quartic, Rotate_Quartic,
  100.      Scale_Quartic, Invert_Quartic};
  101.  
  102. extern long Ray_Quartic_Tests, Ray_Quartic_Tests_Succeeded;
  103.  
  104. int
  105. All_Quartic_Intersections(Object, Ray, Depth_Queue)
  106.    OBJECT *Object;
  107.    RAY *Ray;
  108.    PRIOQ *Depth_Queue;
  109. {
  110.    QUARTIC *Shape = (QUARTIC *) Object;
  111.    DBL Depths [4];
  112.    VECTOR Intersection_Point;
  113.    INTERSECTION Local_Element;
  114.    int cnt, i, j, Intersection_Found;
  115.  
  116.    Intersection_Found = FALSE;
  117.    Ray_Quartic_Tests++;
  118.    cnt = Intersect_Quartic(Ray, Shape, &Depths [0]);
  119.    if (cnt > 0) Ray_Quartic_Tests_Succeeded++;
  120.    for (i=0;i<cnt;i++) {
  121.       for (j=0;j<i;j++)
  122.      if (Depths [i] == Depths [j]) goto l0;
  123.  
  124.       if (Depths [i] < Small_Tolerance)
  125.          continue;
  126.  
  127.       Local_Element.Depth = Depths [i];
  128.       Local_Element.Object = Shape -> Parent_Object;
  129.       VScale (Intersection_Point, Ray -> Direction, Depths [i]);
  130.       VAdd (Intersection_Point, Intersection_Point, Ray -> Initial);
  131.       Local_Element.Point = Intersection_Point;
  132.       Local_Element.Shape = (SHAPE *)Shape;
  133.       pq_add (Depth_Queue, &Local_Element);
  134.       Intersection_Found = TRUE;
  135.       l0:;
  136.       }
  137.    return (Intersection_Found);
  138. }
  139.  
  140. /* Intersection of a ray and a quartic */
  141. int
  142. Intersect_Quartic(Ray, Shape, Depths)
  143.    RAY *Ray;
  144.    QUARTIC *Shape;
  145.    DBL *Depths;
  146. {
  147.    DBL x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  148.    DBL xx,yy,zz,xx2,yy2,zz2,xx3,yy3,zz3,xx4,yy4,zz4;
  149.    DBL *a, t [5];
  150.  
  151.    x = Ray->Initial.x;
  152.    y = Ray->Initial.y;
  153.    z = Ray->Initial.z;
  154.    xx = Ray->Direction.x;
  155.    yy = Ray->Direction.y;
  156.    zz = Ray->Direction.z;
  157.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  158.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  159.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  160.    xx2 = xx * xx;  yy2 = yy * yy;  zz2 = zz * zz;
  161.    xx3 = xx * xx2; yy3 = yy * yy2; zz3 = zz * zz2;
  162.    xx4 = xx * xx3; yy4 = yy * yy3; zz4 = zz * zz3;
  163.    a = Shape->Coeffs;
  164.  
  165.    /*
  166.       Determine the coefficients of t^n, where the line is represented
  167.       as (x,y,z) + (xx,yy,zz) * t.
  168.    */
  169.    t [0] = a [0] * xx4 + a [1] * xx3 * yy + a [2] * xx3 * zz
  170.             + a [4] * xx2 * yy2 + a [5] * xx2 * yy * zz + a [7] * xx2 * zz2
  171.             + a [10] * xx * yy3 + a [11] * xx * yy2 * zz + a [13] * xx * yy * zz2
  172.             + a [16] * xx * zz3 + a [20] * yy4 + a [21] * yy3 * zz
  173.             + a [23] * yy2 * zz2 + a [26] * yy * zz3 + a [30] * zz4;
  174.  
  175.    t [1] = 4 * a [0] * x * xx3 + a [1] * (3 * x * xx2 * yy + xx3 * y)
  176.             + a [2] * (3 * x * xx2 * zz + xx3 * z) + a [3] * xx3
  177.             + a [4] * (2 * x * xx * yy2 + 2 * xx2 * y * yy)
  178.             + a [5] * (xx2 * (y * zz + yy * z) + 2 * x * xx * yy * zz);
  179.     t [1] += a [6] * xx2 * yy + a [7] * (2 * x * xx * zz2 + 2 * xx2 * z * zz)
  180.             + a [8] * xx2 * zz + a [10] * (x * yy3 + 3 * xx * y * yy2)
  181.             + a [11] * (xx * (2 * y * yy * zz + yy2 * z) + x * yy2 * zz)
  182.             + a [12] * xx * yy2;
  183.     t [1] += a [13] * (xx * (y * zz2 + 2 * yy * z * zz) + x * yy * zz2)
  184.             + a [14] * xx * yy * zz + a [16] * (x * zz3 + 3 * xx * z * zz2)
  185.             + a [17] * xx * zz2 + 4 * a [20] * y * yy3;
  186.     t [1] += a [21] * (3 * y * yy2 * zz + yy3 * z)
  187.             + a [22] * yy3 + zz * (2 * a [23] * yy * (y * zz + yy * z)
  188.             + a [24] * yy2 + zz * (a [26] * (y * zz + 3 * yy * z)
  189.             + a [27] * yy + zz * (4 * a [30] * z + a [31])));
  190.  
  191.    t [2] = 6 * a [0] * x2 * xx2 + 3 * a [1] * x * xx * (x * yy + xx * y)
  192.             + 3 * a [2] * x * xx * (x * zz + xx * z) + 3 * a [3] * x * xx2
  193.             + a [4] * (x2 * yy2 + 4 * x * xx * y * yy + xx2 * y2)
  194.             + a [5] * (x2 * yy * zz + 2 * x * xx * (y * zz + yy * z)
  195.             + xx2 * y * z) + a [6] * (2 * x * xx * yy + xx2 * y);
  196.     t [2] += a [7] * (x2 * zz2 + 4 * x * xx * z * zz + xx2 * z2)
  197.             + a [8] * (2 * x * xx * zz + xx2 * z) + a [9] * xx2
  198.             + a [10] * (3 * x * y * yy2 + 3 * xx * y2 * yy)
  199.             + a [11] * (x * yy * (2 * y * zz + yy * z)
  200.             + xx * (y2 * zz + 2 * y * yy * z))
  201.             + a [12] * (x * yy2 + 2 * xx * y * yy);
  202.     t [2] += a [13] * (x * zz * (y * zz + 2 * yy * z)
  203.             + xx * (2 * y * z * zz + yy * z2))
  204.             + a [14] * (x * yy * zz + xx * (y * zz + yy * z)) + a [15] * xx * yy
  205.             + a [16] * (3 * x * z * zz2 + 3 * xx * z2 * zz)
  206.             + a [17] * (x * zz2 + 2 * xx * z * zz) + a [18] * xx * zz;
  207.     t [2] += 6 * a [20] * y2 * yy2 + 3 * a [21] * y * yy * (y * zz + yy * z)
  208.             + 3 * a [22] * y * yy2 + a [23] * (y2 * zz2 + 4 * y * yy * z * zz
  209.             + yy2 * z2) + a [24] * (2 * y * yy * zz + yy2 * z) + a [25] * yy2
  210.             + zz * (3 * a [26] * z * (y * zz + yy * z)
  211.             + a [27] * (y * zz + 2 * yy * z) + a [28] * yy
  212.             + 6 * a [30] * z2 * zz + 3 * a [31] * z * zz + a [32] * zz);
  213.  
  214.    t [3] = 4 * a [0] * x3 * xx + a [1] * x2 * (x * yy + 3 * xx * y)
  215.             + a [2] * x2 * (x * zz + 3 * xx * z) + 3 * a [3] * x2 * xx
  216.             + 2 * a [4] * x * y * (x * yy + xx * y)
  217.             + a [5] * x * (x * (y * zz + yy * z) + 2 * xx * y * z)
  218.             + a [6] * x * (x * yy + 2 * xx * y);
  219.     t [3] += 2 * a [7] * x * z * (x * zz + xx * z)
  220.             + a [8] * x * (x * zz + 2 * xx * z)
  221.             + 2 * a [9] * x * xx + a [10] * (3 * x * y2 * yy-xx * y3)
  222.             + a [11] * (x * y * (y * zz + 2 * yy * z) + xx * y2 * z)
  223.             + a [12] * (2 * x * y * yy + xx * y2)
  224.             + a [13] * (x * z * (2 * y * zz + yy * z) + xx * y * z2);
  225.     t [3] += a [14] * (x * (y * zz + yy * z) + xx * y * z)
  226.             + a [15] * (x * yy + xx * y) + a [16] * (3 * x * z2 * zz + xx * z3)
  227.             + a [17] * (2 * x * z * zz + xx * z2) + a [18] * (x * zz + xx * z)
  228.             + a [19] * xx + 4 * a [20] * y3 * yy
  229.             + a [21] * y2 * (y * zz + 3 * yy * z) + 3 * a [22] * y2 * yy;
  230.     t [3] += 2 * a [23] * y * z * (y * zz + yy * z)
  231.             + a [24] * y * (y * zz + 2 * yy * z) + 2 * a [25] * y * yy
  232.             + a [26] * (3 * y * z2 * zz + yy * z3)
  233.             + a [27] * (2 * y * z * zz + yy * z2) + a [28] * (y * zz + yy * z)
  234.             + a [29] * yy + zz * (4 * a [30] * z3 + 3 * a [31] * z2
  235.             + 2 * a [32] * z + a [33]);
  236.  
  237.    t [4] = a [0] * x4 + a [1] * x3 * y + a [2] * x3 * z + a [3] * x3
  238.             + a [4] * x2 * y2 + a [5] * x2 * y * z + a [6] * x2 * y
  239.             + a [7] * x2 * z2 + a [8] * x2 * z + a [9] * x2 + a [10] * x * y3
  240.             + a [11] * x * y2 * z + a [12] * x * y2 + a [13] * x * y * z2
  241.             + a [14] * x * y * z + a [15] * x * y + a [16] * x * z3;
  242.     t [4] += a [17] * x * z2 + a [18] * x * z + a [19] * x + a [20] * y4
  243.             + a [21] * y3 * z + a [22] * y3 + a [23] * y2 * z2 + a [24] * y2 * z
  244.             + a [25] * y2 + a [26] * y * z3 + a [27] * y * z2 + a [28] * y * z
  245.             + a [29] * y + a [30] * z4 + a [31] * z3 + a [32] * z2 + a [33] * z
  246.             + a [34];
  247.    return solve_quartic (t, Depths);
  248. }
  249.  
  250. int
  251. Inside_Quartic (point, Object)
  252.    VECTOR *point;
  253.    OBJECT *Object;
  254. {
  255.    QUARTIC *Shape = (QUARTIC *) Object;
  256.    DBL Result,*a,x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  257.  
  258.    x = point->x; y = point->y; z = point->z;
  259.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  260.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  261.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  262.    a = Shape->Coeffs;
  263.  
  264.    Result = a [0]*x4 + a [1]*x3*y + a [2]*x3*z + a [3]*x3 + a [4]*x2*y2 +
  265.         a [5]*x2*y*z + a [6]*x2*y + a [7]*x2*z2 + a [8]*x2*z + a [9]*x2 +
  266.         a [10]*x*y3 + a [11]*x*y2*z + a [12]*x*y2 + a [13]*x*y*z2 + a [14]*x*y*z +
  267.         a [15]*x*y + a [16]*x*z3 + a [17]*x*z2 + a [18]*x*z + a [19]*x + a [20]*y4 +
  268.         a [21]*y3*z + a [22]*y3 + a [23]*y2*z2 + a [24]*y2*z + a [25]*y2 +
  269.         a [26]*y*z3 + a [27]*y*z2 + a [28]*y*z + a [29]*y + a [30]*z4 + a [31]*z3 +
  270.         a [32]*z2 + a [33]*z + a [34];
  271.  
  272.    if (Result < Small_Tolerance)
  273.       return (TRUE);
  274.    else
  275.       return (FALSE);
  276. }
  277.  
  278. /* Normal to a quartic */
  279. void
  280. Quartic_Normal(Result, Object, Intersection_Point)
  281.    VECTOR *Result, *Intersection_Point;
  282.    OBJECT *Object;
  283. {
  284.    QUARTIC *Shape = (QUARTIC *) Object;
  285.    DBL *a,x,y,z,x2,y2,z2,x3,y3,z3,Len;
  286.  
  287.    a = Shape->Coeffs;
  288.    x = Intersection_Point->x;
  289.    y = Intersection_Point->y;
  290.    z = Intersection_Point->z;
  291.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  292.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  293.  
  294.    Result->x = 4*a [0]*x3 + 3*x2*(a [1]*y + a [2]*z + a [3]) +
  295.            2*x*(a [4]*y2 + y*(a [5]*z + a [6]) + a [7]*z2 + a [8]*z + a [9]) +
  296.            a [10]*y3 + y2*(a [11]*z + a [12]) + y*(a [13]*z2 + a [14]*z + a [15]) +
  297.            a [16]*z3 + a [17]*z2 + a [18]*z + a [19];
  298.  
  299.    Result->y = a [1]*x3 + x2*(2*a [4]*y + a [5]*z + a [6]) +
  300.            x*(3*a [10]*y2 + 2*y*(a [11]*z + a [12]) + a [13]*z2 + a [14]*z + a [15]) +
  301.            4*a [20]*y3 + 3*y2*(a [21]*z + a [22]) + 2*y*(a [23]*z2 + a [24]*z + a [25]) +
  302.            a [26]*z3 + a [27]*z2 + a [28]*z + a [29];
  303.  
  304.    Result->z = a [2]*x3 + x2*(a [5]*y + 2*a [7]*z + a [8]) +
  305.            x*(a [11]*y2 + y*(2*a [13]*z + a [14]) + 3*a [16]*z2 + 2*a [17]*z + a [18]) +
  306.            a [21]*y3 + y2*(2*a [23]*z + a [24]) + y*(3*a [26]*z2 + 2*a [27]*z + a [28]) +
  307.            4*a [30]*z3 + 3*a [31]*z2 + 2*a [32]*z + a [33];
  308.  
  309.    Len = Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
  310.    Len = sqrt(Len);
  311.    if (Len == 0.0) {
  312.       /* The normal is not defined at this point of the surface.  Set it
  313.          to any arbitrary direction. */
  314.       Result->x = 1.0;
  315.       Result->y = 0.0;
  316.       Result->z = 0.0;
  317.       }
  318.    else {
  319.       Result->x /= Len;        /* normalize the normal */
  320.       Result->y /= Len;
  321.       Result->z /= Len;
  322.       }
  323.    }
  324.  
  325. /* Make a copy of a quartic object */
  326. void *
  327. Copy_Quartic(Object)
  328.    OBJECT *Object;
  329. {
  330.    QUARTIC *New_Shape;
  331.  
  332.    New_Shape = Get_Quartic_Shape ();
  333.    *New_Shape = *((QUARTIC *) Object);
  334.    New_Shape -> Next_Object = NULL;
  335.  
  336.    if (New_Shape->Shape_Texture != NULL)
  337.       New_Shape->Shape_Texture = Copy_Texture (New_Shape->Shape_Texture);
  338.  
  339.    return (New_Shape);
  340. }
  341.  
  342. /* Move a quartic */
  343. void
  344. Translate_Quartic(Object, Vector)
  345.    OBJECT *Object;
  346.    VECTOR *Vector;
  347. {
  348.    TRANSFORMATION Transformation;
  349.    Get_Translation_Transformation(&Transformation, Vector);
  350.    Transform_Quartic((QUARTIC *)Object,
  351.              (MATRIX *)&(Transformation.inverse [0] [0]));
  352.    Translate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  353. }
  354.  
  355. /* Rotation of a quartic */
  356. void
  357. Rotate_Quartic(Object, Vector)
  358.    OBJECT *Object;
  359.    VECTOR *Vector;
  360. {
  361.    TRANSFORMATION Transformation;
  362.    Get_Rotation_Transformation(&Transformation, Vector);
  363.    Transform_Quartic((QUARTIC *)Object,
  364.              (MATRIX *)&(Transformation.inverse [0] [0]));
  365.    Rotate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  366. }
  367.  
  368. /* Resize a quartic */
  369. void
  370. Scale_Quartic(Object, Vector)
  371.    OBJECT *Object;
  372.    VECTOR *Vector;
  373. {
  374.    TRANSFORMATION Transformation;
  375.    Get_Scaling_Transformation (&Transformation, Vector);
  376.    Transform_Quartic((QUARTIC *)Object,
  377.              (MATRIX *)&(Transformation.inverse [0] [0]));
  378.    Scale_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  379. }
  380.  
  381. void
  382. Invert_Quartic(Object)
  383.    OBJECT *Object;
  384. {
  385.    QUARTIC *Shape = (QUARTIC *) Object;
  386.    int i;
  387.    for (i=0;i<35;i ++)
  388.       Shape->Coeffs [i] *= -1.0;
  389. }
  390.  
  391. /* Given a set of power values from a term in the general quartic, return
  392.    the index of the term. */
  393. int
  394. roll_term_indices(i, j, k)
  395.    int i, j, k;
  396. {
  397.    int l, t;
  398.    l = 4 - (i + j + k);
  399.    if (i < 0 || i > 4 || j < 0 || j > 4 ||
  400.        k < 0 || k > 4 || l < 0 || l > 4) {
  401.       PRINT("Error in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n", i, j, k, l);
  402.       return (0);
  403.       }
  404.    for (t=0;t<term_counts [4];t ++)
  405.       if (i == terms4 [t] [0] && j == terms4 [t] [1] &&
  406.           k == terms4 [t] [2] && l == terms4 [t] [3])
  407.      return t;
  408.    PRINT("Term not in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n", i, j, k, l);
  409.    return (0);
  410. }
  411.  
  412. /* Given the total power and the index, return the individual powers */
  413. void
  414. unroll_term_indices( power, index, i, j, k, l)
  415.    int power, index, *i, *j, *k, *l;
  416. {
  417.    if (power == 4) {
  418.       *i = terms4 [index] [0];
  419.       *j = terms4 [index] [1];
  420.       *k = terms4 [index] [2];
  421.       *l = terms4 [index] [3];
  422.       }
  423.    else if (power == 3) {
  424.       *i = terms3 [index] [0];
  425.       *j = terms3 [index] [1];
  426.       *k = terms3 [index] [2];
  427.       *l = terms3 [index] [3];
  428.       }
  429.    else if (power == 2) {
  430.       *i = terms2 [index] [0];
  431.       *j = terms2 [index] [1];
  432.       *k = terms2 [index] [2];
  433.       *l = terms2 [index] [3];
  434.       }
  435.    else if (power == 1) {
  436.       *i = terms1 [index] [0];
  437.       *j = terms1 [index] [1];
  438.       *k = terms1 [index] [2];
  439.       *l = terms1 [index] [3];
  440.       }
  441.    else {
  442.       *i = 0;
  443.       *j = 0;
  444.       *k = 0;
  445.       *l = 0;
  446.       }
  447. }
  448.  
  449. DBL
  450. do_partial_term(q, row, power, i, j, k, l)
  451.    MATRIX *q;
  452.    int row, power, i, j, k, l;
  453. {
  454.    DBL result;
  455.    result = ((DBL) factorials [power] /
  456.           (factorials [i]*factorials [j]*factorials [k]*factorials [l]));
  457.    if (i>0) result *= POW((*q) [0] [row], (DBL) i);
  458.    if (j>0) result *= POW((*q) [1] [row], (DBL) j);
  459.    if (k>0) result *= POW((*q) [2] [row], (DBL) k);
  460.    if (l>0) result *= POW((*q) [3] [row], (DBL) l);
  461.    return result;
  462. }
  463.  
  464. /* Using the transformation matrix q, transform the general quartic equation
  465.    given by a. */
  466. void
  467. Transform_Quartic(Shape, q)
  468.    QUARTIC *Shape;
  469.    MATRIX *q;
  470. {
  471.    int term_index, partial_index;
  472.    int ip, i, i0, i1, i2, i3;
  473.    int jp, j, j0, j1, j2, j3;
  474.    int kp, k, k0, k1, k2, k3;
  475.    DBL b [35], partial_term;
  476.    DBL tempx, tempy, tempz;
  477.  
  478.    /* Trim out any really low values from the transform. */
  479.    for (i=0;i<4;i ++) for (j=0;j<4;j ++)
  480.       if ((*q) [i] [j] > -EPSILON && (*q) [i] [j] < EPSILON) (*q) [i] [j] = 0.0;
  481.  
  482.    for (i=0;i<35;i ++) b [i] = 0.0;
  483.  
  484.    for (term_index=0;term_index<term_counts [4];term_index ++) {
  485.       if (Shape->Coeffs [term_index] == 0.0) continue;
  486.       ip = terms4 [term_index] [0];    /* Power of original x term */
  487.       jp = terms4 [term_index] [1];    /* Power of original y term */
  488.       kp = terms4 [term_index] [2];    /* Power of original z term */
  489.       /* Step through terms in: (q [0] [0]*x + q [0] [1]*y + q [0] [2]*z + q [0] [3])^i */
  490.       for (i=0;i<term_counts [ip];i ++) {
  491.      unroll_term_indices(ip, i, &i0, &i1, &i2, &i3);
  492.      tempx = do_partial_term(q, 0, ip, i0, i1, i2, i3);
  493.      if (tempx == 0.0) continue;
  494.  
  495.       /* Step through terms in: (q [1] [0]*x + q [1] [1]*y + q [1] [2]*z + q [1] [3])^j */
  496.      for (j=0;j<term_counts [jp];j ++) {
  497.         unroll_term_indices(jp, j, &j0, &j1, &j2, &j3);
  498.         tempy = do_partial_term(q, 1, jp, j0, j1, j2, j3);
  499.         if (tempy == 0.0) continue;
  500.  
  501.       /* Step through terms in: (q [2] [0]*x + q [2] [1]*y + q [2] [2]*z + q [2] [3])^k */
  502.         for (k=0;k<term_counts [kp];k ++) {
  503.            unroll_term_indices(kp, k, &k0, &k1, &k2, &k3);
  504.            tempz = do_partial_term(q, 2, kp, k0, k1, k2, k3);
  505.             if (tempz == 0.0) continue;
  506.  
  507.            /* Figure out it's index, and add it into the result */
  508.            partial_index = roll_term_indices(i0 + j0 + k0,i1 + j1 + k1,i2 + j2 + k2);
  509.            partial_term = Shape->Coeffs [term_index] * tempx * tempy * tempz;
  510.            b [partial_index]  += partial_term;
  511.            }
  512.         }
  513.      }
  514.       }
  515.    for (i=0;i<35;i ++) {
  516.       if (b [i] > -EPSILON && b [i] < EPSILON)
  517.      Shape->Coeffs [i] = 0.0;
  518.       else
  519.      Shape->Coeffs [i] = b [i];
  520.       }
  521. }
  522.