Newton.frag

#version 330
#info 3D Newton Fractal (z^p - 1)
#define providesInit
#include "MathUtils.frag"
#include "DE-Raytracer.frag"

void init() { }

// Newton fractal power
const int power = 3;

#group Newton
uniform int iterations;  slider[10,23,50]
uniform int plane;  slider[0,1,1]
uniform float symmetry1;  slider[1,1,8]
uniform float symmetry2; slider[1,1,8]

vec2 complex_mul(vec2 a, vec2 b) { return vec2(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
vec2 complex_div(vec2 a, vec2 b) {
    float denom = b.x * b.x + b.y * b.y;
    return vec2(
        (a.x * b.x + a.y * b.y) / denom,
        (a.y * b.x - a.x * b.y) / denom
    );
}

float newton_sdf2d(vec2 z0) {
	vec2 zeros[power];
	vec2 zerop[power];

	for (int i = 0; i < power; i++) {
    	float angle = float(i) * 2.0 * PI / float(power);
    	zeros[i] = vec2(cos(angle), sin(angle));
		zerop[i] = vec2(1.0, 0.0);
	}

	vec2 z = z0;
	vec2 dz = vec2(1.0, 0.0);
	float eps = 0.001;

	for (int k = 0; k < iterations; k++) {
		for (int i = 0; i < power; i++) {
			float e = length(z - zeros[i]);
			if (e < eps) {
				float de = e * -log(e) / length(dz);
				return de;
			}
		}

		vec2 sz = vec2(0.0, 0.0);
		vec2 sz2 = vec2(0.0, 0.0);

		for (int i = 0; i < power; ++i) {
			vec2 d = z - zeros[i];
			sz += complex_div(zerop[i], d);
			sz2 += complex_div(zerop[i], complex_mul(d, d));
		}

		vec2 sp = vec2(0.0, 0.0);
		vec2 sp2 = vec2(0.0, 0.0);

		vec2 d = sz;
		z -= complex_div(vec2(1.0, 0.0), d);

		vec2 temp = complex_div(-sz2, complex_mul(d, d));
		dz = complex_mul(dz, temp + vec2(1.0, 0.0));
	}

	return 0.0;
}


float DE(vec3 p) {
	int plane_combo = plane;
	float angle1 = 3.14159265 / symmetry1;
	float angle2 = 3.14159265 / symmetry2;
	float theta1 = 0.0;
	float theta2 = 0.0;
	float r = 0.0;
	vec3 p_rotated1 = p;
	vec3 p_rotated2 = p;

	if (plane_combo == 1) {
		r = length(p.xy);

		theta1 = atan(p.x, p.y);
		theta1 = mod(theta1 + angle1, 2.0 * angle1) - angle1;
		p_rotated1.xy = vec2(cos(theta1), sin(theta1)) * r;

		theta2 = atan(p.x, p.y);
		theta2 = mod(theta2 + angle2, 2.0 * angle2) - angle2;
		p_rotated2.xy = vec2(cos(theta2), sin(theta2)) * r;

		float a = newton_sdf2d(p_rotated1.zx);
		float b = newton_sdf2d(p_rotated2.xy);
		return max(a, b);
	} else {
		r = length(p.xz);

		theta1 = atan(p.z, p.x);
		theta1 = mod(theta1 + angle1, 2.0 * angle1) - angle1;
		p_rotated1.xz = vec2(cos(theta1), sin(theta1)) * r;

		theta2 = atan(p.z, p.x);
		theta2 = mod(theta2 + angle2, 2.0 * angle2) - angle2;
		p_rotated2.xz = vec2(cos(theta2), sin(theta2)) * r;

		float a = newton_sdf2d(p_rotated1.xy);
		float b = newton_sdf2d(p_rotated2.xz);
		return max(a, b);
	}
}

#preset Default
FOV = 0.4
Eye = 3.1020507,0.558315865,-2.89073821
Target = 2.39299919,0.378348499,-2.20893347
Up = -0.217950261,0.968333682,0.028938637
EquiRectangular = false
AutoFocus = false
FocalPlane = 1
Aperture = 0
Gamma = 1.8636364
ToneMapping = 1
Exposure = 0.77830188
Brightness = 1.00490195
Contrast = 0.76744185
AvgLumin = 0.5,0.5,0.5
Saturation = 1.6262136
LumCoeff = 0.2125,0.7154,0.0721
Hue = 0
GaussianWeight = 1
AntiAliasScale = 2
DepthToAlpha = false
ShowDepth = false
DepthMagnitude = 1
Detail = -3.03333331
DetailAO = -1.22421516
FudgeFactor = 1
MaxDistance = 100
MaxRaySteps = 166
Dither = 0
NormalBackStep = 1 NotLocked
AO = 0,0,0,1
Specular = 0.67592593
SpecularExp = 5.641026
SpecularMax = 6.806283
SpotLight = 1,1,1,0.05442177
SpotLightDir = 0.11340208,0.1
CamLight = 1,1,1,0.97297298
CamLightMin = 0
Glow = 1,1,1,1
GlowMax = 100
Fog = 0
HardShadow = 0 NotLocked
ShadowSoft = 2
QualityShadows = false
Reflection = 0 NotLocked
DebugSun = false NotLocked
BaseColor = 1,0.792156863,0.0470588235
OrbitStrength = 0
X = 0.5,0.6,0.6,0.699999988
Y = 1,0.6,0,0.400000006
Z = 0.8,0.78,1,0.5
R = 0.4,0.7,1,0.119999997
BackgroundColor = 0.588235294,0.619607843,0.654901961
GradientBackground = 0.3
CycleColors = false
Cycles = 1.1
EnableFloor = false NotLocked
FloorNormal = 0,0,1
FloorHeight = 0
FloorColor = 1,1,1
iterations = 25
plane = 1
symmetry1 = 1
symmetry2 = 1
#endpreset