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)