In [63]:
def mjcore(z, c, niter, func, bound=2, bound_func=None):
    if bound_func is not None:
        print('Warning: bound_func not None. Ignoring r bound')

    M = zeros(z.shape)
    
    Mf = M.flat
    zf = z.flat
    cf = c.flat
    
    for _ in range(niter):
        if bound_func is not None:
            mask = find(bound_func(z))
        else:
            mask = find(abs(z) < bound)
        
        Mf[mask] += 1
        zf[mask]  = func(zf[mask], cf[mask])
    
    Mf[mask] = 0
    
    return M
In [77]:
def mandelbrot(cmin, cmax, hpx, niter, z0=0, func=lambda z, c: z**2 + c, bound=2, bound_func=None):
    vpx = round(hpx * abs((cmax.imag - cmin.imag) / (cmax.real - cmin.real)))
    
    cRe, cIm = meshgrid(linspace(cmin.real, cmax.real, hpx),
                        linspace(cmin.imag, cmax.imag, vpx))
    c = cRe + cIm * 1j
  
    if z0 == 'c':
        z = c
    else:
        z = zeros((vpx, hpx), dtype=complex128) + z0
        
    return mjcore(z, c, niter, func, bound, bound_func)
In [65]:
def julia(zmin, zmax, hpx, niter, c, func=lambda z, c: z**2 + c, bound=2, bound_func=None):
    vpx=round(hpx * abs((zmax-zmin).imag / (zmax-zmin).real))
    
    zRe, zIm = meshgrid(linspace(zmin.real, zmax.real, hpx), 
                        linspace(zmin.imag, zmax.imag, vpx))
    
    z = zRe + zIm * 1j
    cc = empty(z.shape, dtype=complex128)
    cc.flat[:] = c
    
    return mjcore(z, cc, niter, func, bound, bound_func)
In [31]:
Msin = mandelbrot(-1.2 * pi + 0.9 * pi * 1j, 1.2 * pi - 0.9 * pi * 1j, 320, 64, 
    z0='c', func=lambda z,c: c * sin(z), bound_func=lambda z: abs(z.imag)<=50)
Warning: bound_func not None. Ignoring r bound
In [32]:
imshow(Msin)
Out[32]:
<matplotlib.image.AxesImage at 0xb016806c>
In [59]:
Jsin = julia(-2.4 *pi + 1.8 * pi * 1j, 2.4 * pi - 1.8 * pi * 1j, 320, 64, 1 + 0.5j,
    func=lambda z,c: c * sin(z), bound_func=lambda z: abs(z.imag) <= 50)
Warning: bound_func not None. Ignoring r bound
In [60]:
imshow(Jsin)
Out[60]:
<matplotlib.image.AxesImage at 0xb0323ecc>
In [68]:
Mcos = mandelbrot(-1.2 * pi + 0.9 * pi * 1j, 1.2 * pi - 0.9 * pi * 1j, 320, 64, z0='c', 
    func=lambda z, c: 1j * c * cos(z), bound_func=lambda z: abs(z.imag)<=50)
Warning: bound_func not None. Ignoring r bound
In [69]:
imshow(Mcos)
Out[69]:
<matplotlib.image.AxesImage at 0xb02a308c>
In [70]:
Jcos = julia(-2.4 * pi + 1.8 * pi * 1j, 2.4 * pi - 1.8 * pi * 1j, 320, 64, 0.55 + 1.195j,
    func=lambda z, c: 1j * c * cos(z), bound_func=lambda z: abs(z.imag) <= 50)
