home *** CD-ROM | disk | FTP | other *** search
/ back2roots/padua / padua.7z / padua / gfx / dkb212sr.lzh / quartics.c < prev    next >
C/C++ Source or Header  |  1991-05-03  |  17KB  |  546 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+a[ 4]*xx2*yy2+
  170.            a[ 5]*xx2*yy*zz+a[ 7]*xx2*zz2+a[10]*xx*yy3+
  171.        a[11]*xx*yy2*zz+a[13]*xx*yy*zz2+a[16]*xx*zz3+
  172.        a[20]*yy4+a[21]*yy3*zz+a[23]*yy2*zz2+a[26]*yy*zz3+a[30]*zz4);
  173.  
  174.    t[1] = (4*a[ 0]*x*xx3+
  175.            a[ 1]*(3*x*xx2*yy+xx3*y)+
  176.        a[ 2]*(3*x*xx2*zz+xx3*z)+
  177.        a[ 3]*xx3+
  178.        a[ 4]*(2*x*xx*yy2+2*xx2*y*yy)+
  179.        a[ 5]*(xx2*(y*zz+yy*z)+2*x*xx*yy*zz)+
  180.        a[ 6]*xx2*yy+
  181.        a[ 7]*(2*x*xx*zz2+2*xx2*z*zz)+
  182.        a[ 8]*xx2*zz+
  183.        a[10]*(x*yy3+3*xx*y*yy2)+
  184.        a[11]*(xx*(2*y*yy*zz+yy2*z)+x*yy2*zz)+
  185.        a[12]*xx*yy2+
  186.        a[13]*(xx*(y*zz2+2*yy*z*zz)+x*yy*zz2)+
  187.        a[14]*xx*yy*zz+
  188.        a[16]*(x*zz3+3*xx*z*zz2)+
  189.        a[17]*xx*zz2+
  190.        4*a[20]*y*yy3+
  191.        a[21]*(3*y*yy2*zz+yy3*z)+
  192.        a[22]*yy3+
  193.        zz*(2*a[23]*yy*(y*zz+yy*z)+
  194.        a[24]*yy2+
  195.        zz*(a[26]*(y*zz+3*yy*z)+
  196.        a[27]*yy+zz*(4*a[30]*z+a[31]))));
  197.  
  198.    t[2] = (6*a[ 0]*x2*xx2+
  199.            3*a[ 1]*x*xx*(x*yy+xx*y)+
  200.            3*a[ 2]*x*xx*(x*zz+xx*z)+
  201.        3*a[ 3]*x*xx2+
  202.        a[ 4]*(x2*yy2+4*x*xx*y*yy+xx2*y2)+
  203.        a[ 5]*(x2*yy*zz+2*x*xx*(y*zz+yy*z)+xx2*y*z)+
  204.        a[ 6]*(2*x*xx*yy+xx2*y)+
  205.        a[ 7]*(x2*zz2+4*x*xx*z*zz+xx2*z2)+
  206.        a[ 8]*(2*x*xx*zz+xx2*z)+
  207.        a[ 9]*xx2+
  208.        a[10]*(3*x*y*yy2+3*xx*y2*yy)+
  209.        a[11]*(x*yy*(2*y*zz+yy*z)+xx*(y2*zz+2*y*yy*z))+
  210.        a[12]*(x*yy2+2*xx*y*yy)+
  211.        a[13]*(x*zz*(y*zz+2*yy*z)+xx*(2*y*z*zz+yy*z2))+
  212.        a[14]*(x*yy*zz+xx*(y*zz+yy*z))+
  213.        a[15]*xx*yy+
  214.        a[16]*(3*x*z*zz2+3*xx*z2*zz)+
  215.        a[17]*(x*zz2+2*xx*z*zz)+
  216.        a[18]*xx*zz+
  217.        6*a[20]*y2*yy2+
  218.        3*a[21]*y*yy*(y*zz+yy*z)+
  219.        3*a[22]*y*yy2+
  220.        a[23]*(y2*zz2+4*y*yy*z*zz+yy2*z2)+
  221.        a[24]*(2*y*yy*zz+yy2*z)+
  222.        a[25]*yy2+
  223.        zz*(3*a[26]*z*(y*zz+yy*z)+
  224.        a[27]*(y*zz+2*yy*z)+
  225.        a[28]*yy+
  226.        6*a[30]*z2*zz+
  227.        3*a[31]*z*zz+
  228.        a[32]*zz));
  229.  
  230.    t[3] = (4*a[ 0]*x3*xx+
  231.            a[ 1]*x2*(x*yy+3*xx*y)+
  232.        a[ 2]*x2*(x*zz+3*xx*z)+
  233.        3*a[ 3]*x2*xx+
  234.        2*a[ 4]*x*y*(x*yy+xx*y)+
  235.        a[ 5]*x*(x*(y*zz+yy*z)+2*xx*y*z)+
  236.        a[ 6]*x*(x*yy+2*xx*y)+
  237.        2*a[ 7]*x*z*(x*zz+xx*z)+
  238.        a[ 8]*x*(x*zz+2*xx*z)+
  239.        2*a[ 9]*x*xx+
  240.        a[10]*(3*x*y2*yy-xx*y3)+
  241.        a[11]*(x*y*(y*zz+2*yy*z)+xx*y2*z)+
  242.        a[12]*(2*x*y*yy+xx*y2)+
  243.        a[13]*(x*z*(2*y*zz+yy*z)+xx*y*z2)+
  244.        a[14]*(x*(y*zz+yy*z)+xx*y*z)+
  245.        a[15]*(x*yy+xx*y)+
  246.        a[16]*(3*x*z2*zz+xx*z3)+
  247.        a[17]*(2*x*z*zz+xx*z2)+
  248.        a[18]*(x*zz+xx*z)+
  249.        a[19]*xx+
  250.        4*a[20]*y3*yy+
  251.        a[21]*y2*(y*zz+3*yy*z)+
  252.        3*a[22]*y2*yy+
  253.        2*a[23]*y*z*(y*zz+yy*z)+
  254.        a[24]*y*(y*zz+2*yy*z)+
  255.        2*a[25]*y*yy+
  256.        a[26]*(3*y*z2*zz+yy*z3)+
  257.        a[27]*(2*y*z*zz+yy*z2)+
  258.        a[28]*(y*zz+yy*z)+
  259.        a[29]*yy+
  260.        zz*(4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33]));
  261.  
  262.    t[4] = a[ 0]*x4+a[ 1]*x3*y+a[ 2]*x3*z+a[ 3]*x3+a[ 4]*x2*y2+
  263.       a[ 5]*x2*y*z+a[ 6]*x2*y+a[ 7]*x2*z2+a[ 8]*x2*z+a[ 9]*x2+
  264.       a[10]*x*y3+a[11]*x*y2*z+a[12]*x*y2+a[13]*x*y*z2+a[14]*x*y*z+
  265.       a[15]*x*y+a[16]*x*z3+a[17]*x*z2+a[18]*x*z+a[19]*x+a[20]*y4+
  266.       a[21]*y3*z+a[22]*y3+a[23]*y2*z2+a[24]*y2*z+a[25]*y2+
  267.       a[26]*y*z3+a[27]*y*z2+a[28]*y*z+a[29]*y+a[30]*z4+a[31]*z3+
  268.       a[32]*z2+a[33]*z+a[34];
  269.    return solve_quartic(t, Depths);
  270. }
  271.  
  272. int
  273. Inside_Quartic (point, Object)
  274.    VECTOR *point;
  275.    OBJECT *Object;
  276. {
  277.    QUARTIC *Shape = (QUARTIC *) Object;
  278.    DBL Result,*a,x,y,z,x2,y2,z2,x3,y3,z3,x4,y4,z4;
  279.  
  280.    x = point->x; y = point->y; z = point->z;
  281.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  282.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  283.    x4 = x * x3; y4 = y * y3; z4 = z * z3;
  284.    a = Shape->Coeffs;
  285.  
  286.    Result = a[ 0]*x4+a[ 1]*x3*y+a[ 2]*x3*z+a[ 3]*x3+a[ 4]*x2*y2+
  287.         a[ 5]*x2*y*z+a[ 6]*x2*y+a[ 7]*x2*z2+a[ 8]*x2*z+a[ 9]*x2+
  288.         a[10]*x*y3+a[11]*x*y2*z+a[12]*x*y2+a[13]*x*y*z2+a[14]*x*y*z+
  289.         a[15]*x*y+a[16]*x*z3+a[17]*x*z2+a[18]*x*z+a[19]*x+a[20]*y4+
  290.         a[21]*y3*z+a[22]*y3+a[23]*y2*z2+a[24]*y2*z+a[25]*y2+
  291.         a[26]*y*z3+a[27]*y*z2+a[28]*y*z+a[29]*y+a[30]*z4+a[31]*z3+
  292.         a[32]*z2+a[33]*z+a[34];
  293.  
  294.    if (Result < Small_Tolerance)
  295.       return (TRUE);
  296.    else
  297.       return (FALSE);
  298. }
  299.  
  300. /* Normal to a quartic */
  301. void
  302. Quartic_Normal(Result, Object, Intersection_Point)
  303.    VECTOR *Result, *Intersection_Point;
  304.    OBJECT *Object;
  305. {
  306.    QUARTIC *Shape = (QUARTIC *) Object;
  307.    DBL *a,x,y,z,x2,y2,z2,x3,y3,z3,Len;
  308.  
  309.    a = Shape->Coeffs;
  310.    x = Intersection_Point->x;
  311.    y = Intersection_Point->y;
  312.    z = Intersection_Point->z;
  313.    x2 = x * x;  y2 = y * y;  z2 = z * z;
  314.    x3 = x * x2; y3 = y * y2; z3 = z * z2;
  315.  
  316.    Result->x = 4*a[ 0]*x3+3*x2*(a[ 1]*y+a[ 2]*z+a[ 3])+
  317.            2*x*(a[ 4]*y2+y*(a[ 5]*z+a[ 6])+a[ 7]*z2+a[ 8]*z+a[ 9])+
  318.            a[10]*y3+y2*(a[11]*z+a[12])+y*(a[13]*z2+a[14]*z+a[15])+
  319.            a[16]*z3+a[17]*z2+a[18]*z+a[19];
  320.  
  321.    Result->y = a[ 1]*x3+x2*(2*a[ 4]*y+a[ 5]*z+a[ 6])+
  322.            x*(3*a[10]*y2+2*y*(a[11]*z+a[12])+a[13]*z2+a[14]*z+a[15])+
  323.            4*a[20]*y3+3*y2*(a[21]*z+a[22])+2*y*(a[23]*z2+a[24]*z+a[25])+
  324.            a[26]*z3+a[27]*z2+a[28]*z+a[29];
  325.  
  326.    Result->z = a[ 2]*x3+x2*(a[ 5]*y+2*a[ 7]*z+a[ 8])+
  327.            x*(a[11]*y2+y*(2*a[13]*z+a[14])+3*a[16]*z2+2*a[17]*z+a[18])+
  328.            a[21]*y3+y2*(2*a[23]*z+a[24])+y*(3*a[26]*z2+2*a[27]*z+a[28])+
  329.            4*a[30]*z3+3*a[31]*z2+2*a[32]*z+a[33];
  330.  
  331.    Len = Result->x * Result->x + Result->y * Result->y + Result->z * Result->z;
  332.    Len = sqrt(Len);
  333.    if (Len == 0.0) {
  334.       /* The normal is not defined at this point of the surface.  Set it
  335.          to any arbitrary direction. */
  336.       Result->x = 1.0;
  337.       Result->y = 0.0;
  338.       Result->z = 0.0;
  339.       }
  340.    else {
  341.       Result->x /= Len;        /* normalize the normal */
  342.       Result->y /= Len;
  343.       Result->z /= Len;
  344.       }
  345.    }
  346.  
  347. /* Make a copy of a quartic object */
  348. void *
  349. Copy_Quartic(Object)
  350.    OBJECT *Object;
  351. {
  352.    QUARTIC *New_Shape;
  353.  
  354.    New_Shape = Get_Quartic_Shape ();
  355.    *New_Shape = *((QUARTIC *) Object);
  356.    New_Shape -> Next_Object = NULL;
  357.  
  358.    if (New_Shape->Shape_Texture != NULL)
  359.       New_Shape->Shape_Texture = Copy_Texture (New_Shape->Shape_Texture);
  360.  
  361.    return (New_Shape);
  362. }
  363.  
  364. /* Move a quartic */
  365. void
  366. Translate_Quartic(Object, Vector)
  367.    OBJECT *Object;
  368.    VECTOR *Vector;
  369. {
  370.    TRANSFORMATION Transformation;
  371.    Get_Translation_Transformation(&Transformation, Vector);
  372.    Transform_Quartic((QUARTIC *)Object,
  373.              (MATRIX *)&(Transformation.inverse[0][0]));
  374.    Translate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  375. }
  376.  
  377. /* Rotation of a quartic */
  378. void
  379. Rotate_Quartic(Object, Vector)
  380.    OBJECT *Object;
  381.    VECTOR *Vector;
  382. {
  383.    TRANSFORMATION Transformation;
  384.    Get_Rotation_Transformation(&Transformation, Vector);
  385.    Transform_Quartic((QUARTIC *)Object,
  386.              (MATRIX *)&(Transformation.inverse[0][0]));
  387.    Rotate_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  388. }
  389.  
  390. /* Resize a quartic */
  391. void
  392. Scale_Quartic(Object, Vector)
  393.    OBJECT *Object;
  394.    VECTOR *Vector;
  395. {
  396.    TRANSFORMATION Transformation;
  397.    Get_Scaling_Transformation (&Transformation, Vector);
  398.    Transform_Quartic((QUARTIC *)Object,
  399.              (MATRIX *)&(Transformation.inverse[0][0]));
  400.    Scale_Texture (&((QUARTIC *) Object)->Shape_Texture, Vector);
  401. }
  402.  
  403. void
  404. Invert_Quartic(Object)
  405.    OBJECT *Object;
  406. {
  407.    QUARTIC *Shape = (QUARTIC *) Object;
  408.    int i;
  409.    for (i=0;i<35;i++)
  410.       Shape->Coeffs[i] *= -1.0;
  411. }
  412.  
  413. /* Given a set of power values from a term in the general quartic, return
  414.    the index of the term. */
  415. int
  416. roll_term_indices(i, j, k)
  417.    int i, j, k;
  418. {
  419.    int l, t;
  420.    l = 4 - (i + j + k);
  421.    if (i < 0 || i > 4 || j < 0 || j > 4 ||
  422.        k < 0 || k > 4 || l < 0 || l > 4) {
  423.       printf("Error in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n",
  424.          i, j, k, l);
  425.       exit(1);
  426.       }
  427.    for (t=0;t<term_counts[4];t++)
  428.       if (i == terms4[t][0] && j == terms4[t][1] &&
  429.           k == terms4[t][2] && l == terms4[t][3])
  430.      return t;
  431.    printf("Term not  in 'roll_term_indices': i = %d, j = %d, k = %d, l = %d\n",
  432.          i, j, k, l);
  433.    exit(1);
  434. }
  435.  
  436. /* Given the total power and the index, return the individual powers */
  437. void
  438. unroll_term_indices( power, index, i, j, k, l)
  439.    int power, index, *i, *j, *k, *l;
  440. {
  441.    if (power == 4) {
  442.       *i = terms4[index][0];
  443.       *j = terms4[index][1];
  444.       *k = terms4[index][2];
  445.       *l = terms4[index][3];
  446.       }
  447.    else if (power == 3) {
  448.       *i = terms3[index][0];
  449.       *j = terms3[index][1];
  450.       *k = terms3[index][2];
  451.       *l = terms3[index][3];
  452.       }
  453.    else if (power == 2) {
  454.       *i = terms2[index][0];
  455.       *j = terms2[index][1];
  456.       *k = terms2[index][2];
  457.       *l = terms2[index][3];
  458.       }
  459.    else if (power == 1) {
  460.       *i = terms1[index][0];
  461.       *j = terms1[index][1];
  462.       *k = terms1[index][2];
  463.       *l = terms1[index][3];
  464.       }
  465.    else {
  466.       *i = 0;
  467.       *j = 0;
  468.       *k = 0;
  469.       *l = 0;
  470.       }
  471. }
  472.  
  473. DBL
  474. do_partial_term(q, row, power, i, j, k, l)
  475.    MATRIX *q;
  476.    int row, power, i, j, k, l;
  477. {
  478.    DBL result;
  479.    result = ((DBL) factorials[power] /
  480.           (factorials[i]*factorials[j]*factorials[k]*factorials[l]));
  481.    if (i>0) result *= POW((*q)[0][row], (DBL) i);
  482.    if (j>0) result *= POW((*q)[1][row], (DBL) j);
  483.    if (k>0) result *= POW((*q)[2][row], (DBL) k);
  484.    if (l>0) result *= POW((*q)[3][row], (DBL) l);
  485.    return result;
  486. }
  487.  
  488. /* Using the transformation matrix q, transform the general quartic equation
  489.    given by a. */
  490. void
  491. Transform_Quartic(Shape, q)
  492.    QUARTIC *Shape;
  493.    MATRIX *q;
  494. {
  495.    int term_index, partial_index;
  496.    int ip, i, i0, i1, i2, i3;
  497.    int jp, j, j0, j1, j2, j3;
  498.    int kp, k, k0, k1, k2, k3;
  499.    DBL b[35], partial_term;
  500.    DBL tempx, tempy, tempz;
  501.  
  502.    /* Trim out any really low values from the transform. */
  503.    for (i=0;i<4;i++) for (j=0;j<4;j++)
  504.       if ((*q)[i][j] > -EPSILON && (*q)[i][j] < EPSILON) (*q)[i][j] = 0.0;
  505.  
  506.    for (i=0;i<35;i++) b[i] = 0.0;
  507.  
  508.    for (term_index=0;term_index<term_counts[4];term_index++) {
  509.       if (Shape->Coeffs[term_index] == 0.0) continue;
  510.       ip = terms4[term_index][0];    /* Power of original x term */
  511.       jp = terms4[term_index][1];    /* Power of original y term */
  512.       kp = terms4[term_index][2];    /* Power of original z term */
  513.       /* Step through terms in: (q[0][0]*x+q[0][1]*y+q[0][2]*z+q[0][3])^i */
  514.       for (i=0;i<term_counts[ip];i++) {
  515.      unroll_term_indices(ip, i, &i0, &i1, &i2, &i3);
  516.      tempx = do_partial_term(q, 0, ip, i0, i1, i2, i3);
  517.      if (tempx == 0.0) continue;
  518.  
  519.       /* Step through terms in: (q[1][0]*x+q[1][1]*y+q[1][2]*z+q[1][3])^j */
  520.      for (j=0;j<term_counts[jp];j++) {
  521.         unroll_term_indices(jp, j, &j0, &j1, &j2, &j3);
  522.         tempy = do_partial_term(q, 1, jp, j0, j1, j2, j3);
  523.         if (tempy == 0.0) continue;
  524.  
  525.       /* Step through terms in: (q[2][0]*x+q[2][1]*y+q[2][2]*z+q[2][3])^k */
  526.         for (k=0;k<term_counts[kp];k++) {
  527.            unroll_term_indices(kp, k, &k0, &k1, &k2, &k3);
  528.            tempz = do_partial_term(q, 2, kp, k0, k1, k2, k3);
  529.             if (tempz == 0.0) continue;
  530.  
  531.            /* Figure out it's index, and add it into the result */
  532.            partial_index = roll_term_indices(i0+j0+k0,i1+j1+k1,i2+j2+k2);
  533.            partial_term = Shape->Coeffs[term_index] * tempx * tempy * tempz;
  534.            b[partial_index] += partial_term;
  535.            }
  536.         }
  537.      }
  538.       }
  539.    for (i=0;i<35;i++) {
  540.       if (b[i] > -EPSILON && b[i] < EPSILON)
  541.      Shape->Coeffs[i] = 0.0;
  542.       else
  543.      Shape->Coeffs[i] = b[i];
  544.       }
  545. }
  546.