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
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)
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)
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
imshow(Msin)
<matplotlib.image.AxesImage at 0xb016806c>
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
imshow(Jsin)
<matplotlib.image.AxesImage at 0xb0323ecc>
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
imshow(Mcos)
<matplotlib.image.AxesImage at 0xb02a308c>
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
imshow(Jcos)
<matplotlib.image.AxesImage at 0xb024dfac>
Mtan = mandelbrot(-1.4 + 1.4j, 1.4 - 1.4j, 320, 64, z0='c', func=lambda z, c: c * tan(z))
imshow(Mtan)
<matplotlib.image.AxesImage at 0xb01f6dac>
Jtan = julia(-1.4 + 1.4j, 1.4 - 1.4j, 320, 64, (1 + 1j) * sqrt(2) / 2, func=lambda z, c: c * tan(z))
imshow(Jtan)
<matplotlib.image.AxesImage at 0xb022534c>
Mexp1 = mandelbrot(-1.4 + 2j, 1.4 - 2j, 240, 64, func=lambda z, c: c * exp(z**2))
imshow(Mexp1)
<matplotlib.image.AxesImage at 0xb01d628c>
Mexp2=mandelbrot(-2.2 + 2j, 1 - 2j, 240, 64, func=lambda z, c: exp(z**2) + c)
imshow(Mexp2)
<matplotlib.image.AxesImage at 0xb04c520c>
Mcactus = mandelbrot(-0.8 + 0.6j, 0.8 - 0.6j, 320, 64, z0='c', func=lambda z, c: z**3 + (c - 1) * z - c)
imshow(Mcactus)
<matplotlib.image.AxesImage at 0xb04ed6ac>
The next fractal is called the burning ship.
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)
imshow(Mbs)
<matplotlib.image.AxesImage at 0xb04985ec>
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)
imshow(Mbsz)
<matplotlib.image.AxesImage at 0xb04415ec>
And this next fractal is called the bird of prey.
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)
imshow(Mbs3)
<matplotlib.image.AxesImage at 0xb046d04c>
All the previous fractals where taken from http://brunogirin.blogspot.mx/2009/08/fractals-with-octave-burning-ship.html.
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)
imshow(Mj)
<matplotlib.image.AxesImage at 0xb041a90c>
this is a zoom on the bottom right corner of that image (haven't checked that yet!)
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)
imshow(Mjz)
<matplotlib.image.AxesImage at 0xb017e86c>
another zoom:
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)
imshow(Mjzz)
<matplotlib.image.AxesImage at 0xb01a878c>
One last zoom:
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)
imshow(Mjzzz)
<matplotlib.image.AxesImage at 0xb00796cc>
julia_c = 0.3873 - 0.7314j
julia_polynomial = lambda z,c: z**4 + z**3 / (z - 1) + z**2 / (z**3 + 4 * z**2 + 5) + c
julia_c is the point around which we were zooming in the Mandelbrot set (should I check this?)
Jj=julia(-1.2 + 1.6j, 1.2 - 1.6j, 240, 64, c=julia_c, func=julia_polynomial)
imshow(Jj)
<matplotlib.image.AxesImage at 0xafda76cc>
Jjz = julia(0.2 + 0.1j, 0.4 - 0.05j, 320, 64, c=julia_c, func=julia_polynomial)
imshow(Jjz)
<matplotlib.image.AxesImage at 0xafdd260c>
Jjzz = julia(0.366 + 0.09j, 0.386 + 0.075j, 320, 64, c=julia_c, func=julia_polynomial)
imshow(Jjzz)
<matplotlib.image.AxesImage at 0xafd7d54c>
Jjzzz = julia(0.366 + 0.09j, 0.386 + 0.075j, 320, 64, c=julia_c, func=julia_polynomial)
imshow(Jjzzz)
<matplotlib.image.AxesImage at 0xafc8356c>