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  * Menger4D MOD1   from Syntopia & DarkBeam's Menger4 code from M3D
10  * @reference
11  * http://www.fractalforums.com/mandelbulb-3d/custom-formulas-and-transforms-release-t17106/
12  */
13 
14 #include "all_fractal_definitions.h"
15 
cFractalMenger4dMod2()16 cFractalMenger4dMod2::cFractalMenger4dMod2() : cAbstractFractal()
17 {
18 	nameInComboBox = "Menger 4D Mod2";
19 	internalName = "menger4d_mod2";
20 	internalID = fractal::menger4dMod2;
21 	DEType = analyticDEType;
22 	DEFunctionType = linearDEFunction;
23 	cpixelAddition = cpixelDisabledByDefault;
24 	defaultBailout = 100.0;
25 	DEAnalyticFunction = analyticFunctionIFS;
26 	coloringFunction = coloringFunctionDefault;
27 }
28 
FormulaCode(CVector4 & z,const sFractal * fractal,sExtendedAux & aux)29 void cFractalMenger4dMod2::FormulaCode(CVector4 &z, const sFractal *fractal, sExtendedAux &aux)
30 {
31 	if (fractal->transformCommon.functionEnabledFalse
32 			&& aux.i >= fractal->transformCommon.startIterationsD
33 			&& aux.i < fractal->transformCommon.stopIterationsD1)
34 	{
35 		if (fractal->transformCommon.functionEnabledAFalse)
36 		{
37 			z = CVector4(z.x, z.y, z.z, z.Length());
38 			aux.DE += 0.5;
39 		}
40 		if (fractal->transformCommon.functionEnabledBFalse)
41 		{
42 			z = CVector4(z.x + z.y + z.z, -z.x - z.y + z.z, -z.x + z.y - z.z, z.x - z.y - z.z);
43 			aux.DE *= z.Length() / aux.r;
44 		}
45 		z = fabs(z - fractal->transformCommon.offsetA0000);
46 	}
47 
48 		z = fabs(z);
49 		if (aux.i >= fractal->transformCommon.startIterationsC
50 				&& aux.i < fractal->transformCommon.stopIterationsC)
51 		{
52 			z += fractal->transformCommon.additionConstant0000; // offset
53 		}
54 
55 
56 		if (z.x - z.y < 0.0) swap(z.y, z.x);
57 		if (z.x - z.z < 0.0) swap(z.z, z.x);
58 		if (z.y - z.z < 0.0) swap(z.z, z.y);
59 		if (z.x - z.w < 0.0) swap(z.w, z.x);
60 		if (z.y - z.w < 0.0) swap(z.w, z.y);
61 		if (z.z - z.w < 0.0) swap(z.w, z.z);
62 
63 		// 6 plane rotation
64 		if (fractal->transformCommon.functionEnabledRFalse
65 				&& aux.i >= fractal->transformCommon.startIterationsR
66 				&& aux.i < fractal->transformCommon.stopIterationsR)
67 		{
68 			CVector4 tp;
69 			if (fractal->transformCommon.rotation44a.x != 0)
70 			{
71 				tp = z;
72 				double alpha = fractal->transformCommon.rotation44a.x * M_PI_180;
73 				z.x = tp.x * cos(alpha) + tp.y * sin(alpha);
74 				z.y = tp.x * -sin(alpha) + tp.y * cos(alpha);
75 			}
76 			if (fractal->transformCommon.rotation44a.y != 0)
77 			{
78 				tp = z;
79 				double beta = fractal->transformCommon.rotation44a.y * M_PI_180;
80 				z.y = tp.y * cos(beta) + tp.z * sin(beta);
81 				z.z = tp.y * -sin(beta) + tp.z * cos(beta);
82 			}
83 			if (fractal->transformCommon.rotation44a.z != 0)
84 			{
85 				tp = z;
86 				double gamma = fractal->transformCommon.rotation44a.z * M_PI_180;
87 				z.x = tp.x * cos(gamma) + tp.z * sin(gamma);
88 				z.z = tp.x * -sin(gamma) + tp.z * cos(gamma);
89 			}
90 			if (fractal->transformCommon.rotation44b.x != 0)
91 			{
92 				tp = z;
93 				double delta = fractal->transformCommon.rotation44b.x * M_PI_180;
94 				z.x = tp.x * cos(delta) + tp.w * sin(delta);
95 				z.w = tp.x * -sin(delta) + tp.w * cos(delta);
96 			}
97 			if (fractal->transformCommon.rotation44b.y != 0)
98 			{
99 				tp = z;
100 				double epsilon = fractal->transformCommon.rotation44b.y * M_PI_180;
101 				z.y = tp.y * cos(epsilon) + tp.w * sin(epsilon);
102 				z.w = tp.y * -sin(epsilon) + tp.w * cos(epsilon);
103 			}
104 			if (fractal->transformCommon.rotation44b.z != 0)
105 			{
106 				tp = z;
107 				double zeta = fractal->transformCommon.rotation44b.z * M_PI_180;
108 				z.z = tp.z * cos(zeta) + tp.w * sin(zeta);
109 				z.w = tp.z * -sin(zeta) + tp.w * cos(zeta);
110 			}
111 		}
112 
113 		CVector4 scaleM = fractal->transformCommon.scale3 * fractal->transformCommon.scale1111;
114 		CVector4 offsetM = fractal->transformCommon.additionConstant111d5;
115 		z.x = scaleM.x * z.x - offsetM.x;
116 		z.y = scaleM.y * z.y - offsetM.y;
117 		z.w = scaleM.w * z.w - offsetM.w;
118 		z.z *= fractal->transformCommon.scale1;
119 		double t = 0.5 * offsetM.z / scaleM.z;
120 		z.z -= t;
121 		z.z = -fabs(z.z);
122 		z.z += t;
123 		z.z *= fractal->transformCommon.scaleA3;
124 		aux.DE *= fractal->transformCommon.scale3;
125 
126 		if (fractal->transformCommon.functionEnabledSFalse
127 				&& aux.i >= fractal->transformCommon.startIterationsS
128 				&& aux.i < fractal->transformCommon.stopIterationsS)
129 		{
130 			double rr = z.Dot(z);
131 			// if (r2 < 1e-21 && r2 > -1e-21) r2 = (r2 > 0) ? 1e-21 : -1e-21;
132 
133 			if (rr < fractal->transformCommon.minR2p25)
134 			{
135 				z *= fractal->transformCommon.maxMinR2factor;
136 				aux.DE *= fractal->transformCommon.maxMinR2factor;
137 			}
138 			else if (rr < fractal->transformCommon.maxR2d1)
139 			{
140 				double tglad_factor2 = fractal->transformCommon.maxR2d1 / rr;
141 				z *= tglad_factor2;
142 				aux.DE *= tglad_factor2;
143 			}
144 		}
145 
146 		aux.DE = aux.DE * fractal->analyticDE.scale1 + fractal->analyticDE.offset0;
147 	}
148