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 = -1.748764520194788535  # coordinates by Claude Heiland-Allen
b = 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)

def iteration(C):
    Z, dZ = np.zeros_like(C), np.zeros_like(C)

    def iterate2(C, Z, dZ):
        Z, dZ = Z * Z + C, 2 * Z * dZ + 1
        Z, dZ = Z * Z + C, 2 * Z * dZ + 1
        return Z, dZ

    for i in range(0, n, 2):
        M = abs(Z) < r
        Z[M], dZ[M] = iterate2(C[M], Z[M], dZ[M])

    return Z, dZ

Z, dZ = iteration(C)
D = np.zeros(C.shape)

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")

M = 50 * (2 / d) * np.pi * np.exp(- B)  # adjust marker size 50 as needed
k, l = min(d, h) + 1, max(0, h - d) // 8  # adjust zoom level 8 as needed

for i in range(8):
    X, Y = C[i*l:i*l+k, 0:d].real, C[i*l:i*l+k, 0:d].imag
    S, T = M[0:k, 0:d] ** 2, D[i*l:i*l+k, 0:d] ** 0.5
    ax = fig.add_subplot(3, 4, 5 + i)
    ax.scatter(X, Y, s=S, c=T, cmap=plt.cm.nipy_spectral)
    ax.set_xticks([])
    ax.set_yticks([])
    ax.axis('equal')

fig.savefig("Mandelbrot_numpy_set_3.png", dpi=200)