In [16]:
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
In [17]:
Jc1 = julia(-1.6 + 1.2j, 1.6 - 1.2j, 640, 64, -0.75 + 0.2j)
In [18]:
gray()
In [19]:
Jc1
Out[19]:
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.]])
In [20]:
imshow(Jc1)
Out[20]:
<matplotlib.image.AxesImage at 0xaff511ac>
In [21]:
imshow(Jc1, cmap=cm.hot)
Out[21]:
<matplotlib.image.AxesImage at 0xafefc10c>
In [22]:
Jn3 = julia(-1.6 + 1.2j, 1.6 - 1.2j, 320, 64, 0.4, func=lambda z,c: z**3 + c)
In [23]:
imshow(Jn3)
Out[23]:
<matplotlib.image.AxesImage at 0xaff2814c>
In [24]:
Jn4=julia(-1.6 + 1.2j, 1.6 - 1.2j, 320, 64, -1, func= lambda z, c: z**4 + c)
In [25]:
imshow(Jn4)
Out[25]:
<matplotlib.image.AxesImage at 0xafed108c>
In [ ]: