home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / quatern.c < prev    next >
C/C++ Source or Header  |  1999-10-20  |  12KB  |  557 lines
  1. /*****************************************************************************
  2. *                quatern.c
  3. *
  4. *  This module implements Quaternion algebra julia fractals.
  5. *
  6. *  This file was written by Pascal Massimino.
  7. *  Revised and updated for POV-Ray 3.x by Tim Wegner
  8. *
  9. *  from Persistence of Vision(tm) Ray Tracer
  10. *  Copyright 1996,1999 Persistence of Vision Team
  11. *---------------------------------------------------------------------------
  12. *  NOTICE: This source code file is provided so that users may experiment
  13. *  with enhancements to POV-Ray and to port the software to platforms other
  14. *  than those supported by the POV-Ray Team.  There are strict rules under
  15. *  which you are permitted to use this file.  The rules are in the file
  16. *  named POVLEGAL.DOC which should be distributed with this file.
  17. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  18. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  19. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  20. *
  21. * This program is based on the popular DKB raytracer version 2.12.
  22. * DKBTrace was originally written by David K. Buck.
  23. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  24. *
  25. *****************************************************************************/
  26.  
  27. #include "frame.h"
  28. #include "povray.h"
  29. #include "vector.h"
  30. #include "povproto.h"
  31. #include "fractal.h"
  32. #include "quatern.h"
  33. #include "spheres.h"
  34.  
  35.  
  36.  
  37. /*****************************************************************************
  38. * Local preprocessor defines
  39. ******************************************************************************/
  40.  
  41. #ifndef Fractal_Tolerance
  42. #define Fractal_Tolerance 1e-8
  43. #endif
  44.  
  45. #define Deriv_z2(n1,n2,n3,n4)               \
  46. {                                           \
  47.   tmp = (n1)*x - (n2)*y - (n3)*z - (n4)*w;  \
  48.   (n2) = (n1)*y + x*(n2);                   \
  49.   (n3) = (n1)*z + x*(n3);                   \
  50.   (n4) = (n1)*w + x*(n4);                   \
  51.   (n1) = tmp;                               \
  52. }
  53.  
  54. #define Deriv_z3(n1,n2,n3,n4)              \
  55. {                                          \
  56.   dtmp = 2.0*((n2)*y + (n3)*z + (n4)*w);   \
  57.   dtmp2 = 6.0*x*(n1) - dtmp;               \
  58.   (n1) = ( (n1)*x3 - x*dtmp )*3.0;         \
  59.   (n2) = (n2)*x4 + y*dtmp2;                \
  60.   (n3) = (n3)*x4 + z*dtmp2;                \
  61.   (n4) = (n4)*x4 + w*dtmp2;                \
  62. }
  63.  
  64.  
  65. /*****************************************************************************
  66. * Local typedefs
  67. ******************************************************************************/
  68.  
  69.  
  70.  
  71. /*****************************************************************************
  72. * Static functions
  73. ******************************************************************************/
  74.  
  75.  
  76.  
  77. /*****************************************************************************
  78. * Local variables
  79. ******************************************************************************/
  80.  
  81.  
  82.  
  83. /*****************************************************************************
  84. *
  85. * FUNCTION
  86. *
  87. * INPUT
  88. *
  89. * OUTPUT
  90. *
  91. * RETURNS
  92. *
  93. * AUTHOR
  94. *
  95. *   Pascal Massimino
  96. *
  97. * DESCRIPTION
  98. *
  99. *   -
  100. *
  101. * CHANGES
  102. *
  103. *   Dec 1994 : Creation.
  104. *
  105. ******************************************************************************/
  106.  
  107. int Iteration_z3(VECTOR point, FRACTAL *Julia)
  108. {
  109.   int i;
  110.   DBL x, y, z, w;
  111.   DBL d, x2, tmp;
  112.   DBL Exit_Value;
  113.  
  114.   Sx[0] = x = point[X];
  115.   Sy[0] = y = point[Y];
  116.   Sz[0] = z = point[Z];
  117.   Sw[0] = w = (Julia->SliceDist  
  118.                   - Julia->Slice[X]*x 
  119.                   - Julia->Slice[Y]*y 
  120.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  121.   
  122.   Exit_Value = Julia->Exit_Value;
  123.  
  124.   for (i = 1; i <= Julia->n; ++i)
  125.   {
  126.     d = y * y + z * z + w * w;
  127.  
  128.     x2 = x * x;
  129.  
  130.     if ((d + x2) > Exit_Value)
  131.     {
  132.       return (FALSE);
  133.     }
  134.  
  135.     tmp = 3.0 * x2 - d;
  136.  
  137.     Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
  138.     Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
  139.     Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
  140.     Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
  141.   }
  142.  
  143.   return (TRUE);
  144. }
  145.  
  146.  
  147.  
  148. /*****************************************************************************
  149. *
  150. * FUNCTION
  151. *
  152. * INPUT
  153. *
  154. * OUTPUT
  155. *
  156. * RETURNS
  157. *
  158. * AUTHOR
  159. *
  160. *   Pascal Massimino
  161. *
  162. * DESCRIPTION
  163. *
  164. *   -
  165. *
  166. * CHANGES
  167. *
  168. *   Dec 1994 : Creation.
  169. *
  170. ******************************************************************************/
  171.  
  172. int Iteration_Julia(VECTOR point, FRACTAL *Julia)
  173. {
  174.   int i;
  175.   DBL x, y, z, w;
  176.   DBL d, x2;
  177.   DBL Exit_Value;
  178.  
  179.   Sx[0] = x = point[X];
  180.   Sy[0] = y = point[Y];
  181.   Sz[0] = z = point[Z];
  182.   Sw[0] = w = (Julia->SliceDist  
  183.                   - Julia->Slice[X]*x 
  184.                   - Julia->Slice[Y]*y 
  185.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  186.  
  187.   Exit_Value = Julia->Exit_Value;
  188.  
  189.   for (i = 1; i <= Julia->n; ++i)
  190.   {
  191.     d = y * y + z * z + w * w;
  192.  
  193.     x2 = x * x;
  194.  
  195.     if ((d + x2) > Exit_Value)
  196.     {
  197.       return (FALSE);
  198.     }
  199.  
  200.     x *= 2.0;
  201.  
  202.     Sy[i] = y = x * y + Julia->Julia_Parm[Y];
  203.     Sz[i] = z = x * z + Julia->Julia_Parm[Z];
  204.     Sw[i] = w = x * w + Julia->Julia_Parm[T];
  205.     Sx[i] = x = x2 - d + Julia->Julia_Parm[X];;
  206.  
  207.   }
  208.  
  209.   return (TRUE);
  210. }
  211.  
  212.  
  213.  
  214. /*****************************************************************************
  215. *
  216. * FUNCTION
  217. *
  218. * INPUT
  219. *
  220. * OUTPUT
  221. *
  222. * RETURNS
  223. *
  224. * AUTHOR
  225. *
  226. *   Pascal Massimino
  227. *
  228. * DESCRIPTION
  229. *
  230. * D_Iteration puts in *Dist a lower bound for the distance from *point to the
  231. * set
  232. *
  233. * CHANGES
  234. *
  235. *   Dec 1994 : Creation.
  236. *
  237. ******************************************************************************/
  238.  
  239. /*----------- Distance estimator + iterations ------------*/
  240.  
  241. int D_Iteration_z3(VECTOR point, FRACTAL *Julia, DBL *Dist)
  242. {
  243.   int i, j;
  244.   DBL Norm, d;
  245.   DBL xx, yy, zz;
  246.   DBL x, y, z, w;
  247.   DBL tmp, x2;
  248.   DBL Exit_Value;
  249.   DBL Pow;
  250.  
  251.   x = Sx[0] = point[X];
  252.   y = Sy[0] = point[Y];
  253.   z = Sz[0] = point[Z];
  254.   w = Sw[0] = (Julia->SliceDist  
  255.                   - Julia->Slice[X]*x 
  256.                   - Julia->Slice[Y]*y 
  257.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  258.  
  259.   Exit_Value = Julia->Exit_Value;
  260.  
  261.   for (i = 1; i <= Julia->n; i++)
  262.   {
  263.     d = y * y + z * z + w * w;
  264.  
  265.     x2 = x * x;
  266.  
  267.     if ((Norm = d + x2) > Exit_Value)
  268.     {
  269.       /* Distance estimator */
  270.  
  271.       x = Sx[0];
  272.       y = Sy[0];
  273.       z = Sz[0];
  274.       w = Sw[0];
  275.  
  276.       Pow = 1.0 / 3.0;
  277.  
  278.       for (j = 1; j < i; ++j)
  279.       {
  280.         xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
  281.         yy = x * Sy[j] + y * Sx[j] - z * Sw[j] + w * Sz[j];
  282.         zz = x * Sz[j] + y * Sw[j] + z * Sx[j] - w * Sy[j];
  283.         w  = x * Sw[j] - y * Sz[j] + z * Sy[j] + w * Sx[j];
  284.  
  285.         x = xx;
  286.         y = yy;
  287.         z = zz;
  288.  
  289.         Pow /= 3.0;
  290.       }
  291.  
  292.       *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
  293.  
  294.       return (FALSE);
  295.     }
  296.  
  297.     tmp = 3.0 * x2 - d;
  298.  
  299.     Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
  300.     Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
  301.     Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
  302.     Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
  303.   }
  304.  
  305.   *Dist = Precision;
  306.  
  307.   return (TRUE);
  308. }
  309.  
  310.  
  311.  
  312. /*****************************************************************************
  313. *
  314. * FUNCTION
  315. *
  316. * INPUT
  317. *
  318. * OUTPUT
  319. *
  320. * RETURNS
  321. *
  322. * AUTHOR
  323. *
  324. *   Pascal Massimino
  325. *
  326. * DESCRIPTION
  327. *
  328. *   -
  329. *
  330. * CHANGES
  331. *
  332. *   Dec 1994 : Creation.
  333. *
  334. ******************************************************************************/
  335.  
  336. int D_Iteration_Julia(VECTOR point, FRACTAL *Julia, DBL *Dist)
  337. {
  338.   int i, j;
  339.   DBL Norm, d;
  340.   DBL Exit_Value;
  341.   DBL x, y, z, w;
  342.   DBL xx, yy, zz, x2;
  343.   DBL Pow;
  344.  
  345.   x = Sx[0] = point[X];
  346.   y = Sy[0] = point[Y];
  347.   z = Sz[0] = point[Z];
  348.   w = Sw[0] = (Julia->SliceDist  
  349.                   - Julia->Slice[X]*x 
  350.                   - Julia->Slice[Y]*y 
  351.                   - Julia->Slice[Z]*z)/Julia->Slice[T]; 
  352.  
  353.   Exit_Value = Julia->Exit_Value;
  354.  
  355.   for (i = 1; i <= Julia->n; i++)
  356.   {
  357.     d = y * y + z * z + w * w;
  358.  
  359.     x2 = x * x;
  360.  
  361.     if ((Norm = d + x2) > Exit_Value)
  362.     {
  363.       /* Distance estimator */
  364.  
  365.       x = Sx[0];
  366.       y = Sy[0];
  367.       z = Sz[0];
  368.       w = Sw[0];
  369.  
  370.       Pow = 1.0 / 2.0;
  371.  
  372.       for (j = 1; j < i; ++j)
  373.       {
  374.         xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
  375.         yy = x * Sy[j] + y * Sx[j] + w * Sz[j] - z * Sw[j];
  376.         zz = x * Sz[j] + z * Sx[j] + y * Sw[j] - w * Sy[j];
  377.         w  = x * Sw[j] + w * Sx[j] + z * Sy[j] - y * Sz[j];
  378.  
  379.         x = xx;
  380.         y = yy;
  381.         z = zz;
  382.  
  383.         Pow /= 2.0;
  384.       }
  385.  
  386.       *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
  387.  
  388.       return (FALSE);
  389.     }
  390.  
  391.     x *= 2.0;
  392.  
  393.     Sy[i] = y = x * y + Julia->Julia_Parm[Y];
  394.     Sz[i] = z = x * z + Julia->Julia_Parm[Z];
  395.     Sw[i] = w = x * w + Julia->Julia_Parm[T];
  396.     Sx[i] = x = x2 - d + Julia->Julia_Parm[X];
  397.  
  398.   }
  399.  
  400.   *Dist = Precision;
  401.  
  402.   return (TRUE);
  403. }
  404.  
  405. /*****************************************************************************
  406. *
  407. * FUNCTION
  408. *
  409. * INPUT
  410. *   
  411. * OUTPUT
  412. *   
  413. * RETURNS
  414. *
  415. * AUTHOR
  416. *
  417. *   Pascal Massimino
  418. *
  419. * DESCRIPTION
  420. *
  421. * Provided the iterations sequence has been built, perform the computation of
  422. * the Normal
  423. *
  424. * CHANGES
  425. *
  426. *   Dec 1994 : Creation.
  427. *
  428. ******************************************************************************/
  429.  
  430. void Normal_Calc_z3(VECTOR Result, int N_Max, FRACTAL *fractal)
  431. {
  432.   DBL
  433.   n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
  434.   n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
  435.   n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
  436.  
  437.   DBL x, y, z, w;
  438.   int i;
  439.   DBL tmp, dtmp, dtmp2, x2, x3, x4;
  440.  
  441.   x = Sx[0];
  442.   y = Sy[0];
  443.   z = Sz[0];
  444.   w = Sw[0];
  445.  
  446.   for (i = 1; i <= N_Max; i++)
  447.   {
  448.     tmp = y * y + z * z + w * w;
  449.  
  450.     x2 = x * x;
  451.     x3 = x2 - tmp;
  452.     x4 = 3.0 * x2 - tmp;
  453.  
  454.     Deriv_z3(n11, n12, n13, n14);
  455.     Deriv_z3(n21, n22, n23, n24);
  456.     Deriv_z3(n31, n32, n33, n34);
  457.  
  458.     x = Sx[i];
  459.     y = Sy[i];
  460.     z = Sz[i];
  461.     w = Sw[i];
  462.   }
  463.  
  464.   Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
  465.   Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
  466.   Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
  467. }
  468.  
  469.  
  470.  
  471. /*****************************************************************************
  472. *
  473. * FUNCTION
  474. *
  475. * INPUT
  476. *
  477. * OUTPUT
  478. *
  479. * RETURNS
  480. *
  481. * AUTHOR
  482. *
  483. *   Pascal Massimino
  484. *
  485. * DESCRIPTION
  486. *
  487. *   -
  488. *
  489. * CHANGES
  490. *
  491. *   Dec 1994 : Creation.
  492. *
  493. ******************************************************************************/
  494.  
  495. void Normal_Calc_Julia(VECTOR Result, int N_Max, FRACTAL *fractal)
  496. {
  497.   DBL
  498.   n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
  499.   n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
  500.   n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
  501.   DBL tmp;
  502.   DBL x, y, z, w;
  503.   int i;
  504.  
  505.   x = Sx[0];
  506.   y = Sy[0];
  507.   z = Sz[0];
  508.   w = Sw[0];
  509.  
  510.   for (i = 1; i <= N_Max; i++)
  511.   {
  512.     Deriv_z2(n11, n12, n13, n14);
  513.     Deriv_z2(n21, n22, n23, n24);
  514.     Deriv_z2(n31, n32, n33, n34);
  515.  
  516.     x = Sx[i];
  517.     y = Sy[i];
  518.     z = Sz[i];
  519.     w = Sw[i];
  520.   }
  521.  
  522.   Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
  523.   Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
  524.   Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
  525. }
  526.  
  527.  
  528.  
  529. /*****************************************************************************
  530. *
  531. * FUNCTION
  532. *
  533. * INPUT
  534. *
  535. * OUTPUT
  536. *   
  537. * RETURNS
  538. *   
  539. * AUTHOR
  540. *
  541. *   Pascal Massimino
  542. *   
  543. * DESCRIPTION
  544. *
  545. *   -
  546. *
  547. * CHANGES
  548. *
  549. *   Dec 1994 : Creation.
  550. *
  551. ******************************************************************************/
  552.  
  553. int F_Bound_Julia(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL *Depth_Max)
  554. {
  555.   return (Intersect_Sphere(Ray, Fractal->Center, Fractal->Radius_Squared, Depth_Min, Depth_Max));
  556. }
  557.