gfx/povray/fractal.c

   1 /*****************************************************************************
   2 *                fractal.c
   3 *
   4 *  This module implements the fractal sets primitive.
   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 "bbox.h"
  31 #include "matrices.h"
  32 #include "objects.h"
  33 #include "spheres.h"
  34 #include "fractal.h"
  35 #include "quatern.h"
  36 #include "hcmplx.h"
  37
  38
  39
  40 /*****************************************************************************
  41 * Local preprocessor defines
  42 ******************************************************************************/
  43
  44 #ifndef Fractal_Tolerance
  45 #define Fractal_Tolerance 1e-7
  46 #endif
  47
  48
  49
  50 /*****************************************************************************
  51 * Local typedefs
  52 ******************************************************************************/
  53
  54
  55
  56 /*****************************************************************************
  57 * Static functions
  58 ******************************************************************************/
  59
  60 static int All_Fractal_Intersections (OBJECT * Object, RAY * Ray, ISTACK * Depth_Stack);
  61 static int Inside_Fractal (VECTOR IPoint, OBJECT * Object);
  62 static void Fractal_Normal (VECTOR Result, OBJECT * Object, INTERSECTION * Intersect);
  63 static FRACTAL *Copy_Fractal (OBJECT * Object);
  64 static void Translate_Fractal (OBJECT * Object, VECTOR Vector, TRANSFORM *Trans);
  65 static void Rotate_Fractal (OBJECT * Object, VECTOR Vector, TRANSFORM *Trans);
  66 static void Scale_Fractal (OBJECT * Object, VECTOR Vector, TRANSFORM *Trans);
  67 static void Transform_Fractal (OBJECT * Object, TRANSFORM * Trans);
  68 static void Invert_Fractal (OBJECT * Object);
  69 static void Destroy_Fractal (OBJECT * Object);
  70 static void Compute_Fractal_BBox (FRACTAL * Fractal);
  71
  72 /*****************************************************************************
  73 * Local variables
  74 ******************************************************************************/
  75
  76 static METHODS Fractal_Methods =
  77 {
  78   All_Fractal_Intersections,
  79   Inside_Fractal, Fractal_Normal,
  80   (COPY_METHOD)Copy_Fractal,
  81   Translate_Fractal, Rotate_Fractal,
  82   Scale_Fractal, Transform_Fractal, Invert_Fractal,
  83   Destroy_Fractal
  84 };
  85
  86 static int Allocated_Iteration_Stack_Length = 0;
  87
  88 DBL *Sx = NULL, *Sy = NULL, *Sz = NULL, *Sw = NULL;
  89 DBL Precision;
  90 VECTOR Direction;
  91
  92 static COMPLEX_FUNCTION_METHOD Complex_Function_List[] =
  93 {
  94   /* must match STYPE list in fractal.h */
  95   Complex_Exp,
  96   Complex_Log,
  97   Complex_Sin,
  98   Complex_ASin,
  99   Complex_Cos,
 100   Complex_ACos,
 101   Complex_Tan,
 102   Complex_ATan,
 103   Complex_Sinh,
 104   Complex_ASinh,
 105   Complex_Cosh,
 106   Complex_ACosh,
 107   Complex_Tanh,
 108   Complex_ATanh,
 109   Complex_Pwr
 110 };
 111
 112 /*****************************************************************************
 113 *
 114 * FUNCTION
 115 *
 116 * INPUT
 117 *
 118 * OUTPUT
 119 *
 120 * RETURNS
 121 *
 122 * AUTHOR
 123 *
 124 *   Pascal Massimino
 125 *
 126 * DESCRIPTION
 127 *
 128 *   -
 129 *
 130 * CHANGES
 131 *
 132 *   Dec 1994 : Creation.
 133 *
 134 ******************************************************************************/
 135
 136 static int All_Fractal_Intersections(OBJECT *Object, RAY *Ray, ISTACK *Depth_Stack)
 137 {
 138   int Intersection_Found;
 139   int Last = 0;
 140   int CURRENT, NEXT;
 141   DBL Depth, Depth_Max;
 142   DBL Dist, Dist_Next, Len;
 143
 144   VECTOR IPoint, Mid_Point, Next_Point, Real_Pt;
 145   VECTOR Real_Normal, F_Normal;
 146   RAY New_Ray;
 147   FRACTAL *Fractal = (FRACTAL *) Object;
 148
 149   Increase_Counter(stats[Ray_Fractal_Tests]);
 150
 151   Intersection_Found = FALSE;
 152   Precision = Fractal->Precision;
 153
 154   /* Get into Fractal's world. */
 155
 156   if (Fractal->Trans != NULL)
 157   {
 158     MInvTransDirection(Direction, Ray->Direction, Fractal->Trans);
 159     VDot(Len, Direction, Direction);
 160
 161     if (Len == 0.0)
 162     {
 163       return (FALSE);
 164     }
 165
 166     if (Len != 1.0)
 167     {
 168       Len = 1.0 / sqrt(Len);
 169       VScaleEq(Direction, Len);
 170     }
 171
 172     Assign_Vector(New_Ray.Direction, Direction);
 173     MInvTransPoint(New_Ray.Initial, Ray->Initial, Fractal->Trans);
 174   }
 175   else
 176   {
 177     Assign_Vector(Direction, Ray->Direction);
 178     New_Ray = *Ray;
 179     Len = 1.0;
 180   }
 181
 182   /* Bound fractal. */
 183
 184   if (!F_Bound(&New_Ray, Fractal, &Depth, &Depth_Max))
 185   {
 186     return (FALSE);
 187   }
 188
 189   if (Depth_Max < Fractal_Tolerance)
 190   {
 191     return (FALSE);
 192   }
 193
 194   if (Depth < Fractal_Tolerance)
 195   {
 196     Depth = Fractal_Tolerance;
 197   }
 198
 199   /* Jump to starting point */
 200
 201   VScale(Next_Point, Direction, Depth);
 202   VAddEq(Next_Point, New_Ray.Initial);
 203
 204   CURRENT = D_Iteration(Next_Point, Fractal, &Dist);
 205
 206   /* Light ray starting inside ? */
 207
 208   if (CURRENT)
 209   {
 210     VAddScaledEq(Next_Point, 2.0 * Fractal_Tolerance, Direction);
 211
 212     Depth += 2.0 * Fractal_Tolerance;
 213
 214     if (Depth > Depth_Max)
 215     {
 216       return (FALSE);
 217     }
 218
 219     CURRENT = D_Iteration(Next_Point, Fractal, &Dist);
 220   }
 221
 222   /* Ok. Trace it */
 223
 224   while (Depth < Depth_Max)
 225   {
 226     /*
 227      * Get close to the root: Advance with Next_Point, keeping track of last
 228      * position in IPoint...
 229      */
 230
 231     while (1)
 232     {
 233       if (Dist < Precision)
 234       {
 235         Dist = Precision;
 236       }
 237
 238       Depth += Dist;
 239
 240       if (Depth > Depth_Max)
 241       {
 242         if (Intersection_Found)
 243         {
 244           Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
 245         }
 246
 247         return (Intersection_Found);
 248       }
 249
 250       Assign_Vector(IPoint, Next_Point);
 251       VAddScaledEq(Next_Point, Dist, Direction);
 252
 253       NEXT = D_Iteration(Next_Point, Fractal, &Dist_Next);
 254
 255       if (NEXT != CURRENT)
 256       {
 257         /* Set surface was crossed... */
 258
 259         Depth -= Dist;
 260         break;
 261       }
 262       else
 263       {
 264         Dist = Dist_Next; /* not reached */
 265       }
 266     }
 267
 268     /* then, polish the root via bisection method... */
 269
 270     while (Dist > Fractal_Tolerance)
 271     {
 272       Dist *= 0.5;
 273       VAddScaled(Mid_Point, IPoint, Dist, Direction);
 274
 275       Last = Iteration(Mid_Point, Fractal);
 276
 277       if (Last == CURRENT)
 278       {
 279         Assign_Vector(IPoint, Mid_Point);
 280
 281         Depth += Dist;
 282
 283         if (Depth > Depth_Max)
 284         {
 285           if (Intersection_Found)
 286           {
 287             Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
 288           }
 289
 290           return (Intersection_Found);
 291         }
 292       }
 293     }
 294
 295     if (CURRENT == FALSE) /* Mid_Point isn't inside the set */
 296     {
 297       VAddScaledEq(IPoint, Dist, Direction);
 298
 299       Depth += Dist;
 300
 301       Iteration(IPoint, Fractal);
 302     }
 303     else
 304     {
 305       if (Last != CURRENT)
 306       {
 307         Iteration(IPoint, Fractal);
 308       }
 309     }
 310
 311     if (Fractal->Trans != NULL)
 312     {
 313       MTransPoint(Real_Pt, IPoint, Fractal->Trans);
 314       Normal_Calc(Fractal, F_Normal);
 315       MTransNormal(Real_Normal, F_Normal, Fractal->Trans);
 316     }
 317     else
 318     {
 319       Assign_Vector(Real_Pt, IPoint);
 320       Normal_Calc(Fractal, Real_Normal);
 321     }
 322
 323     if (Point_In_Clip(Real_Pt, Object->Clip))
 324     {
 325       VNormalize(Real_Normal, Real_Normal);
 326       push_normal_entry(Depth * Len, Real_Pt, Real_Normal, Object, Depth_Stack);
 327       Intersection_Found = TRUE;
 328
 329       /* If fractal isn't used with CSG we can exit now. */
 330
 331       if (!(Fractal->Type & IS_CHILD_OBJECT))
 332       {
 333         break;
 334       }
 335     }
 336
 337     /* Start over where work was left */
 338
 339     Assign_Vector(IPoint, Next_Point);
 340     Dist = Dist_Next;
 341     CURRENT = NEXT;
 342
 343   }
 344
 345   if (Intersection_Found)
 346   {
 347     Increase_Counter(stats[Ray_Fractal_Tests_Succeeded]);
 348   }
 349   return (Intersection_Found);
 350 }
 351
 352 /*****************************************************************************
 353 *
 354 * FUNCTION
 355 *
 356 * INPUT
 357 *
 358 * OUTPUT
 359 *
 360 * RETURNS
 361 *
 362 * AUTHOR
 363 *
 364 *   Pascal Massimino
 365 *
 366 * DESCRIPTION
 367 *
 368 *   -
 369 *
 370 * CHANGES
 371 *
 372 *   Dec 1994 : Creation.
 373 *
 374 ******************************************************************************/
 375
 376 static int Inside_Fractal(VECTOR IPoint, OBJECT *Object)
 377 {
 378   FRACTAL *Fractal = (FRACTAL *) Object;
 379   int Result;
 380   static VECTOR New_Point;
 381
 382   if (Fractal->Trans != NULL)
 383   {
 384     MInvTransPoint(New_Point, IPoint, Fractal->Trans);
 385
 386     Result = Iteration(New_Point, Fractal);
 387   }
 388   else
 389   {
 390     Result = Iteration(IPoint, Fractal);
 391   }
 392
 393   if (Fractal->Inverted)
 394   {
 395     return (!Result);
 396   }
 397   else
 398   {
 399     return (Result);
 400   }
 401 }
 402
 403 /*****************************************************************************
 404 *
 405 * FUNCTION
 406 *
 407 * INPUT
 408 *
 409 * OUTPUT
 410 *
 411 * RETURNS
 412 *
 413 * AUTHOR
 414 *
 415 *   Pascal Massimino
 416 *
 417 * DESCRIPTION
 418 *
 419 *   -
 420 *
 421 * CHANGES
 422 *
 423 *   Dec 1994 : Creation.
 424 *
 425 ******************************************************************************/
 426
 427 static void Fractal_Normal(VECTOR Result, OBJECT *Object, INTERSECTION *Intersect)
 428 {
 429   Assign_Vector(Result, Intersect->INormal);
 430 }
 431
 432 /*****************************************************************************
 433 *
 434 * FUNCTION
 435 *
 436 * INPUT
 437 *
 438 * OUTPUT
 439 *
 440 * RETURNS
 441 *
 442 * AUTHOR
 443 *
 444 *   Pascal Massimino
 445 *
 446 * DESCRIPTION
 447 *
 448 *   -
 449 *
 450 * CHANGES
 451 *
 452 *   Dec 1994 : Creation.
 453 *
 454 ******************************************************************************/
 455
 456 static void Translate_Fractal(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 457 {
 458   Transform_Fractal(Object, Trans);
 459 }
 460
 461 /*****************************************************************************
 462 *
 463 * FUNCTION
 464 *
 465 * INPUT
 466 *
 467 * OUTPUT
 468 *
 469 * RETURNS
 470 *
 471 * AUTHOR
 472 *
 473 *   Pascal Massimino
 474 *
 475 * DESCRIPTION
 476 *
 477 *   -
 478 *
 479 * CHANGES
 480 *
 481 *   Dec 1994 : Creation.
 482 *
 483 ******************************************************************************/
 484
 485 static void Rotate_Fractal(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 486 {
 487   Transform_Fractal(Object, Trans);
 488 }
 489
 490 /*****************************************************************************
 491 *
 492 * FUNCTION
 493 *
 494 * INPUT
 495 *
 496 * OUTPUT
 497 *
 498 * RETURNS
 499 *
 500 * AUTHOR
 501 *
 502 *   Pascal Massimino
 503 *
 504 * DESCRIPTION
 505 *
 506 *   -
 507 *
 508 * CHANGES
 509 *
 510 *   Dec 1994 : Creation.
 511 *
 512 ******************************************************************************/
 513
 514 static void Scale_Fractal(OBJECT *Object, VECTOR Vector, TRANSFORM *Trans)
 515 {
 516   Transform_Fractal(Object, Trans);
 517 }
 518
 519 /*****************************************************************************
 520 *
 521 * FUNCTION
 522 *
 523 * INPUT
 524 *
 525 * OUTPUT
 526 *
 527 * RETURNS
 528 *
 529 * AUTHOR
 530 *
 531 *   Pascal Massimino
 532 *
 533 * DESCRIPTION
 534 *
 535 *   -
 536 *
 537 * CHANGES
 538 *
 539 *   Dec 1994 : Creation.
 540 *   Mar 1996 : Moved call to Recompute_BBox to Compute_Fractal_BBox() (TW)
 541 *
 542 ******************************************************************************/
 543
 544 static void Transform_Fractal(OBJECT *Object, TRANSFORM *Trans)
 545 {
 546   FRACTAL *Fractal = (FRACTAL *) Object;
 547
 548   if (((FRACTAL *) Object)->Trans == NULL)
 549   {
 550     ((FRACTAL *) Object)->Trans = Create_Transform();
 551   }
 552
 553   Compose_Transforms(Fractal->Trans, Trans);
 554
 555   Compute_Fractal_BBox((FRACTAL *) Object);
 556 }
 557
 558 /*****************************************************************************
 559 *
 560 * FUNCTION
 561 *
 562 * INPUT
 563 *
 564 * OUTPUT
 565 *
 566 * RETURNS
 567 *
 568 * AUTHOR
 569 *
 570 *   Pascal Massimino
 571 *
 572 * DESCRIPTION
 573 *
 574 *   -
 575 *
 576 * CHANGES
 577 *
 578 *   Dec 1994 : Creation.
 579 *
 580 ******************************************************************************/
 581
 582 static void Invert_Fractal(OBJECT *Object)
 583 {
 584   ((FRACTAL *) Object)->Inverted ^= TRUE;
 585 }
 586
 587 /*****************************************************************************
 588 *
 589 * FUNCTION
 590 *
 591 * INPUT
 592 *
 593 * OUTPUT
 594 *
 595 * RETURNS
 596 *
 597 * AUTHOR
 598 *
 599 *   Pascal Massimino
 600 *
 601 * DESCRIPTION
 602 *
 603 *   -
 604 *
 605 * CHANGES
 606 *
 607 *   Dec 1994 : Creation.
 608 *   Mar 1996 : Added call to recompute_BBox() to bottom (TW)
 609 *
 610 ******************************************************************************/
 611
 612 static void Compute_Fractal_BBox(FRACTAL *Fractal)
 613 {
 614   DBL R;
 615
 616   switch (Fractal->Algebra)
 617   {
 618     case QUATERNION_TYPE:
 619
 620       R = 1.0 + sqrt(Sqr(Fractal->Julia_Parm[X]) + Sqr(Fractal->Julia_Parm[Y]) + Sqr(Fractal->Julia_Parm[Z]) + Sqr(Fractal->Julia_Parm[T]));
 621       R += Fractal_Tolerance; /* fix bug when Julia_Parameter exactly 0 */
 622
 623       if (R > 2.0)
 624       {
 625         R = 2.0;
 626       }
 627
 628       Fractal->Exit_Value = Sqr(R);
 629
 630       break;
 631
 632     case HYPERCOMPLEX_TYPE:
 633     default:
 634
 635       R = 4.0;
 636
 637       Fractal->Exit_Value = 16.0;
 638
 639       break;
 640   }
 641
 642   Fractal->Radius_Squared = Sqr(R);
 643
 644   Fractal->Inverted = FALSE;
 645
 646   Make_BBox(Fractal->BBox, -R, -R, -R, 2.0 * R, 2.0 * R, 2.0 * R);
 647
 648   Recompute_BBox(&Fractal->BBox, Fractal->Trans);
 649 }
 650
 651 /*****************************************************************************
 652 *
 653 * FUNCTION
 654 *
 655 * INPUT
 656 *
 657 * OUTPUT
 658 *
 659 * RETURNS
 660 *
 661 * AUTHOR
 662 *
 663 *   Pascal Massimino
 664 *
 665 * DESCRIPTION
 666 *
 667 *   -
 668 *
 669 * CHANGES
 670 *
 671 *   Dec 1994 : Creation.
 672 *
 673 ******************************************************************************/
 674
 675 FRACTAL *Create_Fractal()
 676 {
 677   FRACTAL *New;
 678
 679   New = (FRACTAL *) POV_MALLOC(sizeof(FRACTAL), "Fractal Set");
 680
 681   INIT_OBJECT_FIELDS(New, BASIC_OBJECT, &Fractal_Methods)
 682
 683   New->Trans = NULL;
 684
 685   Make_Vector(New->Center, 0.0, 0.0, 0.0);
 686
 687   New->Julia_Parm[X] = 1.0;
 688   New->Julia_Parm[Y] = 0.0;
 689   New->Julia_Parm[Z] = 0.0;
 690   New->Julia_Parm[T] = 0.0;
 691
 692   New->Slice[X] = 0.0;
 693   New->Slice[Y] = 0.0;
 694   New->Slice[Z] = 0.0;
 695   New->Slice[T] = 1.0;
 696   New->SliceDist = 0.0;
 697
 698   New->Exit_Value = 4.0;
 699
 700   New->n = 20;
 701
 702   New->Precision = 1.0 / 20.0;
 703
 704   New->Inverted = FALSE;
 705
 706   New->Algebra = QUATERNION_TYPE;
 707
 708   New->Sub_Type = SQR_STYPE;
 709
 710   New->Bound = NULL;
 711
 712   New->Clip = NULL;
 713
 714   New->Normal_Calc_Method = NULL;
 715   New->Iteration_Method   = NULL;
 716   New->D_Iteration_Method = NULL;
 717   New->F_Bound_Method     = NULL;
 718   New->Complex_Function_Method = NULL;
 719
 720   New->Radius_Squared = 0.0;
 721   New->exponent.x = 0.0;
 722   New->exponent.y = 0.0;
 723
 724   return (New);
 725 }
 726
 727 /*****************************************************************************
 728 *
 729 * FUNCTION
 730 *
 731 * INPUT
 732 *
 733 * OUTPUT
 734 *
 735 * RETURNS
 736 *
 737 * AUTHOR
 738 *
 739 *   Pascal Massimino
 740 *
 741 * DESCRIPTION
 742 *
 743 *   -
 744 *
 745 * CHANGES
 746 *
 747 *   Dec 1994 : Creation.
 748 *
 749 ******************************************************************************/
 750
 751 static FRACTAL *Copy_Fractal(OBJECT *Object)
 752 {
 753   FRACTAL *New;
 754
 755   New = Create_Fractal();
 756
 757   *New = *((FRACTAL *) Object);
 758
 759   New->Trans = Copy_Transform(((FRACTAL *) Object)->Trans);
 760
 761   return (New);
 762 }
 763
 764 /*****************************************************************************
 765 *
 766 * FUNCTION
 767 *
 768 * INPUT
 769 *
 770 * OUTPUT
 771 *
 772 * RETURNS
 773 *
 774 * AUTHOR
 775 *
 776 *   Pascal Massimino
 777 *
 778 * DESCRIPTION
 779 *
 780 *   -
 781 *
 782 * CHANGES
 783 *
 784 *   Dec 1994 : Creation.
 785 *
 786 ******************************************************************************/
 787
 788 static void Destroy_Fractal(OBJECT *Object)
 789 {
 790   Destroy_Transform(((FRACTAL *) Object)->Trans);
 791   POV_FREE(Object);
 792 }
 793
 794 /*****************************************************************************
 795 *
 796 * FUNCTION
 797 *
 798 * INPUT
 799 *
 800 * OUTPUT
 801 *
 802 * RETURNS
 803 *
 804 * AUTHOR
 805 *
 806 *   Pascal Massimino
 807 *
 808 * DESCRIPTION
 809 *
 810 *   -
 811 *
 812 * CHANGES
 813 *
 814 *   Dec 1994 : Creation.
 815 *
 816 ******************************************************************************/
 817
 818 void SetUp_Fractal(FRACTAL *Fractal)
 819 {
 820   switch (Fractal->Algebra)
 821   {
 822     case QUATERNION_TYPE:
 823
 824       switch(Fractal->Sub_Type)
 825       {
 826         case CUBE_STYPE:
 827           Fractal->Iteration_Method = Iteration_z3;
 828           Fractal->Normal_Calc_Method = Normal_Calc_z3;
 829           Fractal->D_Iteration_Method = D_Iteration_z3;
 830           break;
 831         case SQR_STYPE:
 832           Fractal->Iteration_Method = Iteration_Julia;
 833           Fractal->D_Iteration_Method = D_Iteration_Julia;
 834           Fractal->Normal_Calc_Method = Normal_Calc_Julia;
 835           break;
 836         default:
 837           Error("illegal function: quaternion only supports sqr and cube");
 838       }
 839       Fractal->F_Bound_Method = F_Bound_Julia;
 840
 841       break;
 842
 843     case HYPERCOMPLEX_TYPE:
 844
 845       switch (Fractal->Sub_Type)
 846       {
 847         case RECIPROCAL_STYPE:
 848
 849           Fractal->Iteration_Method = Iteration_HCompl_Reciprocal;
 850           Fractal->Normal_Calc_Method = Normal_Calc_HCompl_Reciprocal;
 851           Fractal->D_Iteration_Method = D_Iteration_HCompl_Reciprocal;
 852           Fractal->F_Bound_Method = F_Bound_HCompl_Reciprocal;
 853
 854           break;
 855
 856         case EXP_STYPE:
 857         case LOG_STYPE:
 858         case SIN_STYPE:
 859         case ASIN_STYPE:
 860         case COS_STYPE:
 861         case ACOS_STYPE:
 862         case TAN_STYPE:
 863         case ATAN_STYPE:
 864         case SINH_STYPE:
 865         case ASINH_STYPE:
 866         case COSH_STYPE:
 867         case ACOSH_STYPE:
 868         case TANH_STYPE:
 869         case ATANH_STYPE:
 870         case PWR_STYPE:
 871
 872           Fractal->Iteration_Method = Iteration_HCompl_Func;
 873           Fractal->Normal_Calc_Method = Normal_Calc_HCompl_Func;
 874           Fractal->D_Iteration_Method = D_Iteration_HCompl_Func;
 875           Fractal->F_Bound_Method = F_Bound_HCompl_Func;
 876           Fractal->Complex_Function_Method = Complex_Function_List[Fractal->Sub_Type];
 877
 878           break;
 879
 880         case CUBE_STYPE:
 881
 882           Fractal->Iteration_Method = Iteration_HCompl_z3;
 883           Fractal->Normal_Calc_Method = Normal_Calc_HCompl_z3;
 884           Fractal->D_Iteration_Method = D_Iteration_HCompl_z3;
 885           Fractal->F_Bound_Method = F_Bound_HCompl_z3;
 886
 887           break;
 888
 889         default:  /* SQR_STYPE or else... */
 890
 891           Fractal->Iteration_Method = Iteration_HCompl;
 892           Fractal->Normal_Calc_Method = Normal_Calc_HCompl;
 893           Fractal->D_Iteration_Method = D_Iteration_HCompl;
 894           Fractal->F_Bound_Method = F_Bound_HCompl;
 895
 896           break;
 897       }
 898
 899       break;
 900
 901     default:
 902
 903       Error("Algebra unknown in fractal.");
 904   }
 905
 906   Allocate_Iteration_Stack(Fractal->n);
 907
 908   Compute_Fractal_BBox(Fractal);
 909 }
 910
 911 /*****************************************************************************
 912 *
 913 * FUNCTION
 914 *
 915 * INPUT
 916 *
 917 * OUTPUT
 918 *
 919 * RETURNS
 920 *
 921 * AUTHOR
 922 *
 923 *   Pascal Massimino
 924 *
 925 * DESCRIPTION
 926 *
 927 *   -
 928 *
 929 * CHANGES
 930 *
 931 *   Dec 1994 : Creation.
 932 *
 933 ******************************************************************************/
 934
 935 void Allocate_Iteration_Stack(int n)
 936 {
 937   if (n > Allocated_Iteration_Stack_Length)
 938   {
 939     Sx = (DBL *) POV_REALLOC(Sx, (n + 1) * sizeof(DBL), "x iteration stack");
 940     Sy = (DBL *) POV_REALLOC(Sy, (n + 1) * sizeof(DBL), "y iteration stack");
 941     Sz = (DBL *) POV_REALLOC(Sz, (n + 1) * sizeof(DBL), "w iteration stack");
 942     Sw = (DBL *) POV_REALLOC(Sw, (n + 1) * sizeof(DBL), "z iteration stack");
 943
 944     Allocated_Iteration_Stack_Length = n;
 945   }
 946 }
 947
 948 /*****************************************************************************
 949 *
 950 * FUNCTION
 951 *
 952 * INPUT
 953 *
 954 * OUTPUT
 955 *
 956 * RETURNS
 957 *
 958 * AUTHOR
 959 *
 960 *   Pascal Massimino
 961 *
 962 * DESCRIPTION
 963 *
 964 *   -
 965 *
 966 * CHANGES
 967 *
 968 *   Dec 1994 : Creation.
 969 *
 970 ******************************************************************************/
 971
 972 void Free_Iteration_Stack()
 973 {
 974   if (Sx != NULL)
 975   {
 976     POV_FREE(Sx);
 977     POV_FREE(Sy);
 978     POV_FREE(Sz);
 979     POV_FREE(Sw);
 980   }
 981
 982   Sx = NULL;
 983   Sy = NULL;
 984   Sz = NULL;
 985   Sw = NULL;
 986
 987   Allocated_Iteration_Stack_Length = 0;
 988 }
 989
 990