import numpy as np
import matplotlib.pyplot as plt
d, h = 200, 1000 # pixel density (= image width) and image height
n, r = 800, 5000 # number of iterations and escape radius (r > 2)
a, b = -1.748764520194788535, 3e-13 # https://mathr.co.uk/web/m-location-analysis.html
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 = 4.0 * np.exp((A + B * 1j) * 1j) + (a + b * 1j)
Z, dZ = np.zeros_like(C), np.zeros_like(C)
D = 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
fig = plt.figure(figsize=(12.8, 9.6))
fig.subplots_adjust(left=0.05, right=0.95, bottom=0.05, top=0.95)
N = abs(Z) > 2 # exterior distance estimation
D[N] = np.log(abs(Z[N])) * abs(Z[N]) / abs(dZ[N])
ax1 = fig.add_subplot(3, 1, 1)
ax1.imshow(D.T ** 0.05, cmap=plt.cm.nipy_spectral, origin="lower")
X, Y = C.real, C.imag # zoom images (adjust circle size 50 and zoom level 8 as needed)
R, c, z = 50 * (2 / d) * np.pi * np.exp(- B), min(d, h) + 1, max(0, h - d) // 8
for i in range(8):
axs = fig.add_subplot(3, 4, 5 + i)
axs.scatter(X[i*z:i*z+c,0:d], Y[i*z:i*z+c,0:d], s=R[0:c,0:d]**2, c=D[i*z:i*z+c,0:d]**.5, cmap=plt.cm.nipy_spectral)
axs.set_xticks([])
axs.set_yticks([])
axs.axis('equal')
fig.savefig("Mandelbrot_numpy_set_3.png", dpi=200)