Data update

This commit is contained in:
Ingy döt Net 2023-07-22 21:42:12 -07:00
parent 9817d9b99b
commit 07c7092a52
140 changed files with 2712 additions and 241 deletions

View file

@ -4,39 +4,32 @@ gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
d, h = 800, 500 # pixel density (= image width) and image height
n, r = 200, 500 # number of iterations and escape radius (r > 2)
direction, height = 45, 1.5 # direction and height of the incoming light
stripes, damping = 4, 2.0 # stripe density and damping parameter
x = range(0, 2, length=d+1)
y = range(0, 2 * h / d, length=h+1)
A, B = collect(x) .- 1, collect(y) .- h / d
C = (2.0 + 1.0im) .* (A' .+ B .* im) .- 0.5
C = 2.0 .* (A' .+ B .* im) .- 0.5
Z, dZ, ddZ = zero(C), zero(C), zero(C)
Z, dZ = zero(C), zero(C)
D, S, T = zeros(size(C)), zeros(size(C)), zeros(size(C))
for k in 1:n
M = abs.(Z) .< r
S[M], T[M] = S[M] .+ cos.(stripes .* angle.(Z[M])), T[M] .+ 1
Z[M], dZ[M], ddZ[M] = Z[M] .^ 2 .+ C[M], 2 .* Z[M] .* dZ[M] .+ 1, 2 .* (dZ[M] .^ 2 .+ Z[M] .* ddZ[M])
S[M], T[M] = S[M] .+ exp.(.- abs.(Z[M])), T[M] .+ 1
Z[M], dZ[M] = Z[M] .^ 2 .+ C[M], 2 .* Z[M] .* dZ[M] .+ 1
end
N = abs.(Z) .>= r # normal map effect 1 (potential function)
P, Q = S[N] ./ T[N], (S[N] .+ cos.(stripes .* angle.(Z[N]))) ./ (T[N] .+ 1)
F = log2.(log.(abs.(Z[N])) ./ log(r)) # fraction between 0 and 1 (for interpolation)
R = Q .+ (P .- Q) .* F .* F .* (3 .- 2 .* F) # hermite interpolation (r is between q and p)
U, H = Z[N] ./ dZ[N], 1 .+ R ./ damping # normal vectors to the equipotential lines and height perturbation
U, v = U ./ abs.(U), exp(direction / 180 * pi * im) # unit normal vectors and vector in light direction
D[N] = max.((real.(U) .* real(v) .+ imag.(U) .* imag(v) .+ H .* height) ./ (1 + height), 0)
heatmap(S .^ 0.1, c=:balance)
savefig("Mandelbrot_set_1.png")
heatmap(D .^ 1.0, c=:bone_1)
savefig("Mandelbrot_normal_map_1.png")
N = abs.(Z) .>= r # normalized iteration count
T[N] = T[N] .- log2.(log.(abs.(Z[N])) ./ log(r))
N = abs.(Z) .>= r # normal map effect 2 (distance estimation)
U = Z[N] .* dZ[N] .* ((1 .+ log.(abs.(Z[N]))) .* conj.(dZ[N] .^ 2) .- log.(abs.(Z[N])) .* conj.(Z[N] .* ddZ[N]))
U, v = U ./ abs.(U), exp(direction / 180 * pi * im) # unit normal vectors and vector in light direction
D[N] = max.((real.(U) .* real(v) .+ imag.(U) .* imag(v) .+ height) ./ (1 + height), 0)
heatmap(T .^ 0.1, c=:balance)
savefig("Mandelbrot_set_2.png")
heatmap(D .^ 1.0, c=:afmhot)
savefig("Mandelbrot_normal_map_2.png")
N = abs.(Z) .> 2 # exterior distance estimation
D[N] = 0.5 .* log.(abs.(Z[N])) .* abs.(Z[N]) ./ abs.(dZ[N])
heatmap(D .^ 0.1, c=:balance)
savefig("Mandelbrot_set_3.png")

View file

@ -1,39 +1,41 @@
using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
d, h = 200, 1200 # pixel density (= image width) and image height
n, r = 8000, 10000 # number of iterations and escape radius (r > 2)
d, h = 800, 500 # pixel density (= image width) and image height
n, r = 200, 500 # number of iterations and escape radius (r > 2)
a = -.743643887037158704752191506114774 # real coordinate of the zoom center
b = 0.131825904205311970493132056385139 # imaginary coordinate of the center
direction, height = 45, 1.5 # direction and height of the incoming light
stripes, damping = 4.0, 2.0 # stripe density and damping parameter
x = range(0, 2, length=d+1)
y = range(0, 2 * h / d, length=h+1)
A, B = collect(x) .* pi, collect(y) .* pi
C = 8.0 .* exp.((A' .+ B .* im) .* im) .+ (a + b * im)
A, B = collect(x) .- 1, collect(y) .- h / d
C = (2.0 + 1.0im) .* (A' .+ B .* im) .- 0.5
Z, dZ = zero(C), zero(C)
D = zeros(size(C))
Z, dZ, ddZ = zero(C), zero(C), zero(C)
D, S, T = zeros(size(C)), zeros(size(C)), zeros(size(C))
for k in 1:n
M = abs2.(Z) .< abs2(r)
Z[M], dZ[M] = Z[M] .^ 2 .+ C[M], 2 .* Z[M] .* dZ[M] .+ 1
M = abs.(Z) .< r
S[M], T[M] = S[M] .+ cos.(stripes .* angle.(Z[M])), T[M] .+ 1
Z[M], dZ[M], ddZ[M] = Z[M] .^ 2 .+ C[M], 2 .* Z[M] .* dZ[M] .+ 1, 2 .* (dZ[M] .^ 2 .+ Z[M] .* ddZ[M])
end
N = abs.(Z) .> 2 # exterior distance estimation
D[N] = log.(abs.(Z[N])) .* abs.(Z[N]) ./ abs.(dZ[N])
N = abs.(Z) .>= r # normal map effect 1 (equipotential lines)
P, Q = S[N] ./ T[N], (S[N] .+ cos.(stripes .* angle.(Z[N]))) ./ (T[N] .+ 1)
R = Q .+ (P .- Q) .* log2.(log.(abs.(Z[N])) ./ log(r)) # linear interpolation
U, V = Z[N] ./ dZ[N], 1 .+ R ./ damping # normal vectors and variations in inclination
U, v = U ./ abs.(U), exp(direction / 180 * pi * im) # unit vectors
D[N] = max.((real.(U) .* real(v) .+ imag.(U) .* imag(v) .+ V .* height) ./ (1 + height), 0)
heatmap(D' .^ 0.05, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_map.png")
heatmap(D .^ 1.0, c=:bone_1)
savefig("Mandelbrot_normal_map_1.png")
X, Y = real(C), imag(C) # zoom images (adjust circle size 50 and zoom level 20 as needed)
R, c, z = 50 * (2 / d) * pi .* exp.(.- B), min(d, h) + 1, max(0, h - d) ÷ 20
N = abs.(Z) .> 2 # normal map effect 2 (equidistant lines)
U = Z[N] .* dZ[N] .* ((1 .+ log.(abs.(Z[N]))) .* conj.(dZ[N] .^ 2) .- log.(abs.(Z[N])) .* conj.(Z[N] .* ddZ[N]))
U, v = U ./ abs.(U), exp(direction / 180 * pi * im) # unit vectors
D[N] = max.((real.(U) .* real(v) .+ imag.(U) .* imag(v) .+ height) ./ (1 + height), 0)
gr(c=:nipy_spectral, axis=true, ticks=true, legend=false, markerstrokewidth=0)
p1 = scatter(X[1z+1:1z+c,1:d], Y[1z+1:1z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[1z+1:1z+c,1:d].^0.5)
p2 = scatter(X[2z+1:2z+c,1:d], Y[2z+1:2z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[2z+1:2z+c,1:d].^0.4)
p3 = scatter(X[3z+1:3z+c,1:d], Y[3z+1:3z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[3z+1:3z+c,1:d].^0.3)
p4 = scatter(X[4z+1:4z+c,1:d], Y[4z+1:4z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[4z+1:4z+c,1:d].^0.2)
plot(p1, p2, p3, p4, layout=(2, 2))
savefig("Mercator_Mandelbrot_zoom.png")
heatmap(D .^ 1.0, c=:afmhot)
savefig("Mandelbrot_normal_map_2.png")

View file

@ -1,46 +1,39 @@
using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
setprecision(BigFloat, 256) # set precision to 256 bits (default)
setrounding(BigFloat, RoundNearest) # set rounding mode (default)
d, h = 200, 1200 # pixel density (= image width) and image height
n, r = 8000, 10000 # number of iterations and escape radius (r > 2)
d, h = 50, 1000 # pixel density (= image width) and image height
n, r = 80000, 100000 # number of iterations and escape radius (r > 2)
a = BigFloat("-1.256827152259138864846434197797294538253477389787308085590211144291")
b = BigFloat(".37933802890364143684096784819544060002129071484943239316486643285025")
S = zeros(Complex{Float64}, n+1)
let c = a + b * im, z = zero(c)
for k in 1:n+1
S[k] = z
if abs2(z) < abs2(r)
z = z ^ 2 + c
else
println("The reference sequence diverges within $(k-1) iterations.")
break
end
end
end
a = -.743643887037158704752191506114774 # real coordinate of the zoom center
b = 0.131825904205311970493132056385139 # imaginary coordinate of the center
x = range(0, 2, length=d+1)
y = range(0, 2 * h / d, length=h+1)
A, B = collect(x) .* pi, collect(y) .* pi
C = 8.0 .* exp.((A' .+ B .* im) .* im)
C = 8.0 .* exp.((A' .+ B .* im) .* im) .+ (a + b * im)
E, Z, dZ = zero(C), zero(C), zero(C)
D, I, J = zeros(size(C)), ones(Int64, size(C)), ones(Int64, size(C))
Z, dZ = zero(C), zero(C)
D = zeros(size(C))
for k in 1:n
M, R = abs2.(Z) .< abs2(r), abs2.(Z) .< abs2.(E)
E[R], I[R] = Z[R], J[R] # rebase when z is closer to zero
E[M], I[M] = (2 .* S[I[M]] .+ E[M]) .* E[M] .+ C[M], I[M] .+ 1
Z[M], dZ[M] = S[I[M]] .+ E[M], 2 .* Z[M] .* dZ[M] .+ 1
M = abs2.(Z) .< abs2(r)
Z[M], dZ[M] = Z[M] .^ 2 .+ C[M], 2 .* Z[M] .* dZ[M] .+ 1
end
N = abs.(Z) .> 2 # exterior distance estimation
D[N] = log.(abs.(Z[N])) .* abs.(Z[N]) ./ abs.(dZ[N])
heatmap(D' .^ 0.015, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_deep_map.png")
heatmap(D' .^ 0.05, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_map.png")
X, Y = real(C), imag(C) # zoom images (adjust circle size 50 and zoom level 20 as needed)
R, c, z = 50 * (2 / d) * pi .* exp.(.- B), min(d, h) + 1, max(0, h - d) ÷ 20
gr(c=:nipy_spectral, axis=true, ticks=true, legend=false, markerstrokewidth=0)
p1 = scatter(X[1z+1:1z+c,1:d], Y[1z+1:1z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[1z+1:1z+c,1:d].^0.5)
p2 = scatter(X[2z+1:2z+c,1:d], Y[2z+1:2z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[2z+1:2z+c,1:d].^0.4)
p3 = scatter(X[3z+1:3z+c,1:d], Y[3z+1:3z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[3z+1:3z+c,1:d].^0.3)
p4 = scatter(X[4z+1:4z+c,1:d], Y[4z+1:4z+c,1:d], markersize=R[1:c].^0.5, marker_z=D[4z+1:4z+c,1:d].^0.2)
plot(p1, p2, p3, p4, layout=(2, 2))
savefig("Mercator_Mandelbrot_zoom.png")

View file

@ -0,0 +1,46 @@
using Plots
gr(aspect_ratio=:equal, axis=true, ticks=true, legend=false, dpi=200)
setprecision(BigFloat, 256) # set precision to 256 bits (default)
setrounding(BigFloat, RoundNearest) # set rounding mode (default)
d, h = 50, 1000 # pixel density (= image width) and image height
n, r = 80000, 100000 # number of iterations and escape radius (r > 2)
a = BigFloat("-1.256827152259138864846434197797294538253477389787308085590211144291")
b = BigFloat(".37933802890364143684096784819544060002129071484943239316486643285025")
S = zeros(Complex{Float64}, n+1)
let c = a + b * im, z = zero(c)
for k in 1:n+1
S[k] = z
if abs2(z) < abs2(r)
z = z ^ 2 + c
else
println("The reference sequence diverges within $(k-1) iterations.")
break
end
end
end
x = range(0, 2, length=d+1)
y = range(0, 2 * h / d, length=h+1)
A, B = collect(x) .* pi, collect(y) .* pi
C = 8.0 .* exp.((A' .+ B .* im) .* im)
E, Z, dZ = zero(C), zero(C), zero(C)
D, I, J = zeros(size(C)), ones(Int64, size(C)), ones(Int64, size(C))
for k in 1:n
M, R = abs2.(Z) .< abs2(r), abs2.(Z) .< abs2.(E)
E[R], I[R] = Z[R], J[R] # rebase when z is closer to zero
E[M], I[M] = (2 .* S[I[M]] .+ E[M]) .* E[M] .+ C[M], I[M] .+ 1
Z[M], dZ[M] = S[I[M]] .+ E[M], 2 .* Z[M] .* dZ[M] .+ 1
end
N = abs.(Z) .> 2 # exterior distance estimation
D[N] = log.(abs.(Z[N])) .* abs.(Z[N]) ./ abs.(dZ[N])
heatmap(D' .^ 0.015, c=:nipy_spectral)
savefig("Mercator_Mandelbrot_deep_map.png")

View file

@ -0,0 +1,5 @@
from functools import reduce
def mandelbrot(x, y, c): return ' *'[abs(reduce(lambda z, _: z*z + c, range(99), 0)) < 2]
print('\n'.join(''.join(mandelbrot(x, y, x/50 + y/50j) for x in range(-100, 25)) for y in range(-50, 50)))

View file

@ -4,38 +4,31 @@ import matplotlib.pyplot as plt
d, h = 800, 500 # pixel density (= image width) and image height
n, r = 200, 500 # number of iterations and escape radius (r > 2)
direction, height = 45, 1.5 # direction and height of the incoming light
stripes, damping = 4, 2.0 # stripe density and damping parameter
x = np.linspace(0, 2, num=d+1)
y = np.linspace(0, 2 * h / d, num=h+1)
A, B = np.meshgrid(x - 1, y - h / d)
C = (2.0 + 1.0j) * (A + B * 1j) - 0.5
C = 2.0 * (A + B * 1j) - 0.5
Z, dZ, ddZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
Z, dZ = np.zeros_like(C), np.zeros_like(C)
D, S, T = np.zeros(C.shape), np.zeros(C.shape), np.zeros(C.shape)
for k in range(n):
M = abs(Z) < r
S[M], T[M] = S[M] + np.cos(stripes * np.angle(Z[M])), T[M] + 1
Z[M], dZ[M], ddZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1, 2 * (dZ[M] ** 2 + Z[M] * ddZ[M])
S[M], T[M] = S[M] + np.exp(- abs(Z[M])), T[M] + 1
Z[M], dZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1
N = abs(Z) >= r # normal map effect 1 (potential function)
P, Q = S[N] / T[N], (S[N] + np.cos(stripes * np.angle(Z[N]))) / (T[N] + 1)
F = np.log2(np.log(np.abs(Z[N])) / np.log(r)) # fraction between 0 and 1 (for interpolation)
R = Q + (P - Q) * F * F * (3 - 2 * F) # hermite interpolation (r is between q and p)
U, H = Z[N] / dZ[N], 1 + R / damping # normal vectors to the equipotential lines and height perturbation
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and vector in light direction
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + H * height) / (1 + height), 0)
plt.imshow(S ** 0.1, cmap=plt.cm.twilight_shifted, origin="lower")
plt.savefig("Mandelbrot_set_1.png", dpi=200)
plt.imshow(D ** 1.0, cmap=plt.cm.bone, origin="lower")
plt.savefig("Mandelbrot_normal_map_1.png", dpi=200)
N = abs(Z) >= r # normalized iteration count
T[N] = T[N] - np.log2(np.log(np.abs(Z[N])) / np.log(r))
N = abs(Z) >= r # normal map effect 2 (distance estimation)
U = Z[N] * dZ[N] * ((1 + np.log(abs(Z[N]))) * np.conj(dZ[N] ** 2) - np.log(abs(Z[N])) * np.conj(Z[N] * ddZ[N]))
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit normal vectors and vector in light direction
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + height) / (1 + height), 0)
plt.imshow(T ** 0.1, cmap=plt.cm.twilight_shifted, origin="lower")
plt.savefig("Mandelbrot_set_2.png", dpi=200)
plt.imshow(D ** 1.0, cmap=plt.cm.afmhot, origin="lower")
plt.savefig("Mandelbrot_normal_map_2.png", dpi=200)
N = abs(Z) > 2 # exterior distance estimation
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
plt.imshow(D ** 0.1, cmap=plt.cm.twilight_shifted, origin="lower")
plt.savefig("Mandelbrot_set_3.png", dpi=200)

View file

@ -1,37 +1,40 @@
import numpy as np
import matplotlib.pyplot as plt
d, h = 200, 1200 # pixel density (= image width) and image height
n, r = 8000, 10000 # number of iterations and escape radius (r > 2)
d, h = 800, 500 # pixel density (= image width) and image height
n, r = 200, 500 # number of iterations and escape radius (r > 2)
a = -.743643887037158704752191506114774 # real coordinate of the zoom center
b = 0.131825904205311970493132056385139 # imaginary coordinate of the center
direction, height = 45, 1.5 # direction and height of the incoming light
stripes, damping = 4.0, 2.0 # stripe density and damping parameter
x = np.linspace(0, 2, num=d+1)
y = np.linspace(0, 2 * h / d, num=h+1)
A, B = np.meshgrid(x * np.pi, y * np.pi)
C = 8.0 * np.exp((A + B * 1j) * 1j) + (a + b * 1j)
A, B = np.meshgrid(x - 1, y - h / d)
C = (2.0 + 1.0j) * (A + B * 1j) - 0.5
Z, dZ = np.zeros_like(C), np.zeros_like(C)
D = np.zeros(C.shape)
Z, dZ, ddZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
D, S, T = np.zeros(C.shape), np.zeros(C.shape), np.zeros(C.shape)
for k in range(n):
M = Z.real ** 2 + Z.imag ** 2 < r ** 2
Z[M], dZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1
M = abs(Z) < r
S[M], T[M] = S[M] + np.cos(stripes * np.angle(Z[M])), T[M] + 1
Z[M], dZ[M], ddZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1, 2 * (dZ[M] ** 2 + Z[M] * ddZ[M])
N = abs(Z) > 2 # exterior distance estimation
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
N = abs(Z) >= r # normal map effect 1 (equipotential lines)
P, Q = S[N] / T[N], (S[N] + np.cos(stripes * np.angle(Z[N]))) / (T[N] + 1)
R = Q + (P - Q) * np.log2(np.log(np.abs(Z[N])) / np.log(r)) # linear interpolation
U, V = Z[N] / dZ[N], 1 + R / damping # normal vectors and variations in inclination
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit vectors
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + V * height) / (1 + height), 0)
plt.imshow(D.T ** 0.05, cmap=plt.cm.nipy_spectral, origin="lower")
plt.savefig("Mercator_Mandelbrot_map.png", dpi=200)
plt.imshow(D ** 1.0, cmap=plt.cm.bone, origin="lower")
plt.savefig("Mandelbrot_normal_map_1.png", dpi=200)
X, Y = C.real, C.imag # zoom images (adjust circle size 100 and zoom level 20 as needed)
R, c, z = 100 * (2 / d) * np.pi * np.exp(- B), min(d, h) + 1, max(0, h - d) // 20
N = abs(Z) > 2 # normal map effect 2 (equidistant lines)
U = Z[N] * dZ[N] * ((1 + np.log(abs(Z[N]))) * np.conj(dZ[N] ** 2) - np.log(abs(Z[N])) * np.conj(Z[N] * ddZ[N]))
U, v = U / abs(U), np.exp(direction / 180 * np.pi * 1j) # unit vectors
D[N] = np.maximum((U.real * v.real + U.imag * v.imag + height) / (1 + height), 0)
fig, ax = plt.subplots(2, 2, figsize=(12, 12))
ax[0, 0].scatter(X[1*z:1*z+c,0:d], Y[1*z:1*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[1*z:1*z+c,0:d]**0.5, cmap=plt.cm.nipy_spectral)
ax[0, 1].scatter(X[2*z:2*z+c,0:d], Y[2*z:2*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[2*z:2*z+c,0:d]**0.4, cmap=plt.cm.nipy_spectral)
ax[1, 0].scatter(X[3*z:3*z+c,0:d], Y[3*z:3*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[3*z:3*z+c,0:d]**0.3, cmap=plt.cm.nipy_spectral)
ax[1, 1].scatter(X[4*z:4*z+c,0:d], Y[4*z:4*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[4*z:4*z+c,0:d]**0.2, cmap=plt.cm.nipy_spectral)
plt.savefig("Mercator_Mandelbrot_zoom.png", dpi=100)
plt.imshow(D ** 1.0, cmap=plt.cm.afmhot, origin="lower")
plt.savefig("Mandelbrot_normal_map_2.png", dpi=200)

View file

@ -1,44 +1,37 @@
import numpy as np
import matplotlib.pyplot as plt
import decimal as dc # decimal floating point arithmetic with arbitrary precision
dc.getcontext().prec = 80 # set precision to 80 digits (about 256 bits)
d, h = 200, 1200 # pixel density (= image width) and image height
n, r = 8000, 10000 # number of iterations and escape radius (r > 2)
d, h = 50, 1000 # pixel density (= image width) and image height
n, r = 80000, 100000 # number of iterations and escape radius (r > 2)
a = dc.Decimal("-1.256827152259138864846434197797294538253477389787308085590211144291")
b = dc.Decimal(".37933802890364143684096784819544060002129071484943239316486643285025")
S = np.zeros(n+1, dtype=np.complex128)
u, v = dc.Decimal(0), dc.Decimal(0)
for k in range(n+1):
S[k] = float(u) + float(v) * 1j
if u ** 2 + v ** 2 < r ** 2:
u, v = u ** 2 - v ** 2 + a, 2 * u * v + b
else:
print("The reference sequence diverges within %s iterations." % k)
break
a = -.743643887037158704752191506114774 # real coordinate of the zoom center
b = 0.131825904205311970493132056385139 # imaginary coordinate of the center
x = np.linspace(0, 2, num=d+1)
y = np.linspace(0, 2 * h / d, num=h+1)
A, B = np.meshgrid(x * np.pi, y * np.pi)
C = 8.0 * np.exp((A + B * 1j) * 1j)
C = 8.0 * np.exp((A + B * 1j) * 1j) + (a + b * 1j)
E, Z, dZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
D, I, J = np.zeros(C.shape), np.zeros(C.shape, dtype=np.int64), np.zeros(C.shape, dtype=np.int64)
Z, dZ = np.zeros_like(C), np.zeros_like(C)
D = np.zeros(C.shape)
for k in range(n):
Z2 = Z.real ** 2 + Z.imag ** 2
M, R = Z2 < r ** 2, Z2 < E.real ** 2 + E.imag ** 2
E[R], I[R] = Z[R], J[R] # rebase when z is closer to zero
E[M], I[M] = (2 * S[I[M]] + E[M]) * E[M] + C[M], I[M] + 1
Z[M], dZ[M] = S[I[M]] + E[M], 2 * Z[M] * dZ[M] + 1
M = Z.real ** 2 + Z.imag ** 2 < r ** 2
Z[M], dZ[M] = Z[M] ** 2 + C[M], 2 * Z[M] * dZ[M] + 1
N = abs(Z) > 2 # exterior distance estimation
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
plt.imshow(D.T ** 0.015, cmap=plt.cm.nipy_spectral, origin="lower")
plt.savefig("Mercator_Mandelbrot_deep_map.png", dpi=200)
plt.imshow(D.T ** 0.05, cmap=plt.cm.nipy_spectral, origin="lower")
plt.savefig("Mercator_Mandelbrot_map.png", dpi=200)
X, Y = C.real, C.imag # zoom images (adjust circle size 100 and zoom level 20 as needed)
R, c, z = 100 * (2 / d) * np.pi * np.exp(- B), min(d, h) + 1, max(0, h - d) // 20
fig, ax = plt.subplots(2, 2, figsize=(12, 12))
ax[0, 0].scatter(X[1*z:1*z+c,0:d], Y[1*z:1*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[1*z:1*z+c,0:d]**0.5, cmap=plt.cm.nipy_spectral)
ax[0, 1].scatter(X[2*z:2*z+c,0:d], Y[2*z:2*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[2*z:2*z+c,0:d]**0.4, cmap=plt.cm.nipy_spectral)
ax[1, 0].scatter(X[3*z:3*z+c,0:d], Y[3*z:3*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[3*z:3*z+c,0:d]**0.3, cmap=plt.cm.nipy_spectral)
ax[1, 1].scatter(X[4*z:4*z+c,0:d], Y[4*z:4*z+c,0:d], s=R[0:c,0:d]**2.0, c=D[4*z:4*z+c,0:d]**0.2, cmap=plt.cm.nipy_spectral)
plt.savefig("Mercator_Mandelbrot_zoom.png", dpi=100)

View file

@ -1,8 +1,44 @@
print(
'\n'.join(
''.join(
' *'[(z:=0, c:=x/50+y/50j, [z:=z*z+c for _ in range(99)], abs(z))[-1]<2]
for x in range(-100,25)
)
for y in range(-50,50)
))
import numpy as np
import matplotlib.pyplot as plt
import decimal as dc # decimal floating point arithmetic with arbitrary precision
dc.getcontext().prec = 80 # set precision to 80 digits (about 256 bits)
d, h = 50, 1000 # pixel density (= image width) and image height
n, r = 80000, 100000 # number of iterations and escape radius (r > 2)
a = dc.Decimal("-1.256827152259138864846434197797294538253477389787308085590211144291")
b = dc.Decimal(".37933802890364143684096784819544060002129071484943239316486643285025")
S = np.zeros(n+1, dtype=np.complex128)
u, v = dc.Decimal(0), dc.Decimal(0)
for k in range(n+1):
S[k] = float(u) + float(v) * 1j
if u ** 2 + v ** 2 < r ** 2:
u, v = u ** 2 - v ** 2 + a, 2 * u * v + b
else:
print("The reference sequence diverges within %s iterations." % k)
break
x = np.linspace(0, 2, num=d+1)
y = np.linspace(0, 2 * h / d, num=h+1)
A, B = np.meshgrid(x * np.pi, y * np.pi)
C = 8.0 * np.exp((A + B * 1j) * 1j)
E, Z, dZ = np.zeros_like(C), np.zeros_like(C), np.zeros_like(C)
D, I, J = np.zeros(C.shape), np.zeros(C.shape, dtype=np.int64), np.zeros(C.shape, dtype=np.int64)
for k in range(n):
Z2 = Z.real ** 2 + Z.imag ** 2
M, R = Z2 < r ** 2, Z2 < E.real ** 2 + E.imag ** 2
E[R], I[R] = Z[R], J[R] # rebase when z is closer to zero
E[M], I[M] = (2 * S[I[M]] + E[M]) * E[M] + C[M], I[M] + 1
Z[M], dZ[M] = S[I[M]] + E[M], 2 * Z[M] * dZ[M] + 1
N = abs(Z) > 2 # exterior distance estimation
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
plt.imshow(D.T ** 0.015, cmap=plt.cm.nipy_spectral, origin="lower")
plt.savefig("Mercator_Mandelbrot_deep_map.png", dpi=200)

View file

@ -1,5 +1,8 @@
from functools import reduce
def mandelbrot(x, y, c): return ' *'[abs(reduce(lambda z, _: z*z + c, range(99), 0)) < 2]
print('\n'.join(''.join(mandelbrot(x, y, x/50 + y/50j) for x in range(-100, 25)) for y in range(-50, 50)))
print(
'\n'.join(
''.join(
' *'[(z:=0, c:=x/50+y/50j, [z:=z*z+c for _ in range(99)], abs(z))[-1]<2]
for x in range(-100,25)
)
for y in range(-50,50)
))