1 /*****************************************************************************
2 * quatern.cpp
3 *
4 * This module implements Quaternion algebra julia fractals.
5 *
6 * This file was written by Pascal Massimino.
7 * Revised and updated for POV-Ray 3.x by Tim Wegner
8 *
9 * from Persistence of Vision(tm) Ray Tracer version 3.6.
10 * Copyright 1991-2003 Persistence of Vision Team
11 * Copyright 2003-2004 Persistence of Vision Raytracer Pty. Ltd.
12 *---------------------------------------------------------------------------
13 * NOTICE: This source code file is provided so that users may experiment
14 * with enhancements to POV-Ray and to port the software to platforms other
15 * than those supported by the POV-Ray developers. There are strict rules
16 * regarding how you are permitted to use this file. These rules are contained
17 * in the distribution and derivative versions licenses which should have been
18 * provided with this file.
19 *
20 * These licences may be found online, linked from the end-user license
21 * agreement that is located at http://www.povray.org/povlegal.html
22 *---------------------------------------------------------------------------
23 * This program is based on the popular DKB raytracer version 2.12.
24 * DKBTrace was originally written by David K. Buck.
25 * DKBTrace Ver 2.0-2.12 were written by David K. Buck & Aaron A. Collins.
26 *---------------------------------------------------------------------------
27 * $File: //depot/povray/3.6-release/source/quatern.cpp $
28 * $Revision: #3 $
29 * $Change: 3032 $
30 * $DateTime: 2004/08/02 18:43:41 $
31 * $Author: chrisc $
32 * $Log$
33 *****************************************************************************/
34
35 #include "frame.h"
36 #include "povray.h"
37 #include "vector.h"
38 #include "fractal.h"
39 #include "quatern.h"
40 #include "spheres.h"
41
42 BEGIN_POV_NAMESPACE
43
44 /*****************************************************************************
45 * Local preprocessor defines
46 ******************************************************************************/
47
48 #define Deriv_z2(n1,n2,n3,n4) \
49 { \
50 tmp = (n1)*x - (n2)*y - (n3)*z - (n4)*w; \
51 (n2) = (n1)*y + x*(n2); \
52 (n3) = (n1)*z + x*(n3); \
53 (n4) = (n1)*w + x*(n4); \
54 (n1) = tmp; \
55 }
56
57 #define Deriv_z3(n1,n2,n3,n4) \
58 { \
59 dtmp = 2.0*((n2)*y + (n3)*z + (n4)*w); \
60 dtmp2 = 6.0*x*(n1) - dtmp; \
61 (n1) = ( (n1)*x3 - x*dtmp )*3.0; \
62 (n2) = (n2)*x4 + y*dtmp2; \
63 (n3) = (n3)*x4 + z*dtmp2; \
64 (n4) = (n4)*x4 + w*dtmp2; \
65 }
66
67
68 /*****************************************************************************
69 * Local typedefs
70 ******************************************************************************/
71
72
73
74 /*****************************************************************************
75 * Static functions
76 ******************************************************************************/
77
78
79
80 /*****************************************************************************
81 * Local variables
82 ******************************************************************************/
83
84
85
86 /*****************************************************************************
87 *
88 * FUNCTION
89 *
90 * INPUT
91 *
92 * OUTPUT
93 *
94 * RETURNS
95 *
96 * AUTHOR
97 *
98 * Pascal Massimino
99 *
100 * DESCRIPTION
101 *
102 * -
103 *
104 * CHANGES
105 *
106 * Dec 1994 : Creation.
107 *
108 ******************************************************************************/
109
Iteration_z3(VECTOR point,FRACTAL * Julia)110 int Iteration_z3(VECTOR point, FRACTAL *Julia)
111 {
112 int i;
113 DBL x, y, z, w;
114 DBL d, x2, tmp;
115 DBL Exit_Value;
116
117 Sx[0] = x = point[X];
118 Sy[0] = y = point[Y];
119 Sz[0] = z = point[Z];
120 Sw[0] = w = (Julia->SliceDist
121 - Julia->Slice[X]*x
122 - Julia->Slice[Y]*y
123 - Julia->Slice[Z]*z)/Julia->Slice[T];
124
125 Exit_Value = Julia->Exit_Value;
126
127 for (i = 1; i <= Julia->n; ++i)
128 {
129 d = y * y + z * z + w * w;
130
131 x2 = x * x;
132
133 if ((d + x2) > Exit_Value)
134 {
135 return (false);
136 }
137
138 tmp = 3.0 * x2 - d;
139
140 Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
141 Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
142 Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
143 Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
144 }
145
146 return (true);
147 }
148
149
150
151 /*****************************************************************************
152 *
153 * FUNCTION
154 *
155 * INPUT
156 *
157 * OUTPUT
158 *
159 * RETURNS
160 *
161 * AUTHOR
162 *
163 * Pascal Massimino
164 *
165 * DESCRIPTION
166 *
167 * -
168 *
169 * CHANGES
170 *
171 * Dec 1994 : Creation.
172 *
173 ******************************************************************************/
174
Iteration_Julia(VECTOR point,FRACTAL * Julia)175 int Iteration_Julia(VECTOR point, FRACTAL *Julia)
176 {
177 int i;
178 DBL x, y, z, w;
179 DBL d, x2;
180 DBL Exit_Value;
181
182 Sx[0] = x = point[X];
183 Sy[0] = y = point[Y];
184 Sz[0] = z = point[Z];
185 Sw[0] = w = (Julia->SliceDist
186 - Julia->Slice[X]*x
187 - Julia->Slice[Y]*y
188 - Julia->Slice[Z]*z)/Julia->Slice[T];
189
190 Exit_Value = Julia->Exit_Value;
191
192 for (i = 1; i <= Julia->n; ++i)
193 {
194 d = y * y + z * z + w * w;
195
196 x2 = x * x;
197
198 if ((d + x2) > Exit_Value)
199 {
200 return (false);
201 }
202
203 x *= 2.0;
204
205 Sy[i] = y = x * y + Julia->Julia_Parm[Y];
206 Sz[i] = z = x * z + Julia->Julia_Parm[Z];
207 Sw[i] = w = x * w + Julia->Julia_Parm[T];
208 Sx[i] = x = x2 - d + Julia->Julia_Parm[X];;
209
210 }
211
212 return (true);
213 }
214
215
216
217 /*****************************************************************************
218 *
219 * FUNCTION
220 *
221 * INPUT
222 *
223 * OUTPUT
224 *
225 * RETURNS
226 *
227 * AUTHOR
228 *
229 * Pascal Massimino
230 *
231 * DESCRIPTION
232 *
233 * D_Iteration puts in *Dist a lower bound for the distance from *point to the
234 * set
235 *
236 * CHANGES
237 *
238 * Dec 1994 : Creation.
239 *
240 ******************************************************************************/
241
242 /*----------- Distance estimator + iterations ------------*/
243
D_Iteration_z3(VECTOR point,FRACTAL * Julia,DBL * Dist)244 int D_Iteration_z3(VECTOR point, FRACTAL *Julia, DBL *Dist)
245 {
246 int i, j;
247 DBL Norm, d;
248 DBL xx, yy, zz;
249 DBL x, y, z, w;
250 DBL tmp, x2;
251 DBL Exit_Value;
252 DBL Pow;
253
254 x = Sx[0] = point[X];
255 y = Sy[0] = point[Y];
256 z = Sz[0] = point[Z];
257 w = Sw[0] = (Julia->SliceDist
258 - Julia->Slice[X]*x
259 - Julia->Slice[Y]*y
260 - Julia->Slice[Z]*z)/Julia->Slice[T];
261
262 Exit_Value = Julia->Exit_Value;
263
264 for (i = 1; i <= Julia->n; i++)
265 {
266 d = y * y + z * z + w * w;
267
268 x2 = x * x;
269
270 if ((Norm = d + x2) > Exit_Value)
271 {
272 /* Distance estimator */
273
274 x = Sx[0];
275 y = Sy[0];
276 z = Sz[0];
277 w = Sw[0];
278
279 Pow = 1.0 / 3.0;
280
281 for (j = 1; j < i; ++j)
282 {
283 xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
284 yy = x * Sy[j] + y * Sx[j] - z * Sw[j] + w * Sz[j];
285 zz = x * Sz[j] + y * Sw[j] + z * Sx[j] - w * Sy[j];
286 w = x * Sw[j] - y * Sz[j] + z * Sy[j] + w * Sx[j];
287
288 x = xx;
289 y = yy;
290 z = zz;
291
292 Pow /= 3.0;
293 }
294
295 *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
296
297 return (false);
298 }
299
300 tmp = 3.0 * x2 - d;
301
302 Sx[i] = x = x * (x2 - 3.0 * d) + Julia->Julia_Parm[X];
303 Sy[i] = y = y * tmp + Julia->Julia_Parm[Y];
304 Sz[i] = z = z * tmp + Julia->Julia_Parm[Z];
305 Sw[i] = w = w * tmp + Julia->Julia_Parm[T];
306 }
307
308 *Dist = Precision;
309
310 return (true);
311 }
312
313
314
315 /*****************************************************************************
316 *
317 * FUNCTION
318 *
319 * INPUT
320 *
321 * OUTPUT
322 *
323 * RETURNS
324 *
325 * AUTHOR
326 *
327 * Pascal Massimino
328 *
329 * DESCRIPTION
330 *
331 * -
332 *
333 * CHANGES
334 *
335 * Dec 1994 : Creation.
336 *
337 ******************************************************************************/
338
D_Iteration_Julia(VECTOR point,FRACTAL * Julia,DBL * Dist)339 int D_Iteration_Julia(VECTOR point, FRACTAL *Julia, DBL *Dist)
340 {
341 int i, j;
342 DBL Norm, d;
343 DBL Exit_Value;
344 DBL x, y, z, w;
345 DBL xx, yy, zz, x2;
346 DBL Pow;
347
348 x = Sx[0] = point[X];
349 y = Sy[0] = point[Y];
350 z = Sz[0] = point[Z];
351 w = Sw[0] = (Julia->SliceDist
352 - Julia->Slice[X]*x
353 - Julia->Slice[Y]*y
354 - Julia->Slice[Z]*z)/Julia->Slice[T];
355
356 Exit_Value = Julia->Exit_Value;
357
358 for (i = 1; i <= Julia->n; i++)
359 {
360 d = y * y + z * z + w * w;
361
362 x2 = x * x;
363
364 if ((Norm = d + x2) > Exit_Value)
365 {
366 /* Distance estimator */
367
368 x = Sx[0];
369 y = Sy[0];
370 z = Sz[0];
371 w = Sw[0];
372
373 Pow = 1.0 / 2.0;
374
375 for (j = 1; j < i; ++j)
376 {
377 xx = x * Sx[j] - y * Sy[j] - z * Sz[j] - w * Sw[j];
378 yy = x * Sy[j] + y * Sx[j] + w * Sz[j] - z * Sw[j];
379 zz = x * Sz[j] + z * Sx[j] + y * Sw[j] - w * Sy[j];
380 w = x * Sw[j] + w * Sx[j] + z * Sy[j] - y * Sz[j];
381
382 x = xx;
383 y = yy;
384 z = zz;
385
386 Pow /= 2.0;
387 }
388
389 *Dist = Pow * sqrt(Norm / (x * x + y * y + z * z + w * w)) * log(Norm);
390
391 return (false);
392 }
393
394 x *= 2.0;
395
396 Sy[i] = y = x * y + Julia->Julia_Parm[Y];
397 Sz[i] = z = x * z + Julia->Julia_Parm[Z];
398 Sw[i] = w = x * w + Julia->Julia_Parm[T];
399 Sx[i] = x = x2 - d + Julia->Julia_Parm[X];
400
401 }
402
403 *Dist = Precision;
404
405 return (true);
406 }
407
408 /*****************************************************************************
409 *
410 * FUNCTION
411 *
412 * INPUT
413 *
414 * OUTPUT
415 *
416 * RETURNS
417 *
418 * AUTHOR
419 *
420 * Pascal Massimino
421 *
422 * DESCRIPTION
423 *
424 * Provided the iterations sequence has been built, perform the computation of
425 * the Normal
426 *
427 * CHANGES
428 *
429 * Dec 1994 : Creation.
430 *
431 ******************************************************************************/
432
Normal_Calc_z3(VECTOR Result,int N_Max,FRACTAL *)433 void Normal_Calc_z3(VECTOR Result, int N_Max, FRACTAL *)
434 {
435 DBL
436 n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
437 n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
438 n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
439
440 DBL x, y, z, w;
441 int i;
442 DBL tmp, dtmp, dtmp2, x2, x3, x4;
443
444 x = Sx[0];
445 y = Sy[0];
446 z = Sz[0];
447 w = Sw[0];
448
449 for (i = 1; i <= N_Max; i++)
450 {
451 tmp = y * y + z * z + w * w;
452
453 x2 = x * x;
454 x3 = x2 - tmp;
455 x4 = 3.0 * x2 - tmp;
456
457 Deriv_z3(n11, n12, n13, n14);
458 Deriv_z3(n21, n22, n23, n24);
459 Deriv_z3(n31, n32, n33, n34);
460
461 x = Sx[i];
462 y = Sy[i];
463 z = Sz[i];
464 w = Sw[i];
465 }
466
467 Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
468 Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
469 Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
470 }
471
472
473
474 /*****************************************************************************
475 *
476 * FUNCTION
477 *
478 * INPUT
479 *
480 * OUTPUT
481 *
482 * RETURNS
483 *
484 * AUTHOR
485 *
486 * Pascal Massimino
487 *
488 * DESCRIPTION
489 *
490 * -
491 *
492 * CHANGES
493 *
494 * Dec 1994 : Creation.
495 *
496 ******************************************************************************/
497
Normal_Calc_Julia(VECTOR Result,int N_Max,FRACTAL *)498 void Normal_Calc_Julia(VECTOR Result, int N_Max, FRACTAL *)
499 {
500 DBL
501 n11 = 1.0, n12 = 0.0, n13 = 0.0, n14 = 0.0,
502 n21 = 0.0, n22 = 1.0, n23 = 0.0, n24 = 0.0,
503 n31 = 0.0, n32 = 0.0, n33 = 1.0, n34 = 0.0;
504 DBL tmp;
505 DBL x, y, z, w;
506 int i;
507
508 x = Sx[0];
509 y = Sy[0];
510 z = Sz[0];
511 w = Sw[0];
512
513 for (i = 1; i <= N_Max; i++)
514 {
515 Deriv_z2(n11, n12, n13, n14);
516 Deriv_z2(n21, n22, n23, n24);
517 Deriv_z2(n31, n32, n33, n34);
518
519 x = Sx[i];
520 y = Sy[i];
521 z = Sz[i];
522 w = Sw[i];
523 }
524
525 Result[X] = n11 * x + n12 * y + n13 * z + n14 * w;
526 Result[Y] = n21 * x + n22 * y + n23 * z + n24 * w;
527 Result[Z] = n31 * x + n32 * y + n33 * z + n34 * w;
528 }
529
530
531
532 /*****************************************************************************
533 *
534 * FUNCTION
535 *
536 * INPUT
537 *
538 * OUTPUT
539 *
540 * RETURNS
541 *
542 * AUTHOR
543 *
544 * Pascal Massimino
545 *
546 * DESCRIPTION
547 *
548 * -
549 *
550 * CHANGES
551 *
552 * Dec 1994 : Creation.
553 *
554 ******************************************************************************/
555
F_Bound_Julia(RAY * Ray,FRACTAL * Fractal,DBL * Depth_Min,DBL * Depth_Max)556 int F_Bound_Julia(RAY *Ray, FRACTAL *Fractal, DBL *Depth_Min, DBL *Depth_Max)
557 {
558 return (Intersect_Sphere(Ray, Fractal->Center, Fractal->Radius_Squared, Depth_Min, Depth_Max));
559 }
560
561 END_POV_NAMESPACE
562