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