home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / Pattern.c < prev    next >
C/C++ Source or Header  |  2001-01-27  |  112KB  |  4,788 lines
  1. /**************************************************************************
  2. *                pattern.c
  3. *
  4. *  -------------------------------------------------------
  5. *  ATTENTION:
  6. *  This is an unofficial version of pattern.c modified by
  7. *  Ryoichi Suzuki, rsuzuki@etl.go.jp, for use  with 
  8. *  "isosurface" shape type.
  9. *  -------------------------------------------------------
  10. *  This module implements texturing functions that return a value to be
  11. *  used in a pigment or normal.
  12. *
  13. *  from Persistence of Vision(tm) Ray Tracer
  14. *  Copyright 1996,1999 Persistence of Vision Team
  15. *---------------------------------------------------------------------------
  16. *  NOTICE: This source code file is provided so that users may experiment
  17. *  with enhancements to POV-Ray and to port the software to platforms other
  18. *  than those supported by the POV-Ray Team.  There are strict rules under
  19. *  which you are permitted to use this file.  The rules are in the file
  20. *  named POVLEGAL.DOC which should be distributed with this file.
  21. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  22. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  23. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  24. *
  25. * This program is based on the popular DKB raytracer version 2.12.
  26. * DKBTrace was originally written by David K. Buck.
  27. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  28. *
  29. * Modifications by Hans-Detlev Fink, January 1999, used with permission
  30. * Modifications by Thomas Willhalm, March 1999, used with permission
  31. *
  32. *****************************************************************************/
  33.  
  34. /*
  35.  * Some texture ideas garnered from SIGGRAPH '85 Volume 19 Number 3,
  36.  * "An Image Synthesizer" By Ken Perlin.
  37.  * Further Ideas Garnered from "The RenderMan Companion" (Addison Wesley).
  38.  */
  39.  
  40. #include "frame.h"
  41. #include "vector.h"
  42. #include "povproto.h"
  43. #include "matrices.h"
  44. #include "pattern.h"
  45. #include "povray.h"
  46. #include "texture.h"
  47. #include "image.h"
  48. #include "txttest.h"
  49. #include "colour.h"
  50. /** poviso: Jul.14 '96 R.S. **/
  51. #ifdef POVISO
  52. #include "isosrf.h"
  53. #endif
  54. /** --- **/
  55. #include "objects.h"
  56. #ifdef BlobPatternPatch
  57. #include "pigment.h"
  58. #include "blob.h"
  59. #endif
  60. #ifdef ProximityPatch
  61. #include "ray.h"
  62. #endif
  63. #ifdef SplineWavePatch
  64. #include "splines.h"
  65. #endif
  66. /*****************************************************************************
  67. * Local preprocessor defines
  68. ******************************************************************************/
  69.  
  70.  
  71. /*****************************************************************************
  72. * Static functions
  73. ******************************************************************************/
  74.  
  75. #ifdef CellsPatch
  76. static DBL cells (VECTOR EPoint);
  77. #endif
  78. #ifdef VanSicklePatternPatch
  79. static DBL blotches (VECTOR EPoint);
  80. static DBL banded (VECTOR EPoint, TPATTERN *TPat);
  81. static DBL sheet (VECTOR EPoint, TPATTERN *TPat);
  82. #endif
  83. static DBL agate (VECTOR EPoint, TPATTERN *TPat);
  84. static DBL brick (VECTOR EPoint, TPATTERN *TPat);
  85. static DBL checker (VECTOR EPoint);
  86. #ifdef TrianglulairSquarePatch 
  87. static DBL square (VECTOR EPoint);
  88. static DBL ternaire (VECTOR EPoint);
  89. #endif
  90.  
  91. #ifdef CracklePatch
  92. static DBL crackle (VECTOR EPoint, TPATTERN *TPat);
  93. #else
  94. static DBL crackle (VECTOR EPoint);
  95. #endif
  96. static DBL gradient (VECTOR EPoint, TPATTERN *TPat);
  97. static DBL granite (VECTOR EPoint);
  98. static DBL leopard (VECTOR EPoint);
  99. static DBL magnet1m (VECTOR EPoint, TPATTERN *TPat);
  100. static DBL magnet1j (VECTOR EPoint, TPATTERN *TPat);
  101. static DBL magnet2m (VECTOR EPoint, TPATTERN *TPat);
  102. static DBL magnet2j (VECTOR EPoint, TPATTERN *TPat);
  103. static DBL mandel3 (VECTOR EPoint, TPATTERN *TPat);
  104. static DBL mandel4 (VECTOR EPoint, TPATTERN *TPat);
  105. static DBL julia (VECTOR EPoint, TPATTERN *TPat);
  106. static DBL julia3 (VECTOR EPoint, TPATTERN *TPat);
  107. static DBL julia4 (VECTOR EPoint, TPATTERN *TPat);
  108. static DBL marble (VECTOR EPoint, TPATTERN *TPat);
  109. static DBL onion (VECTOR EPoint);
  110. static DBL radial (VECTOR EPoint);
  111. #ifdef PolaricalPatch
  112. static DBL polarical (VECTOR EPoint);
  113. #endif
  114. static DBL spiral1 (VECTOR EPoint, TPATTERN *TPat);
  115. static DBL spiral2 (VECTOR EPoint, TPATTERN *TPat);
  116. static DBL wood (VECTOR EPoint, TPATTERN *TPat);
  117. static DBL hexagon (VECTOR EPoint);
  118. static DBL planar_pattern (VECTOR EPoint);
  119. static DBL spherical (VECTOR EPoint);
  120. static DBL boxed (VECTOR EPoint);
  121. static DBL cylindrical (VECTOR EPoint);
  122. static DBL density_file (VECTOR EPoint, TPATTERN *TPat);
  123. #ifdef SolidPatternPatch
  124. static DBL SolidPat(VECTOR EPoint, TPATTERN *TPat);  /*Chris Huff 7/20/99 solid pattern*/
  125. #endif
  126. #ifdef ClothPatternPatch
  127.   static DBL ClothPat(VECTOR EPoint); /*Chris Huff cloth pattern*/
  128.   static DBL Cloth2Pat(VECTOR EPoint); /*Chris Huff cloth2 pattern*/
  129. #endif
  130. #ifdef TorodialPatch
  131.   static DBL ToroidalSpiral(VECTOR EPoint, TPATTERN *TPat); /*Chris Huff torodil pattern */
  132. #endif
  133. /*YS*/ /*moved this to pattern.h*/
  134. /* removed static */
  135. /*long PickInCube (VECTOR tv, VECTOR p1);*/
  136. /*YS*/
  137.  
  138. static DBL ripples_pigm (VECTOR EPoint, TPATTERN *TPat);
  139. static DBL waves_pigm  (VECTOR EPoint, TPATTERN *TPat);
  140. static DBL dents_pigm  (VECTOR EPoint);
  141. static DBL wrinkles_pigm (VECTOR EPoint);
  142. static DBL quilted_pigm (VECTOR EPoint, TPATTERN *TPat);
  143. static TURB *Search_For_Turb (WARP *Warps);
  144. /* static TURB *Copy_Turb (TURB *Old);   Unused function [AED] */
  145. static unsigned short readushort(FILE *infile);
  146. #ifdef ProximityPatch
  147.   static DBL proximity(VECTOR EPoint, TPATTERN *TPat);/*Chris Huff proximity pattern*/
  148. #endif
  149. #ifdef ObjectPatternPatch
  150. static DBL object(VECTOR EPoint, TPATTERN *TPat);/*Chris Huff proximity pattern*/
  151. #endif
  152. #ifdef BlobPatternPatch
  153. static DBL blob_pattern (VECTOR EPoint, TPATTERN *TPat);/*Chris Huff-blob pattern*/
  154. static DBL blob_pigment(VECTOR EPoint, TPATTERN *TPat, INTERSECTION *Intersection);/*Chris Huff-blob pigment*/
  155. #endif
  156. /** poviso: Jul. 14, 96 R.S. **/
  157. #ifdef POVISO
  158. static DBL func_pat (VECTOR EPoint, TPATTERN *TPat);
  159. #endif
  160. /** --- **/
  161. /* -hdf Apr 98 */
  162. static DBL slope (VECTOR EPoint, TPATTERN *TPat, INTERSECTION *Intersection);
  163.  
  164. #ifdef PigmentPatternPatch
  165. static DBL pigment_pattern (VECTOR EPoint, TPATTERN *TPat, INTERSECTION *Inter);
  166. #endif
  167.  
  168. /*****************************************************************************
  169. *
  170. * FUNCTION
  171. *
  172. *   agate
  173. *
  174. * INPUT
  175. *
  176. *   EPoint -- The point in 3d space at which the pattern is evaluated.
  177. *   TPat   -- Texture pattern struct
  178. *
  179. * OUTPUT
  180. *
  181. * RETURNS
  182. *
  183. *   DBL value in the range 0.0 to 1.0
  184. *
  185. * AUTHOR
  186. *
  187. *   POV-Ray Team
  188. *
  189. * DESCRIPTION
  190. *
  191. * CHANGES
  192. *   Oct 1994    : adapted from agate pigment by [CY]
  193. *
  194. ******************************************************************************/
  195.  
  196. static DBL agate (VECTOR EPoint, TPATTERN *TPat)
  197. {
  198.   register DBL noise, turb_val;
  199.   TURB* Turb;
  200.  
  201.   Turb=Search_For_Turb(TPat->Warps);
  202.  
  203.   turb_val = TPat->Vals.Agate_Turb_Scale * Turbulence(EPoint,Turb);
  204.  
  205.   noise = 0.5 * (cycloidal(1.3 * turb_val + 1.1 * EPoint[Z]) + 1.0);
  206.  
  207.   if (noise < 0.0)
  208.   {
  209.     noise = 0.0;
  210.   }
  211.   else
  212.   {
  213.     noise = min(1.0, noise);
  214.     noise = pow(noise, 0.77);
  215.   }
  216.  
  217.   return(noise);
  218. }
  219.  
  220.  
  221. /*****************************************************************************
  222. *
  223. * FUNCTION
  224. *
  225. *   brick
  226. *
  227. * INPUT
  228. *
  229. *   EPoint -- The point in 3d space at which the pattern
  230. *   is evaluated.
  231. *   TPat   -- Texture pattern struct
  232. *   
  233. * OUTPUT
  234. *   
  235. * RETURNS
  236. *
  237. *   DBL value exactly 0.0 or 1.0
  238. *   
  239. * AUTHOR
  240. *
  241. *   Dan Farmer
  242. *   
  243. * DESCRIPTION
  244. *
  245. * CHANGES
  246. *   Oct 1994    : adapted from pigment by [CY]
  247. *
  248. ******************************************************************************/
  249.  
  250. static DBL brick (VECTOR EPoint, TPATTERN *TPat)
  251. {
  252.   int ibrickx, ibricky, ibrickz;
  253.   DBL brickheight, brickwidth, brickdepth;
  254.   DBL brickmortar, mortarheight, mortarwidth, mortardepth;
  255.   DBL brickx, bricky, brickz;
  256.   DBL x, y, z, fudgit;
  257.  
  258.   fudgit=Small_Tolerance+TPat->Vals.Brick.Mortar;
  259.  
  260.   x =  EPoint[X]+fudgit;
  261.   y =  EPoint[Y]+fudgit;
  262.   z =  EPoint[Z]+fudgit;
  263.  
  264.   brickwidth  = TPat->Vals.Brick.Size[X];
  265.   brickheight = TPat->Vals.Brick.Size[Y];
  266.   brickdepth  = TPat->Vals.Brick.Size[Z];
  267.   brickmortar = (DBL)TPat->Vals.Brick.Mortar;
  268.  
  269.   mortarwidth  = brickmortar / brickwidth;
  270.   mortarheight = brickmortar / brickheight;
  271.   mortardepth  = brickmortar / brickdepth;
  272.  
  273.   /* 1) Check mortar layers in the X-Z plane (ie: top view) */
  274.  
  275.   bricky = y / brickheight;
  276.   ibricky = (int) bricky;
  277.   bricky -= (DBL) ibricky;
  278.  
  279.   if (bricky < 0.0)
  280.   {
  281.     bricky += 1.0;
  282.   }
  283.  
  284.   if (bricky <= mortarheight)
  285.   {
  286.     return(0.0);
  287.   }
  288.  
  289.   bricky = (y / brickheight) * 0.5;
  290.   ibricky = (int) bricky;
  291.   bricky -= (DBL) ibricky;
  292.  
  293.   if (bricky < 0.0)
  294.   {
  295.     bricky += 1.0;
  296.   }
  297.  
  298.  
  299.   /* 2) Check ODD mortar layers in the Y-Z plane (ends) */
  300.  
  301.   brickx = (x / brickwidth);
  302.   ibrickx = (int) brickx;
  303.   brickx -= (DBL) ibrickx;
  304.  
  305.   if (brickx < 0.0)
  306.   {
  307.     brickx += 1.0;
  308.   }
  309.  
  310.   if ((brickx <= mortarwidth) && (bricky <= 0.5))
  311.   {
  312.     return(0.0);
  313.   }
  314.  
  315.   /* 3) Check EVEN mortar layers in the Y-Z plane (ends) */
  316.  
  317.   brickx = (x / brickwidth) + 0.5;
  318.   ibrickx = (int) brickx;
  319.   brickx -= (DBL) ibrickx;
  320.  
  321.   if (brickx < 0.0)
  322.   {
  323.     brickx += 1.0;
  324.   }
  325.  
  326.   if ((brickx <= mortarwidth) && (bricky > 0.5))
  327.   {
  328.     return(0.0);
  329.   }
  330.  
  331.   /* 4) Check ODD mortar layers in the Y-X plane (facing) */
  332.  
  333.   brickz = (z / brickdepth);
  334.   ibrickz = (int) brickz;
  335.   brickz -= (DBL) ibrickz;
  336.  
  337.   if (brickz < 0.0)
  338.   {
  339.     brickz += 1.0;
  340.   }
  341.  
  342.   if ((brickz <= mortardepth) && (bricky > 0.5))
  343.   {
  344.     return(0.0);
  345.   }
  346.  
  347.   /* 5) Check EVEN mortar layers in the X-Y plane (facing) */
  348.  
  349.   brickz = (z / brickdepth) + 0.5;
  350.   ibrickz = (int) brickz;
  351.   brickz -= (DBL) ibrickz;
  352.  
  353.   if (brickz < 0.0)
  354.   {
  355.     brickz += 1.0;
  356.   }
  357.  
  358.   if ((brickz <= mortardepth) && (bricky <= 0.5))
  359.   {
  360.     return(0.0);
  361.   }
  362.  
  363.   /* If we've gotten this far, color me brick. */
  364.  
  365.   return(1.0);
  366. }
  367.  
  368.  
  369. /*****************************************************************************
  370. *
  371. * FUNCTION
  372. *
  373. *   checker
  374. *
  375. * INPUT
  376. *
  377. *   EPoint -- The point in 3d space at which the pattern
  378. *   is evaluated.
  379. *
  380. * OUTPUT
  381. *
  382. * RETURNS
  383. *
  384. *   DBL value exactly 0.0 or 1.0
  385. *
  386. * AUTHOR
  387. *
  388. *   POV-Team
  389. *
  390. * DESCRIPTION
  391. *
  392. * CHANGES
  393. *   Oct 1994    : adapted from pigment by [CY]
  394. *
  395. ******************************************************************************/
  396.  
  397. static DBL checker (VECTOR EPoint)
  398. {
  399.   int value;
  400.  
  401.   value = (int)(floor(EPoint[X]+Small_Tolerance) +
  402.                 floor(EPoint[Y]+Small_Tolerance) +
  403.                 floor(EPoint[Z]+Small_Tolerance));
  404.  
  405.   if (value & 1)
  406.   {
  407.     return (1.0);
  408.   }
  409.   else
  410.   {
  411.     return (0.0);
  412.   }
  413. }
  414.  
  415.  
  416. #ifdef TrianglulairSquarePatch 
  417.  
  418. /*****************************************************************************
  419. *
  420. * FUNCTION
  421. *
  422. *   square 
  423. *
  424. * INPUT
  425. *
  426. *   EPoint -- The point in 3d space at which the pattern
  427. *   is evaluated.
  428. *
  429. * OUTPUT
  430. *
  431. * RETURNS
  432. *
  433. *   DBL value exactly 0.0, 1.0, 2.0 or 3.0
  434. *
  435. * AUTHOR
  436. *
  437. *   J. Grimbert
  438. *
  439. * DESCRIPTION
  440. *   Paving the XZ plan with 4 'colours', in square
  441. *
  442. * CHANGES
  443. *
  444. ******************************************************************************/
  445.  
  446. static DBL square (VECTOR EPoint)
  447. {
  448.   int valueX,valueZ;
  449.  
  450.   valueX = (int)(floor(EPoint[X])); 
  451.   valueZ = (int)(floor(EPoint[Z])); 
  452.  
  453.   if (valueX & 1)
  454.   {
  455.     if (valueZ & 1)
  456.     {
  457.        return (2.0);
  458.     }
  459.     else
  460.     {
  461.        return (3.0);
  462.     }
  463.   }
  464.   else
  465.   {
  466.     if (valueZ & 1)
  467.     {
  468.        return (1.0);
  469.     }
  470.     else
  471.     {
  472.        return (0.0);
  473.     }
  474.   }
  475. }
  476.  
  477. /*****************************************************************************
  478. *
  479. * FUNCTION
  480. *
  481. *   ternaire
  482. *
  483. * INPUT
  484. *
  485. *   EPoint -- The point in 3d space at which the pattern
  486. *   is evaluated.
  487. *
  488. * OUTPUT
  489. *
  490. * RETURNS
  491. *
  492. *   DBL value exactly 0.0, 1.0, 2.0, 3.0, 4.0 or 5.0
  493. *
  494. * AUTHOR
  495. *
  496. *   J. Grimbert
  497. *
  498. * DESCRIPTION
  499. *   Paving the XZ plan with 6 'colours', in triangle around the origin
  500. *
  501. * CHANGES
  502. *
  503. ******************************************************************************/
  504.  
  505. /* 1.73205080756887729352 is sqrt(3) */
  506. /* .86602540378443864676 is sqrt(3)/2 */
  507. #define SQR3_2 .86602540378443864676
  508. #define SQR3    1.73205080756887729352
  509. static DBL ternaire (VECTOR EPoint)
  510. {
  511. DBL answer;
  512. DBL x,z;
  513. DBL xs,zs;
  514. int a,b;
  515. DBL k,slop1,slop2;
  516. int mask;
  517.  
  518. x=EPoint[X];
  519. z=EPoint[Z];
  520.  
  521.  xs = x-3.0*floor(x/3.0);
  522.  zs = z-SQR3*floor(z/SQR3);
  523.  
  524.  /* xs,zs is in { [0.0, 3.0 [, [0.0, SQR3 [ } 
  525.  ** but there is some symetry to simplify the testing
  526.  */
  527.  
  528.  a = (int)floor(xs);
  529.  xs -= a;
  530.  b = (zs <SQR3_2 ? 0: 1);
  531.  if (b)
  532.  {
  533.    zs = SQR3 - zs; /* mirror */
  534.  } 
  535.  
  536.  k = 1.0 - xs;
  537.  if ((xs != 0.0)&&( k != 0.0 )) /* second condition should never occurs */
  538.  {
  539.     slop1 = zs/xs;
  540.     slop2 = zs/k; /* just in case */
  541.    switch( (slop1<SQR3?1:0)+(slop2<SQR3?2:0))
  542.     { 
  543.      case 3:
  544.       answer = 0.0;
  545.       break;
  546.      case 2:
  547.       answer = 1.0;
  548.        break;
  549.      case 1: 
  550.       answer = 3.0;
  551.       break;
  552.     }
  553.  }
  554.  else
  555.  {
  556.    answer = 1.0;
  557.  }
  558.  mask = (int) answer;
  559.  answer = (mask & 1) ? fmod(answer+2.0*a,6.0): fmod(6.0+answer-2.0*a,6.0);
  560.  if (b)
  561.  {
  562.    answer = 5.0 - answer;
  563.  }
  564.  
  565.  return answer;
  566. }
  567. #endif
  568.  
  569.  
  570. /*****************************************************************************
  571. *
  572. * FUNCTION
  573. *
  574. *   crackle
  575. *
  576. * INPUT
  577. *
  578. *   EPoint -- The point in 3d space at which the pattern
  579. *   is evaluated.
  580. * OUTPUT
  581. *
  582. * RETURNS
  583. *
  584. *   DBL value in the range 0.0 to 1.0
  585. *
  586. * AUTHOR
  587. *
  588. *   Jim McElhiney
  589. *
  590. * DESCRIPTION
  591. *
  592. *   "crackle":
  593. *
  594. *   New colour function by Jim McElhiney,
  595. *     CompuServe 71201,1326, aka mcelhiney@acm.org
  596. *
  597. *   Large scale, without turbulence, makes a pretty good stone wall.
  598. *   Small scale, without turbulence, makes a pretty good crackle ceramic glaze.
  599. *   Highly turbulent (with moderate displacement) makes a good marble, solving
  600. *   the problem of apparent parallel layers in Perlin's method.
  601. *   2 octaves of full-displacement turbulence make a great "drizzled paint"
  602. *   pattern, like a 1950's counter top.
  603. *   Rule of thumb:  put a single colour transition near 0 in your colour map.
  604. *
  605. *   Mathematically, the set crackle(p)=0 is a 3D Voronoi diagram of a field of
  606. *   semirandom points, and crackle(p)>0 is distance from set along shortest path.
  607. *   (A Voronoi diagram is the locus of points equidistant from their 2 nearest
  608. *   neighbours from a set of disjoint points, like the membranes in suds are
  609. *   to the centres of the bubbles).
  610. *
  611. *   All "crackle" specific source code and examples are in the public domain.
  612. *
  613. * CHANGES
  614. *   Oct 1994    : adapted from pigment by [CY]
  615. *
  616. ******************************************************************************/
  617. #ifdef CracklePatch
  618. static long IntPickInCube(int tvx, int tvy, int tvz, VECTOR  p1);
  619.  
  620. static DBL crackle (VECTOR EPoint, TPATTERN *TPat ) {
  621.   int    i;
  622.   long   thisseed;
  623.   DBL    sum, minsum, minsum2, minsum3, tf;
  624.   VECTOR minvec;
  625.   VECTOR tv, dv, t1;
  626.   int addx,addy,addz;
  627.  
  628.   VECTOR flo;
  629.   int cvc;
  630.   static int vali=0, vals[3];
  631.   static int valid[125];
  632.  
  633.   DBL Metric;
  634.   DBL Offset;
  635.  
  636.   int UseSquare;
  637.   int UseUnity;
  638.   int flox,floy,floz;
  639.   /*int seed,temp;*/
  640.  
  641.   Metric = TPat->Vals.Crackle.Metric[X];
  642.   Offset = TPat->Vals.Crackle.Offset;
  643.  
  644.   UseSquare = ( Metric == 2);
  645.   UseUnity = ( Metric == 1);
  646.  
  647.   Assign_Vector(tv,EPoint);
  648.  
  649.   /*
  650.    * Check to see if the input point is in the same unit cube as the last
  651.    * call to this function, to use cache of cubelets for speed.
  652.    */
  653.  
  654.   thisseed = PickInCube(tv, t1);
  655.  
  656.   if (thisseed != TPat->Vals.Crackle.lastseed)
  657.   {
  658.     /*
  659.      * No, not same unit cube.  Calculate the random points for this new
  660.      * cube and its 80 neighbours which differ in any axis by 1 or 2.
  661.      * Why distance of 2?  If there is 1 point in each cube, located
  662.      * randomly, it is possible for the closest random point to be in the
  663.      * cube 2 over, or the one two over and one up.  It is NOT possible
  664.      * for it to be two over and two up.  Picture a 3x3x3 cube with 9 more
  665.      * cubes glued onto each face.
  666.      */
  667.  
  668.  
  669.     flo[X] = floor(tv[X] - EPSILON);
  670.     flo[Y] = floor(tv[Y] - EPSILON);
  671.     flo[Z] = floor(tv[Z] - EPSILON);
  672.     
  673.     Assign_Vector( (TPat->Vals.Crackle.lastcenter), flo );
  674.  
  675.     /* Now store a points for this cube and each of the 80 neighbour cubes. */
  676.  
  677.     vals[0]=25*(2+(-2))+5*(2+(-1))+2+(-1);
  678.     vals[1]=25*(2+(-2))+5*(2+(-1))+2+(0);
  679.     vals[2]=25*(2+(-2))+5*(2+(-1))+2+(1);
  680.  
  681.     flox = (int)flo[X];
  682.     floy = (int)flo[Y];
  683.     floz = (int)flo[Z];
  684.  
  685.     for (addx = -2; addx <= 2; addx++)
  686.     {
  687.       for (addy = -2; addy <= 2; addy++)
  688.       {
  689.           for (addz = -2; addz <= 2; addz++)
  690.           {
  691.             /* For each cubelet in a 5x5 cube. */
  692.           cvc = 25*(2+addx)+5*(2+addy)+2+addz;
  693.  
  694.           if ((abs(addx)==2)+(abs(addy)==2)+(abs(addz)==2) <= 1)
  695.             {
  696.               /* Yes, it's within a 3d knight move away. */
  697.  
  698.               /*VAdd(sv, tv, add);*/
  699.             /*
  700.             sv[X] = tv[X] + (DBL)addx;
  701.             sv[Y] = tv[Y] + (DBL)addy;
  702.             sv[Z] = tv[Z] + (DBL)addz;
  703.               PickInCube(sv, t1);
  704.  
  705.               TPat->Vals.Crackle.cv[cvc][X] = t1[X];
  706.               TPat->Vals.Crackle.cv[cvc][Y] = t1[Y];
  707.               TPat->Vals.Crackle.cv[cvc][Z] = t1[Z];
  708.             */
  709. #define INLINE_PICK_IN_CUBE 0
  710. #if INLINE_PICK_IN_CUBE
  711.             /* do our own PickInCube and use as much integer math as possible */
  712. #ifdef NoiseTranslateFixPatch
  713.             seed = Hash3d((flox+addx)&0xFFF,(floy+addy)&0xFFF,(floz+addz)&0xFFF);
  714. #else
  715.             seed = Hash3d(flox+addx,floy+addy,floz+addz);
  716. #endif
  717.             temp = POV_GET_OLD_RAND(); /* save current seed */
  718.             POV_SRAND(seed);
  719.             TPat->Vals.Crackle.cv[cvc][X] = flox+addx + FRAND();
  720.             TPat->Vals.Crackle.cv[cvc][Y] = floy+addy + FRAND();
  721.             TPat->Vals.Crackle.cv[cvc][Z] = floz+addz + FRAND();
  722.             POV_SRAND(temp);  /* restore */
  723. #else
  724.               IntPickInCube(flox+addx,floy+addy,floz+addz, t1);
  725.  
  726.               TPat->Vals.Crackle.cv[cvc][X] = t1[X];
  727.               TPat->Vals.Crackle.cv[cvc][Y] = t1[Y];
  728.               TPat->Vals.Crackle.cv[cvc][Z] = t1[Z];
  729. #endif
  730.               valid[cvc]=1;
  731.             }
  732.             else {
  733.             valid[cvc]=0;
  734.             }
  735.           }
  736.       }
  737.     }
  738.  
  739.     TPat->Vals.Crackle.lastseed = thisseed;
  740.   }
  741.  
  742.   cvc=125;
  743.   /*
  744.    * Find the 2 points with the 2 shortest distances from the input point.
  745.    * Loop invariant:  minsum is shortest dist, minsum2 is 2nd shortest
  746.    */
  747.  
  748.   /* Set up the loop so the invariant is true:  minsum <= minsum2 */
  749.  
  750.   VSub(dv, TPat->Vals.Crackle.cv[vals[0]], tv);  
  751.   if ( UseSquare ) {
  752.       minsum  = VSumSqr(dv);
  753.       if ( Offset ) minsum += Offset*Offset;
  754.   }
  755.   else if ( UseUnity ) {
  756.       minsum = dv[X] + dv[Y] + dv[Z];
  757.       if ( Offset ) minsum += Offset;
  758.   }
  759.   else {
  760.       minsum = pow( fabs( dv[X] ), Metric ) +
  761.                pow( fabs( dv[Y] ), Metric ) +
  762.                pow( fabs( dv[Z] ), Metric );
  763.       if ( Offset ) minsum += pow( Offset, Metric );
  764.   }
  765.   Assign_Vector( &minvec, TPat->Vals.Crackle.cv+vals[0] );
  766.   VSub(dv, TPat->Vals.Crackle.cv[vals[1]], tv);  
  767.   if ( UseSquare ) {
  768.       minsum2  = VSumSqr(dv);
  769.       if ( Offset ) minsum2 += Offset*Offset;
  770.   }
  771.   else if ( UseUnity ) {
  772.       minsum2 = dv[X] + dv[Y] + dv[Z];
  773.       if ( Offset ) minsum2 += Offset;
  774.   }
  775.   else {
  776.       minsum2 = pow( fabs( dv[X] ), Metric ) +
  777.                pow( fabs( dv[Y] ), Metric ) +
  778.                pow( fabs( dv[Z] ), Metric );
  779.       if ( Offset ) minsum2 += pow( Offset, Metric );
  780.   }
  781.   VSub(dv, TPat->Vals.Crackle.cv[vals[2]], tv);  
  782.   if ( UseSquare ) {
  783.       minsum3  = VSumSqr(dv);
  784.       if ( Offset ) minsum3 += Offset*Offset;
  785.   }
  786.   else if ( UseUnity ) {
  787.       minsum3 = dv[X] + dv[Y] + dv[Z];
  788.       if ( Offset ) minsum3 += Offset;
  789.   }
  790.   else {
  791.       minsum3 = pow( fabs( dv[X] ), Metric ) +
  792.                pow( fabs( dv[Y] ), Metric ) +
  793.                pow( fabs( dv[Z] ), Metric );
  794.       if ( Offset ) minsum3 += pow( Offset, Metric );
  795.   }
  796.  
  797.   if (minsum2 < minsum)
  798.   {
  799.     tf = minsum; minsum = minsum2; minsum2 = tf;
  800.     Assign_Vector( &minvec, TPat->Vals.Crackle.cv+vals[1] );
  801.   }
  802.   if (minsum3 < minsum)
  803.   {
  804.     tf = minsum; minsum = minsum3; minsum3 = tf;
  805.     Assign_Vector( &minvec, TPat->Vals.Crackle.cv+vals[2] );
  806.   }
  807.   if ( minsum3 < minsum2 ) {
  808.       tf = minsum2; minsum2=minsum3; minsum3= tf;
  809.   }
  810.  
  811.   /* Loop for the 81 cubelets to find closest and 2nd closest. */
  812.  
  813.   for (i = vals[2]+1; i < cvc; i++) if (valid[i])
  814.   {
  815.     VSub(dv, TPat->Vals.Crackle.cv[i], tv);
  816.  
  817.     if ( UseSquare ) {
  818.         sum  = VSumSqr(dv);
  819.         if ( Offset ) sum += Offset*Offset;
  820.     }
  821.     else if ( UseUnity ) {
  822.         sum = dv[X] + dv[Y] + dv[Z];
  823.         if ( Offset ) sum += Offset;
  824.     }
  825.     else {
  826.         sum = pow( fabs( dv[X] ), Metric ) +
  827.               pow( fabs( dv[Y] ), Metric ) +
  828.               pow( fabs( dv[Z] ), Metric );
  829.         if ( Offset ) sum += pow( Offset, Metric );
  830.     }
  831.  
  832.     if (sum < minsum)
  833.     {
  834.       minsum3 = minsum2;
  835.       minsum2 = minsum;
  836.       minsum = sum;
  837.       Assign_Vector( &minvec, TPat->Vals.Crackle.cv+i );
  838.     }
  839.     else if (sum < minsum2) 
  840.     {
  841.         minsum3 = minsum2;
  842.     minsum2 = sum;
  843.     }
  844.     else if ( sum < minsum3 ) {
  845.         minsum3 = sum;
  846.     }
  847.   }
  848.  
  849.   if ( TPat->Vals.Crackle.IsSolid ) {
  850.       tf = Noise( minvec );
  851.   }
  852.   else if (UseSquare) {
  853.       tf = TPat->Vals.Crackle.Form[X]*sqrt(minsum) + 
  854.            TPat->Vals.Crackle.Form[Y]*sqrt(minsum2) +
  855.            TPat->Vals.Crackle.Form[Z]*sqrt(minsum3); 
  856.   }
  857.   else if ( UseUnity ) {
  858.       tf = TPat->Vals.Crackle.Form[X]*minsum + 
  859.            TPat->Vals.Crackle.Form[Y]*minsum2 +
  860.            TPat->Vals.Crackle.Form[Z]*minsum3; 
  861.   }
  862.   else {
  863.       tf = TPat->Vals.Crackle.Form[X]*pow(minsum, 1.0/Metric) + 
  864.            TPat->Vals.Crackle.Form[Y]*pow(minsum2, 1.0/Metric) + 
  865.            TPat->Vals.Crackle.Form[Y]*pow(minsum3, 1.0/Metric); 
  866.   }
  867.  
  868.  
  869.   return max(min(tf, 1.), 0.);
  870. }
  871.  
  872. /*****************************************************************************
  873. *
  874. * FUNCTION
  875. *
  876. *   IntPickInCube(tvx,tvy,tvz, p1)
  877. *    a version of PickInCube that takes integers for input
  878. *
  879. * INPUT
  880. *
  881. *   ?
  882. *
  883. * OUTPUT
  884. *   
  885. * RETURNS
  886. *
  887. *   long integer hash function used, to speed up cacheing.
  888. *   
  889. * AUTHOR
  890. *
  891. *   original PickInCube by Jim McElhiney
  892. *   this integer one modified by Nathan Kopp
  893. *   
  894. * DESCRIPTION
  895. *
  896. *   A subroutine to go with crackle.
  897. *
  898. *   Pick a random point in the same unit-sized cube as tv, in a
  899. *   predictable way, so that when called again with another point in
  900. *   the same unit cube, p1 is picked to be the same.
  901. *
  902. * CHANGES
  903. *
  904. ******************************************************************************/
  905.  
  906. static long IntPickInCube(int tvx, int tvy, int tvz, VECTOR  p1)
  907. {
  908.   int seed, temp;
  909.  
  910. #ifdef NoiseTranslateFixPatch
  911.   seed = Hash3d(tvx&0xFFF,tvy&0xFFF,tvz&0xFFF);
  912. #else
  913.   seed = Hash3d(tvx,tvy,tvz);
  914. #endif
  915.   temp = POV_GET_OLD_RAND(); /* save current seed */
  916.   POV_SRAND(seed);
  917.  
  918.   p1[X] = tvx + FRAND();
  919.   p1[Y] = tvy + FRAND();
  920.   p1[Z] = tvz + FRAND();
  921.  
  922.   POV_SRAND(temp);  /* restore */
  923.  
  924.   return((long)seed);
  925. }
  926.  
  927.  
  928.  
  929.  
  930. #else
  931. static DBL crackle (VECTOR EPoint)
  932. {
  933.   int    i;
  934.   long   thisseed;
  935.   DBL    sum, minsum, minsum2, tf;
  936.   VECTOR sv, tv, dv, t1, add;
  937.  
  938.   static int cvc;
  939.   static long lastseed = 0x80000000;
  940.   static VECTOR cv[81];
  941.  
  942.   Assign_Vector(tv,EPoint);
  943.  
  944.   /*
  945.    * Check to see if the input point is in the same unit cube as the last
  946.    * call to this function, to use cache of cubelets for speed.
  947.    */
  948.  
  949.   thisseed = PickInCube(tv, t1);
  950.  
  951.   if (thisseed != lastseed)
  952.   {
  953.     /*
  954.      * No, not same unit cube.  Calculate the random points for this new
  955.      * cube and its 80 neighbours which differ in any axis by 1 or 2.
  956.      * Why distance of 2?  If there is 1 point in each cube, located
  957.      * randomly, it is possible for the closest random point to be in the
  958.      * cube 2 over, or the one two over and one up.  It is NOT possible
  959.      * for it to be two over and two up.  Picture a 3x3x3 cube with 9 more
  960.      * cubes glued onto each face.
  961.      */
  962.  
  963.     /* Now store a points for this cube and each of the 80 neighbour cubes. */
  964.  
  965.     cvc = 0;
  966.  
  967.     for (add[X] = -2.0; add[X] < 2.5; add[X] +=1.0)
  968.     {
  969.       for (add[Y] = -2.0; add[Y] < 2.5; add[Y] += 1.0)
  970.       {
  971.         for (add[Z] = -2.0; add[Z] < 2.5; add[Z] += 1.0)
  972.         {
  973.           /* For each cubelet in a 5x5 cube. */
  974.  
  975.           if ((fabs(add[X])>1.5)+(fabs(add[Y])>1.5)+(fabs(add[Z])>1.5) <= 1.0)
  976.           {
  977.             /* Yes, it's within a 3d knight move away. */
  978.  
  979.             VAdd(sv, tv, add);
  980.  
  981.             PickInCube(sv, t1);
  982.  
  983.             cv[cvc][X] = t1[X];
  984.             cv[cvc][Y] = t1[Y];
  985.             cv[cvc][Z] = t1[Z];
  986.             cvc++;
  987.           }
  988.         }
  989.       }
  990.     }
  991.  
  992.     lastseed = thisseed;
  993.   }
  994.  
  995.   /*
  996.    * Find the 2 points with the 2 shortest distances from the input point.
  997.    * Loop invariant:  minsum is shortest dist, minsum2 is 2nd shortest
  998.    */
  999.  
  1000.   /* Set up the loop so the invariant is true:  minsum <= minsum2 */
  1001.  
  1002.   VSub(dv, cv[0], tv);  minsum  = VSumSqr(dv);
  1003.   VSub(dv, cv[1], tv);  minsum2 = VSumSqr(dv);
  1004.  
  1005.   if (minsum2 < minsum)
  1006.   {
  1007.     tf = minsum; minsum = minsum2; minsum2 = tf;
  1008.   }
  1009.  
  1010.   /* Loop for the 81 cubelets to find closest and 2nd closest. */
  1011.  
  1012.   for (i = 2; i < cvc; i++)
  1013.   {
  1014.     VSub(dv, cv[i], tv);
  1015.  
  1016.     sum = VSumSqr(dv);
  1017.  
  1018.     if (sum < minsum)
  1019.     {
  1020.       minsum2 = minsum;
  1021.       minsum = sum;
  1022.     }
  1023.     else
  1024.     {
  1025.       if (sum < minsum2)
  1026.       {
  1027.         minsum2 = sum;
  1028.       }
  1029.     }
  1030.   }
  1031.  
  1032.   /* Crackle value is absolute value of diff in dist to closest 2 points. */
  1033.  
  1034.   tf = sqrt(minsum2) - sqrt(minsum);      /* minsum is known <= minsum2 */
  1035.  
  1036.   /*
  1037.    * Note that the theoretical range of this function is 0 to root 3.
  1038.    * In practice, it rarely exceeds 0.9, and only very rarely 1.0
  1039.    */
  1040.  
  1041.   return min(tf, 1.);
  1042. }
  1043. #endif
  1044.  
  1045. /*****************************************************************************
  1046. *
  1047. * FUNCTION
  1048. *
  1049. *   gradient
  1050. *
  1051. * INPUT
  1052. *
  1053. *   EPoint -- The point in 3d space at which the pattern
  1054. *   is evaluated.
  1055. *   
  1056. * OUTPUT
  1057. *   
  1058. * RETURNS
  1059. *
  1060. *   DBL value in the range 0.0 to 1.0
  1061. *   
  1062. * AUTHOR
  1063. *
  1064. *   POV-Ray Team
  1065. *   
  1066. * DESCRIPTION
  1067. *
  1068. *   Gradient Pattern - gradient based on the fractional values of
  1069. *   x, y or z, based on whether or not the given directional vector is
  1070. *   a 1.0 or a 0.0.
  1071. *   The basic concept of this is from DBW Render, but Dave Wecker's
  1072. *   only supports simple Y axis gradients.
  1073. *
  1074. * CHANGES
  1075. *   Oct 1994    : adapted from pigment by [CY]
  1076. *
  1077. ******************************************************************************/
  1078.  
  1079. static DBL gradient (VECTOR EPoint, TPATTERN *TPat)
  1080. {
  1081.   register int i;
  1082.   register DBL temp;
  1083.   DBL value = 0.0;
  1084.  
  1085.   for (i=X; i<=Z; i++)
  1086.   {
  1087.     if (TPat->Vals.Gradient[i] != 0.0)
  1088.     {
  1089.       temp = fabs(EPoint[i]);
  1090.  
  1091.       value += fmod(temp,1.0);
  1092.     }
  1093.   }
  1094.  
  1095.   /* Clamp to 1.0. */
  1096.  
  1097.   value = ((value > 1.0) ? fmod(value, 1.0) : value);
  1098.  
  1099.   return(value);
  1100. }
  1101.  
  1102.  
  1103.  
  1104. /*****************************************************************************
  1105. *
  1106. * FUNCTION
  1107. *
  1108. *   granite
  1109. *
  1110. * INPUT
  1111. *
  1112. *   EPoint -- The point in 3d space at which the pattern
  1113. *   is evaluated.
  1114. *   
  1115. * OUTPUT
  1116. *   
  1117. * RETURNS
  1118. *
  1119. *   DBL value in the range 0.0 to 1.0
  1120. *   
  1121. * AUTHOR
  1122. *
  1123. *   POV-Ray Team
  1124. *   
  1125. * DESCRIPTION
  1126. *
  1127. *   Granite - kind of a union of the "spotted" and the "dented" textures,
  1128. *   using a 1/f fractal noise function for color values. Typically used
  1129. *   with small scaling values. Should work with colour maps for pink granite.
  1130. *
  1131. * CHANGES
  1132. *   Oct 1994    : adapted from pigment by [CY]
  1133. *
  1134. ******************************************************************************/
  1135.  
  1136. static DBL granite (VECTOR EPoint)
  1137. {
  1138.   register int i;
  1139.   register DBL temp, noise = 0.0, freq = 1.0;
  1140.   VECTOR tv1,tv2;
  1141.  
  1142.   VScale(tv1,EPoint,4.0);
  1143.  
  1144.   for (i = 0; i < 6 ; freq *= 2.0, i++)
  1145.   {
  1146.     VScale(tv2,tv1,freq);
  1147.     temp = 0.5 - Noise (tv2);
  1148.  
  1149.     temp = fabs(temp);
  1150.  
  1151.     noise += temp / freq;
  1152.   }
  1153.  
  1154.   return(noise);
  1155. }
  1156.  
  1157.  
  1158.  
  1159. /*****************************************************************************
  1160. *
  1161. * FUNCTION
  1162. *
  1163. *   leopard
  1164. *
  1165. * INPUT
  1166. *
  1167. *   EPoint -- The point in 3d space at which the pattern
  1168. *   is evaluated.
  1169. *
  1170. * OUTPUT
  1171. *
  1172. * RETURNS
  1173. *
  1174. *   DBL value in the range 0.0 to 1.0
  1175. *
  1176. * AUTHOR
  1177. *
  1178. *   Scott Taylor
  1179. *
  1180. * DESCRIPTION
  1181. *
  1182. * CHANGES
  1183. *   Jul 1991 : Creation.
  1184. *   Oct 1994 : adapted from pigment by [CY]
  1185. *
  1186. ******************************************************************************/
  1187.  
  1188. static DBL leopard (VECTOR EPoint)
  1189. {
  1190.   register DBL value, temp1, temp2, temp3;
  1191.  
  1192.   /* This form didn't work with Zortech 386 compiler */
  1193.   /* value = Sqr((sin(x)+sin(y)+sin(z))/3); */
  1194.   /* So we break it down. */
  1195.  
  1196.   temp1 = sin(EPoint[X]);
  1197.   temp2 = sin(EPoint[Y]);
  1198.   temp3 = sin(EPoint[Z]);
  1199.  
  1200.   value = Sqr((temp1 + temp2 + temp3) / 3.0);
  1201.  
  1202.   return(value);
  1203. }
  1204.  
  1205.  
  1206.  
  1207. /*****************************************************************************
  1208. *
  1209. * FUNCTION
  1210. *
  1211. *   mandel
  1212. *
  1213. * INPUT
  1214. *
  1215. *   EPoint -- The point in 3d space at which the pattern
  1216. *   is evaluated.
  1217. *
  1218. * OUTPUT
  1219. *
  1220. * RETURNS
  1221. *
  1222. *   DBL value in the range 0.0 to 1.0
  1223. *
  1224. * AUTHOR
  1225. *
  1226. *   submitted by user, name lost (sorry)
  1227. *
  1228. * DESCRIPTION
  1229. * The mandel pattern computes the standard Mandelbrot fractal pattern and
  1230. * projects it onto the X-Y plane.  It uses the X and Y coordinates to compute
  1231. * the Mandelbrot set.
  1232. *
  1233. * CHANGES
  1234. *   Oct 1994 : adapted from pigment by [CY]
  1235. *
  1236. ******************************************************************************/
  1237.  
  1238. /* NK fractal - various changes */
  1239.  
  1240. static DBL fractal_exterior_color(TPATTERN *TPat, int iters,
  1241.                                   DBL a, DBL b)
  1242. {
  1243.     switch(TPat->Vals.Fractal.exterior_type)
  1244.     {
  1245.       case 0:
  1246.           return  (DBL)TPat->Vals.Fractal.efactor;
  1247.       case 1:
  1248.           return (DBL)iters / (DBL)TPat->Vals.Fractal.Iterations;
  1249.       case 2:
  1250.           return a * (DBL)TPat->Vals.Fractal.efactor;
  1251.       case 3:
  1252.           return b * (DBL)TPat->Vals.Fractal.efactor;
  1253.       case 4:
  1254.           return a*a * (DBL)TPat->Vals.Fractal.efactor;
  1255.       case 5:
  1256.           return b*b * (DBL)TPat->Vals.Fractal.efactor;
  1257.       case 6:
  1258.           return sqrt(a*a+b*b) * (DBL)TPat->Vals.Fractal.efactor;
  1259.     }
  1260.     return 0;
  1261. }
  1262.  
  1263. static DBL fractal_interior_color(TPATTERN *TPat, int iters,
  1264.                                   DBL a, DBL b, DBL mindist2)
  1265. {
  1266.     switch(TPat->Vals.Fractal.interior_type)
  1267.     {
  1268.       case 0:
  1269.           return  (DBL)TPat->Vals.Fractal.ifactor;
  1270.       case 1:
  1271.           return sqrt(mindist2) * (DBL)TPat->Vals.Fractal.ifactor;
  1272.       case 2:
  1273.           return a * (DBL)TPat->Vals.Fractal.ifactor;
  1274.       case 3:
  1275.           return b * (DBL)TPat->Vals.Fractal.ifactor;
  1276.       case 4:
  1277.           return a*a * (DBL)TPat->Vals.Fractal.ifactor;
  1278.       case 5:
  1279.           return b*b * (DBL)TPat->Vals.Fractal.ifactor;
  1280.       case 6:
  1281.           return a*a+b*b * (DBL)TPat->Vals.Fractal.ifactor;
  1282.     }
  1283.     return 0;
  1284. }
  1285.  
  1286. static DBL mandel (VECTOR EPoint, TPATTERN *TPat)
  1287. {
  1288.   int it_max, col;
  1289.   DBL a, b, cf, a2, b2, x, y, dist2, mindist2;
  1290.  
  1291.   a = x = EPoint[X]; a2 = Sqr(a);
  1292.   b = y = EPoint[Y]; b2 = Sqr(b);
  1293.   mindist2 = a2+b2;
  1294.  
  1295.   it_max = TPat->Vals.Fractal.Iterations;
  1296.  
  1297.   for (col = 0; col < it_max; col++)
  1298.   {
  1299.     b  = 2.0 * a * b + y;
  1300.     a  = a2 - b2 + x;
  1301.  
  1302.     a2 = Sqr(a);
  1303.     b2 = Sqr(b);
  1304.     dist2 = a2+b2;
  1305.  
  1306.     if(dist2 < mindist2) mindist2 = dist2;
  1307.     if(dist2 > 4.0)
  1308.     {
  1309.         cf = fractal_exterior_color(TPat, col, a, b);
  1310.         break;
  1311.     }
  1312.   }
  1313.  
  1314.   if(col == it_max)
  1315.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1316.  
  1317.   return(cf);
  1318. }
  1319.  
  1320. static DBL mandel3 (VECTOR EPoint, TPATTERN *TPat)
  1321. {
  1322.   int it_max, col;
  1323.   DBL a, b, cf, a2, b2, x, y, dist2, mindist2;
  1324.  
  1325.   a = x = EPoint[X]; a2 = Sqr(a);
  1326.   b = y = EPoint[Y]; b2 = Sqr(b);
  1327.   mindist2 = a2+b2;
  1328.  
  1329.   it_max = TPat->Vals.Fractal.Iterations;
  1330.  
  1331.   for (col = 0; col < it_max; col++)
  1332.   {
  1333.       b = 3.0*a2*b - b2*b + y;
  1334.       a = a2*a - 3.0*a*b2 + x;
  1335.  
  1336.       a2 = Sqr(a);
  1337.       b2 = Sqr(b);
  1338.       dist2 = a2+b2;
  1339.  
  1340.       if(dist2 < mindist2) mindist2 = dist2;
  1341.       if(dist2 > 4.0)
  1342.       {
  1343.           cf = fractal_exterior_color(TPat, col, a, b);
  1344.           break;
  1345.       }
  1346.   }
  1347.  
  1348.   if(col == it_max)
  1349.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1350.  
  1351.   return(cf);
  1352. }
  1353.  
  1354.  
  1355. static DBL mandel4 (VECTOR EPoint, TPATTERN *TPat)
  1356. {
  1357.   int it_max, col;
  1358.   DBL a, b, cf, a2, b2, x, y, dist2, mindist2;
  1359.  
  1360.   a = x = EPoint[X]; a2 = Sqr(a);
  1361.   b = y = EPoint[Y]; b2 = Sqr(b);
  1362.   mindist2 = a2+b2;
  1363.  
  1364.   it_max = TPat->Vals.Fractal.Iterations;
  1365.  
  1366.   for (col = 0; col < it_max; col++)
  1367.   {
  1368.       b = 4.0 * (a2*a*b - a*b2*b) + y;
  1369.       a = a2*a2 - 6.0*a2*b2 + b2*b2 + x;
  1370.  
  1371.       a2 = Sqr(a);
  1372.       b2 = Sqr(b);
  1373.       dist2 = a2+b2;
  1374.  
  1375.       if(dist2 < mindist2) mindist2 = dist2;
  1376.       if(dist2 > 4.0)
  1377.       {
  1378.           cf = fractal_exterior_color(TPat, col, a, b);
  1379.           break;
  1380.       }
  1381.   }
  1382.  
  1383.   if(col == it_max)
  1384.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1385.  
  1386.   return(cf);
  1387. }
  1388.  
  1389. static DBL julia (VECTOR EPoint, TPATTERN *TPat)
  1390. {
  1391.   int it_max, col;
  1392.   DBL a, b, cf, a2, b2, dist2, mindist2,
  1393.       cr = TPat->Vals.Fractal.Coord[U], ci = TPat->Vals.Fractal.Coord[V];
  1394.  
  1395.   a = EPoint[X]; a2 = Sqr(a);
  1396.   b = EPoint[Y]; b2 = Sqr(b);
  1397.   mindist2 = a2+b2;
  1398.  
  1399.   it_max = TPat->Vals.Fractal.Iterations;
  1400.  
  1401.   for (col = 0; col < it_max; col++)
  1402.   {
  1403.     b  = 2.0 * a * b + ci;
  1404.     a  = a2 - b2 + cr;
  1405.  
  1406.     a2 = Sqr(a);
  1407.     b2 = Sqr(b);
  1408.     dist2 = a2+b2;
  1409.  
  1410.     if(dist2 < mindist2) mindist2 = dist2;
  1411.     if(dist2 > 4.0)
  1412.     {
  1413.         cf = fractal_exterior_color(TPat, col, a, b);
  1414.         break;
  1415.     }
  1416.   }
  1417.  
  1418.   if(col == it_max)
  1419.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1420.  
  1421.   return(cf);
  1422. }
  1423.  
  1424. static DBL julia3 (VECTOR EPoint, TPATTERN *TPat)
  1425. {
  1426.   int it_max, col;
  1427.   DBL a, b, cf, a2, b2, dist2, mindist2,
  1428.       cr = TPat->Vals.Fractal.Coord[U], ci = TPat->Vals.Fractal.Coord[V];
  1429.  
  1430.   a = EPoint[X]; a2 = Sqr(a);
  1431.   b = EPoint[Y]; b2 = Sqr(b);
  1432.   mindist2 = a2+b2;
  1433.  
  1434.   it_max = TPat->Vals.Fractal.Iterations;
  1435.  
  1436.   for (col = 0; col < it_max; col++)
  1437.   {
  1438.     b = 3.0*a2*b - b2*b + ci;
  1439.     a = a2*a - 3.0*a*b2 + cr;
  1440.  
  1441.     a2 = Sqr(a);
  1442.     b2 = Sqr(b);
  1443.     dist2 = a2+b2;
  1444.  
  1445.     if(dist2 < mindist2) mindist2 = dist2;
  1446.     if(dist2 > 4.0)
  1447.     {
  1448.         cf = fractal_exterior_color(TPat, col, a, b);
  1449.       break;
  1450.     }
  1451.   }
  1452.  
  1453.   if(col == it_max)
  1454.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1455.  
  1456.   return(cf);
  1457. }
  1458.  
  1459. static DBL julia4 (VECTOR EPoint, TPATTERN *TPat)
  1460. {
  1461.   int it_max, col;
  1462.   DBL a, b, cf, a2, b2, dist2, mindist2,
  1463.       cr = TPat->Vals.Fractal.Coord[U], ci = TPat->Vals.Fractal.Coord[V];
  1464.  
  1465.   a = EPoint[X]; a2 = Sqr(a);
  1466.   b = EPoint[Y]; b2 = Sqr(b);
  1467.   mindist2 = a2+b2;
  1468.  
  1469.   it_max = TPat->Vals.Fractal.Iterations;
  1470.  
  1471.   for (col = 0; col < it_max; col++)
  1472.   {
  1473.     b = 4.0 * (a2*a*b - a*b2*b) + ci;
  1474.     a = a2*a2 - 6.0*a2*b2 + b2*b2 + cr;
  1475.  
  1476.     a2 = Sqr(a);
  1477.     b2 = Sqr(b);
  1478.     dist2 = a2+b2;
  1479.  
  1480.     if(dist2 < mindist2) mindist2 = dist2;
  1481.     if(dist2 > 4.0)
  1482.     {
  1483.         cf = fractal_exterior_color(TPat, col, a, b);
  1484.         break;
  1485.     }
  1486.   }
  1487.  
  1488.   if(col == it_max)
  1489.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1490.  
  1491.   return(cf);
  1492. }
  1493.  
  1494.  
  1495. static DBL magnet1m (VECTOR EPoint, TPATTERN *TPat)
  1496. {
  1497.   int it_max, col;
  1498.   DBL a, b, cf, a2, b2, x, y, tmp, tmp1r, tmp1i, tmp2r, tmp2i, dist2, mindist2;
  1499.  
  1500.   x = EPoint[X];
  1501.   y = EPoint[Y];
  1502.   a = a2 = 0;
  1503.   b = b2 = 0;
  1504.   mindist2 = 10000;
  1505.  
  1506.   it_max = TPat->Vals.Fractal.Iterations;
  1507.  
  1508.   for (col = 0; col < it_max; col++)
  1509.   {
  1510.       tmp1r = a2-b2 + x-1;
  1511.       tmp1i = 2*a*b + y;
  1512.       tmp2r = 2*a + x-2;
  1513.       tmp2i = 2*b + y;
  1514.       tmp = tmp2r*tmp2r + tmp2i*tmp2i;
  1515.       a = (tmp1r*tmp2r + tmp1i*tmp2i) / tmp;
  1516.       b = (tmp1i*tmp2r - tmp1r*tmp2i) / tmp;
  1517.       b2 = b*b;
  1518.       b = 2*a*b;
  1519.       a = a*a-b2;
  1520.  
  1521.       a2 = Sqr(a);
  1522.       b2 = Sqr(b);
  1523.       dist2 = a2+b2;
  1524.  
  1525.       if(dist2 < mindist2) mindist2 = dist2;
  1526.       tmp1r = a-1;
  1527.       if(dist2 > 10000.0 || tmp1r*tmp1r+b2 < 1/10000.0)
  1528.       {
  1529.           cf = fractal_exterior_color(TPat, col, a, b);
  1530.           break;
  1531.       }
  1532.   }
  1533.  
  1534.   if(col == it_max)
  1535.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1536.  
  1537.   return(cf);
  1538. }
  1539.  
  1540.  
  1541. static DBL magnet1j (VECTOR EPoint, TPATTERN *TPat)
  1542. {
  1543.   int it_max, col;
  1544.   DBL a, b, cf, a2, b2, tmp, tmp1r, tmp1i, tmp2r, tmp2i, dist2, mindist2,
  1545.       cr = TPat->Vals.Fractal.Coord[U], ci = TPat->Vals.Fractal.Coord[V];
  1546.  
  1547.   a = EPoint[X]; a2 = Sqr(a);
  1548.   b = EPoint[Y]; b2 = Sqr(b);
  1549.   mindist2 = a2+b2;
  1550.  
  1551.   it_max = TPat->Vals.Fractal.Iterations;
  1552.  
  1553.   for (col = 0; col < it_max; col++)
  1554.   {
  1555.       tmp1r = a2-b2 + cr-1;
  1556.       tmp1i = 2*a*b + ci;
  1557.       tmp2r = 2*a + cr-2;
  1558.       tmp2i = 2*b + ci;
  1559.       tmp = tmp2r*tmp2r + tmp2i*tmp2i;
  1560.       a = (tmp1r*tmp2r + tmp1i*tmp2i) / tmp;
  1561.       b = (tmp1i*tmp2r - tmp1r*tmp2i) / tmp;
  1562.       b2 = b*b;
  1563.       b = 2*a*b;
  1564.       a = a*a-b2;
  1565.  
  1566.       a2 = Sqr(a);
  1567.       b2 = Sqr(b);
  1568.       dist2 = a2+b2;
  1569.  
  1570.       if(dist2 < mindist2) mindist2 = dist2;
  1571.       tmp1r = a-1;
  1572.       if(dist2 > 10000.0 || tmp1r*tmp1r+b2 < 1/10000.0)
  1573.       {
  1574.           cf = fractal_exterior_color(TPat, col, a, b);
  1575.           break;
  1576.       }
  1577.   }
  1578.  
  1579.   if(col == it_max)
  1580.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1581.  
  1582.   return(cf);
  1583. }
  1584.  
  1585.  
  1586. static DBL magnet2m (VECTOR EPoint, TPATTERN *TPat)
  1587. {
  1588.   int it_max, col;
  1589.   DBL a, b, cf, a2, b2, x, y, tmp, tmp1r, tmp1i, tmp2r, tmp2i,
  1590.       c1r, c2r, c1c2r, c1c2i, dist2, mindist2;
  1591.  
  1592.   x = EPoint[X];
  1593.   y = EPoint[Y];
  1594.   a = a2 = 0;
  1595.   b = b2 = 0;
  1596.   mindist2 = 10000;
  1597.  
  1598.   c1r = x-1; c2r = x-2;
  1599.   c1c2r = c1r*c2r-y*y;
  1600.   c1c2i = (c1r+c2r)*y;
  1601.  
  1602.   it_max = TPat->Vals.Fractal.Iterations;
  1603.  
  1604.   for (col = 0; col < it_max; col++)
  1605.   {
  1606.       tmp1r = a2*a-3*a*b2 + 3*(a*c1r-b*y) + c1c2r;
  1607.       tmp1i = 3*a2*b-b2*b + 3*(a*y+b*c1r) + c1c2i;
  1608.       tmp2r = 3*(a2-b2) + 3*(a*c2r-b*y) + c1c2r + 1;
  1609.       tmp2i = 6*a*b + 3*(a*y+b*c2r) + c1c2i;
  1610.       tmp = tmp2r*tmp2r + tmp2i*tmp2i;
  1611.       a = (tmp1r*tmp2r + tmp1i*tmp2i) / tmp;
  1612.       b = (tmp1i*tmp2r - tmp1r*tmp2i) / tmp;
  1613.       b2 = b*b;
  1614.       b = 2*a*b;
  1615.       a = a*a-b2;
  1616.  
  1617.       a2 = Sqr(a);
  1618.       b2 = Sqr(b);
  1619.       dist2 = a2+b2;
  1620.  
  1621.       if(dist2 < mindist2) mindist2 = dist2;
  1622.       tmp1r = a-1;
  1623.       if(dist2 > 10000.0 || tmp1r*tmp1r+b2 < 1/10000.0)
  1624.       {
  1625.           cf = fractal_exterior_color(TPat, col, a, b);
  1626.           break;
  1627.       }
  1628.   }
  1629.  
  1630.   if(col == it_max)
  1631.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1632.  
  1633.   return(cf);
  1634. }
  1635.  
  1636.  
  1637. static DBL magnet2j (VECTOR EPoint, TPATTERN *TPat)
  1638. {
  1639.   int it_max, col;
  1640.   DBL a, b, cf, a2, b2, tmp, tmp1r, tmp1i, tmp2r, tmp2i, c1r,c2r,c1c2r,c1c2i,
  1641.       cr = TPat->Vals.Fractal.Coord[U], ci = TPat->Vals.Fractal.Coord[V],
  1642.       dist2, mindist2;
  1643.  
  1644.   a = EPoint[X]; a2 = Sqr(a);
  1645.   b = EPoint[Y]; b2 = Sqr(b);
  1646.   mindist2 = a2+b2;
  1647.  
  1648.   c1r = cr-1, c2r = cr-2;
  1649.   c1c2r = c1r*c2r-ci*ci;
  1650.   c1c2i = (c1r+c2r)*ci;
  1651.  
  1652.   it_max = TPat->Vals.Fractal.Iterations;
  1653.  
  1654.   for (col = 0; col < it_max; col++)
  1655.   {
  1656.       tmp1r = a2*a-3*a*b2 + 3*(a*c1r-b*ci) + c1c2r;
  1657.       tmp1i = 3*a2*b-b2*b + 3*(a*ci+b*c1r) + c1c2i;
  1658.       tmp2r = 3*(a2-b2) + 3*(a*c2r-b*ci) + c1c2r + 1;
  1659.       tmp2i = 6*a*b + 3*(a*ci+b*c2r) + c1c2i;
  1660.       tmp = tmp2r*tmp2r + tmp2i*tmp2i;
  1661.       a = (tmp1r*tmp2r + tmp1i*tmp2i) / tmp;
  1662.       b = (tmp1i*tmp2r - tmp1r*tmp2i) / tmp;
  1663.       b2 = b*b;
  1664.       b = 2*a*b;
  1665.       a = a*a-b2;
  1666.  
  1667.       a2 = Sqr(a);
  1668.       b2 = Sqr(b);
  1669.       dist2 = a2+b2;
  1670.  
  1671.       if(dist2 < mindist2) mindist2 = dist2;
  1672.       tmp1r = a-1;
  1673.       if(dist2 > 10000.0 || tmp1r*tmp1r+b2 < 1/10000.0)
  1674.       {
  1675.           cf = fractal_exterior_color(TPat, col, a, b);
  1676.           break;
  1677.       }
  1678.   }
  1679.  
  1680.   if(col == it_max)
  1681.       cf = fractal_interior_color(TPat, col, a, b, mindist2);
  1682.  
  1683.   return(cf);
  1684. }
  1685.  
  1686. /* NK ---- */
  1687.  
  1688. /*****************************************************************************
  1689. *
  1690. * FUNCTION
  1691. *
  1692. *   marble
  1693. *
  1694. * INPUT
  1695. *
  1696. *   EPoint -- The point in 3d space at which the pattern
  1697. *   is evaluated.
  1698. *   TPat   -- Texture pattern struct
  1699. *
  1700. * OUTPUT
  1701. *
  1702. * RETURNS
  1703. *
  1704. *   DBL value in the range 0.0 to 1.0
  1705. *
  1706. * AUTHOR
  1707. *
  1708. *   POV-Ray Team
  1709. *
  1710. * DESCRIPTION
  1711. *
  1712. * CHANGES
  1713. *   Oct 1994 : adapted from pigment by [CY]
  1714. *
  1715. ******************************************************************************/
  1716.  
  1717. static DBL marble (VECTOR EPoint, TPATTERN *TPat)
  1718. {
  1719.   register DBL turb_val;
  1720.   TURB *Turb;
  1721.  
  1722.   if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  1723.   {
  1724.     turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  1725.   }
  1726.   else
  1727.   {
  1728.     turb_val = 0.0;
  1729.   }
  1730.  
  1731.   return(EPoint[X] + turb_val);
  1732. }
  1733.  
  1734.  
  1735.  
  1736. /*****************************************************************************
  1737. *
  1738. * FUNCTION
  1739. *
  1740. *   onion
  1741. *
  1742. * INPUT
  1743. *
  1744. *   EPoint -- The point in 3d space at which the pattern
  1745. *   is evaluated.
  1746. *
  1747. * OUTPUT
  1748. *
  1749. * RETURNS
  1750. *
  1751. *   DBL value in the range 0.0 to 1.0
  1752. *
  1753. * AUTHOR
  1754. *
  1755. *   Scott Taylor
  1756. *
  1757. * DESCRIPTION
  1758. *
  1759. * CHANGES
  1760. *   Jul 1991 : Creation.
  1761. *   Oct 1994 : adapted from pigment by [CY]
  1762. *
  1763. ******************************************************************************/
  1764.  
  1765. static DBL onion (VECTOR EPoint)
  1766. {
  1767.   /* The variable noise is not used as noise in this function */
  1768.  
  1769.   register DBL noise;
  1770.  
  1771. /*
  1772.    This ramp goes 0-1,1-0,0-1,1-0...
  1773.  
  1774.    noise = (fmod(sqrt(Sqr(x)+Sqr(y)+Sqr(z)),2.0)-1.0);
  1775.  
  1776.    if (noise<0.0) {noise = 0.0-noise;}
  1777. */
  1778.  
  1779.   /* This ramp goes 0-1, 0-1, 0-1, 0-1 ... */
  1780.  
  1781.   noise = (fmod(sqrt(Sqr(EPoint[X])+Sqr(EPoint[Y])+Sqr(EPoint[Z])), 1.0));
  1782.  
  1783.   return(noise);
  1784. }
  1785.  
  1786.  
  1787.  
  1788. /*****************************************************************************
  1789. *
  1790. * FUNCTION
  1791. *
  1792. *   radial
  1793. *
  1794. * INPUT
  1795. *
  1796. *   EPoint -- The point in 3d space at which the pattern
  1797. *   is evaluated.
  1798. *
  1799. * OUTPUT
  1800. *
  1801. * RETURNS
  1802. *
  1803. *   DBL value in the range 0.0 to 1.0
  1804. *
  1805. * AUTHOR
  1806. *
  1807. *   Chris Young -- new in vers 2.0
  1808. *
  1809. * DESCRIPTION
  1810. *
  1811. * CHANGES
  1812. *   Oct 1994 : adapted from pigment by [CY]
  1813. *
  1814. ******************************************************************************/
  1815.  
  1816. static DBL radial (VECTOR EPoint)
  1817. {
  1818.   register DBL value;
  1819.  
  1820.   if ((fabs(EPoint[X])<0.001) && (fabs(EPoint[Z])<0.001))
  1821.   {
  1822.     value = 0.25;
  1823.   }
  1824.   else
  1825.   {
  1826.     value = 0.25 + (atan2(EPoint[X],EPoint[Z]) + M_PI) / TWO_M_PI;
  1827.   }
  1828.  
  1829.   return(value);
  1830. }
  1831.  
  1832.  
  1833. #ifdef PolaricalPatch
  1834. /*****************************************************************************
  1835. *
  1836. * FUNCTION
  1837. *
  1838. *   polarical
  1839. *
  1840. * INPUT
  1841. *
  1842. *   EPoint -- The point in 3d space at which the pattern
  1843. *   is evaluated.
  1844. *
  1845. * OUTPUT
  1846. *
  1847. * RETURNS
  1848. *
  1849. *   DBL value in the range 0.0 to 1.0
  1850. *
  1851. * AUTHOR
  1852. *
  1853. *   Krzysztof Garus
  1854. *
  1855. * DESCRIPTION
  1856. *
  1857. * CHANGES
  1858. *
  1859. *
  1860. ******************************************************************************/
  1861.  
  1862. static DBL polarical (VECTOR EPoint)
  1863. {
  1864.   register DBL value;
  1865.   DBL radius = sqrt(Sqr(EPoint[X])+Sqr(EPoint[Z]));
  1866.  
  1867.   if ((radius<0.001) && (fabs(EPoint[Y])<0.001))
  1868.   {
  1869.     value = 0;
  1870.   }
  1871.   else
  1872.   {
  1873.     value = atan2(radius,-EPoint[Y]) / M_PI;
  1874.     
  1875.   }
  1876.  
  1877.   return(value);
  1878. }
  1879. #endif
  1880. /*****************************************************************************
  1881. *
  1882. * FUNCTION
  1883. *
  1884. *   spiral1
  1885. *
  1886. * INPUT
  1887. *
  1888. *   EPoint -- The point in 3d space at which the pattern
  1889. *   is evaluated.
  1890. *   TPat   -- Texture pattern struct
  1891. *
  1892. * OUTPUT
  1893. *
  1894. * RETURNS
  1895. *
  1896. *   DBL value in the range 0.0 to 1.0
  1897. *
  1898. * AUTHOR
  1899. *
  1900. *   Dieter Bayer
  1901. *
  1902. * DESCRIPTION
  1903. *   Spiral whirles around z-axis.
  1904. *   The number of "arms" is defined in the TPat.
  1905. *
  1906. * CHANGES
  1907. *   Aug 1994 : Creation.
  1908. *   Oct 1994 : adapted from pigment by [CY]
  1909. *
  1910. ******************************************************************************/
  1911.  
  1912. static DBL spiral1(VECTOR EPoint, TPATTERN *TPat)
  1913. {
  1914.   DBL rad, phi, turb_val;
  1915.   DBL x = EPoint[X];
  1916.   DBL y = EPoint[Y];
  1917.   DBL z = EPoint[Z];
  1918.   TURB *Turb;
  1919.  
  1920.   if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  1921.   {
  1922.     turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  1923.   }
  1924.   else
  1925.   {
  1926.     turb_val = 0.0;
  1927.   }
  1928.  
  1929.   /* Get distance from z-axis. */
  1930.  
  1931.   rad = sqrt(x * x + y * y);
  1932.  
  1933.   /* Get angle in x,y-plane (0...2 PI). */
  1934.  
  1935.   if (rad == 0.0)
  1936.   {
  1937.     phi = 0.0;
  1938.   }
  1939.   else
  1940.   {
  1941.     if (x < 0.0)
  1942.     {
  1943.       phi = 3.0 * M_PI_2 - asin(y / rad);
  1944.     }
  1945.     else
  1946.     {
  1947.       phi = M_PI_2 + asin(y / rad);
  1948.     }
  1949.   }
  1950.  
  1951.   return(z + rad + (DBL)TPat->Vals.Arms * phi / TWO_M_PI + turb_val);
  1952. }
  1953.  
  1954.  
  1955.  
  1956. /*****************************************************************************
  1957. *
  1958. * FUNCTION
  1959. *
  1960. *   spiral2
  1961. *
  1962. * INPUT
  1963. *
  1964. *   EPoint -- The point in 3d space at which the pattern
  1965. *   is evaluated.
  1966. *   TPat   -- Texture pattern struct
  1967. *
  1968. * OUTPUT
  1969. *
  1970. * RETURNS
  1971. *
  1972. *   DBL value in the range 0.0 to 1.0
  1973. *
  1974. * AUTHOR
  1975. *
  1976. *   Dieter Bayer
  1977. *
  1978. * DESCRIPTION
  1979. *   Spiral whirles around z-axis.
  1980. *   The number of "arms" is defined in the TPat.
  1981. *
  1982. * CHANGES
  1983. *   Aug 1994 : Creation.
  1984. *   Oct 1994 : adapted from pigment by [CY]
  1985. *
  1986. ******************************************************************************/
  1987.  
  1988.  
  1989. static DBL spiral2(VECTOR EPoint, TPATTERN *TPat)
  1990. {
  1991.   DBL rad, phi, turb_val;
  1992.   DBL x = EPoint[X];
  1993.   DBL y = EPoint[Y];
  1994.   DBL z = EPoint[Z];
  1995.   TURB *Turb;
  1996.  
  1997.   if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  1998.   {
  1999.     turb_val = Turb->Turbulence[X] * Turbulence(EPoint,Turb);
  2000.   }
  2001.   else
  2002.   {
  2003.     turb_val = 0.0;
  2004.   }
  2005.  
  2006.   /* Get distance from z-axis. */
  2007.  
  2008.   rad = sqrt(x * x + y * y);
  2009.  
  2010.   /* Get angle in x,y-plane (0...2 PI) */
  2011.  
  2012.   if (rad == 0.0)
  2013.   {
  2014.     phi = 0.0;
  2015.   }
  2016.   else
  2017.   {
  2018.     if (x < 0.0)
  2019.     {
  2020.       phi = 3.0 * M_PI_2 - asin(y / rad);
  2021.     }
  2022.     else
  2023.     {
  2024.       phi = M_PI_2 + asin(y / rad);
  2025.     }
  2026.   }
  2027.  
  2028.   turb_val = Triangle_Wave(z + rad + (DBL)TPat->Vals.Arms * phi / TWO_M_PI +
  2029.                            turb_val);
  2030.  
  2031.   return(Triangle_Wave(rad) + turb_val);
  2032. }
  2033.  
  2034.  
  2035.  
  2036. /*****************************************************************************
  2037. *
  2038. * FUNCTION
  2039. *
  2040. *   wood
  2041. *
  2042. * INPUT
  2043. *
  2044. *   EPoint -- The point in 3d space at which the pattern
  2045. *   is evaluated.
  2046. *   TPat   -- Texture pattern struct
  2047. *
  2048. * OUTPUT
  2049. *
  2050. * RETURNS
  2051. *
  2052. *   DBL value in the range 0.0 to 1.0
  2053. *
  2054. * AUTHOR
  2055. *
  2056. *   POV-Ray Team
  2057. *
  2058. * DESCRIPTION
  2059. *
  2060. * CHANGES
  2061. *   Oct 1994 : adapted from pigment by [CY]
  2062. *
  2063. ******************************************************************************/
  2064.  
  2065.  
  2066. static DBL wood (VECTOR EPoint, TPATTERN *TPat)
  2067. {
  2068.   register DBL length;
  2069.   VECTOR WoodTurbulence;
  2070.   VECTOR point;
  2071.   DBL x=EPoint[X];
  2072.   DBL y=EPoint[Y];
  2073.   TURB *Turb;
  2074.  
  2075.   if ((Turb=Search_For_Turb(TPat->Warps)) != NULL)
  2076.   {
  2077.     DTurbulence (WoodTurbulence, EPoint,Turb);
  2078.     point[X] = cycloidal((x + WoodTurbulence[X]) * Turb->Turbulence[X]);
  2079.     point[Y] = cycloidal((y + WoodTurbulence[Y]) * Turb->Turbulence[Y]);
  2080.   }
  2081.   else
  2082.   {
  2083.     point[X] = 0.0;
  2084.     point[Y] = 0.0;
  2085.   }
  2086.   point[Z] = 0.0;
  2087.  
  2088.   point[X] += x;
  2089.   point[Y] += y;
  2090.  
  2091.   /* point[Z] += z; Deleted per David Buck --  BP 7/91 */
  2092.  
  2093.   VLength (length, point);
  2094.  
  2095.   return(length);
  2096. }
  2097.  
  2098.  
  2099. /*****************************************************************************
  2100. *
  2101. * FUNCTION
  2102. *
  2103. *   hexagon
  2104. *
  2105. * INPUT
  2106. *
  2107. *   EPoint -- The point in 3d space at which the pattern
  2108. *   is evaluated.
  2109. *
  2110. * OUTPUT
  2111. *
  2112. * RETURNS
  2113. *
  2114. *   DBL value exactly 0.0, 1.0 or 2.0
  2115. *
  2116. * AUTHOR
  2117. *
  2118. *   Ernest MacDougal Campbell III
  2119. *   
  2120. * DESCRIPTION
  2121. *
  2122. *   TriHex pattern -- Ernest MacDougal Campbell III (EMC3) 11/23/92
  2123. *
  2124. *   Creates a hexagon pattern in the XZ plane.
  2125. *
  2126. *   This algorithm is hard to explain.  First it scales the point to make
  2127. *   a few of the later calculations easier, then maps some points to be
  2128. *   closer to the Origin.  A small area in the first quadrant is subdivided
  2129. *   into a 6 x 6 grid.  The position of the point mapped into that grid
  2130. *   determines its color.  For some points, just the grid location is enough,
  2131. *   but for others, we have to calculate which half of the block it's in
  2132. *   (this is where the atan2() function comes in handy).
  2133. *
  2134. * CHANGES
  2135. *   Nov 1992 : Creation.
  2136. *   Oct 1994 : adapted from pigment by [CY]
  2137. *
  2138. ******************************************************************************/
  2139.  
  2140. #define xfactor 0.5;         /* each triangle is split in half for the grid */
  2141. #define zfactor 0.866025404; /* sqrt(3)/2 -- Height of an equilateral triangle */
  2142.  
  2143. static DBL hexagon (VECTOR EPoint)
  2144. {
  2145.   int xm, zm;
  2146.   int brkindx;
  2147.   DBL xs, zs, xl, zl, value = 0.0;
  2148.   DBL x=EPoint[X];
  2149.   DBL z=EPoint[Z];
  2150.  
  2151.  
  2152.   /* Keep all numbers positive.  Also, if z is negative, map it in such a
  2153.    * way as to avoid mirroring across the x-axis.  The value 5.196152424
  2154.    * is (sqrt(3)/2) * 6 (because the grid is 6 blocks high)
  2155.    */
  2156.  
  2157.   x = fabs(x);
  2158.  
  2159.   /* Avoid mirroring across x-axis. */
  2160.  
  2161.   z = z < 0.0 ? 5.196152424 - fabs(z) : z;
  2162.  
  2163.   /* Scale point to make calcs easier. */
  2164.  
  2165.   xs = x/xfactor;
  2166.   zs = z/zfactor;
  2167.  
  2168.   /* Map points into the 6 x 6 grid where the basic formula works. */
  2169.  
  2170.   xs -= floor(xs/6.0) * 6.0;
  2171.   zs -= floor(zs/6.0) * 6.0;
  2172.  
  2173.   /* Get a block in the 6 x 6 grid. */
  2174.  
  2175.   xm = (int) FLOOR(xs) % 6;
  2176.   zm = (int) FLOOR(zs) % 6;
  2177.  
  2178.   switch (xm)
  2179.   {
  2180.     /* These are easy cases: Color depends only on xm and zm. */
  2181.  
  2182.     case 0:
  2183.     case 5:
  2184.  
  2185.       switch (zm)
  2186.       {
  2187.         case 0:
  2188.         case 5: value = 0; break;
  2189.  
  2190.         case 1:
  2191.         case 2: value = 1; break;
  2192.  
  2193.         case 3:
  2194.         case 4: value = 2; break;
  2195.       }
  2196.  
  2197.       break;
  2198.  
  2199.     case 2:
  2200.     case 3:
  2201.  
  2202.       switch (zm)
  2203.       {
  2204.         case 0:
  2205.         case 1: value = 2; break;
  2206.  
  2207.         case 2:
  2208.         case 3: value = 0; break;
  2209.  
  2210.         case 4:
  2211.         case 5: value = 1; break;
  2212.       }
  2213.  
  2214.       break;
  2215.  
  2216.     /* These cases are harder.  These blocks are divided diagonally
  2217.      * by the angled edges of the hexagons.  Some slope positive, and
  2218.      * others negative.  We flip the x value of the negatively sloped
  2219.      * pieces.  Then we check to see if the point in question falls
  2220.      * in the upper or lower half of the block.  That info, plus the
  2221.      * z status of the block determines the color.
  2222.      */
  2223.  
  2224.     case 1:
  2225.     case 4:
  2226.  
  2227.       /* Map the point into the block at the origin. */
  2228.  
  2229.       xl = xs-xm;
  2230.       zl = zs-zm;
  2231.  
  2232.       /* These blocks have negative slopes so we flip it horizontally. */
  2233.  
  2234.       if (((xm + zm) % 2) == 1)
  2235.       {
  2236.         xl = 1.0 - xl;
  2237.       }
  2238.  
  2239.       /* Avoid a divide-by-zero error. */
  2240.  
  2241.       if (xl == 0.0)
  2242.       {
  2243.         xl = 0.0001;
  2244.       }
  2245.  
  2246.       /* Is the angle less-than or greater-than 45 degrees? */
  2247.  
  2248.       brkindx = (zl / xl) < 1.0;
  2249.  
  2250.       /* was...
  2251.        * brkindx = (atan2(zl,xl) < (45 * M_PI_180));
  2252.        * ...but because of the mapping, it's easier and cheaper,
  2253.        * CPU-wise, to just use a good ol' slope.
  2254.        */
  2255.  
  2256.       switch (brkindx)
  2257.       {
  2258.         case TRUE:
  2259.  
  2260.           switch (zm)
  2261.           {
  2262.             case 0:
  2263.             case 3: value = 0; break;
  2264.  
  2265.             case 2:
  2266.             case 5: value = 1; break;
  2267.  
  2268.             case 1:
  2269.             case 4: value = 2; break;
  2270.           }
  2271.  
  2272.           break;
  2273.  
  2274.         case FALSE:
  2275.  
  2276.           switch (zm)
  2277.           {
  2278.             case 0:
  2279.             case 3: value = 2; break;
  2280.  
  2281.             case 2:
  2282.             case 5: value = 0; break;
  2283.  
  2284.             case 1:
  2285.             case 4: value = 1; break;
  2286.           }
  2287.  
  2288.           break;
  2289.       }
  2290.   }
  2291.  
  2292.   value = fmod(value, 3.0);
  2293.  
  2294.   return(value);
  2295. }
  2296.  
  2297. /* In addition to clipping the value to 
  2298.    lie between 0.0 to 1.0, it also fudges 1.0-value.
  2299.  */
  2300.  
  2301. #define CLIP_DENSITY(r) if((r)<0.0){(r)=1.0;}else{if((r)>1.0){(r)=0.0;}else{(r)=1.0-(r);}}
  2302.  
  2303. static DBL planar_pattern (VECTOR EPoint)
  2304. {
  2305.   register DBL value;
  2306.  
  2307.   value = fabs(EPoint[Y]);
  2308.   CLIP_DENSITY(value);
  2309.  
  2310.   return(value);
  2311. }
  2312.  
  2313. static DBL spherical (VECTOR EPoint)
  2314. {
  2315.   register DBL value;
  2316.  
  2317.   VLength(value, EPoint);
  2318.   CLIP_DENSITY(value);
  2319.  
  2320.   return(value);
  2321. }
  2322.  
  2323. static DBL boxed (VECTOR EPoint)
  2324. {
  2325.   register DBL value;
  2326.  
  2327.   value = max(fabs(EPoint[X]), max(fabs(EPoint[Y]), fabs(EPoint[Z])));
  2328.   CLIP_DENSITY(value);
  2329.  
  2330.   return(value);
  2331. }
  2332.  
  2333. static DBL cylindrical (VECTOR EPoint)
  2334. {
  2335.   register DBL value;
  2336.  
  2337.   value = sqrt(Sqr(EPoint[X]) + Sqr(EPoint[Z]));
  2338.   CLIP_DENSITY(value);
  2339.  
  2340.   return(value);
  2341. }
  2342.  
  2343.  
  2344.  
  2345. /*****************************************************************************
  2346. *
  2347. * FUNCTION
  2348. *
  2349. *   PickInCube(tv, p1)
  2350. *
  2351. * INPUT
  2352. *
  2353. *   ?
  2354. *
  2355. * OUTPUT
  2356. *   
  2357. * RETURNS
  2358. *
  2359. *   long integer hash function used, to speed up cacheing.
  2360. *   
  2361. * AUTHOR
  2362. *
  2363. *   Jim McElhiney
  2364. *   
  2365. * DESCRIPTION
  2366. *
  2367. *   A subroutine to go with crackle.
  2368. *
  2369. *   Pick a random point in the same unit-sized cube as tv, in a
  2370. *   predictable way, so that when called again with another point in
  2371. *   the same unit cube, p1 is picked to be the same.
  2372. *
  2373. * CHANGES
  2374. *
  2375. ******************************************************************************/
  2376.  
  2377. #ifdef CracklePatch
  2378. long PickInCube(VECTOR tv, VECTOR  p1)
  2379. #else
  2380. static long PickInCube(VECTOR tv, VECTOR  p1)
  2381. #endif
  2382. {
  2383.   int seed, temp;
  2384.   VECTOR flo;
  2385.  
  2386.   /*
  2387.    * This uses floor() not FLOOR, so it will not be a mirror
  2388.    * image about zero in the range -1.0 to 1.0. The viewer
  2389.    * won't see an artefact around the origin.
  2390.    */
  2391.  
  2392.   flo[X] = floor(tv[X] - EPSILON);
  2393.   flo[Y] = floor(tv[Y] - EPSILON);
  2394.   flo[Z] = floor(tv[Z] - EPSILON);
  2395.  
  2396. #ifdef NoiseTranslateFixPatch
  2397.   seed = Hash3d((int)flo[X]&0xFFF, (int)flo[Y]&0xFFF, (int)flo[Z]&0xFFF);
  2398. #else
  2399.   seed = Hash3d((int)flo[X], (int)flo[Y], (int)flo[Z]);
  2400. #endif
  2401.   temp = POV_GET_OLD_RAND(); /* save current seed */
  2402.   
  2403.   POV_SRAND(seed);
  2404.  
  2405.   p1[X] = flo[X] + FRAND();
  2406.   p1[Y] = flo[Y] + FRAND();
  2407.   p1[Z] = flo[Z] + FRAND();
  2408.  
  2409.   POV_SRAND(temp);  /* restore */
  2410.  
  2411.   return((long)seed);
  2412. }
  2413.  
  2414.  
  2415.  
  2416. /*****************************************************************************
  2417. *
  2418. * FUNCTION
  2419. *
  2420. *   Evaluate_Pattern
  2421. *
  2422. * INPUT
  2423. *
  2424. *   EPoint -- The point in 3d space at which the pattern
  2425. *   is evaluated.
  2426. *   TPat   -- Texture pattern struct
  2427. *   Intersection - intersection structure
  2428. *   
  2429. * OUTPUT
  2430. *   
  2431. * RETURNS
  2432. *
  2433. *   DBL result usual 0.0 to 1.0 but may be 2.0 in hexagon
  2434. *   
  2435. * AUTHOR
  2436. *
  2437. *   adapted from Add_Pigment by Chris Young
  2438. *   
  2439. * DESCRIPTION
  2440. *
  2441. * CHANGES
  2442. *   added parameter Intersection   -hdf- May 98
  2443. *  Removed Warp_EPoint call - moved it outside
  2444. *
  2445. ******************************************************************************/
  2446.  
  2447. DBL Evaluate_TPat (TPATTERN *TPat, VECTOR EPoint, INTERSECTION *Intersection)
  2448. {
  2449.   DBL value = 0.0;
  2450.  
  2451.   /* NK 19 Nov 1999 removed Warp_EPoint call */
  2452.  
  2453.   switch (TPat->Type)
  2454.   {
  2455.     case AGATE_PATTERN:    value = agate    (EPoint, TPat);   break;
  2456.  
  2457.     case BOZO_PATTERN:
  2458.     case SPOTTED_PATTERN:
  2459.     case BUMPS_PATTERN:    value = Noise    (EPoint);         break;
  2460.  
  2461.     case BRICK_PATTERN:    value = brick    (EPoint, TPat);   break;
  2462.     case CHECKER_PATTERN:  value = checker  (EPoint);         break;
  2463. #ifdef TrianglulairSquarePatch 
  2464.     case SQUARE_PATTERN:   value = square   (EPoint);         break;
  2465.     case TERNAIRE_PATTERN: value = ternaire (EPoint);         break;
  2466. #endif
  2467. #ifdef CracklePatch
  2468.     case CRACKLE_PATTERN:  value = crackle  (EPoint, TPat);   break;
  2469. #else
  2470.     case CRACKLE_PATTERN:  value = crackle  (EPoint);         break;
  2471. #endif
  2472.     case GRADIENT_PATTERN: value = gradient (EPoint, TPat);   break;
  2473.     case GRANITE_PATTERN:  value = granite  (EPoint);         break;
  2474.     case HEXAGON_PATTERN:  value = hexagon  (EPoint);         break;
  2475.     case LEOPARD_PATTERN:  value = leopard  (EPoint);         break;
  2476.     case MAGNET1M_PATTERN: value = magnet1m (EPoint, TPat);   break;
  2477.     case MAGNET1J_PATTERN: value = magnet1j (EPoint, TPat);   break;
  2478.     case MAGNET2M_PATTERN: value = magnet2m (EPoint, TPat);   break;
  2479.     case MAGNET2J_PATTERN: value = magnet2j (EPoint, TPat);   break;
  2480.     case MANDEL_PATTERN:   value = mandel   (EPoint, TPat);   break;
  2481.     case MANDEL3_PATTERN:  value = mandel3  (EPoint, TPat);   break;
  2482.     case MANDEL4_PATTERN:  value = mandel4  (EPoint, TPat);   break;
  2483.     case JULIA_PATTERN:    value = julia    (EPoint, TPat);   break;
  2484.     case JULIA3_PATTERN:   value = julia3   (EPoint, TPat);   break;
  2485.     case JULIA4_PATTERN:   value = julia4   (EPoint, TPat);   break;
  2486.     case MARBLE_PATTERN:   value = marble   (EPoint, TPat);   break;
  2487.     case ONION_PATTERN:    value = onion    (EPoint);         break;
  2488.     case RADIAL_PATTERN:   value = radial   (EPoint);         break;
  2489. #ifdef PolaricalPatch
  2490.     case POLARICAL_PATTERN:value = polarical(EPoint);         break;
  2491. #endif    
  2492.       case SPIRAL1_PATTERN:  value = spiral1  (EPoint, TPat);   break;
  2493.     case SPIRAL2_PATTERN:  value = spiral2  (EPoint, TPat);   break;
  2494.     case WOOD_PATTERN:     value = wood     (EPoint, TPat);   break;
  2495. #ifdef CellsPatch
  2496.     case CELLS_PATTERN:value = cells    (EPoint);         break;
  2497. #endif
  2498. #ifdef VanSicklePatternPatch
  2499.     case BLOTCHES_PATTERN: value = blotches    (EPoint);         break;
  2500.     case BANDS_PATTERN:value = banded    (EPoint, TPat);   break;
  2501.     case SHEET_PATTERN:value = sheet    (EPoint, TPat);   break;
  2502. #endif
  2503.  
  2504.     case WAVES_PATTERN:    value = waves_pigm    (EPoint, TPat);   break;
  2505.     case RIPPLES_PATTERN:  value = ripples_pigm  (EPoint, TPat);   break;
  2506.     case WRINKLES_PATTERN: value = wrinkles_pigm (EPoint);   break;
  2507.     case DENTS_PATTERN:    value = dents_pigm    (EPoint);   break;
  2508.     case QUILTED_PATTERN:  value = quilted_pigm  (EPoint, TPat);   break;
  2509.  
  2510.     case PLANAR_PATTERN:      value = planar_pattern (EPoint);      break;
  2511.     case BOXED_PATTERN:       value = boxed          (EPoint);      break;
  2512.     case SPHERICAL_PATTERN:   value = spherical      (EPoint);      break;
  2513.     case CYLINDRICAL_PATTERN: value = cylindrical    (EPoint);      break;
  2514.     case DENSITY_FILE_PATTERN:value = density_file (EPoint, TPat);  break;
  2515. #ifdef SolidPatternPatch
  2516.     case SOLID_PATTERN       :value = SolidPat     (EPoint, TPat);  break; /*Chris Huff solid pattern*/
  2517. #endif
  2518. #ifdef ClothPatternPatch
  2519.     case CLOTH_PATTERN       :value = ClothPat(EPoint);        break;/*Chris Huff cloth pattern*/
  2520.     case CLOTH2_PATTERN       :value = Cloth2Pat(EPoint);        break;/*Chris Huff cloth2 pattern*/
  2521. #endif
  2522. #ifdef TorodialPatch
  2523.     case TOROIDAL_SPIRAL_PATTERN :value = ToroidalSpiral(EPoint, TPat);   break;/*Chris Huff toroidalSpiral pattern*/
  2524. #endif
  2525. #ifdef BlobPatternPatch
  2526.     case BLOB_PATTERN:        value = blob_pattern (EPoint, TPat);  break;/*Chris Huff blob pattern*/
  2527. #endif
  2528. #ifdef ObjectPatternPatch
  2529.     case OBJECT_PATTERN:      value = object(EPoint, TPat);         break;/*Chris Huff object pattern*/
  2530. #endif
  2531. #ifdef ProximityPatch
  2532.     case PROXIMITY_PATTERN:   value = proximity(EPoint, TPat);      break;/*Chris Huff proximity pattern*/
  2533. #endif
  2534. /** poviso: July 96  R.S. **/
  2535. #ifdef POVISO
  2536.     case FUNCTION_PATTERN:  value = func_pat  (EPoint, TPat);   break;
  2537. #endif
  2538. /** --- **/
  2539.     /* -hdf- Apr 98 */
  2540.     case SLOPE_PATTERN:    
  2541.         #if defined (EvalPatternPatch) || defined(EvalPigmentPatch)
  2542.       if (!Intersection)
  2543.         Error("Cannot use slope pattern for eval_pigment, eval_pattern, or sky_sphere.\n");
  2544.           #endif
  2545.       value = slope (EPoint, TPat, Intersection);   
  2546.       break;
  2547.  
  2548.     
  2549.     /* NK 1998 */
  2550.     case IMAGE_PATTERN:    value = image_pattern(EPoint, TPat); break;
  2551.     /* NK ---- */
  2552.  
  2553. #ifdef PigmentPatternPatch
  2554.     case PIGMENT_PATTERN: value = pigment_pattern(EPoint, TPat, Intersection); break;
  2555. #endif
  2556.  
  2557.     default: Error("Problem in Evaluate_TPat.");
  2558.   }
  2559.  
  2560.   if (TPat->Frequency !=0.0)
  2561.   {
  2562.     value = fmod(value * TPat->Frequency + TPat->Phase, 1.00001);
  2563.   }
  2564.  
  2565.   /* allow negative Frequency */
  2566.  
  2567.   if (value < 0.0)
  2568.   {
  2569.     value -= floor(value);
  2570.   }
  2571.  
  2572.   switch (TPat->Wave_Type)
  2573.   {
  2574.     case RAMP_WAVE:
  2575.       break;
  2576.  
  2577.     case SINE_WAVE:
  2578.       value = (1.0+cycloidal(value))*0.5;
  2579.       break;
  2580.  
  2581.     case TRIANGLE_WAVE:
  2582.       value = Triangle_Wave(value);
  2583.       break;
  2584.  
  2585.     case SCALLOP_WAVE:
  2586.       value = fabs(cycloidal(value*0.5));
  2587.       break;
  2588.  
  2589.     case CUBIC_WAVE:
  2590.       value = Sqr(value)*((-2.0 * value) + 3.0);
  2591.       break;
  2592.  
  2593.     case POLY_WAVE:
  2594.       value = pow(value, TPat->Exponent);
  2595.       break;
  2596.  
  2597. #ifdef AtanWavePatch
  2598.     case ATAN_WAVE:
  2599.       value = (atan(value*M_PI_2)/M_PI_2)*0.5+0.5;
  2600.       break;
  2601. #endif
  2602.  
  2603. #ifdef SplineWavePatch
  2604.     case SPLINE_WAVE:
  2605.     {
  2606.         EXPRESS tmp;
  2607.         value = Get_Spline_Val(TPat->spline_wave, value, tmp);
  2608.     }
  2609.     break;
  2610. #endif
  2611.     default: Error("Unknown Wave Type %d.",TPat->Wave_Type);
  2612.    }
  2613.  
  2614.   return(value);
  2615. }
  2616.  
  2617.  
  2618. /*****************************************************************************
  2619. *
  2620. * FUNCTION
  2621. *
  2622. * INPUT
  2623. *
  2624. * OUTPUT
  2625. *
  2626. * RETURNS
  2627. *
  2628. * AUTHOR
  2629. *
  2630. * DESCRIPTION
  2631. *
  2632. * CHANGES
  2633. *
  2634. ******************************************************************************/
  2635.  
  2636. void Init_TPat_Fields (TPATTERN *Tpat)
  2637. {
  2638.   Tpat->Type       = NO_PATTERN;
  2639.   Tpat->Wave_Type  = RAMP_WAVE;
  2640.   Tpat->Flags      = NO_FLAGS;
  2641.   Tpat->References = 1;
  2642.   Tpat->Exponent   = 1.0;
  2643.   Tpat->Frequency  = 1.0;
  2644.   Tpat->Phase      = 0.0;
  2645.   Tpat->Warps      = NULL;
  2646.   Tpat->Next       = NULL;
  2647.   Tpat->Blend_Map  = NULL;
  2648. #ifdef SplineWavePatch
  2649.   Tpat->spline_wave  = NULL;
  2650. #endif  
  2651. }
  2652.  
  2653.  
  2654.  
  2655. /*****************************************************************************
  2656. *
  2657. * FUNCTION
  2658. *
  2659. * INPUT
  2660. *
  2661. * OUTPUT
  2662. *
  2663. * RETURNS
  2664. *
  2665. * AUTHOR
  2666. *
  2667. * DESCRIPTION
  2668. *
  2669. * CHANGES
  2670. *
  2671. ******************************************************************************/
  2672.  
  2673. void Copy_TPat_Fields (TPATTERN *New, TPATTERN  *Old)
  2674. {
  2675.   /* is this necessary, or even wanted? */
  2676.   /*
  2677.   if ( New->Type == CRACKLE_PATTERN )
  2678.   {
  2679.     if (New->Vals.Crackle.cv)
  2680.       POV_FREE(New->Vals.Crackle.cv);
  2681.   }
  2682.   */
  2683.   *New = *Old;
  2684.   
  2685.   /* Copy warp chain */
  2686.   New->Warps = Copy_Warps(Old->Warps);
  2687.  
  2688.   New->Blend_Map = Copy_Blend_Map(Old->Blend_Map);
  2689.  
  2690.   /* Note, cannot copy Old->Next because we don't know what kind of
  2691.      thing this is.  It must be copied by Copy_Pigment, Copy_Tnormal etc.
  2692.   */
  2693.  
  2694.   /* NK 1998 - added IMAGE_PATTERN */
  2695.   if ((Old->Type == BITMAP_PATTERN) || (Old->Type == IMAGE_PATTERN))
  2696.   {
  2697.      New->Vals.Image = Copy_Image(Old->Vals.Image);
  2698.   }
  2699.  
  2700. #ifdef SplineWavePatch
  2701.   if(Old->spline_wave != NULL)
  2702.       New->spline_wave = Copy_Spline(Old->spline_wave);
  2703. #endif
  2704.   
  2705. #ifdef ObjectPatternPatch
  2706.   if (Old->Type == OBJECT_PATTERN)
  2707.   {
  2708.     if(Old->Vals.Object != NULL)
  2709.     {
  2710.       New->Vals.Object = (OBJECT*)Copy_Object(Old->Vals.Object);
  2711.     }
  2712.   }
  2713. #endif
  2714.  
  2715. #ifdef ProximityPatch
  2716.   if (Old->Type == PROXIMITY_PATTERN)
  2717.   {
  2718.     if(Old->Vals.Proximity.proxObject != NULL)
  2719.     {
  2720.       New->Vals.Proximity.proxObject = (OBJECT*)Copy_Object(Old->Vals.Proximity.proxObject);
  2721.     }
  2722.   }
  2723. #endif
  2724.  
  2725. #ifdef BlobPatternPatch
  2726.   if(Old->Type == BLOB_PATTERN ||
  2727.      Old->Type == BLOB_PIGMENT)
  2728.   {
  2729.     if(Old->Vals.Blob.blob_dat != NULL)
  2730.     {
  2731.         BLOB_PATTERN_DATA * currentComponent = Old->Vals.Blob.blob_dat;
  2732.         BLOB_PATTERN_DATA * nextComponent = currentComponent->next;/*tPat->blob_dat->next;*/
  2733.         BLOB_PATTERN_DATA * destComponent = NULL;
  2734.       
  2735.         while(currentComponent != NULL)/*Chris Huff blob pattern*/
  2736.         {
  2737.             BLOB_PATTERN_DATA * componentCopy = (BLOB_PATTERN_DATA *)POV_MALLOC(sizeof(BLOB_PATTERN_DATA), "spherical blob pattern component");
  2738.             nextComponent = currentComponent->next;
  2739.           
  2740.             /*START: copy component*/
  2741.             *componentCopy = *currentComponent;
  2742.             if(currentComponent->transform)
  2743.                 componentCopy->transform = Copy_Transform(currentComponent->transform);
  2744.                 
  2745.             if(currentComponent->pigment)
  2746.                 componentCopy->pigment = Copy_Pigment(currentComponent->pigment);
  2747. /*            if(currentComponent->pattern)
  2748.                 Copy_TPat_Fields(componentCopy->pattern, currentComponent->pattern);*/
  2749.             if(currentComponent->blob)
  2750.                 componentCopy->blob = (BLOB *)Copy_Object((OBJECT *)currentComponent->blob);
  2751.             componentCopy->next = NULL;
  2752.             /*END: copy component*/
  2753.           
  2754.             /*add in new component*/
  2755.             if(destComponent == NULL)/*This is the first copy*/
  2756.             {
  2757.                 New->Vals.Blob.blob_dat = componentCopy;
  2758.                 destComponent = New->Vals.Blob.blob_dat;
  2759.             }
  2760.             else
  2761.             {
  2762.                 destComponent->next = componentCopy;
  2763.                 destComponent = destComponent->next;
  2764.             }
  2765.           
  2766.             /*move to next component*/
  2767.             currentComponent = nextComponent;
  2768.         }
  2769.     }
  2770.   }
  2771. #endif
  2772.   if (Old->Type == DENSITY_FILE_PATTERN)
  2773.   {
  2774.      New->Vals.Density_File = Copy_Density_File(Old->Vals.Density_File);
  2775.   }
  2776.  
  2777. #ifdef CracklePatch
  2778.   if (Old->Type == CRACKLE_PATTERN)
  2779.   {
  2780.     New->Vals.Crackle.cv =(VECTOR*) POV_MALLOC( 125*sizeof(VECTOR), "crackle cache");
  2781.   }
  2782.  
  2783. #endif
  2784. #ifdef PigmentPatternPatch
  2785.   if (Old->Type == PIGMENT_PATTERN )
  2786.   {
  2787.     New->Vals.Pigment = Copy_Pigment(Old->Vals.Pigment);
  2788.   }
  2789. #endif
  2790. }
  2791.  
  2792.  
  2793.  
  2794. /*****************************************************************************
  2795. *
  2796. * FUNCTION
  2797. *
  2798. * INPUT
  2799. *
  2800. * OUTPUT
  2801. *
  2802. * RETURNS
  2803. *
  2804. * AUTHOR
  2805. *
  2806. * DESCRIPTION
  2807. *
  2808. * CHANGES
  2809. *
  2810. ******************************************************************************/
  2811.  
  2812. void Destroy_TPat_Fields(TPATTERN *Tpat)
  2813. {
  2814.   Destroy_Warps(Tpat->Warps);
  2815.   Destroy_Blend_Map(Tpat->Blend_Map);
  2816.   /* Note, cannot destroy Tpat->Next nor pattern itself because we don't
  2817.      know what kind of thing this is.  It must be destroied by Destroy_Pigment, etc.
  2818.   */
  2819.  
  2820.   /* NK 1998 - added IMAGE_PATTERN */
  2821.   if ((Tpat->Type == BITMAP_PATTERN))/* || (Tpat->Type == IMAGE_PATTERN))*/
  2822.   {
  2823.      Destroy_Image(Tpat->Vals.Image);
  2824.   }
  2825.   
  2826. #ifdef SplineWavePatch
  2827.   if(Tpat->spline_wave != NULL)
  2828.       Destroy_Spline(Tpat->spline_wave);
  2829. #endif
  2830.   if (Tpat->Type == DENSITY_FILE_PATTERN)
  2831.   {
  2832.      Destroy_Density_File(Tpat->Vals.Density_File);
  2833.   }
  2834. #ifdef ObjectPatternPatch
  2835.   if (Tpat->Type == OBJECT_PATTERN)
  2836.   {
  2837.     if(Tpat->Vals.Object != NULL)
  2838.     {
  2839.         Destroy_Object((OBJECT *)Tpat->Vals.Object);
  2840.     }
  2841.   }
  2842. #endif
  2843. #ifdef ProximityPatch
  2844.   if (Tpat->Type == PROXIMITY_PATTERN)
  2845.   {
  2846.     if(Tpat->Vals.Proximity.proxObject != NULL)
  2847.     {
  2848.         Destroy_Object((OBJECT *)Tpat->Vals.Proximity.proxObject);
  2849.     }
  2850.   }
  2851. #endif
  2852. #ifdef BlobPatternPatch
  2853.   if (Tpat->Type == BLOB_PATTERN ||
  2854.       Tpat->Type == BLOB_PIGMENT)
  2855.   {
  2856.     if(Tpat->Vals.Blob.blob_dat != NULL)/*Chris Huff blob pattern*/
  2857.     {
  2858.         BLOB_PATTERN_DATA * currentComponent = Tpat->Vals.Blob.blob_dat;
  2859.         BLOB_PATTERN_DATA * nextComponent = NULL;/*tPat->blob_dat->next;*/
  2860.         while(currentComponent != NULL)
  2861.         {
  2862.           nextComponent = currentComponent->next;
  2863.         
  2864.           if(currentComponent->pigment != NULL)
  2865.               Destroy_Pigment(currentComponent->pigment);
  2866. /*          if(currentComponent->pattern != NULL)
  2867.               Destroy_TPat_Fields(currentComponent->pattern);*/
  2868.           if(currentComponent->blob != NULL)
  2869.               Destroy_Object((OBJECT *)currentComponent->blob);
  2870.           POV_FREE(currentComponent->transform);
  2871.           POV_FREE(currentComponent);
  2872.         
  2873.           currentComponent = nextComponent;
  2874.         }
  2875.     }
  2876.   }
  2877. #endif
  2878. #ifdef CracklePatch
  2879.   if ( Tpat->Type == CRACKLE_PATTERN )
  2880.   {
  2881.     if ( Tpat->Vals.Crackle.cv )
  2882.     POV_FREE( Tpat->Vals.Crackle.cv );
  2883.     Tpat->Vals.Crackle.cv = NULL;
  2884.   }
  2885. #endif
  2886. #ifdef PigmentPatternPatch
  2887.   if (Tpat->Type == PIGMENT_PATTERN )
  2888.   {
  2889.     Destroy_Pigment(Tpat->Vals.Pigment);
  2890.     Tpat->Vals.Pigment = NULL;
  2891.   }
  2892. #endif
  2893. }
  2894.  
  2895.  
  2896.  
  2897.  
  2898. /*****************************************************************************
  2899. *
  2900. * FUNCTION
  2901. *
  2902. * INPUT
  2903. *
  2904. * OUTPUT
  2905. *
  2906. * RETURNS
  2907. *
  2908. * AUTHOR
  2909. *
  2910. * DESCRIPTION
  2911. *
  2912. * CHANGES
  2913. *
  2914. ******************************************************************************/
  2915.  
  2916. TURB *Create_Turb()
  2917. {
  2918.   TURB *New;
  2919.  
  2920.   New = (TURB *)POV_MALLOC(sizeof(TURB),"turbulence struct");
  2921.  
  2922.   Make_Vector(New->Turbulence, 0.0, 0.0, 0.0);
  2923.  
  2924.   New->Octaves = 6;
  2925.   New->Omega = 0.5;
  2926.   New->Lambda = 2.0;
  2927.  
  2928.   return(New);
  2929. }
  2930.  
  2931.  
  2932.  
  2933. /*****************************************************************************
  2934. *
  2935. * FUNCTION
  2936. *
  2937. * INPUT
  2938. *
  2939. * OUTPUT
  2940. *
  2941. * RETURNS
  2942. *
  2943. * AUTHOR
  2944. *
  2945. * DESCRIPTION
  2946. *
  2947. * CHANGES
  2948. *
  2949. ******************************************************************************/
  2950.  
  2951. #if 0   /* Unused function [AED] */
  2952. static TURB *Copy_Turb(TURB *Old)
  2953. {
  2954.   TURB *New;
  2955.  
  2956.   if (Old != NULL)
  2957.   {
  2958.     New = Create_Turb();
  2959.  
  2960.     *New = *Old;
  2961.   }
  2962.   else
  2963.   {
  2964.     New=NULL;
  2965.   }
  2966.  
  2967.   return(New);
  2968. }
  2969. #endif
  2970.  
  2971.  
  2972. /*****************************************************************************
  2973. *
  2974. * FUNCTION
  2975. *
  2976. *   Translate_Tpattern
  2977. *
  2978. * INPUT
  2979. *
  2980. * OUTPUT
  2981. *
  2982. * RETURNS
  2983. *
  2984. * AUTHOR
  2985. *
  2986. *   POV-Ray Team
  2987. *
  2988. * DESCRIPTION
  2989. *
  2990. * CHANGES
  2991. *
  2992. ******************************************************************************/
  2993.  
  2994. void Translate_Tpattern(TPATTERN *Tpattern,VECTOR Vector)
  2995. {
  2996.   TRANSFORM Trans;
  2997.  
  2998.   if (Tpattern != NULL)
  2999.   {
  3000.     Compute_Translation_Transform (&Trans, Vector);
  3001.  
  3002.     Transform_Tpattern (Tpattern, &Trans);
  3003.   }
  3004. }
  3005.  
  3006.  
  3007.  
  3008. /*****************************************************************************
  3009. *
  3010. * FUNCTION
  3011. *
  3012. *   Rotate_Tpattern
  3013. *
  3014. * INPUT
  3015. *
  3016. * OUTPUT
  3017. *
  3018. * RETURNS
  3019. *
  3020. * AUTHOR
  3021. *
  3022. *   POV-Ray Team
  3023. *
  3024. * DESCRIPTION
  3025. *
  3026. * CHANGES
  3027. *
  3028. ******************************************************************************/
  3029.  
  3030. void Rotate_Tpattern(TPATTERN *Tpattern,VECTOR Vector)
  3031. {
  3032.   TRANSFORM Trans;
  3033.  
  3034.   if (Tpattern != NULL)
  3035.   {
  3036.     Compute_Rotation_Transform (&Trans, Vector);
  3037.  
  3038.     Transform_Tpattern (Tpattern, &Trans);
  3039.   }
  3040. }
  3041.  
  3042.  
  3043.  
  3044. /*****************************************************************************
  3045. *
  3046. * FUNCTION
  3047. *
  3048. *   Scale_Tpattern
  3049. *
  3050. * INPUT
  3051. *
  3052. * OUTPUT
  3053. *
  3054. * RETURNS
  3055. *
  3056. * AUTHOR
  3057. *
  3058. *   POV-Ray Team
  3059. *
  3060. * DESCRIPTION
  3061. *
  3062. * CHANGES
  3063. *
  3064. ******************************************************************************/
  3065.  
  3066. void Scale_Tpattern(TPATTERN *Tpattern,VECTOR Vector)
  3067. {
  3068.   TRANSFORM Trans;
  3069.  
  3070.   if (Tpattern != NULL)
  3071.   {
  3072.     Compute_Scaling_Transform (&Trans, Vector);
  3073.  
  3074.     Transform_Tpattern (Tpattern, &Trans);
  3075.   }
  3076. }
  3077.  
  3078.  
  3079.  
  3080. /*****************************************************************************
  3081. *
  3082. * FUNCTION
  3083. *
  3084. *   Transform_Tpattern
  3085. *
  3086. * INPUT
  3087. *
  3088. * OUTPUT
  3089. *
  3090. * RETURNS
  3091. *
  3092. * AUTHOR
  3093. *
  3094. *   POV-Ray Team
  3095. *
  3096. * DESCRIPTION
  3097. *
  3098. * CHANGES
  3099. *
  3100. ******************************************************************************/
  3101.  
  3102. void Transform_Tpattern(TPATTERN *Tpattern,TRANSFORM *Trans)
  3103. {
  3104.   WARP *Temp;
  3105.  
  3106.   if (Tpattern != NULL)
  3107.   {
  3108.     if (Tpattern->Warps == NULL)
  3109.     {
  3110.       Tpattern->Warps=Create_Warp(TRANSFORM_WARP);
  3111.     }
  3112.     else
  3113.     {
  3114.       if (Tpattern->Warps->Warp_Type != TRANSFORM_WARP)
  3115.       {
  3116.         Temp=Tpattern->Warps;
  3117.  
  3118.         Tpattern->Warps=Create_Warp(TRANSFORM_WARP);
  3119.  
  3120.         Tpattern->Warps->Next_Warp=Temp;
  3121.       }
  3122.     }
  3123.  
  3124.     Compose_Transforms (&( ((TRANS *)(Tpattern->Warps))->Trans), Trans);
  3125.   }
  3126. }
  3127.  
  3128.  
  3129. /*****************************************************************************
  3130. *
  3131. * FUNCTION
  3132. *
  3133. *   ripples_pigm
  3134. *
  3135. * INPUT
  3136. *
  3137. *   EPoint -- The point in 3d space at which the pattern
  3138. *   is evaluated.
  3139. *   TPat   -- Texture pattern struct
  3140. *
  3141. * OUTPUT
  3142. *   
  3143. * RETURNS
  3144. *
  3145. *   DBL value in the range 0.0 to 1.0
  3146. *   
  3147. * AUTHOR
  3148. *
  3149. *   POV-Ray Team
  3150. *
  3151. * DESCRIPTION   : Note this pattern is only used for pigments and textures.
  3152. *                 Normals have a specialized pattern for this.
  3153. *
  3154. * CHANGES
  3155. *   Nov 1994 : adapted from normal by [CY]
  3156. *
  3157. ******************************************************************************/
  3158.  
  3159. static DBL ripples_pigm (VECTOR EPoint, TPATTERN *TPat)
  3160. {
  3161.   register unsigned int i;
  3162.   register DBL length, index;
  3163.   DBL scalar =0.0;
  3164.   VECTOR point;
  3165.  
  3166.   for (i = 0 ; i < Number_Of_Waves ; i++)
  3167.   {
  3168.     VSub (point, EPoint, Wave_Sources[i]);
  3169.     VLength (length, point);
  3170.  
  3171.     if (length == 0.0)
  3172.       length = 1.0;
  3173.  
  3174.     index = length * TPat->Frequency + TPat->Phase;
  3175.  
  3176.     scalar += cycloidal(index);
  3177.   }
  3178.  
  3179.   scalar = 0.5*(1.0+(scalar / (DBL)Number_Of_Waves));
  3180.  
  3181.   return(scalar);
  3182. }
  3183.  
  3184.  
  3185. /*****************************************************************************
  3186. *
  3187. * FUNCTION
  3188. *
  3189. *   waves_pigm
  3190. *
  3191. * INPUT
  3192. *
  3193. *   EPoint -- The point in 3d space at which the pattern
  3194. *   is evaluated.
  3195. *   TPat   -- Texture pattern struct
  3196. *
  3197. * OUTPUT
  3198. *
  3199. * RETURNS
  3200. *
  3201. *   DBL value in the range 0.0 to 1.0
  3202. *
  3203. * AUTHOR
  3204. *
  3205. *   POV-Ray Team
  3206. *
  3207. * DESCRIPTION   : Note this pattern is only used for pigments and textures.
  3208. *                 Normals have a specialized pattern for this.
  3209. *
  3210. * CHANGES
  3211. *   Nov 1994 : adapted from normal by [CY]
  3212. *
  3213. ******************************************************************************/
  3214.  
  3215. static DBL waves_pigm (VECTOR EPoint, TPATTERN *TPat)
  3216. {
  3217.   register unsigned int i;
  3218.   register DBL length, index;
  3219.   DBL scalar = 0.0;
  3220.   VECTOR point;
  3221.  
  3222.   for (i = 0 ; i < Number_Of_Waves ; i++)
  3223.   {
  3224.     VSub (point, EPoint, Wave_Sources[i]);
  3225.     VLength (length, point);
  3226.  
  3227.     if (length == 0.0)
  3228.     {
  3229.       length = 1.0;
  3230.     }
  3231.  
  3232.     index = length * TPat->Frequency * frequency[i] + TPat->Phase;
  3233.  
  3234.     scalar += cycloidal(index)/frequency[i];
  3235.   }
  3236.  
  3237.   scalar = 0.2*(2.5+(scalar / (DBL)Number_Of_Waves));
  3238.  
  3239.   return(scalar);
  3240. }
  3241.  
  3242.  
  3243. /*****************************************************************************
  3244. *
  3245. * FUNCTION
  3246. *
  3247. *   dents_pigm
  3248. *
  3249. * INPUT
  3250. *
  3251. *   EPoint -- The point in 3d space at which the pattern
  3252. *   is evaluated.
  3253. *
  3254. * OUTPUT
  3255. *
  3256. * RETURNS
  3257. *
  3258. *   DBL value in the range 0.0 to 1.0
  3259. *   
  3260. * AUTHOR
  3261. *
  3262. *   POV-Ray Team
  3263. *   
  3264. * DESCRIPTION   : Note this pattern is only used for pigments and textures.
  3265. *                 Normals have a specialized pattern for this.
  3266. *
  3267. * CHANGES
  3268. *   Nov 1994 : adapted from normal by [CY]
  3269. *
  3270. ******************************************************************************/
  3271.  
  3272. static DBL dents_pigm (VECTOR EPoint)
  3273. {
  3274.   DBL noise;
  3275.  
  3276.   noise = Noise (EPoint);
  3277.  
  3278.   return(noise * noise * noise);
  3279. }
  3280.  
  3281.  
  3282.  
  3283. /*****************************************************************************
  3284. *
  3285. * FUNCTION
  3286. *
  3287. *   wrinkles_pigm
  3288. *
  3289. * INPUT
  3290. *
  3291. *   EPoint -- The point in 3d space at which the pattern
  3292. *   is evaluated.
  3293. *   
  3294. * OUTPUT
  3295. *   
  3296. * RETURNS
  3297. *
  3298. *   DBL value in the range 0.0 to 1.0
  3299. *   
  3300. * AUTHOR
  3301. *
  3302. *   POV-Ray Team
  3303. *   
  3304. * DESCRIPTION   : Note this pattern is only used for pigments and textures.
  3305. *                 Normals have a specialized pattern for this.
  3306. *
  3307. * CHANGES
  3308. *   Nov 1994 : adapted from normal by [CY]
  3309. *
  3310. ******************************************************************************/
  3311.  
  3312.  
  3313. static DBL wrinkles_pigm (VECTOR EPoint)
  3314. {
  3315.   register int i;
  3316.   DBL lambda = 2.0;
  3317.   DBL omega = 0.5;
  3318.   DBL value;
  3319.   VECTOR temp;
  3320.  
  3321.   value = Noise(EPoint);
  3322.  
  3323.   for (i = 1; i < 10; i++)
  3324.   {
  3325.     VScale(temp,EPoint,lambda);
  3326.  
  3327.     value += omega * Noise(temp);
  3328.  
  3329.     lambda *= 2.0;
  3330.  
  3331.     omega *= 0.5;
  3332.   }
  3333.  
  3334.   return(value/2.0);
  3335. }
  3336.  
  3337.  
  3338.  
  3339. /*****************************************************************************
  3340. *
  3341. * FUNCTION
  3342. *
  3343. *   quilted_pigm
  3344. *
  3345. * INPUT
  3346. *   
  3347. * OUTPUT
  3348. *   
  3349. * RETURNS
  3350. *   
  3351. * AUTHOR
  3352. *
  3353. *   Dan Farmer & Chris Young
  3354. *   
  3355. * DESCRIPTION
  3356. *
  3357. * CHANGES
  3358. *
  3359. ******************************************************************************/
  3360.  
  3361. static DBL quilted_pigm (VECTOR EPoint, TPATTERN *TPat)
  3362. {
  3363.   VECTOR value;
  3364.   DBL t;
  3365.  
  3366.   value[X] = EPoint[X]-FLOOR(EPoint[X])-0.5;
  3367.   value[Y] = EPoint[Y]-FLOOR(EPoint[Y])-0.5;
  3368.   value[Z] = EPoint[Z]-FLOOR(EPoint[Z])-0.5;
  3369.  
  3370.   t = sqrt(value[X]*value[X]+value[Y]*value[Y]+value[Z]*value[Z]);
  3371.  
  3372.   t = quilt_cubic(t, TPat->Vals.Quilted.Control0, TPat->Vals.Quilted.Control1);
  3373.  
  3374.   value[X] *= t;
  3375.   value[Y] *= t;
  3376.   value[Z] *= t;
  3377.  
  3378.   return((fabs(value[X])+fabs(value[Y])+fabs(value[Z]))/3.0);
  3379. }
  3380.  
  3381.  
  3382.  
  3383. /*****************************************************************************
  3384. *
  3385. * FUNCTION
  3386. *
  3387. * INPUT
  3388. *
  3389. * OUTPUT
  3390. *
  3391. * RETURNS
  3392. *
  3393. * AUTHOR
  3394. *
  3395. * DESCRIPTION
  3396. *
  3397. * CHANGES
  3398. *
  3399. ******************************************************************************/
  3400.  
  3401. #define INV_SQRT_3_4 1.154700538
  3402. DBL quilt_cubic(DBL t,DBL p1,DBL p2)
  3403. {
  3404.  DBL it=(1-t);
  3405.  DBL itsqrd=it*it;
  3406. /* DBL itcubed=it*itsqrd; */
  3407.  DBL tsqrd=t*t;
  3408.  DBL tcubed=t*tsqrd;
  3409.  DBL val;
  3410.  
  3411. /* Originally coded as...
  3412.  
  3413.  val= (DBL)(itcubed*n1+(tcubed)*n2+3*t*(itsqrd)*p1+3*(tsqrd)*(it)*p2);
  3414.  
  3415.  re-written by CEY to optimise because n1=0 n2=1 always.
  3416.  
  3417. */
  3418.  
  3419.  val = (tcubed + 3.0*t*itsqrd*p1 + 3.0*tsqrd*it*p2) * INV_SQRT_3_4;
  3420.  
  3421.  return(val);
  3422. }
  3423.  
  3424.  
  3425.  
  3426. /*****************************************************************************
  3427. *
  3428. * FUNCTION
  3429. *
  3430. * INPUT
  3431. *
  3432. * OUTPUT
  3433. *
  3434. * RETURNS
  3435. *
  3436. * AUTHOR
  3437. *
  3438. * DESCRIPTION
  3439. *
  3440. * CHANGES
  3441. *
  3442. ******************************************************************************/
  3443.  
  3444. void Search_Blend_Map (DBL value,BLEND_MAP *Blend_Map,BLEND_MAP_ENTRY **Prev,BLEND_MAP_ENTRY  **Cur)
  3445. {
  3446.   BLEND_MAP_ENTRY *P, *C;
  3447.   int Max_Ent=Blend_Map->Number_Of_Entries-1;
  3448.  
  3449.   /* if greater than last, use last. */
  3450.  
  3451.   if (value >= Blend_Map->Blend_Map_Entries[Max_Ent].value)
  3452.   {
  3453.     P = C = &(Blend_Map->Blend_Map_Entries[Max_Ent]);
  3454.   }
  3455.   else
  3456.   {
  3457.     P = C = &(Blend_Map->Blend_Map_Entries[0]);
  3458.  
  3459.     while (value > C->value)
  3460.     {
  3461.       P = C++;
  3462.     }
  3463.   }
  3464.  
  3465.   if (value == C->value)
  3466.   {
  3467.     P = C;
  3468.   }
  3469.  
  3470.   *Prev = P;
  3471.   *Cur  = C;
  3472. }
  3473.  
  3474.  
  3475.  
  3476. /*****************************************************************************
  3477. *
  3478. * FUNCTION
  3479. *
  3480. * INPUT
  3481. *
  3482. * OUTPUT
  3483. *
  3484. * RETURNS
  3485. *
  3486. * AUTHOR
  3487. *
  3488. * DESCRIPTION
  3489. *
  3490. * CHANGES
  3491. *
  3492. ******************************************************************************/
  3493.  
  3494. static TURB *Search_For_Turb(WARP *Warps)
  3495. {
  3496.   WARP* Temp=Warps;
  3497.  
  3498.   if (Temp!=NULL)
  3499.   {
  3500.     while (Temp->Next_Warp != NULL)
  3501.     {
  3502.       Temp=Temp->Next_Warp;
  3503.     }
  3504.  
  3505.     if (Temp->Warp_Type != CLASSIC_TURB_WARP)
  3506.     {
  3507.        Temp=NULL;
  3508.     }
  3509.   }
  3510.  
  3511.   return ((TURB *)Temp);
  3512. }
  3513.  
  3514.  
  3515. /*****************************************************************************
  3516. *
  3517. * FUNCTION
  3518. *
  3519. *   density_file
  3520. *
  3521. * INPUT
  3522. *
  3523. * OUTPUT
  3524. *
  3525. * RETURNS
  3526. *
  3527. * AUTHOR
  3528. *
  3529. *   Dieter Bayer
  3530. *
  3531. * DESCRIPTION
  3532. *
  3533. * CHANGES
  3534. *
  3535. *   Dec 1996 : Creation.
  3536. *
  3537. ******************************************************************************/
  3538.  
  3539. static DBL density_file(VECTOR EPoint, TPATTERN *TPat)
  3540. {
  3541.   int x, y, z;
  3542.   int x1, y1, z1;
  3543.   int x2, y2, z2;
  3544.   DBL xx, yy, zz;
  3545.   DBL xi, yi, zi;
  3546.   DBL f111, f112, f121, f122, f211, f212, f221, f222;
  3547.   DBL density = 0.0;
  3548.   DENSITY_FILE_DATA *Data;
  3549.  
  3550.   if ((TPat->Vals.Density_File != NULL) &&
  3551.       ((Data = TPat->Vals.Density_File->Data) != NULL))
  3552.   {
  3553.     if ((EPoint[X] >= 0.0) && (EPoint[X] < 1.0) &&
  3554.         (EPoint[Y] >= 0.0) && (EPoint[Y] < 1.0) &&
  3555.         (EPoint[Z] >= 0.0) && (EPoint[Z] < 1.0))
  3556.     {
  3557.       switch (TPat->Vals.Density_File->Interpolation)
  3558.       {
  3559.         case NO_INTERPOLATION:
  3560.  
  3561.           x = (int)(EPoint[X] * (DBL)Data->Sx);
  3562.           y = (int)((EPoint[Y] )* (DBL)Data->Sy); 
  3563.           z = (int)(EPoint[Z] * (DBL)Data->Sz);
  3564.  
  3565.           if ((x < 0) || (x >= Data->Sx) ||
  3566.               (y < 0) || (y >= Data->Sy) ||
  3567.               (z < 0) || (z >= Data->Sz))
  3568.           {
  3569.             density = 0.0;
  3570.           }
  3571.           else
  3572.           {
  3573.             density = (DBL)Data->Density[z][y][x] / 255.0;
  3574.           }
  3575.  
  3576.           break;
  3577.  
  3578.         case TRILINEAR_INTERPOLATION:
  3579.  
  3580.           xx = EPoint[X] * (DBL)(Data->Sx - 1);
  3581.           yy = (EPoint[Y]) * (DBL)(Data->Sy - 1);
  3582.           zz = EPoint[Z] * (DBL)(Data->Sz - 1);
  3583.  
  3584.           x1 = (int)xx;
  3585.           y1 = (int)yy;
  3586.           z1 = (int)zz;
  3587.  
  3588.           x2 = x1 + 1;
  3589.           y2 = y1 + 1;
  3590.           z2 = z1 + 1;
  3591.  
  3592.           xx -= floor(xx);
  3593.           yy -= floor(yy);
  3594.           zz -= floor(zz);
  3595.  
  3596.           xi = 1.0 - xx;
  3597.           yi = 1.0 - yy;
  3598.           zi = 1.0 - zz;
  3599.  
  3600.           f111 = (DBL)Data->Density[z1][y1][x1] / 255.0;
  3601.           f112 = (DBL)Data->Density[z1][y1][x2] / 255.0;
  3602.           f121 = (DBL)Data->Density[z1][y2][x1] / 255.0;
  3603.           f122 = (DBL)Data->Density[z1][y2][x2] / 255.0;
  3604.           f211 = (DBL)Data->Density[z2][y1][x1] / 255.0;
  3605.           f212 = (DBL)Data->Density[z2][y1][x2] / 255.0;
  3606.           f221 = (DBL)Data->Density[z2][y2][x1] / 255.0;
  3607.           f222 = (DBL)Data->Density[z2][y2][x2] / 255.0;
  3608.  
  3609.           density = f111 * zi * yi * xi +
  3610.                     f112 * zi * yi * xx +
  3611.                     f121 * zi * yy * xi +
  3612.                     f122 * zi * yy * xx +
  3613.                     f211 * zz * yi * xi +
  3614.                     f212 * zz * yi * xx +
  3615.                     f221 * zz * yy * xi +
  3616.                     f222 * zz * yy * xx;
  3617.  
  3618.           break;
  3619.       }
  3620.     }
  3621.     else
  3622.     {
  3623.       density = 0.0;
  3624.     }
  3625.  
  3626. /*
  3627.     fprintf(stderr, "x = %3d, y = %3d, z = %3d, density = %5.4f\n", x, y, z, density);
  3628. */
  3629.   }
  3630.   return(density);
  3631. }
  3632.  
  3633.  
  3634. /*****************************************************************************
  3635. *
  3636. * FUNCTION
  3637. *
  3638. *   Create_Density_File
  3639. *
  3640. * INPUT
  3641. *
  3642. * OUTPUT
  3643. *
  3644. * RETURNS
  3645. *
  3646. * AUTHOR
  3647. *
  3648. *   Dieter Bayer
  3649. *
  3650. * DESCRIPTION
  3651. *
  3652. *   Create a density file structure.
  3653. *
  3654. * CHANGES
  3655. *
  3656. *   Dec 1996 : Creation.
  3657. *
  3658. ******************************************************************************/
  3659.  
  3660. DENSITY_FILE *Create_Density_File()
  3661. {
  3662.   DENSITY_FILE *New;
  3663.  
  3664.   New = (DENSITY_FILE *)POV_MALLOC(sizeof(DENSITY_FILE), "density file");
  3665.  
  3666.   New->Interpolation = NO_INTERPOLATION;
  3667.  
  3668.   New->Data = (DENSITY_FILE_DATA *)POV_MALLOC(sizeof(DENSITY_FILE_DATA), "density file data");
  3669.  
  3670.   New->Data->References = 1;
  3671.  
  3672.   New->Data->Name = NULL;
  3673.  
  3674.   New->Data->Sx =
  3675.   New->Data->Sy =
  3676.   New->Data->Sz = 0;
  3677.  
  3678.   New->Data->Density = NULL;
  3679.  
  3680.   return (New);
  3681. }
  3682.  
  3683.  
  3684.  
  3685. /*****************************************************************************
  3686. *
  3687. * FUNCTION
  3688. *
  3689. *   Copy_Density_File
  3690. *
  3691. * INPUT
  3692. *
  3693. * OUTPUT
  3694. *
  3695. * RETURNS
  3696. *
  3697. * AUTHOR
  3698. *
  3699. *   Dieter Bayer
  3700. *
  3701. * DESCRIPTION
  3702. *
  3703. *   Copy a density file structure.
  3704. *
  3705. * CHANGES
  3706. *
  3707. *   Dec 1996 : Creation.
  3708. *
  3709. ******************************************************************************/
  3710.  
  3711. DENSITY_FILE *Copy_Density_File(DENSITY_FILE *Old)
  3712. {
  3713.   DENSITY_FILE *New;
  3714.  
  3715.   if (Old != NULL)
  3716.   {
  3717.     New = (DENSITY_FILE *)POV_MALLOC(sizeof(DENSITY_FILE), "density file");
  3718.  
  3719.     *New = *Old;
  3720.  
  3721.     New->Data->References++;
  3722.   }
  3723.   else          /* tw */
  3724.     New = NULL; /* tw */
  3725.  
  3726.   return(New);
  3727. }
  3728.  
  3729.  
  3730.  
  3731. /*****************************************************************************
  3732. *
  3733. * FUNCTION
  3734. *
  3735. *   Destroy_Density_File
  3736. *
  3737. * INPUT
  3738. *
  3739. * OUTPUT
  3740. *
  3741. * RETURNS
  3742. *
  3743. * AUTHOR
  3744. *
  3745. *   Dieter Bayer
  3746. *
  3747. * DESCRIPTION
  3748. *
  3749. *   Destroy a density file structure.
  3750. *
  3751. * CHANGES
  3752. *
  3753. *   Dec 1996 : Creation.
  3754. *
  3755. ******************************************************************************/
  3756.  
  3757. void Destroy_Density_File(DENSITY_FILE *Density_File)
  3758. {
  3759.     int y, z;
  3760.     if (Density_File != NULL)
  3761.     {
  3762.         if ((--(Density_File->Data->References)) == 0)
  3763.         {
  3764.             POV_FREE(Density_File->Data->Name);
  3765.             for (z = 0; z < Density_File->Data->Sz; z++)
  3766.             {
  3767.                 for (y = 0; y < Density_File->Data->Sy; y++)
  3768.                     POV_FREE(Density_File->Data->Density[z][y]);
  3769.                 POV_FREE(Density_File->Data->Density[z]);
  3770.             }
  3771.             if(Density_File->Data->Density!= NULL)  /*YS aug 2000 bug fix*/
  3772.                 POV_FREE(Density_File->Data->Density);
  3773.             POV_FREE(Density_File->Data);
  3774.         }
  3775.         POV_FREE(Density_File);
  3776.     }
  3777. }
  3778.  
  3779.  
  3780. void Read_Density_File(DENSITY_FILE *df)
  3781. {
  3782.   unsigned char ***map;
  3783.   int y, z, sx, sy, sz;
  3784.   FILE *file;
  3785.  
  3786.   if (df == NULL)
  3787.   {
  3788.      return;
  3789.   }
  3790.   
  3791.   /* Allocate and read density file. */
  3792.  
  3793.   if ((df != NULL) && (df->Data->Name != NULL))
  3794.   {
  3795.     if ((file = Locate_File(df->Data->Name, READ_BINFILE_STRING, ".df3", ".DF3",NULL,TRUE)) == NULL)
  3796.     {
  3797.       Error("Cannot read media density file.\n");
  3798.     }
  3799.     
  3800.     sx = df->Data->Sx = readushort(file);
  3801.     sy = df->Data->Sy = readushort(file);
  3802.     sz = df->Data->Sz = readushort(file);
  3803.  
  3804.     map = (unsigned char ***)POV_MALLOC(sz*sizeof(unsigned char **), "media density file data");
  3805.  
  3806.     for (z = 0; z < sz; z++)
  3807.     {
  3808.       map[z] = (unsigned char **)POV_MALLOC(sy*sizeof(unsigned char *), "media density file data");
  3809.  
  3810.       for (y = 0; y < sy; y++)
  3811.       {
  3812.         map[z][y] = (unsigned char *)POV_MALLOC(sx*sizeof(unsigned char), "media density file data");
  3813.  
  3814.         fread(map[z][y], sizeof(unsigned char), (size_t)sx, file);
  3815.       }
  3816.     }
  3817.  
  3818.     df->Data->Density = map;
  3819.  
  3820.     if (file != NULL)        /* -hdf99- */
  3821.     {
  3822.       fclose(file);
  3823.     }
  3824.   }
  3825. }
  3826.  
  3827. #ifdef SolidPatternPatch
  3828.  /*Chris Huff 7/20/99 solid pattern*/
  3829. static DBL SolidPat(VECTOR EPoint, TPATTERN *TPat)
  3830. {
  3831.     return TPat->Vals.SolidVal;
  3832. }
  3833. #endif
  3834. #ifdef ClothPatternPatch
  3835. /*Chris Huff cloth pattern*/
  3836. static DBL ClothPat(VECTOR EPoint)
  3837. {
  3838.     DBL xVal = EPoint[X];
  3839.     DBL zVal = EPoint[Z];
  3840.     xVal = ((xVal > 1.0) ? fmod(xVal, 1.0) : xVal);
  3841.     zVal = ((zVal > 1.0) ? fmod(zVal, 1.0) : zVal);
  3842. /*    DBL xVal = ((EPoint[X] > 1.0) ? fmod(EPoint[X], 1.0) : EPoint[X]);
  3843.     DBL zVal = ((EPoint[Z] > 1.0) ? fmod(EPoint[Z], 1.0) : EPoint[Z]);*/
  3844.  
  3845.     if(checker(EPoint))
  3846.         return fabs(xVal);/*fabs(cycloidal(xVal/2));*/
  3847.     else
  3848.         return fabs(zVal);/*fabs(cycloidal(zVal/2));*/
  3849. }
  3850.  
  3851. /*Chris Huff cloth2 pattern*/
  3852. static DBL Cloth2Pat(VECTOR EPoint)
  3853. {
  3854.     DBL xVal = fabs(EPoint[X]);
  3855.     DBL zVal = fabs(EPoint[Z]);
  3856.     xVal = ((xVal > 1.0) ? fmod(xVal, 1.0) : xVal);
  3857.     zVal = ((zVal > 1.0) ? fmod(zVal, 1.0) : zVal);
  3858.  
  3859.     if(checker(EPoint))
  3860.         return (fabs(xVal)/2);/*fabs(cycloidal(xVal/2));*/
  3861.     else
  3862.         return ((fabs(zVal)/2)+0.5);/*fabs(cycloidal(zVal/2));*/
  3863. }
  3864. #endif
  3865. #ifdef TorodialPatch
  3866. /*Chris Huff torodial pattern*/
  3867. static DBL ToroidalSpiral(VECTOR EPoint, TPATTERN *TPat)
  3868. {
  3869.     DBL x = EPoint[X];
  3870.     DBL y = EPoint[Y];
  3871.     DBL z = EPoint[Z];
  3872.     DBL longitude = atan2(x, y);
  3873.     DBL latitude = atan2(sqrt(x*x + y*y)-1.0, z);
  3874.  
  3875.     return fmod(((latitude + longitude*(TPat->Vals.SolidVal))/2*M_PI), 1.0);
  3876. }
  3877. #endif
  3878. static unsigned short readushort(FILE *infile)
  3879. {
  3880.   int i0, i1 = 0; /* To quiet warnings */
  3881.  
  3882.   if ((i0  = fgetc(infile)) == EOF || (i1  = fgetc(infile)) == EOF)
  3883.   {
  3884.     Error("Error reading density_file\n");
  3885.   }
  3886.  
  3887.   return (unsigned short)((((unsigned short)i0) << 8) | ((unsigned short)i1));
  3888. }
  3889.  
  3890. /** poviso: May 97 R.S. **/
  3891. #ifdef POVISO
  3892. /*****************************************************************************
  3893. *
  3894. * FUNCTION
  3895. *
  3896. * INPUT
  3897. *
  3898. * OUTPUT
  3899. *
  3900. * RETURNS
  3901. *
  3902. * AUTHOR  R.Suzuki
  3903. *
  3904. * DESCRIPTION
  3905. *
  3906. * CHANGES    May 97 
  3907. *
  3908. ******************************************************************************/
  3909.  
  3910.  
  3911.  
  3912. static DBL func_pat (VECTOR EPoint, TPATTERN *TPat)
  3913. {
  3914.   DBL value;
  3915.   VECTOR V1;
  3916.   FUNCTION *TFunc;
  3917.  
  3918.   Assign_Vector(V1,EPoint);
  3919.   TFunc=(FUNCTION*)TPat->Vals.Function;
  3920.   if (TFunc!=NULL) value = (TFunc->iso_func)(TFunc, V1);
  3921.   value = ((value > 1.0) ? fmod(value, 1.0) : value);
  3922.   return(value);
  3923. }
  3924.  
  3925. #endif
  3926. /** --- **/
  3927.  
  3928. #ifdef PigmentPatternPatch
  3929. static DBL pigment_pattern (VECTOR EPoint, TPATTERN *TPat, INTERSECTION *Inter)
  3930. {
  3931.   DBL value;
  3932.   COLOUR Col;
  3933.   int colour_found=FALSE;
  3934.  
  3935.   if (TPat->Vals.Pigment)
  3936.   {
  3937.     colour_found = Compute_Pigment(Col,TPat->Vals.Pigment,EPoint,Inter);
  3938.   }
  3939.  
  3940.   if(!colour_found)
  3941.   {
  3942.     value = 0.0;
  3943.   }
  3944.   else
  3945.   {
  3946.     value = GREY_SCALE(Col);
  3947.   }
  3948.  
  3949.   return(value);
  3950. }
  3951. #endif
  3952.  
  3953. /*****************************************************************************
  3954. *
  3955. * FUNCTION
  3956. *
  3957. *   slope
  3958. *
  3959. * INPUT
  3960. *
  3961. *   EPoint -- The point in 3d space at which the pattern
  3962. *             is evaluated.
  3963. *   TPat   -- Texture pattern struct
  3964. *   Intersection - intersection struct
  3965. *
  3966. * OUTPUT
  3967. *
  3968. * RETURNS
  3969. *
  3970. *   DBL value in the range 0.0 to 1.0, 0.0 if normal is NULL
  3971. *
  3972. * AUTHOR
  3973. *
  3974. *   -hdf-
  3975. *
  3976. * DESCRIPTION   :
  3977. *
  3978. *   calculates the surface slope from surface normal vector
  3979. *
  3980. * CHANGES
  3981. *   Apr 1998 : written by H.-D. Fink
  3982. *   May 1998 : modified by M.C. Andrews - now combines slope and 'gradient'.
  3983. *
  3984. ******************************************************************************/
  3985.  
  3986. static DBL slope (VECTOR EPoint, TPATTERN *TPat, INTERSECTION *Intersection)
  3987. {
  3988.   DBL value, value1, value2;
  3989.  
  3990.   if (Intersection == NULL) return 0.0; /* just in case ... */
  3991.  
  3992.   if (TPat->Vals.Slope.Slope_Base > 0)
  3993.     /* short case 1: slope vector in x, y or z direction */
  3994.     value1 = Intersection->PNormal[TPat->Vals.Slope.Slope_Base - 1];
  3995.   else if (TPat->Vals.Slope.Slope_Base < 0)
  3996.     /* short case 2: slope vector in negative x, y or z direction */
  3997.     value1 = -Intersection->PNormal[-TPat->Vals.Slope.Slope_Base - 1];
  3998.   else
  3999.     /* projection slope onto normal vector */
  4000.     VDot(value1, Intersection->PNormal, TPat->Vals.Slope.Slope_Vector);
  4001.  
  4002.   /* Clamp to 1.0. */
  4003.   /* should never be necessary since both vectors are normalized */
  4004.   if      (value1 >  1.0) value1 =  1.0;
  4005.   else if (value1 < -1.0) value1 = -1.0;
  4006.  
  4007.   value1 = (value1 + 1.0) * 0.5;    /* normalize to [0..1] interval */
  4008.  
  4009.   if (!TPat->Vals.Slope.Altit_Len) return value1; /* no altitude defined */
  4010.  
  4011.   /* Calculate projection of Epoint along altitude vector */
  4012.   if (TPat->Vals.Slope.Altit_Base > 0)
  4013.     /* short case 1: altitude vector in x, y or z direction */
  4014.     value2 = EPoint[TPat->Vals.Slope.Altit_Base - 1];
  4015.   else if (TPat->Vals.Slope.Altit_Base < 0)
  4016.     /* short case 2: altitude vector in negative x, y or z direction */
  4017.     value2 = -EPoint[-TPat->Vals.Slope.Altit_Base - 1];
  4018.   else
  4019.     /* projection of Epoint along altitude vector */
  4020.     VDot(value2, EPoint, TPat->Vals.Slope.Altit_Vector);
  4021.  
  4022.   /* If set, use offset and scalings for slope and altitude. */
  4023.   if (0.0 != TPat->Vals.Slope.Slope_Mod[V])
  4024.   {
  4025.     value1 = (value1 - TPat->Vals.Slope.Slope_Mod[U]) / TPat->Vals.Slope.Slope_Mod[V];
  4026.   }
  4027.   if (0.0 != TPat->Vals.Slope.Altit_Mod[V])
  4028.   {
  4029.     value2 = (value2 - TPat->Vals.Slope.Altit_Mod[U]) / TPat->Vals.Slope.Altit_Mod[V];
  4030.   }
  4031.  
  4032.   value = TPat->Vals.Slope.Slope_Len * value1 + TPat->Vals.Slope.Altit_Len * value2;
  4033.  
  4034.   /* Clamp to 1.0. */
  4035.   value = (value < 0.0) ? 1.0 + fmod(value, 1.0) : fmod(value, 1.0);
  4036.   return value;
  4037.  
  4038. }
  4039.  
  4040. #ifdef CellsPatch
  4041. /*****************************************************************************
  4042. *
  4043. * FUNCTION
  4044. *
  4045. *   cells
  4046. *
  4047. * INPUT
  4048. *
  4049. *   EPoint -- The point in 3d space at which the pattern
  4050. *   is evaluated.
  4051. * OUTPUT
  4052. *
  4053. * RETURNS
  4054. *
  4055. *   DBL value in the range 0.0 to 1.0
  4056. *
  4057. * AUTHOR
  4058. *
  4059. *   John VanSickle
  4060. *
  4061. * DESCRIPTION
  4062. *
  4063. *   "cells":
  4064. *
  4065. *   New colour function by John VanSickle,
  4066. *     vansickl@erols.com
  4067. *
  4068. *   Assigns a pseudorandom value to each unit cube.  The value for the cube in
  4069. *   which the evaluted point lies is returned.
  4070. *
  4071. *   All "cells" specific source code and examples are in the public domain.
  4072. *
  4073. * CHANGES
  4074. *   Oct 1994    : original crackle code adapted by [CY]
  4075. *   Jul 1999    : adapted for a different pattern by [JV]
  4076. *
  4077. ******************************************************************************/
  4078.  
  4079. static DBL cells (VECTOR EPoint)
  4080. {
  4081.   int    temp,seed;
  4082.   DBL    tf;
  4083.  
  4084.   /* select a random value based on the cube from which this came. */
  4085.  
  4086. #ifdef NoiseTranslateFixPatch
  4087.   seed = Hash3d((int)EPoint[X]&0xFFF, (int)EPoint[Y]&0xFFF, (int)EPoint[Z]&0xFFF);
  4088. #else
  4089.   seed = Hash3d((int)EPoint[X], (int)EPoint[Y], (int)EPoint[Z]);
  4090. #endif
  4091.   temp = POV_GET_OLD_RAND(); /* save current seed */
  4092.   
  4093.   POV_SRAND(seed);
  4094.  
  4095.   tf = FRAND();
  4096.  
  4097.   POV_SRAND(temp);  /* restore */
  4098.  
  4099.   return min(tf,1.);
  4100. }
  4101. #endif
  4102.  
  4103. #ifdef VanSicklePatternPatch
  4104. /*****************************************************************************
  4105. *
  4106. * FUNCTION
  4107. *
  4108. *   blotches
  4109. *
  4110. * INPUT
  4111. *
  4112. *   EPoint -- The point in 3d space at which the pattern
  4113. *   is evaluated.
  4114. * OUTPUT
  4115. *
  4116. * RETURNS
  4117. *
  4118. *   DBL value in the range 0.0 to 1.0
  4119. *
  4120. * AUTHOR
  4121. *
  4122. *   John VanSickle, based on code by Jim McElhiney
  4123. *
  4124. * DESCRIPTION
  4125. *
  4126. *   "blotches":
  4127. *
  4128. *   New colour function by John VanSickle,
  4129. *     vansickl@erols.com
  4130. *
  4131. *   Assigns a pseudorandom point to each unit cube, and another float from 0 to
  4132. *   1 to each cube.  The nearest pseudorandom point to the evaluated point is
  4133. *   selected, and its associated pseudorandom value is returned.  It creates
  4134. *   regions of uniform value.
  4135. *
  4136. *   All "blotches" specific source code and examples are in the public domain.
  4137. *
  4138. * CHANGES
  4139. *   Oct 1994    : original crackle code adapted by [CY]
  4140. *   Jul 1999    : adapted for a different pattern by [JV]
  4141. *
  4142. ******************************************************************************/
  4143.  
  4144. static DBL blotches (VECTOR EPoint)
  4145. {
  4146.   int    i,j,temp,seed;
  4147.   long   thisseed;
  4148.   DBL    sum, minsum,  tf;
  4149.   VECTOR sv, tv, dv, t1, add;
  4150.  
  4151.   static int bvc;
  4152.   static long lastbeed = 0x80000000;
  4153.   static VECTOR bv[81];
  4154.  
  4155.   Assign_Vector(tv,EPoint);
  4156.  
  4157.   /*
  4158.    * Check to see if the input point is in the same unit cube as the last
  4159.    * call to this function, to use cache of cubelets for speed.
  4160.    */
  4161.  
  4162.   thisseed = PickInCube(tv, t1);
  4163.  
  4164.   if (thisseed != lastbeed)
  4165.   {
  4166.     /*
  4167.      * No, not same unit cube.  Calculate the random points for this new
  4168.      * cube and its 80 neighbours which differ in any axis by 1 or 2.
  4169.      * Why distance of 2?  If there is 1 point in each cube, located
  4170.      * randomly, it is possible for the closest random point to be in the
  4171.      * cube 2 over, or the one two over and one up.  It is NOT possible
  4172.      * for it to be two over and two up.  Picture a 3x3x3 cube with 9 more
  4173.      * cubes glued onto each face.
  4174.      */
  4175.  
  4176.     /* Now store a points for this cube and each of the 80 neighbour cubes. */
  4177.  
  4178.     bvc = 0;
  4179.  
  4180.     for (add[X] = -2.0; add[X] < 2.5; add[X] +=1.0)
  4181.     {
  4182.       for (add[Y] = -2.0; add[Y] < 2.5; add[Y] += 1.0)
  4183.       {
  4184.         for (add[Z] = -2.0; add[Z] < 2.5; add[Z] += 1.0)
  4185.         {
  4186.           /* For each cubelet in a 5x5 cube. */
  4187.  
  4188.           if ((fabs(add[X])>1.5)+(fabs(add[Y])>1.5)+(fabs(add[Z])>1.5) <= 1.0)
  4189.           {
  4190.             /* Yes, it's within a 3d knight move away. */
  4191.  
  4192.             VAdd(sv, tv, add);
  4193.  
  4194.             PickInCube(sv, t1);
  4195.  
  4196.             bv[bvc][X] = t1[X];
  4197.             bv[bvc][Y] = t1[Y];
  4198.             bv[bvc][Z] = t1[Z];
  4199.             bvc++;
  4200.           }
  4201.         }
  4202.       }
  4203.     }
  4204.  
  4205.     lastbeed = thisseed;
  4206.   }
  4207.  
  4208.   /*
  4209.    * Find the point with the shortest distance from the input point.
  4210.    * Loop invariant:  minsum is shortest dist
  4211.    */
  4212.  
  4213.   VSub(dv, bv[0], tv);  minsum  = VSumSqr(dv);
  4214.   j=0;
  4215.  
  4216.   /* Loop for the 81 cubelets to find closest. */
  4217.  
  4218.   for (i = 1; i < bvc; i++)
  4219.   {
  4220.     VSub(dv, bv[i], tv);
  4221.  
  4222.     sum = VSumSqr(dv);
  4223.  
  4224.     if (sum < minsum) {
  4225.       minsum = sum;
  4226.       j=i;
  4227.     }
  4228.   }
  4229.  
  4230.   /* select a random value based on the cube from which this came. */
  4231.  
  4232. #ifdef NoiseTranslateFixPatch
  4233.   seed = Hash3d((int)bv[j][X]&0xFFF, (int)bv[j][Y]&0xFFF, (int)bv[j][Z]&0xFFF);
  4234. #else
  4235.   seed = Hash3d((int)bv[j][X], (int)bv[j][Y], (int)bv[j][Z]);
  4236. #endif
  4237.   temp = POV_GET_OLD_RAND(); /* save current seed */
  4238.   
  4239.   POV_SRAND(seed);
  4240.  
  4241.   /* The first three values were used for the random point, and shouldn't
  4242.      be used for the value of the region */
  4243.  
  4244.   FRAND();
  4245.   FRAND();
  4246.   FRAND();
  4247.  
  4248.   tf = FRAND();
  4249.  
  4250.   POV_SRAND(temp);  /* restore */
  4251.  
  4252.   return min(tf,1.);
  4253. }
  4254.  
  4255. /*****************************************************************************
  4256. *
  4257. * FUNCTION
  4258. *
  4259. *   bands
  4260. *
  4261. * INPUT
  4262. *
  4263. *   EPoint -- The point in 3d space at which the pattern
  4264. *   is evaluated.
  4265. *   
  4266. * OUTPUT
  4267. *   
  4268. * RETURNS
  4269. *
  4270. *   DBL value in the range 0.0 to 1.0
  4271. *   
  4272. * AUTHOR
  4273. *
  4274. *   John VanSickle, based on the gradient code
  4275. *   
  4276. * DESCRIPTION
  4277. *
  4278. *   banded Pattern - just like gradient, but is not mirrored around the origin,
  4279. *   and the bands are of uniform width, regardless of the direction of the
  4280. *   vector.
  4281. *
  4282. * CHANGES
  4283. *   Oct 1994    : adapted from pigment by [CY]
  4284. *   Jul 1999    : adapted from gradient by [JV]
  4285. *
  4286. ******************************************************************************/
  4287.  
  4288. /* íMUY IMPORTANTE! This function expects the Vals.Gradient vector in TPat to
  4289.  * be normalized (ie, unit length).  If the vector is greater or less than
  4290.  * unit length, there will be a scaling effect on the pattern.  I could have
  4291.  * divided by the magnitude of the vector, but since this function will be
  4292.  * called more than once by the renderer, it is more appropriate for the
  4293.  * parser to normalize the vector (where it need be done only once). */
  4294.  
  4295. static DBL banded (VECTOR EPoint, TPATTERN *TPat)
  4296. {
  4297.   register DBL value;
  4298.  
  4299.   value = EPoint[X]*TPat->Vals.Gradient[X]
  4300.         + EPoint[Y]*TPat->Vals.Gradient[Y]+EPoint[Z]*TPat->Vals.Gradient[Z];
  4301.  
  4302.   return(value-floor(value));
  4303. }
  4304.  
  4305. /*****************************************************************************
  4306. *
  4307. * FUNCTION
  4308. *
  4309. *   sheet
  4310. *
  4311. * INPUT
  4312. *
  4313. *   EPoint -- The point in 3d space at which the pattern
  4314. *   is evaluated.
  4315. *   
  4316. * OUTPUT
  4317. *   
  4318. * RETURNS
  4319. *
  4320. *   DBL value in the range 0.0 to 1.0
  4321. *   
  4322. * AUTHOR
  4323. *
  4324. *   John VanSickle, based on the gradient code
  4325. *   
  4326. * DESCRIPTION
  4327. *
  4328. *   Takes the dot product of the evaluated point and a specified vector.
  4329. *   If the product is <0, 0 is returned.  If the product is >1, 1 is returned.
  4330. *   Otherwise the product is returned.  This is basically a non-repeating,
  4331. *   non-mirrored gradient, but like the banded pattern above, the width of
  4332. *   the band is independent of direction of the vector.
  4333. *
  4334. * CHANGES
  4335. *   Oct 1994    : adapted from pigment by [CY]
  4336. *   Jul 1999    : adapted from gradient by [JV]
  4337. *
  4338. ******************************************************************************/
  4339.  
  4340. /* íMUY IMPORTANTE! This function expects the Vals.Gradient vector in TPat to
  4341.  * be normalized (ie, unit length).  If the vector is greater or less than
  4342.  * unit length, there will be a scaling effect on the pattern.  I could have
  4343.  * divided by the magnitude of the vector, but since this function will be
  4344.  * called more than once by the renderer, it is more appropriate for the
  4345.  * parser to normalize the vector (where it need be done only once). */
  4346.  
  4347. static DBL sheet (VECTOR EPoint, TPATTERN *TPat)
  4348. {
  4349.   register DBL value;
  4350.  
  4351.   value = EPoint[X]*TPat->Vals.Gradient[X]
  4352.         + EPoint[Y]*TPat->Vals.Gradient[Y]+EPoint[Z]*TPat->Vals.Gradient[Z];
  4353.  
  4354.   return(value>1.0 ? 1.0 : (value<0.0? 0.0:value));
  4355. }
  4356. #endif
  4357. #ifdef BlobPatternPatch
  4358. /*Chris Huff-blob pattern*/
  4359. DBL blob_comp_strength(VECTOR EPoint, BLOB_PATTERN_DATA * comp)
  4360. {
  4361.     VECTOR TPoint;
  4362.     DBL distance = 0;
  4363.     DBL strength = comp->strength;
  4364.     DBL radius = comp->radius;
  4365.     DBL falloff = comp->falloff;
  4366.     DBL fieldVal = 0;
  4367.     MInvTransPoint (TPoint, EPoint, comp->transform);
  4368.     /*Calculate distance*/
  4369.     switch(comp->type)
  4370.     {
  4371.         case 0:
  4372.         {/*spherical component*/
  4373.             VDist(distance, TPoint, comp->center);
  4374.             if(distance < radius)/*Clip to the component shape*/
  4375.             {
  4376.                 if(comp->inverse == TRUE)
  4377.                 {
  4378.                     fieldVal = strength*eval_density_func(1-(distance/radius), falloff, comp->function);
  4379.                 }
  4380.                 else
  4381.                 {
  4382.                     fieldVal = strength*eval_density_func(distance/radius, falloff, comp->function);
  4383.                 }
  4384.             }
  4385.         }
  4386.         break;
  4387.         case 1:
  4388.         {/*cylinderical component*/
  4389.             /*calculate closest point on axis of cylinder to evaluation point*/
  4390.             DBL x1 = comp->center[X];
  4391.             DBL x2 = comp->pointB[X];
  4392.             DBL x3 = TPoint[X];
  4393.             DBL y1 = comp->center[Y];
  4394.             DBL y2 = comp->pointB[Y];
  4395.             DBL y3 = TPoint[Y];
  4396.             DBL z1 = comp->center[Z];
  4397.             DBL z2 = comp->pointB[Z];
  4398.             DBL z3 = TPoint[Z];
  4399.             DBL dist = ((x3 - x1)*(x2 - x1) + (y3 - y1)*(y2 - y1) + (z3 - z1)*(z2 - z1));
  4400.             
  4401.             VECTOR nearestPt;
  4402.             VECTOR cylAxis;
  4403.             DBL len;
  4404.             VSub(cylAxis, comp->pointB, comp->center);
  4405.             VLength(len, cylAxis);
  4406.             dist = dist/(len*len);
  4407.                         
  4408.             if(dist <= 0)
  4409.             {
  4410.                 VDist(distance, TPoint, comp->center);
  4411.             }
  4412.             else
  4413.             {
  4414.                 if(dist <= 1)
  4415.                 {
  4416.                     VScale(nearestPt, cylAxis, dist);
  4417.                     VAddEq(nearestPt, comp->center);
  4418.                     VDist(distance, TPoint, nearestPt);
  4419.                 }
  4420.                 else
  4421.                 {
  4422.                     VDist(distance, TPoint, comp->pointB);
  4423.                 }
  4424.             }
  4425.             if(distance < radius)/*Clip to the component shape*/
  4426.             {
  4427.                 if(comp->inverse == TRUE)
  4428.                 {
  4429.                     fieldVal = strength*eval_density_func(1-(distance/radius), falloff, comp->function);
  4430.                 }
  4431.                 else
  4432.                 {
  4433.                     fieldVal = strength*eval_density_func(distance/radius, falloff, comp->function);
  4434.                 }
  4435.             }
  4436.         }
  4437.         break;
  4438.         case 2:
  4439.         {/*box component*/
  4440.             VECTOR halfSize;
  4441.             VECTOR center;
  4442.             VHalf(center, comp->center, comp->pointB);
  4443.             
  4444.             Assign_Vector(halfSize, comp->pointB);
  4445.             VSubEq(halfSize, comp->center);
  4446.             VInverseScaleEq(halfSize, 2);
  4447.             
  4448.             distance = max(fabs(TPoint[X]-center[X])/halfSize[X],
  4449.                        max(fabs(TPoint[Y]-center[Y])/halfSize[Y],
  4450.                            fabs(TPoint[Z]-center[Z])/halfSize[Z]));
  4451.             if(distance < 1)
  4452.             {/*Clip to the component shape*/
  4453.                 if(comp->inverse == TRUE)
  4454.                 {
  4455.                     fieldVal = strength*eval_density_func(1-distance, falloff, comp->function);
  4456.                 }
  4457.                 else
  4458.                 {
  4459.                     fieldVal = strength*eval_density_func(distance, falloff, comp->function);
  4460.                 }
  4461.             }
  4462.         }
  4463.         break;
  4464.         case 3:
  4465.         {/*pigment component*/
  4466.             COLOUR tempCol;
  4467.             int colour_found;
  4468.             
  4469.             if (comp->pigment)
  4470.             {
  4471.                 colour_found = Compute_Pigment (tempCol, comp->pigment, TPoint, NULL);
  4472.             }
  4473.  
  4474.             if(!colour_found)
  4475.             {
  4476.                 distance = 0.0;
  4477.             }
  4478.             else
  4479.             {
  4480.                 distance = GREY_SCALE(tempCol);
  4481.             }
  4482.             if(comp->inverse == TRUE)
  4483.             {
  4484.                 fieldVal = strength*eval_density_func(1-distance, falloff, comp->function);
  4485.             }
  4486.             else
  4487.             {
  4488.                 fieldVal = strength*eval_density_func(distance, falloff, comp->function);
  4489.             }
  4490.         }
  4491.         break;
  4492.         case 4:
  4493.         {/*blob component*/
  4494.             fieldVal = calculate_blob_field_value((BLOB*)comp->blob, TPoint);
  4495. /*            if(comp->inverse == TRUE)
  4496.             {
  4497.                 fieldVal = 1-fieldVal;
  4498.             }*/
  4499.         }
  4500.         break;
  4501.         default:
  4502.         {/*spherical component*/
  4503.             VDist(distance, TPoint, comp->center);
  4504.             if(distance < radius)/*Clip to the component shape*/
  4505.             {
  4506.                 fieldVal = strength*eval_density_func(distance/radius, falloff, comp->function);
  4507.             }
  4508.         }
  4509.     }
  4510.     return fieldVal;
  4511. }
  4512. static DBL blob_pattern (VECTOR EPoint, TPATTERN *TPat)
  4513. {
  4514.     DBL totalVal = 0;
  4515.     BLOB_PATTERN_DATA * currentComponent = NULL;
  4516.     DBL max_density = TPat->Vals.Blob.max_density;
  4517.     DBL threshold = TPat->Vals.Blob.blob_threshold;
  4518.     
  4519.     if(TPat->Vals.Blob.blob_dat == NULL)
  4520.         return 0;/*If the list of components is empty, return 0*/
  4521.     
  4522.     currentComponent = TPat->Vals.Blob.blob_dat;
  4523.     while(currentComponent != NULL)
  4524.     {
  4525.         totalVal += blob_comp_strength(EPoint, currentComponent);
  4526.         currentComponent = currentComponent->next;
  4527.     }
  4528.     
  4529.     
  4530.     if(totalVal < threshold) totalVal = threshold;
  4531.     else if(totalVal > max_density) totalVal = max_density;
  4532.     
  4533.     if((max_density - threshold) != 0 && (totalVal - threshold) != 0)
  4534.         totalVal = (totalVal - threshold)/(max_density - threshold);
  4535.     else
  4536.         totalVal = 0;
  4537.     
  4538.     return totalVal;
  4539. }
  4540. /*Chris Huff proximity and blob patterns*/
  4541. DBL eval_density_func(DBL val, DBL falloff, int func)
  4542. {
  4543.     DBL temp = 0;
  4544.     switch(func)/*decide on density function and calculate the density*/
  4545.     {
  4546.         case 0:
  4547.             if(falloff == 1)
  4548.                 temp = 1 - val;
  4549.             else
  4550.             {
  4551.                 if(falloff == 2)
  4552.                     temp = Sqr(1 - Sqr(val));
  4553.                 else
  4554.                     temp = pow(1 - pow(val, falloff), falloff);
  4555.             }
  4556.             break;
  4557.         case 1:
  4558.             if(falloff == 1)
  4559.                 temp = val;
  4560.             else
  4561.                 if(falloff == 2)
  4562.                     temp = 1/Sqr(val);
  4563.                 else
  4564.                     temp = 1/pow(val, falloff);
  4565.             break;
  4566.         case 2:
  4567.             if(falloff == 1)
  4568.                 temp = 1 - val;
  4569.             else
  4570.                 if(falloff == 2)
  4571.                     temp = 1 - Sqr(val);
  4572.                 else
  4573.                     temp = 1 - pow(val, falloff);
  4574.             break;
  4575.         case 3:
  4576.             if(falloff == 1)
  4577.                 temp = val;
  4578.             else
  4579.                 if(falloff == 2)
  4580.                     temp = Sqr(val);
  4581.                 else
  4582.                     temp = pow(val, falloff);
  4583.             break;
  4584.         case 4:/*No filter*/
  4585.             temp = val;
  4586.             break;
  4587.         default:
  4588.         temp = Sqr(1 - Sqr(val));
  4589.     }
  4590.     return temp;
  4591. }
  4592. #endif
  4593.  
  4594. #ifdef ProximityPatch
  4595. static DBL proximity(VECTOR EPoint, TPATTERN *TPat)
  4596. {
  4597.    INTERSECTION Intersect;
  4598.    RAY Ray;
  4599.    int k=0;
  4600.    int attempts=0;
  4601.    int insideObjBBox=0;
  4602.    DBL temp = 0;
  4603.    VECTOR tempVect;
  4604.    VECTOR resVect;
  4605.    DBL result = 0;
  4606.    DBL x,y,z,r,t;
  4607.    VECTOR minExtent;
  4608.    VECTOR maxExtent;
  4609.    Assign_Vector(Ray.Initial, EPoint);
  4610.    if(Test_Flag(TPat->Vals.Proximity.proxObject, INFINITE_FLAG))
  4611.        return 0;
  4612.  
  4613.    Make_min_max_from_BBox(minExtent, maxExtent, TPat->Vals.Proximity.proxObject->BBox)
  4614.    if( (EPoint[X] > minExtent[X])&&(EPoint[Y] > minExtent[Y])&&(EPoint[Z] > minExtent[Z])&&
  4615.       (EPoint[X] < maxExtent[X])&&(EPoint[Y] < maxExtent[Y])&&(EPoint[Z] < maxExtent[Z]))
  4616.    {
  4617.      insideObjBBox = 1;
  4618.    }
  4619.    
  4620.    Make_Vector(tempVect,0,0,0);
  4621.    Make_Vector(resVect,0,0,0);
  4622.    Initialize_Ray_Containers( &Ray );
  4623.    
  4624.    /*detect whether the point is inside the object or not*/
  4625.    if(Inside_Object(EPoint, TPat->Vals.Proximity.proxObject))
  4626.    {
  4627.            if(!((TPat->Vals.Proximity.sides == 0)||(TPat->Vals.Proximity.sides == 2))) 
  4628.                return 0;
  4629.    }
  4630.    else 
  4631.    {
  4632.            if(!((TPat->Vals.Proximity.sides == 1)||(TPat->Vals.Proximity.sides == 2))) 
  4633.                return 0;
  4634.    }
  4635.    while((k<TPat->Vals.Proximity.samples)&&(attempts < TPat->Vals.Proximity.sample_bailout))
  4636.    {
  4637.        attempts++;
  4638.        switch(TPat->Vals.Proximity.sampleMthd)
  4639.        {
  4640.            case 0:
  4641.            {
  4642.            /*pick a random vector in the bounding box of the sample object*/
  4643.                x = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[X];
  4644.                x += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[X];
  4645.                y = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[Y];
  4646.                y += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[Y];
  4647.                z = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[Z];
  4648.                z += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[Z];
  4649.                Make_Vector( Ray.Direction, x, y, z);
  4650.                VSubEq(Ray.Direction, EPoint);
  4651.            } break;
  4652.            
  4653.            case 1:
  4654.            {
  4655.            /* Pick a random vector with length <= 1 using the "trig method"
  4656.              described in the comp.graphics.algorithms FAQ. 
  4657.              Chris Huff: I got this code from the WyzPov reflection blur*/
  4658.                z = (FRAND() * 2) - 1;
  4659.                t = FRAND() * M_PI * 2;
  4660.                r = sqrt(1 - z*z);
  4661.                x = r * cos(t);
  4662.                y = r * sin(t);
  4663.                Make_Vector( Ray.Direction, x, y, z);
  4664.            } break;
  4665.            
  4666.            
  4667.            case 2:
  4668.            if(insideObjBBox)
  4669.            {
  4670.            /* Pick a random vector with length <= 1 using the "trig method"
  4671.              described in the comp.graphics.algorithms FAQ. 
  4672.              Chris Huff: I got this code from the WyzPov reflection blur*/
  4673.                z = (FRAND() * 2) - 1;
  4674.                t = FRAND() * M_PI * 2;
  4675.                r = sqrt(1 - z*z);
  4676.                x = r * cos(t);
  4677.                y = r * sin(t);
  4678.                Make_Vector( Ray.Direction, x, y, z);
  4679.            }
  4680.            else
  4681.            {
  4682.            /*pick a random vector in the bounding box of the sample object*/
  4683.                x = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[X];
  4684.                x += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[X];
  4685.                y = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[Y];
  4686.                y += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[Y];
  4687.                z = FRAND()*TPat->Vals.Proximity.proxObject->BBox.Lengths[Z];
  4688.                z += TPat->Vals.Proximity.proxObject->BBox.Lower_Left[Z];
  4689.                Make_Vector( Ray.Direction, x, y, z);
  4690.                VSubEq(Ray.Direction, EPoint);
  4691.            } break;
  4692.            case 3:/*planar sampling*/
  4693.            {
  4694.                z = (FRAND() * 2) - 1;
  4695.                t = FRAND() * M_PI * 2;
  4696.                r = sqrt(1 - z*z);
  4697.                x = r * cos(t);
  4698.                y = r * sin(t);
  4699.                Make_Vector( Ray.Direction, x, y, 0);   
  4700.            } break;
  4701.        }
  4702.        
  4703.        VAddEq(Ray.Direction, TPat->Vals.Proximity.sampleWeight);
  4704.        VNormalizeEq(Ray.Direction);
  4705.        
  4706.        if( Intersection( &Intersect, TPat->Vals.Proximity.proxObject, &Ray ) )
  4707.        {
  4708. /*           VDist(temp, EPoint, Intersect.IPoint);*/
  4709.            switch(TPat->Vals.Proximity.proxCalcMthd)
  4710.            {
  4711.                case 0:
  4712.                    VDist(temp, EPoint, Intersect.IPoint);
  4713.                    result += temp;
  4714.                break;
  4715.                case 1:
  4716.                {
  4717.                    VDist(temp, EPoint, Intersect.IPoint);
  4718.                    if((temp < result)||(result == 0))/*the last part is necessary because result == 0 at first*/
  4719.                        result = temp;
  4720.                }
  4721.                break;
  4722.                case 2:
  4723.                {
  4724.                    VSub(tempVect, Intersect.IPoint, EPoint);
  4725.                    VAddEq(resVect, tempVect);
  4726.                }
  4727.                break;
  4728.  
  4729.                
  4730.                case 3:
  4731.                    VDist(temp, EPoint, Intersect.IPoint);
  4732.                    result += 1/(1+pow(temp, TPat->Vals.Proximity.falloff));
  4733.                break;
  4734.                
  4735.                default:
  4736.                result += temp;
  4737.            }
  4738.            k++;
  4739.        }
  4740.    }
  4741.    
  4742.    switch(TPat->Vals.Proximity.proxCalcMthd)
  4743.    {
  4744.        case 0:
  4745.            result /= k;/*divide by total samples to get average*/
  4746.        break;
  4747.        case 1:
  4748.          {;}
  4749.        break;
  4750.        case 2:
  4751.          VLength(result, resVect);
  4752.            break;
  4753.        case 3:
  4754.            result /= k;/*divide by total to get average*/
  4755.            break;
  4756.        
  4757.        default:
  4758.          result /= k;/*divide by total to get average*/
  4759.    }
  4760.  
  4761.    
  4762.    result /= TPat->Vals.Proximity.proxDist;
  4763.    if(result>TPat->Vals.Proximity.max_density)
  4764.      result = TPat->Vals.Proximity.max_density;
  4765.     
  4766.    return result;
  4767. }
  4768.  
  4769. #endif
  4770.  
  4771.  
  4772. #ifdef ObjectPatternPatch
  4773. static DBL object(VECTOR EPoint, TPATTERN *TPat)
  4774. {
  4775.    if(TPat->Vals.Object != NULL)
  4776.    {
  4777.        if(Inside_Object(EPoint, TPat->Vals.Object))
  4778.         return 1.0;
  4779.        else
  4780.         return 0.0;
  4781.    }
  4782.    else
  4783.     return 0.0;
  4784. }
  4785. #pragma peephole off 
  4786.  
  4787. #endif
  4788.