def julia(zmin, zmax, hpx, niter, c, func=lambda z, c: z**2 + c):
vpx=round(hpx * abs((zmax-zmin).imag / (zmax-zmin).real))
x = linspace(zmin.real, zmax.real, hpx)
y = linspace(zmin.imag, zmax.imag, vpx)
zRe, zIm = meshgrid(x, y)
z = zRe + zIm * 1j
M = zeros((vpx,hpx))
for _ in range(niter):
mask = find(abs(z)<2)
M.flat[mask] = M.flat[mask] + 1
z.flat[mask] = func(z.flat[mask], c)
#z.flat[mask] = z.flat[mask]**2 + c
M.flat[mask] = 0
return M
Jc1 = julia(-1.6 + 1.2j, 1.6 - 1.2j, 640, 64, -0.75 + 0.2j)
gray()
Jc1
array([[ 0., 1., 1., ..., 1., 1., 0.], [ 1., 1., 1., ..., 1., 1., 1.], [ 1., 1., 1., ..., 1., 1., 1.], ..., [ 1., 1., 1., ..., 1., 1., 1.], [ 1., 1., 1., ..., 1., 1., 1.], [ 0., 1., 1., ..., 1., 1., 0.]])
imshow(Jc1)
<matplotlib.image.AxesImage at 0xaff511ac>
imshow(Jc1, cmap=cm.hot)
<matplotlib.image.AxesImage at 0xafefc10c>
Jn3 = julia(-1.6 + 1.2j, 1.6 - 1.2j, 320, 64, 0.4, func=lambda z,c: z**3 + c)
imshow(Jn3)
<matplotlib.image.AxesImage at 0xaff2814c>
Jn4=julia(-1.6 + 1.2j, 1.6 - 1.2j, 320, 64, -1, func= lambda z, c: z**4 + c)
imshow(Jn4)
<matplotlib.image.AxesImage at 0xafed108c>