1 /**
2 * Mandelbulber v2, a 3D fractal generator _%}}i*<. ______
3 * Copyright (C) 2020 Mandelbulber Team _>]|=||i=i<, / ____/ __ __
4 * \><||i|=>>%) / / __/ /___/ /_
5 * This file is part of Mandelbulber. )<=i=]=|=i<> / /__ /_ __/_ __/
6 * The project is licensed under GPLv3, -<>>=|><|||` \____/ /_/ /_/
7 * see also COPYING file in this folder. ~+{i%+++
8 *
9 * quadratic iteration in real or imaginary scator algebra
10 * @reference
11 * http://www.fractalforums.com/new-theories-and-research/
12 * ix-possibly-the-holy-grail-fractal-%28in-fff-lore%29
13 * https://luz.izt.uam.mx/drupal/en/fractals/ix
14 * @author Manuel Fernandez-Guasti
15 * This formula contains aux.DE
16 */
17
18 #include "all_fractal_definitions.h"
19
cFractalScatorPower2()20 cFractalScatorPower2::cFractalScatorPower2() : cAbstractFractal()
21 {
22 nameInComboBox = "Scator Power 2";
23 internalName = "scator_power2";
24 internalID = fractal::scatorPower2;
25 DEType = analyticDEType;
26 DEFunctionType = customDEFunction;
27 cpixelAddition = cpixelEnabledByDefault;
28 defaultBailout = 100.0;
29 DEAnalyticFunction = analyticFunctionCustomDE;
30 coloringFunction = coloringFunctionDefault;
31 }
32
FormulaCode(CVector4 & z,const sFractal * fractal,sExtendedAux & aux)33 void cFractalScatorPower2::FormulaCode(CVector4 &z, const sFractal *fractal, sExtendedAux &aux)
34 {
35 // r calc
36 double r;
37 CVector4 zz = z * z;
38 if (fractal->transformCommon.functionEnabledXFalse)
39 {
40 r = aux.r;
41 }
42 else if (!fractal->transformCommon.functionEnabledYFalse)
43 {
44 r = sqrt(zz.x - zz.y - zz.z + (zz.y * zz.z) / zz.x);
45 }
46 else
47 { // this should be used for imaginary scators
48 r = sqrt(zz.x + zz.y + zz.z + (zz.y * zz.z) / zz.x);
49 }
50 // Scator real enabled by default
51 CVector4 newZ = z;
52 // double temp1;
53 if (!fractal->transformCommon.functionEnabledFalse)
54 { // scator real
55 newZ.x = zz.x + zz.y + zz.z;
56 newZ.y = z.x * z.y;
57 newZ.z = z.x * z.z;
58 newZ *= fractal->transformCommon.constantMultiplier122;
59 // temp1 = newZ.Length();
60 newZ.x += (zz.y * zz.z) / zz.x;
61 newZ.y *= (1.0 + zz.z / zz.x);
62 newZ.z *= (1.0 + zz.y / zz.x);
63 // r = sqrt(zz.x - zz.y - zz.z + (zz.y * zz.z) / zz.x);
64 }
65 else
66 { // scator imaginary
67 newZ.x = zz.x - zz.y - zz.z;
68 newZ.y = z.x * z.y;
69 newZ.z = z.x * z.z;
70 newZ *= fractal->transformCommon.constantMultiplier122;
71 // temp1 = newZ.Length();
72 newZ.x += (zz.y * zz.z) / zz.x;
73 newZ.y *= (1.0 - zz.z / zz.x);
74 newZ.z *= (1.0 - zz.y / zz.x);
75 // r = sqrt(zz.x + zz.y + zz.z + (zz.y * zz.z) / zz.x);
76 }
77 z = newZ;
78 // double temp2 = newZ.Length();
79 // temp2 = temp1 / temp2;
80 /* aux.DE = aux.DE * 2.0 * aux.r;
81 double newx = z.x * z.x - z.y * z.y - z.z * z.z;
82 double newy = 2.0 * z.x * z.y;
83 double newz = 2.0 * z.x * z.z;
84 z.x = newx;
85 z.y = newy;
86 z.z = newz; */
87
88 // addCpixel
89 if (fractal->transformCommon.addCpixelEnabledFalse
90 && aux.i >= fractal->transformCommon.startIterationsE
91 && aux.i < fractal->transformCommon.stopIterationsE)
92 {
93 CVector4 c = aux.const_c;
94 CVector4 tempC = c;
95 if (fractal->transformCommon.alternateEnabledFalse) // alternate
96 {
97 tempC = aux.c;
98 switch (fractal->mandelbulbMulti.orderOfXYZ)
99 {
100 case multi_OrderOfXYZ_xyz:
101 default: tempC = CVector4(tempC.x, tempC.y, tempC.z, tempC.w); break;
102 case multi_OrderOfXYZ_xzy: tempC = CVector4(tempC.x, tempC.z, tempC.y, tempC.w); break;
103 case multi_OrderOfXYZ_yxz: tempC = CVector4(tempC.y, tempC.x, tempC.z, tempC.w); break;
104 case multi_OrderOfXYZ_yzx: tempC = CVector4(tempC.y, tempC.z, tempC.x, tempC.w); break;
105 case multi_OrderOfXYZ_zxy: tempC = CVector4(tempC.z, tempC.x, tempC.y, tempC.w); break;
106 case multi_OrderOfXYZ_zyx: tempC = CVector4(tempC.z, tempC.y, tempC.x, tempC.w); break;
107 }
108 aux.c = tempC;
109 }
110 else
111 {
112 switch (fractal->mandelbulbMulti.orderOfXYZ)
113 {
114 case multi_OrderOfXYZ_xyz:
115 default: tempC = CVector4(c.x, c.y, c.z, c.w); break;
116 case multi_OrderOfXYZ_xzy: tempC = CVector4(c.x, c.z, c.y, c.w); break;
117 case multi_OrderOfXYZ_yxz: tempC = CVector4(c.y, c.x, c.z, c.w); break;
118 case multi_OrderOfXYZ_yzx: tempC = CVector4(c.y, c.z, c.x, c.w); break;
119 case multi_OrderOfXYZ_zxy: tempC = CVector4(c.z, c.x, c.y, c.w); break;
120 case multi_OrderOfXYZ_zyx: tempC = CVector4(c.z, c.y, c.x, c.w); break;
121 }
122 }
123 z += tempC * fractal->transformCommon.constantMultiplier111;
124 }
125
126 // analytic DE calc
127 if (fractal->analyticDE.enabled)
128 {
129 if (!fractal->analyticDE.enabledFalse)
130 {
131 aux.DE = 2.0 * r * aux.DE * fractal->analyticDE.scale1 + fractal->analyticDE.offset1;
132 }
133 else
134 { // vec3
135 // rd calc
136 double rd;
137 zz = z * z;
138 if (fractal->transformCommon.functionEnabledXFalse)
139 {
140 rd = z.Length();
141 }
142 else if (!fractal->transformCommon.functionEnabledYFalse)
143 {
144 rd = sqrt(zz.x - zz.y - zz.z + (zz.y * zz.z) / zz.x);
145 }
146 else
147 {
148 rd = sqrt(zz.x + zz.y + zz.z + (zz.y * zz.z) / zz.x);
149 }
150 double vecDE = fractal->transformCommon.scaleA1 * rd / r;
151 aux.DE =
152 max(r * 2.0, vecDE) * aux.DE * fractal->analyticDE.scale1 + fractal->analyticDE.offset1;
153 }
154 aux.dist = 0.5 * r * log(r) / aux.DE;
155 }
156 // force bailout
157 // double tp = fractal->transformCommon.scale1;
158 // double tpz = fabs(z.z);
159 // if (tpz < -tp && tpz > tp) tpz = (tpz > 0) ? z.z += 100000.0 : z.z -= 100000.0;
160 if (r > fractal->transformCommon.scale2) z.z += 100000.0;
161 aux.DE0 = r; // temp for testing
162 }
163