home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: Graphics / Graphics.zip / povsrc31.zip / hcmplx.c < prev    next >
C/C++ Source or Header  |  1999-10-20  |  31KB  |  1,432 lines
  1. /*****************************************************************************
  2. *                hcmplx.c
  3. *
  4. *  This module implements hypercomplex Julia fractals.
  5. *
  6. *  This file was written by Pascal Massimino.
  7. *
  8. *  from Persistence of Vision(tm) Ray Tracer
  9. *  Copyright 1996,1999 Persistence of Vision Team
  10. *---------------------------------------------------------------------------
  11. *  NOTICE: This source code file is provided so that users may experiment
  12. *  with enhancements to POV-Ray and to port the software to platforms other
  13. *  than those supported by the POV-Ray Team.  There are strict rules under
  14. *  which you are permitted to use this file.  The rules are in the file
  15. *  named POVLEGAL.DOC which should be distributed with this file.
  16. *  If POVLEGAL.DOC is not available or for more info please contact the POV-Ray
  17. *  Team Coordinator by email to team-coord@povray.org or visit us on the web at
  18. *  http://www.povray.org. The latest version of POV-Ray may be found at this site.
  19. *
  20. * This program is based on the popular DKB raytracer version 2.12.
  21. * DKBTrace was originally written by David K. Buck.
  22. * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
  23. *
  24. *****************************************************************************/
  25.  
  26. #include "frame.h"
  27. #include "povray.h"
  28. #include "vector.h"
  29. #include "povproto.h"
  30. #include "fractal.h"
  31. #include "spheres.h"
  32. #include "hcmplx.h"
  33.  
  34.  
  35.  
  36. /*****************************************************************************
  37. * Local preprocessor defines
  38. ******************************************************************************/
  39.  
  40. #ifndef Fractal_Tolerance
  41. #define Fractal_Tolerance 1e-8
  42. #endif
  43.  
  44.  
  45. #define HMult(xr,yr,zr,wr,x1,y1,z1,w1,x2,y2,z2,w2)        \
  46.     (xr) = (x1) * (x2) - (y1) * (y2) - (z1) * (z2) + (w1) * (w2);   \
  47.     (yr) = (y1) * (x2) + (x1) * (y2) - (w1) * (z2) - (z1) * (w2);   \
  48.     (zr) = (z1) * (x2) - (w1) * (y2) + (x1) * (z2) - (y1) * (w2);   \
  49.     (wr) = (w1) * (x2) + (z1) * (y2) + (y1) * (z2) + (x1) * (w2);
  50.  
  51. #define HSqr(xr,yr,zr,wr,x,y,z,w)         \
  52.     (xr) = (x) * (x) - (y) * (y) - (z) * (z) + (w) * (w) ;  \
  53.     (yr) = 2.0 * ( (x) * (y) - (z) * (w) );     \
  54.     (zr) = 2.0 * ( (z) * (x) - (w) * (y) );       \
  55.     (wr) = 2.0 * ( (w) * (x) + (z) * (y) );
  56.  
  57.  
  58.  
  59. /*****************************************************************************
  60. * Local typedefs
  61. ******************************************************************************/
  62.  
  63.  
  64.  
  65. /*****************************************************************************
  66. * Static functions
  67. ******************************************************************************/
  68.  
  69. static int HReciprocal (DBL * xr, DBL * yr, DBL * zr, DBL * wr, DBL x, DBL y, DBL z, DBL w);
  70. static void HFunc (DBL * xr, DBL * yr, DBL * zr, DBL * wr, DBL x, DBL y, DBL z, DBL w, FRACTAL *f);
  71.  
  72.  
  73.  
  74. /*****************************************************************************
  75. * Local variables
  76. ******************************************************************************/
  77.  
  78. static CMPLX exponent = {0,0};
  79.  
  80. /******** Computations with Hypercomplexes **********/
  81.  
  82.  
  83.  
  84. /*****************************************************************************
  85. *
  86. * FUNCTION
  87. *
  88. * INPUT
  89. *
  90. * OUTPUT
  91. *
  92. * RETURNS
  93. *
  94. * AUTHOR
  95. *
  96. *   Pascal Massimino
  97. *
  98. * DESCRIPTION
  99. *
  100. * CHANGES
  101. *
  102. ******************************************************************************/
  103.  
  104. static int HReciprocal(DBL *xr, DBL  *yr, DBL  *zr, DBL  *wr, DBL  x, DBL  y, DBL  z, DBL  w)
  105. {
  106.   DBL det, mod, xt_minus_yz;
  107.  
  108.   det = ((x - w) * (x - w) + (y + z) * (y + z)) * ((x + w) * (x + w) + (y - z) * (y - z));
  109.  
  110.   if (det == 0.0)
  111.   {
  112.     return (-1);
  113.   }
  114.  
  115.   mod = (x * x + y * y + z * z + w * w);
  116.  
  117.   xt_minus_yz = x * w - y * z;
  118.  
  119.   *xr = (x * mod - 2 * w * xt_minus_yz) / det;
  120.   *yr = (-y * mod - 2 * z * xt_minus_yz) / det;
  121.   *zr = (-z * mod - 2 * y * xt_minus_yz) / det;
  122.   *wr = (w * mod - 2 * x * xt_minus_yz) / det;
  123.  
  124.   return (0);
  125. }
  126.  
  127.  
  128.  
  129. /*****************************************************************************
  130. *
  131. * FUNCTION Hfunc
  132. *
  133. * INPUT 4D Hypercomplex number, pointer to fractal object
  134. *
  135. * OUTPUT  calculates the 4D generalization of fractal->Complex_Function
  136. *
  137. * RETURNS void
  138. *
  139. * AUTHOR
  140. *
  141. *   Pascal Massimino
  142. *
  143. * DESCRIPTION
  144. *   Hypercomplex numbers allow generalization of any complex->complex
  145. *   function in a uniform way. This function implements a general
  146. *   unary 4D function based on the corresponding complex function. 
  147. *
  148. * CHANGES
  149. *  Generalized to use Complex_Function()   TW 
  150. *
  151. ******************************************************************************/
  152.  
  153. static void HFunc(DBL *xr, DBL  *yr, DBL  *zr, DBL  *wr, DBL  x, DBL  y, DBL  z, DBL  w, FRACTAL *f)
  154. {
  155.   CMPLX a, b, ra, rb;
  156.   
  157.   /* convert to duplex form */
  158.   a.x = x - w;
  159.   a.y = y + z;
  160.   b.x = x + w;
  161.   b.y = y - z;
  162.   
  163.   if(f->Sub_Type == PWR_STYPE)
  164.   {
  165.      exponent = f->exponent;
  166.   }
  167.   
  168.   /* apply function to each part */
  169.   Complex_Function(&ra, &a, f);
  170.   Complex_Function(&rb, &b, f);
  171.  
  172.   /* convert back */
  173.   *xr = .5 * (ra.x + rb.x);
  174.   *yr = .5 * (ra.y + rb.y);
  175.   *zr = .5 * (ra.y - rb.y);
  176.   *wr = .5 * (rb.x - ra.x);
  177. }
  178.  
  179.  
  180.  
  181. /*****************************************************************************
  182. *
  183. * FUNCTION
  184. *
  185. * INPUT
  186. *
  187. * OUTPUT
  188. *
  189. * RETURNS
  190. *
  191. * AUTHOR
  192. *
  193. *   Pascal Massimino
  194. *
  195. * DESCRIPTION
  196. *
  197. * CHANGES
  198. *
  199. ******************************************************************************/
  200.  
  201. /*------------------ z2 Iteration method ------------------*/
  202.  
  203. int Iteration_HCompl(VECTOR IPoint, FRACTAL *HCompl)
  204. {
  205.   int i;
  206.   DBL yz, xw;
  207.   DBL Exit_Value;
  208.   DBL x, y, z, w;
  209.  
  210.   x = Sx[0] = IPoint[X];
  211.   y = Sy[0] = IPoint[Y];
  212.   z = Sz[0] = IPoint[Z];
  213.   w = Sw[0] = (HCompl->SliceDist
  214.                   - HCompl->Slice[X]*x 
  215.                   - HCompl->Slice[Y]*y 
  216.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  217.   
  218.   Exit_Value = HCompl->Exit_Value;
  219.  
  220.   for (i = 1; i <= HCompl->n; ++i)
  221.   {
  222.     yz = y * y + z * z;
  223.     xw = x * x + w * w;
  224.  
  225.     if ((xw + yz) > Exit_Value)
  226.     {
  227.       return (FALSE);
  228.     }
  229.  
  230.     Sx[i] = xw - yz + HCompl->Julia_Parm[X];
  231.     Sy[i] = 2.0 * (x * y - z * w) + HCompl->Julia_Parm[Y];
  232.     Sz[i] = 2.0 * (x * z - w * y) + HCompl->Julia_Parm[Z];
  233.     Sw[i] = 2.0 * (x * w + y * z) + HCompl->Julia_Parm[T];
  234.  
  235.     w = Sw[i];
  236.     x = Sx[i];
  237.  
  238.     z = Sz[i];
  239.     y = Sy[i];
  240.   }
  241.  
  242.   return (TRUE);
  243. }
  244.  
  245.  
  246.  
  247. /*****************************************************************************
  248. *
  249. * FUNCTION
  250. *
  251. * INPUT
  252. *
  253. * OUTPUT
  254. *
  255. * RETURNS
  256. *
  257. * AUTHOR
  258. *
  259. *   Pascal Massimino
  260. *
  261. * DESCRIPTION
  262. *
  263. * CHANGES
  264. *
  265. ******************************************************************************/
  266.  
  267. int D_Iteration_HCompl(VECTOR IPoint, FRACTAL *HCompl, DBL *Dist)
  268. {
  269.   int i;
  270.   DBL yz, xw;
  271.   DBL Exit_Value, F_Value, Step;
  272.   DBL x, y, z, w;
  273.   VECTOR H_Normal;
  274.  
  275.   x = Sx[0] = IPoint[X];
  276.   y = Sy[0] = IPoint[Y];
  277.   z = Sz[0] = IPoint[Z];
  278.   w = Sw[0] = (HCompl->SliceDist
  279.                   - HCompl->Slice[X]*x 
  280.                   - HCompl->Slice[Y]*y 
  281.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  282.  
  283.   Exit_Value = HCompl->Exit_Value;
  284.  
  285.   for (i = 1; i <= HCompl->n; ++i)
  286.   {
  287.     yz = y * y + z * z;
  288.     xw = x * x + w * w;
  289.  
  290.     if ((F_Value = xw + yz) > Exit_Value)
  291.     {
  292.       Normal_Calc_HCompl(H_Normal, i - 1, HCompl);
  293.  
  294.       VDot(Step, H_Normal, Direction);
  295.  
  296.       if (Step < -Fractal_Tolerance)
  297.       {
  298.         Step = -2.0 * Step;
  299.  
  300.         if ((F_Value > Precision * Step) && (F_Value < 30 * Precision * Step))
  301.         {
  302.           *Dist = F_Value / Step;
  303.  
  304.           return (FALSE);
  305.         }
  306.       }
  307.  
  308.       *Dist = Precision;
  309.  
  310.       return (FALSE);
  311.     }
  312.  
  313.     Sx[i] = xw - yz + HCompl->Julia_Parm[X];
  314.     Sy[i] = 2.0 * (x * y - z * w) + HCompl->Julia_Parm[Y];
  315.     Sz[i] = 2.0 * (x * z - w * y) + HCompl->Julia_Parm[Z];
  316.     Sw[i] = 2.0 * (x * w + y * z) + HCompl->Julia_Parm[T];
  317.  
  318.     w = Sw[i];
  319.     x = Sx[i];
  320.  
  321.     z = Sz[i];
  322.     y = Sy[i];
  323.   }
  324.  
  325.   *Dist = Precision;
  326.  
  327.   return (TRUE);
  328. }
  329.  
  330.  
  331.  
  332. /*****************************************************************************
  333. *
  334. * FUNCTION
  335. *
  336. * INPUT
  337. *
  338. * OUTPUT
  339. *
  340. * RETURNS
  341. *
  342. * AUTHOR
  343. *
  344. *   Pascal Massimino
  345. *
  346. * DESCRIPTION
  347. *
  348. * CHANGES
  349. *
  350. ******************************************************************************/
  351.  
  352. void Normal_Calc_HCompl(VECTOR Result, int N_Max, FRACTAL *fractal)
  353. {
  354.   DBL n1, n2, n3, n4;
  355.   int i;
  356.   DBL x, y, z, w;
  357.   DBL xx, yy, zz, ww;
  358.   DBL Pow;
  359.  
  360.   /*
  361.    * Algebraic properties of hypercomplexes allows simplifications in
  362.    * computations...
  363.    */
  364.  
  365.   x = Sx[0];
  366.   y = Sy[0];
  367.   z = Sz[0];
  368.   w = Sw[0];
  369.  
  370.   Pow = 2.0;
  371.  
  372.   for (i = 1; i < N_Max; ++i)
  373.   {
  374.  
  375.     /*
  376.      * For a map z->f(z), f depending on c, one must perform here :
  377.      *
  378.      * (x,y,z,w) * df/dz(Sx[i],Sy[i],Sz[i],Sw[i]) -> (x,y,z,w) ,
  379.      *
  380.      * up to a constant.
  381.      */
  382.  
  383.     /******************* Case z->z^2+c *****************/
  384.  
  385.     /* the df/dz part needs no work */
  386.  
  387.     HMult(xx, yy, zz, ww, Sx[i], Sy[i], Sz[i], Sw[i], x, y, z, w);
  388.  
  389.     w = ww;
  390.     z = zz;
  391.     y = yy;
  392.     x = xx;
  393.  
  394.     Pow *= 2.0;
  395.   }
  396.  
  397.   n1 = Sx[N_Max] * Pow;
  398.   n2 = Sy[N_Max] * Pow;
  399.   n3 = Sz[N_Max] * Pow;
  400.   n4 = Sw[N_Max] * Pow;
  401.  
  402.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  403.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  404.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  405. }
  406.  
  407.  
  408.  
  409. /*****************************************************************************
  410. *
  411. * FUNCTION
  412. *
  413. * INPUT
  414. *
  415. * OUTPUT
  416. *
  417. * RETURNS
  418. *
  419. * AUTHOR
  420. *
  421. *   Pascal Massimino
  422. *
  423. * DESCRIPTION
  424. *
  425. * CHANGES
  426. *
  427. ******************************************************************************/
  428.  
  429. int F_Bound_HCompl(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL  *Depth_Max)
  430. {
  431.   return (Intersect_Sphere(Ray, Fractal->Center, Fractal->Radius_Squared, Depth_Min, Depth_Max));
  432. }
  433.  
  434. /****************************************************************/
  435. /*--------------------------- z3 -------------------------------*/
  436. /****************************************************************/
  437.  
  438.  
  439.  
  440. /*****************************************************************************
  441. *
  442. * FUNCTION
  443. *
  444. * INPUT
  445. *
  446. * OUTPUT
  447. *
  448. * RETURNS
  449. *
  450. * AUTHOR
  451. *
  452. *   Pascal Massimino
  453. *
  454. * DESCRIPTION
  455. *
  456. * CHANGES
  457. *
  458. ******************************************************************************/
  459.  
  460. int Iteration_HCompl_z3(VECTOR IPoint, FRACTAL *HCompl)
  461. {
  462.   int i;
  463.   DBL Norm, xx, yy, zz, ww;
  464.   DBL Exit_Value;
  465.   DBL x, y, z, w;
  466.  
  467.   x = Sx[0] = IPoint[X];
  468.   y = Sy[0] = IPoint[Y];
  469.   z = Sz[0] = IPoint[Z];
  470.   w = Sw[0] = (HCompl->SliceDist
  471.                   - HCompl->Slice[X]*x 
  472.                   - HCompl->Slice[Y]*y 
  473.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  474.  
  475.   Exit_Value = HCompl->Exit_Value;
  476.  
  477.   for (i = 1; i <= HCompl->n; ++i)
  478.   {
  479.     Norm = x * x + y * y + z * z + w * w;
  480.  
  481.     /* is this test correct ? */
  482.     if (Norm > Exit_Value)
  483.     {
  484.       return (FALSE);
  485.     }
  486.  
  487.     /*************** Case: z->z^2+c *********************/
  488.     HSqr(xx, yy, zz, ww, x, y, z, w);
  489.  
  490.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  491.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  492.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  493.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  494.  
  495.   }
  496.  
  497.   return (TRUE);
  498. }
  499.  
  500.  
  501.  
  502. /*****************************************************************************
  503. *
  504. * FUNCTION
  505. *
  506. * INPUT
  507. *
  508. * OUTPUT
  509. *
  510. * RETURNS
  511. *
  512. * AUTHOR
  513. *
  514. *   Pascal Massimino
  515. *
  516. * DESCRIPTION
  517. *
  518. * CHANGES
  519. *
  520. ******************************************************************************/
  521.  
  522. int D_Iteration_HCompl_z3(VECTOR IPoint, FRACTAL *HCompl, DBL *Dist)
  523. {
  524.   int i;
  525.   DBL xx, yy, zz, ww;
  526.   DBL Exit_Value, F_Value, Step;
  527.   DBL x, y, z, w;
  528.   VECTOR H_Normal;
  529.  
  530.   x = Sx[0] = IPoint[X];
  531.   y = Sy[0] = IPoint[Y];
  532.   z = Sz[0] = IPoint[Z];
  533.   w = Sw[0] = (HCompl->SliceDist
  534.                   - HCompl->Slice[X]*x 
  535.                   - HCompl->Slice[Y]*y 
  536.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  537.  
  538.   Exit_Value = HCompl->Exit_Value;
  539.  
  540.   for (i = 1; i <= HCompl->n; ++i)
  541.   {
  542.     F_Value = x * x + y * y + z * z + w * w;
  543.  
  544.     if (F_Value > Exit_Value)
  545.     {
  546.       Normal_Calc_HCompl_z3(H_Normal, i - 1, HCompl);
  547.  
  548.       VDot(Step, H_Normal, Direction);
  549.  
  550.       if (Step < -Fractal_Tolerance)
  551.       {
  552.         Step = -2.0 * Step;
  553.  
  554.         if ((F_Value > Precision * Step) && (F_Value < 30 * Precision * Step))
  555.         {
  556.           *Dist = F_Value / Step;
  557.  
  558.           return (FALSE);
  559.         }
  560.       }
  561.  
  562.       *Dist = Precision;
  563.  
  564.       return (FALSE);
  565.     }
  566.  
  567.     /*************** Case: z->z^2+c *********************/
  568.  
  569.     HSqr(xx, yy, zz, ww, x, y, z, w);
  570.  
  571.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  572.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  573.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  574.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  575.   }
  576.  
  577.   *Dist = Precision;
  578.  
  579.   return (TRUE);
  580. }
  581.  
  582.  
  583.  
  584. /*****************************************************************************
  585. *
  586. * FUNCTION
  587. *
  588. * INPUT
  589. *
  590. * OUTPUT
  591. *
  592. * RETURNS
  593. *
  594. * AUTHOR
  595. *
  596. *   Pascal Massimino
  597. *
  598. * DESCRIPTION
  599. *
  600. * CHANGES
  601. *
  602. ******************************************************************************/
  603.  
  604. void Normal_Calc_HCompl_z3(VECTOR Result, int N_Max, FRACTAL *fractal)
  605. {
  606.   DBL n1, n2, n3, n4;
  607.   int i;
  608.   DBL x, y, z, w;
  609.   DBL xx, yy, zz, ww;
  610.  
  611.   /*
  612.    * Algebraic properties of hypercomplexes allows simplifications in
  613.    * computations...
  614.    */
  615.  
  616.   x = Sx[0];
  617.   y = Sy[0];
  618.   z = Sz[0];
  619.   w = Sw[0];
  620.  
  621.   for (i = 1; i < N_Max; ++i)
  622.   {
  623.     /*
  624.      * For a map z->f(z), f depending on c, one must perform here :
  625.      * 
  626.      * (x,y,z,w) * df/dz(Sx[i],Sy[i],Sz[i],Sw[i]) -> (x,y,z,w) ,
  627.      * 
  628.      * up to a constant.
  629.      */
  630.  
  631.     /******************* Case z->z^2+c *****************/
  632.  
  633.     /* the df/dz part needs no work */
  634.  
  635.     HMult(xx, yy, zz, ww, Sx[i], Sy[i], Sz[i], Sw[i], x, y, z, w);
  636.  
  637.     x = xx;
  638.     y = yy;
  639.     z = zz;
  640.     w = ww;
  641.   }
  642.  
  643.   n1 = Sx[N_Max];
  644.   n2 = Sy[N_Max];
  645.   n3 = Sz[N_Max];
  646.   n4 = Sw[N_Max];
  647.  
  648.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  649.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  650.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  651. }
  652.  
  653.  
  654.  
  655. /*****************************************************************************
  656. *
  657. * FUNCTION
  658. *
  659. * INPUT
  660. *
  661. * OUTPUT
  662. *
  663. * RETURNS
  664. *
  665. * AUTHOR
  666. *
  667. *   Pascal Massimino
  668. *
  669. * DESCRIPTION
  670. *
  671. * CHANGES
  672. *
  673. ******************************************************************************/
  674.  
  675. int F_Bound_HCompl_z3(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL *Depth_Max)
  676. {
  677.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  678. }
  679.  
  680. /*--------------------------- Inv -------------------------------*/
  681.  
  682.  
  683. /*****************************************************************************
  684. *
  685. * FUNCTION
  686. *
  687. * INPUT
  688. *
  689. * OUTPUT
  690. *
  691. * RETURNS
  692. *
  693. * AUTHOR
  694. *
  695. *   Pascal Massimino
  696. *
  697. * DESCRIPTION
  698. *
  699. * CHANGES
  700. *
  701. ******************************************************************************/
  702.  
  703. int Iteration_HCompl_Reciprocal(VECTOR IPoint, FRACTAL *HCompl)
  704. {
  705.   int i;
  706.   DBL Norm, xx, yy, zz, ww;
  707.   DBL Exit_Value;
  708.   DBL x, y, z, w;
  709.  
  710.   x = Sx[0] = IPoint[X];
  711.   y = Sy[0] = IPoint[Y];
  712.   z = Sz[0] = IPoint[Z];
  713.   w = Sw[0] = (HCompl->SliceDist
  714.                   - HCompl->Slice[X]*x 
  715.                   - HCompl->Slice[Y]*y 
  716.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  717.  
  718.   Exit_Value = HCompl->Exit_Value;
  719.  
  720.   for (i = 1; i <= HCompl->n; ++i)
  721.   {
  722.     Norm = x * x + y * y + z * z + w * w;
  723.  
  724.     if (Norm > Exit_Value)
  725.     {
  726.       return (FALSE);
  727.     }
  728.  
  729.     HReciprocal(&xx, &yy, &zz, &ww, x, y, z, w);
  730.  
  731.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  732.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  733.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  734.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  735.  
  736.   }
  737.  
  738.   return (TRUE);
  739. }
  740.  
  741.  
  742.  
  743. /*****************************************************************************
  744. *
  745. * FUNCTION
  746. *
  747. * INPUT
  748. *
  749. * OUTPUT
  750. *
  751. * RETURNS
  752. *
  753. * AUTHOR
  754. *
  755. *   Pascal Massimino
  756. *
  757. * DESCRIPTION
  758. *
  759. * CHANGES
  760. *
  761. ******************************************************************************/
  762.  
  763. int D_Iteration_HCompl_Reciprocal(VECTOR IPoint, FRACTAL *HCompl, DBL *Dist)
  764. {
  765.   int i;
  766.   DBL xx, yy, zz, ww;
  767.   DBL Exit_Value, F_Value, Step;
  768.   DBL x, y, z, w;
  769.   VECTOR H_Normal;
  770.  
  771.   x = Sx[0] = IPoint[X];
  772.   y = Sy[0] = IPoint[Y];
  773.   z = Sz[0] = IPoint[Z];
  774.   w = Sw[0] = (HCompl->SliceDist
  775.                   - HCompl->Slice[X]*x 
  776.                   - HCompl->Slice[Y]*y 
  777.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  778.  
  779.   Exit_Value = HCompl->Exit_Value;
  780.  
  781.   for (i = 1; i <= HCompl->n; ++i)
  782.   {
  783.     F_Value = x * x + y * y + z * z + w * w;
  784.  
  785.     if (F_Value > Exit_Value)
  786.     {
  787.       Normal_Calc_HCompl_Reciprocal(H_Normal, i - 1, HCompl);
  788.  
  789.       VDot(Step, H_Normal, Direction);
  790.  
  791.       if (Step < -Fractal_Tolerance)
  792.       {
  793.         Step = -2.0 * Step;
  794.  
  795.         if ((F_Value > Precision * Step) && F_Value < (30 * Precision * Step))
  796.         {
  797.           *Dist = F_Value / Step;
  798.  
  799.           return (FALSE);
  800.         }
  801.       }
  802.  
  803.       *Dist = Precision;
  804.  
  805.       return (FALSE);
  806.     }
  807.  
  808.     HReciprocal(&xx, &yy, &zz, &ww, x, y, z, w);
  809.  
  810.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  811.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  812.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  813.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  814.  
  815.   }
  816.  
  817.   *Dist = Precision;
  818.  
  819.   return (TRUE);
  820. }
  821.  
  822.  
  823.  
  824. /*****************************************************************************
  825. *
  826. * FUNCTION
  827. *
  828. * INPUT
  829. *
  830. * OUTPUT
  831. *
  832. * RETURNS
  833. *
  834. * AUTHOR
  835. *
  836. *   Pascal Massimino
  837. *
  838. * DESCRIPTION
  839. *
  840. * CHANGES
  841. *
  842. ******************************************************************************/
  843.  
  844. void Normal_Calc_HCompl_Reciprocal(VECTOR Result, int N_Max, FRACTAL *fractal)
  845. {
  846.   DBL n1, n2, n3, n4;
  847.   int i;
  848.   DBL x, y, z, w;
  849.   DBL xx, yy, zz, ww;
  850.   DBL xxx, yyy, zzz, www;
  851.  
  852.   /*
  853.    * Algebraic properties of hypercomplexes allows simplifications in
  854.    * computations...
  855.    */
  856.  
  857.   x = Sx[0];
  858.   y = Sy[0];
  859.   z = Sz[0];
  860.   w = Sw[0];
  861.  
  862.   for (i = 1; i < N_Max; ++i)
  863.   {
  864.     /******************* Case: z->1/z+c *****************/
  865.  
  866.     HReciprocal(&xx, &yy, &zz, &ww, Sx[i], Sy[i], Sz[i], Sw[i]);
  867.  
  868.     HSqr(xxx, yyy, zzz, www, xx, yy, zz, ww);
  869.  
  870.     HMult(xx, yy, zz, ww, x, y, z, w, -xxx, -yyy, -zzz, -www);
  871.  
  872.     x = xx;
  873.     y = yy;
  874.     z = zz;
  875.     w = ww;
  876.   }
  877.  
  878.   n1 = Sx[N_Max];
  879.   n2 = Sy[N_Max];
  880.   n3 = Sz[N_Max];
  881.   n4 = Sw[N_Max];
  882.  
  883.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  884.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  885.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  886. }
  887.  
  888.  
  889.  
  890. /*****************************************************************************
  891. *
  892. * FUNCTION
  893. *
  894. * INPUT
  895. *
  896. * OUTPUT
  897. *
  898. * RETURNS
  899. *
  900. * AUTHOR
  901. *
  902. *   Pascal Massimino
  903. *
  904. * DESCRIPTION
  905. *
  906. * CHANGES
  907. *
  908. ******************************************************************************/
  909.  
  910. int F_Bound_HCompl_Reciprocal(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL  *Depth_Max)
  911. {
  912.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  913. }
  914.  
  915. /*--------------------------- Func -------------------------------*/
  916.  
  917.  
  918. /*****************************************************************************
  919. *
  920. * FUNCTION
  921. *
  922. * INPUT
  923. *
  924. * OUTPUT
  925. *
  926. * RETURNS
  927. *
  928. * AUTHOR
  929. *
  930. *   Pascal Massimino
  931. *
  932. * DESCRIPTION
  933. *
  934. * CHANGES
  935. *
  936. ******************************************************************************/
  937.  
  938. int Iteration_HCompl_Func(VECTOR IPoint, FRACTAL *HCompl)
  939. {
  940.   int i;
  941.   DBL Norm, xx, yy, zz, ww;
  942.   DBL Exit_Value;
  943.   DBL x, y, z, w;
  944.  
  945.   x = Sx[0] = IPoint[X];
  946.   y = Sy[0] = IPoint[Y];
  947.   z = Sz[0] = IPoint[Y];
  948.   w = Sw[0] = (HCompl->SliceDist
  949.                   - HCompl->Slice[X]*x 
  950.                   - HCompl->Slice[Y]*y 
  951.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  952.  
  953.   Exit_Value = HCompl->Exit_Value;
  954.  
  955.   for (i = 1; i <= HCompl->n; ++i)
  956.   {
  957.     Norm = x * x + y * y + z * z + w * w;
  958.  
  959.     if (Norm > Exit_Value)
  960.     {
  961.       return (FALSE);
  962.     }
  963.  
  964.     HFunc(&xx, &yy, &zz, &ww, x, y, z, w, HCompl);
  965.  
  966.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  967.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  968.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  969.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  970.  
  971.   }
  972.  
  973.   return (TRUE);
  974. }
  975.  
  976.  
  977.  
  978. /*****************************************************************************
  979. *
  980. * FUNCTION
  981. *
  982. * INPUT
  983. *
  984. * OUTPUT
  985. *
  986. * RETURNS
  987. *
  988. * AUTHOR
  989. *
  990. *   Pascal Massimino
  991. *
  992. * DESCRIPTION
  993. *
  994. * CHANGES
  995. *
  996. ******************************************************************************/
  997.  
  998. int D_Iteration_HCompl_Func(VECTOR IPoint, FRACTAL *HCompl, DBL *Dist)
  999. {
  1000.   int i;
  1001.   DBL xx, yy, zz, ww;
  1002.   DBL Exit_Value, F_Value, Step;
  1003.   DBL x, y, z, w;
  1004.   VECTOR H_Normal;
  1005.  
  1006.   x = Sx[0] = IPoint[X];
  1007.   y = Sy[0] = IPoint[Y];
  1008.   z = Sz[0] = IPoint[Z];
  1009.   w = Sw[0] = (HCompl->SliceDist
  1010.                   - HCompl->Slice[X]*x 
  1011.                   - HCompl->Slice[Y]*y 
  1012.                   - HCompl->Slice[Z]*z)/HCompl->Slice[T]; 
  1013.  
  1014.   Exit_Value = HCompl->Exit_Value;
  1015.  
  1016.   for (i = 1; i <= HCompl->n; ++i)
  1017.   {
  1018.     F_Value = x * x + y * y + z * z + w * w;
  1019.  
  1020.     if (F_Value > Exit_Value)
  1021.     {
  1022.       Normal_Calc_HCompl_Func(H_Normal, i - 1, HCompl);
  1023.  
  1024.       VDot(Step, H_Normal, Direction);
  1025.  
  1026.       if (Step < -Fractal_Tolerance)
  1027.       {
  1028.         Step = -2.0 * Step;
  1029.  
  1030.         if ((F_Value > Precision * Step) && F_Value < (30 * Precision * Step))
  1031.         {
  1032.           *Dist = F_Value / Step;
  1033.  
  1034.           return (FALSE);
  1035.         }
  1036.       }
  1037.  
  1038.       *Dist = Precision;
  1039.  
  1040.       return (FALSE);
  1041.     }
  1042.  
  1043.     HFunc(&xx, &yy, &zz, &ww, x, y, z, w, HCompl);
  1044.  
  1045.     x = Sx[i] = xx + HCompl->Julia_Parm[X];
  1046.     y = Sy[i] = yy + HCompl->Julia_Parm[Y];
  1047.     z = Sz[i] = zz + HCompl->Julia_Parm[Z];
  1048.     w = Sw[i] = ww + HCompl->Julia_Parm[T];
  1049.  
  1050.   }
  1051.  
  1052.   *Dist = Precision;
  1053.  
  1054.   return (TRUE);
  1055. }
  1056.  
  1057.  
  1058.  
  1059. /*****************************************************************************
  1060. *
  1061. * FUNCTION
  1062. *
  1063. * INPUT
  1064. *
  1065. * OUTPUT
  1066. *
  1067. * RETURNS
  1068. *
  1069. * AUTHOR
  1070. *
  1071. *   Pascal Massimino
  1072. *
  1073. * DESCRIPTION
  1074. *
  1075. * CHANGES
  1076. *
  1077. ******************************************************************************/
  1078.  
  1079. void Normal_Calc_HCompl_Func(VECTOR Result, int N_Max, FRACTAL *fractal)
  1080. {
  1081.   DBL n1, n2, n3, n4;
  1082.   int i;
  1083.   DBL x, y, z, w;
  1084.   DBL xx, yy, zz, ww;
  1085.   DBL xxx, yyy, zzz, www;
  1086.  
  1087.   /*
  1088.    * Algebraic properties of hypercomplexes allows simplifications in
  1089.    * computations...
  1090.    */
  1091.  
  1092.   x = Sx[0];
  1093.   y = Sy[0];
  1094.   z = Sz[0];
  1095.   w = Sw[0];
  1096.  
  1097.   for (i = 1; i < N_Max; ++i)
  1098.   {
  1099.     /**************** Case: z-> f(z)+c ************************/
  1100.  
  1101.     HFunc(&xx, &yy, &zz, &ww, Sx[i], Sy[i], Sz[i], Sw[i], fractal);
  1102.  
  1103.     HMult(xxx, yyy, zzz, www, xx, yy, zz, ww, x, y, z, w);
  1104.  
  1105.     x = xxx;
  1106.     y = yyy;
  1107.     z = zzz;
  1108.     w = www;
  1109.   }
  1110.  
  1111.   n1 = Sx[N_Max];
  1112.   n2 = Sy[N_Max];
  1113.   n3 = Sz[N_Max];
  1114.   n4 = Sw[N_Max];
  1115.  
  1116.   Result[X] = x * n1 + y * n2 + z * n3 + w * n4;
  1117.   Result[Y] = -y * n1 + x * n2 - w * n3 + z * n4;
  1118.   Result[Z] = -z * n1 - w * n2 + x * n3 + y * n4;
  1119. }
  1120.  
  1121.  
  1122.  
  1123. /*****************************************************************************
  1124. *
  1125. * FUNCTION
  1126. *
  1127. * INPUT
  1128. *
  1129. * OUTPUT
  1130. *
  1131. * RETURNS
  1132. *
  1133. * AUTHOR
  1134. *
  1135. *   Pascal Massimino
  1136. *
  1137. * DESCRIPTION
  1138. *
  1139. * CHANGES
  1140. *
  1141. ******************************************************************************/
  1142.  
  1143. int F_Bound_HCompl_Func(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL  *Depth_Max)
  1144. {
  1145.   return F_Bound_HCompl(Ray, Fractal, Depth_Min, Depth_Max);
  1146. }
  1147.  
  1148. /*****************************************************************************
  1149. *
  1150. * FUNCTION  Complex transcental functions
  1151. *
  1152. * INPUT     pointer to source complex number
  1153. *
  1154. * OUTPUT    fn(input)
  1155. *
  1156. * RETURNS   void
  1157. *
  1158. * AUTHOR
  1159. *
  1160. *   Tim Wegner
  1161. *
  1162. * DESCRIPTION  Calculate common functions on complexes
  1163. *   Since our purpose is fractals, error checking is lax.            
  1164. *
  1165. * CHANGES
  1166. *
  1167. ******************************************************************************/
  1168.  
  1169. void Complex_Mult (CMPLX *target, CMPLX *source1, CMPLX *source2)
  1170. {
  1171.   DBL tmpx;
  1172.   tmpx = source1->x * source2->x - source1->y * source2->y;
  1173.   target->y = source1->x * source2->y + source1->y * source2->x;
  1174.   target->x = tmpx;
  1175. }
  1176.  
  1177. void Complex_Div (CMPLX *target, CMPLX *source1, CMPLX *source2)
  1178. {
  1179.   DBL mod,tmpx,yxmod,yymod;
  1180.   mod = Sqr(source2->x) + Sqr(source2->y);
  1181.   if (mod==0)
  1182.      return;
  1183.   yxmod = source2->x/mod;
  1184.   yymod = - source2->y/mod;
  1185.   tmpx = source1->x * yxmod - source1->y * yymod;
  1186.   target->y = source1->x * yymod + source1->y * yxmod;
  1187.   target->x = tmpx;
  1188. } /* End Complex_Mult() */
  1189.  
  1190. void Complex_Exp (CMPLX *target, CMPLX *source)
  1191. {
  1192.   DBL expx;
  1193.   expx = exp(source->x);
  1194.   target->x = expx * cos(source->y);
  1195.   target->y = expx * sin(source->y);
  1196. } /* End Complex_Exp() */
  1197.  
  1198. void Complex_Sin (CMPLX *target, CMPLX *source)
  1199. {
  1200.   target->x = sin(source->x) * cosh(source->y);
  1201.   target->y = cos(source->x) * sinh(source->y);
  1202. } /* End Complex_Sin() */
  1203.  
  1204. void Complex_Sinh (CMPLX *target, CMPLX *source)
  1205. {
  1206.   target->x = sinh(source->x) * cos(source->y);
  1207.   target->y = cosh(source->x) * sin(source->y);
  1208. } /* End Complex_Sinh() */
  1209.  
  1210.  
  1211. void Complex_Cos (CMPLX *target, CMPLX *source)
  1212. {
  1213.   target->x = cos(source->x) * cosh(source->y);
  1214.   target->y = -sin(source->x) * sinh(source->y);
  1215. } /* End Complex_Cos() */
  1216.  
  1217. void Complex_Cosh (CMPLX *target, CMPLX *source)
  1218. {
  1219.   target->x = cosh(source->x) * cos(source->y);
  1220.   target->y = sinh(source->x) * sin(source->y);
  1221. } /* End Complex_Cosh() */
  1222.  
  1223.  
  1224. void Complex_Log (CMPLX *target, CMPLX *source)
  1225. {
  1226.   DBL mod,zx,zy;
  1227.   mod = sqrt(source->x * source->x + source->y * source->y);
  1228.   zx = log(mod);
  1229.   zy = atan2(source->y,source->x);
  1230.  
  1231.   target->x = zx;
  1232.   target->y = zy;
  1233. } /* End Complex_Log() */
  1234.  
  1235. void Complex_Sqrt(CMPLX *target, CMPLX *source)
  1236. {
  1237.   DBL mag;
  1238.   DBL theta;
  1239.  
  1240.   if(source->x == 0.0 && source->y == 0.0)
  1241.   {
  1242.     target->x = target->y = 0.0;
  1243.   }
  1244.   else
  1245.   {   
  1246.     mag   = sqrt(sqrt(Sqr(source->x) + Sqr(source->y)));
  1247.     theta = atan2(source->y, source->x) / 2;
  1248.     target->y = mag * sin(theta);
  1249.     target->x = mag * cos(theta);
  1250.   }
  1251. } /* End Complex_Sqrt() */
  1252.  
  1253. /* rz=Arcsin(z)=-i*Log{i*z+sqrt(1-z*z)} */
  1254. void Complex_ASin(CMPLX *target, CMPLX *source)
  1255. {
  1256.   CMPLX tempz1,tempz2;
  1257.  
  1258.   Complex_Mult(&tempz1, source, source);
  1259.   tempz1.x = 1 - tempz1.x; tempz1.y = -tempz1.y; 
  1260.   Complex_Sqrt( &tempz1, &tempz1);
  1261.  
  1262.   tempz2.x = -source->y; tempz2.y = source->x;
  1263.   tempz1.x += tempz2.x;  tempz1.y += tempz2.y;
  1264.  
  1265.   Complex_Log( &tempz1, &tempz1);
  1266.   target->x = tempz1.y;  target->y = -tempz1.x;   
  1267. }   /* End Complex_ASin() */
  1268.  
  1269.  
  1270. void Complex_ACos(CMPLX *target, CMPLX *source)
  1271. {
  1272.   CMPLX temp;
  1273.  
  1274.   Complex_Mult(&temp, source, source);
  1275.   temp.x -= 1;
  1276.   Complex_Sqrt(&temp, &temp);
  1277.  
  1278.   temp.x += source->x; temp.y += source->y;
  1279.   
  1280.   Complex_Log(&temp, &temp);
  1281.   target->x = temp.y;  target->y = -temp.x; 
  1282. }   /* End Complex_ACos() */
  1283.  
  1284. void Complex_ASinh(CMPLX *target, CMPLX *source)
  1285. {
  1286.   CMPLX temp;
  1287.  
  1288.   Complex_Mult (&temp, source, source);
  1289.   temp.x += 1;   
  1290.   Complex_Sqrt (&temp, &temp);
  1291.   temp.x += source->x; temp.y += source->y; 
  1292.   Complex_Log(target, &temp);
  1293. }  /* End Complex_ASinh */
  1294.  
  1295. /* rz=Arccosh(z)=Log(z+sqrt(z*z-1)} */
  1296. void Complex_ACosh (CMPLX *target, CMPLX *source)
  1297. {
  1298.   CMPLX tempz;
  1299.   Complex_Mult(&tempz, source, source);
  1300.   tempz.x -= 1; 
  1301.   Complex_Sqrt (&tempz, &tempz);
  1302.   tempz.x = source->x + tempz.x; tempz.y = source->y + tempz.y;
  1303.   Complex_Log (target, &tempz);
  1304. }   /* End Complex_ACosh() */
  1305.  
  1306. /* rz=Arctanh(z)=1/2*Log{(1+z)/(1-z)} */
  1307. void Complex_ATanh(CMPLX *target, CMPLX *source)
  1308. {
  1309.   CMPLX temp0,temp1,temp2;
  1310.    
  1311.   if( source->x == 0.0)
  1312.   {
  1313.     target->x = 0;
  1314.     target->y = atan( source->y);
  1315.     return;
  1316.   }
  1317.   else
  1318.   {
  1319.     if( fabs(source->x) == 1.0 && source->y == 0.0)
  1320.     {
  1321.       return;
  1322.     }
  1323.     else if( fabs( source->x) < 1.0 && source->y == 0.0)
  1324.     {
  1325.       target->x = log((1+source->x)/(1-source->x))/2;
  1326.       target->y = 0;
  1327.       return;
  1328.     } 
  1329.     else
  1330.     {
  1331.       temp0.x = 1 + source->x; temp0.y = source->y;
  1332.       temp1.x = 1 - source->x; temp1.y = -source->y; 
  1333.       Complex_Div(&temp2, &temp0, &temp1);
  1334.       Complex_Log(&temp2, &temp2);
  1335.       target->x = .5 * temp2.x; target->y = .5 * temp2.y;
  1336.       return;
  1337.     }
  1338.   }
  1339. }   /* End Complex_ATanh() */
  1340.  
  1341. /* rz=Arctan(z)=i/2*Log{(1-i*z)/(1+i*z)} */
  1342. void Complex_ATan(CMPLX *target, CMPLX *source)
  1343. {
  1344.   CMPLX temp0,temp1,temp2,temp3;
  1345.   if( source->x == 0.0 && source->y == 0.0)
  1346.     target->x = target->y = 0;
  1347.   else if( source->x != 0.0 && source->y == 0.0){
  1348.     target->x = atan(source->x);
  1349.     target->y = 0;
  1350.   }
  1351.   else if( source->x == 0.0 && source->y != 0.0){
  1352.     temp0.x = source->y;  temp0.y = 0.0;
  1353.     Complex_ATanh(&temp0, &temp0);
  1354.     target->x = -temp0.y; target->y = temp0.x;
  1355.   } 
  1356.   else if( source->x != 0.0 && source->y != 0.0)
  1357.   {
  1358.     temp0.x = -source->y; temp0.y = source->x; 
  1359.     temp1.x = 1 - temp0.x; temp1.y = -temp0.y;   
  1360.     temp2.x = 1 + temp0.x; temp2.y = temp0.y; 
  1361.  
  1362.     Complex_Div(&temp3, &temp1, &temp2);
  1363.     Complex_Log(&temp3, &temp3);
  1364.     target->x = -temp3.y * .5; target->y = .5 * temp3.x; 
  1365.   }
  1366. }   /* End Complex_ATanz() */
  1367.  
  1368. void Complex_Tan(CMPLX *target, CMPLX *source)
  1369. {
  1370.   DBL x, y, sinx, cosx, sinhy, coshy, denom;
  1371.   x = 2 * source->x;
  1372.   y = 2 * source->y;
  1373.   sinx = sin(x); cosx = cos(x);
  1374.   sinhy = sinh(y); coshy = cosh(y);
  1375.   denom = cosx + coshy;
  1376.   if(denom == 0)
  1377.     return;
  1378.   target->x = sinx/denom;
  1379.   target->y = sinhy/denom;
  1380. }   /* End Complex_Tan() */
  1381.  
  1382. void Complex_Tanh(CMPLX *target, CMPLX *source)
  1383. {
  1384.   DBL x, y, siny, cosy, sinhx, coshx, denom;
  1385.   x = 2 * source->x;
  1386.   y = 2 * source->y;
  1387.   siny = sin(y); cosy = cos(y);
  1388.   sinhx = sinh(x); coshx = cosh(x);
  1389.   denom = coshx + cosy;
  1390.   if(denom == 0)
  1391.     return;
  1392.   target->x = sinhx/denom;
  1393.   target->y = siny/denom;
  1394. }   /* End Complex_Tanh() */
  1395.  
  1396.  
  1397. void Complex_Power (CMPLX *target, CMPLX *source1, CMPLX *source2)
  1398. {
  1399.   CMPLX cLog, t;
  1400.   DBL e2x;
  1401.  
  1402.   if(source1->x == 0 && source1->y == 0) 
  1403.   {
  1404.     target->x = target->y = 0.0;
  1405.     return;
  1406.   }
  1407.  
  1408.   Complex_Log (&cLog, source1);
  1409.   Complex_Mult (&t, &cLog, source2);
  1410.  
  1411.   if(t.x < -690)
  1412.     e2x = 0;
  1413.   else
  1414.     e2x = exp(t.x);
  1415.   target->x = e2x * cos(t.y);
  1416.   target->y = e2x * sin(t.y);
  1417. }
  1418.  
  1419. #if 1
  1420. void Complex_Pwr (CMPLX *target, CMPLX *source)
  1421. {
  1422.   Complex_Power(target, source, &exponent);
  1423. }
  1424.  
  1425. #else
  1426. void Complex_Pwr (CMPLX *target, CMPLX *source)
  1427. {
  1428.   /* the sqr dunction for testing */
  1429.   Complex_Mult(target, source, source);
  1430. }
  1431. #endif
  1432.