gfx/povray/hcmplx.c

   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