1 /*
2 -----------------------------------------------------------------------------
3 This source file is part of OGRE
4 (Object-oriented Graphics Rendering Engine)
5 For the latest info, see http://www.ogre3d.org/
6
7 Copyright (c) 2000-2013 Torus Knot Software Ltd
8 Also see acknowledgements in Readme.html
9
10 You may use this sample code for anything you like, it is not covered by the
11 same license as the rest of the engine.
12 -----------------------------------------------------------------------------
13 */
14 /* Original author: John C. Hare
15 * Date: June 26, 1996 (although I've had these around for at least 6 years)
16 * Adapted to C++ by W.J. van der Laan 2004
17 */
18 #ifndef H_Q_Julia
19 #define H_Q_Julia
20 #include <cmath>
21
22 /**
23 * Simple, fast, inline quaternion math functions
24 */
25 struct Quat {
26 float r,i,j,k;
27 };
28
qadd(Quat & a,const Quat & b)29 inline void qadd(Quat &a, const Quat &b)
30 {
31 a.r += b.r;
32 a.i += b.i;
33 a.j += b.j;
34 a.k += b.k;
35 }
36
qmult(Quat & c,const Quat & a,const Quat & b)37 inline void qmult(Quat &c, const Quat &a, const Quat &b)
38 {
39 c.r = a.r*b.r - a.i*b.i - a.j*b.j - a.k*b.k;
40 c.i = a.r*b.i + a.i*b.r + a.j*b.k - a.k*b.j;
41 c.j = a.r*b.j + a.j*b.r + a.k*b.i - a.i*b.k;
42 c.k = a.r*b.k + a.k*b.r + a.i*b.j - a.j*b.i;
43 }
44
qsqr(Quat & b,const Quat & a)45 inline void qsqr(Quat &b, const Quat &a)
46 {
47 b.r = a.r*a.r - a.i*a.i - a.j*a.j - a.k*a.k;
48 b.i = 2.0f*a.r*a.i;
49 b.j = 2.0f*a.r*a.j;
50 b.k = 2.0f*a.r*a.k;
51 }
52
53 /**
54 * Implicit function that evaluates the Julia set.
55 */
56 class Julia {
57 private:
58 float global_real, global_imag, global_theta;
59 Quat oc,c,eio,emio;
60 public:
61 Julia(float global_real, float global_imag, float global_theta);
eval(float x,float y,float z)62 inline float eval(float x, float y, float z) {
63 Quat q, temp;
64 int i;
65
66 q.r = x;
67 q.i = y;
68 q.j = z;
69 q.k = 0.0;
70
71 for (i = 30; i > 0; i--) {
72 qsqr(temp, q);
73 qmult(q, emio, temp);
74 qadd(q, c);
75
76 if (q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k > 8.0)
77 break;
78 }
79
80 return((float)i);
81 }
82 };
83
84
Julia(float in_global_real,float in_global_imag,float in_global_theta)85 Julia::Julia(float in_global_real, float in_global_imag, float in_global_theta):
86 global_real(in_global_real), global_imag(in_global_imag), global_theta(in_global_theta) {
87
88 oc.r = global_real;
89 oc.i = global_imag;
90 oc.j = oc.k = 0.0;
91
92 eio.r = cos(global_theta);
93 eio.i = sin(global_theta);
94 eio.j = 0.0;
95 eio.k = 0.0;
96
97 emio.r = cos(-global_theta);
98 emio.i = sin(-global_theta);
99 emio.j = 0.0;
100 emio.k = 0.0;
101
102 /***
103 *** multiply eio*c only once at the beginning of iteration
104 *** (since q |-> sqrt(eio*(q-eio*c)))
105 *** q -> e-io*q^2 - eio*c
106 ***/
107
108 qmult(c, eio,oc);
109
110 }
111
112
113 #endif
114