Warning: bound_func not None. Ignoring r bound
In [71]:
imshow(Jcos)
Out[71]:
<matplotlib.image.AxesImage at 0xb024dfac>
In [72]:
Mtan = mandelbrot(-1.4 + 1.4j, 1.4 - 1.4j, 320, 64, z0='c', func=lambda z, c: c * tan(z))
In [73]:
imshow(Mtan)
Out[73]:
<matplotlib.image.AxesImage at 0xb01f6dac>
In [74]:
Jtan = julia(-1.4 + 1.4j, 1.4 - 1.4j, 320, 64, (1 + 1j) * sqrt(2) / 2, func=lambda z, c: c * tan(z))
In [75]:
imshow(Jtan)
Out[75]:
<matplotlib.image.AxesImage at 0xb022534c>
In [78]:
Mexp1 = mandelbrot(-1.4 + 2j, 1.4 - 2j, 240, 64, func=lambda z, c: c * exp(z**2))
In [79]:
imshow(Mexp1)
Out[79]:
<matplotlib.image.AxesImage at 0xb01d628c>
In [80]:
Mexp2=mandelbrot(-2.2 + 2j, 1 - 2j, 240, 64, func=lambda z, c: exp(z**2) + c)
In [81]:
imshow(Mexp2)
Out[81]:
<matplotlib.image.AxesImage at 0xb04c520c>
In [82]:
Mcactus = mandelbrot(-0.8 + 0.6j, 0.8 - 0.6j, 320, 64, z0='c', func=lambda z, c: z**3 + (c - 1) * z - c)
In [83]:
imshow(Mcactus)
Out[83]:
<matplotlib.image.AxesImage at 0xb04ed6ac>

The next fractal is called the burning ship.

In [84]:
Mbs = mandelbrot(-2.5 - 2j, 1.5 + 1j, 320, 64, func=lambda z, c: (abs(z.real) + 1j * abs(z.imag))**2 + c, bound=200)
In [85]:
imshow(Mbs)
Out[85]:
<matplotlib.image.AxesImage at 0xb04985ec>
In [86]:
Mbsz=mandelbrot(-1.8 - 0.06j, -1.7 + 0.02j, 320, 64, func=lambda z, c: (abs(z.real)+ 1j * abs(z.imag))**2 + c, bound=200)
In [87]:
imshow(Mbsz)
Out[87]:
<matplotlib.image.AxesImage at 0xb04415ec>

And this next fractal is called the bird of prey.

In [88]:
Mbs3 = mandelbrot(-1.5 - 1.5j, 1.5 + 1.5j, 320, 64, func=lambda z, c: (abs(z.real) + 1j * abs(z.imag))**3 + c, bound=200)
In [89]:
imshow(Mbs3)
Out[89]:
<matplotlib.image.AxesImage at 0xb046d04c>

The original mandelbrot and Julia sets

The fractal studied by Mandelbrot

In [91]:
Mj = mandelbrot(-1.4 + 1.05j, 1.4 - 1.05j, 320, 64, func=lambda z, c: z**4 + z**3 / (z - 1) + z**2 / (z**3 + 4 * z**2 + 5) + c)
In [92]:
imshow(Mj)
Out[92]:
<matplotlib.image.AxesImage at 0xb041a90c>

this is a zoom on the bottom right corner of that image (haven't checked that yet!)

In [93]:
Mjz = mandelbrot(0.25 - 0.65j, 0.45 - 0.8j, 320, 64, func=lambda z, c: z**4 + z**3 / (z - 1) + z**2 / (z**3 + 4 * z**2 + 5) + c)
In [94]:
imshow(Mjz)
Out[94]:
<matplotlib.image.AxesImage at 0xb017e86c>

another zoom:

In [95]:
Mjzz = mandelbrot(0.378 - 0.725j, 0.398 - 0.74j, 320, 64, func=lambda z, c: z**4 + z**3 / (z - 1) + z**2 / (z**3 + 4 * z**2 + 5) + c)
In [96]:
imshow(Mjzz)
Out[96]:
<matplotlib.image.AxesImage at 0xb01a878c>

One last zoom:

In [97]:
Mjzzz = mandelbrot(0.3863 - 0.7314j, 0.3883 - 0.7329j, 320, 64, func=lambda z,c: z**4 + z**3 / (z - 1) + z**2 / (z**3 + 4 * z**2 + 5) + c)
In [98]:
imshow(Mjzzz)
Out[98]:
<matplotlib.image.AxesImage at 0xb00796cc>