JuliaSet_IP.py

import numpy as np 
import matplotlib.pyplot as plt 
import matplotlib.animation as animation
from numba import jit

# Julia set
# https://en.wikipedia.org/wiki/Julia_set#Quadratic_polynomials

def motion(event):  
    if (event.xdata is  None) or (event.ydata is  None):
        return
    cx = event.xdata
    cy = event.ydata
    cz = cx + cy*1j
    A  = JuliaSet(N,Nx,Ny,z,cz)
  
    im.set_data(A)
    print('c = %+.3f %+.3fj' % (cz.real,cz.imag))
    plt.draw()

@jit
def JuliaSet(N,Nx,Ny,z,c):
    M = Nx*Ny
    Z = z*np.ones(M,dtype=complex)
    A = np.zeros(M)
    for i in range(M):
        for n in range(N):
            if np.abs(Z[i]) < 2:
                Z[i] = Z[i]**2 + c
            else:
                A[i] = n
                break    
        
    A = A.reshape(Ny,Nx)
    return A

N = 200

h = 0.005

xlm = [-1.6,1.6]
ylm = [-1.2,1.2]

c = -0.747 + 0.210j

x  = np.arange(xlm[0],xlm[1]+h,h)
y  = np.arange(ylm[0],ylm[1]+h,h)
Nx = len(x)
Ny = len(y)

xx, yy = np.meshgrid(x,y)
z  = (xx + yy*1j).flatten()

A = JuliaSet(N,Nx,Ny,z,c)


plt.close('all')
fig = plt.figure()

cmap = 'hot'
interpolation = 'None'
im = plt.imshow(A,extent=[xlm[0],xlm[1],ylm[0],ylm[1]],cmap=cmap,interpolation=interpolation)

plt.clim([0,100])
plt.axis('off')

plt.connect('motion_notify_event', motion)
plt.show()