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)

imshow(Msin)

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)

imshow(Jsin)

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)

imshow(Mcos)

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)

imshow(Jcos)

Mtan = mandelbrot(-1.4 + 1.4j, 1.4 - 1.4j, 320, 64, z0='c', func=lambda z, c: c * tan(z))

imshow(Mtan)

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)

Mexp1 = mandelbrot(-1.4 + 2j, 1.4 - 2j, 240, 64, func=lambda z, c: c * exp(z**2))

imshow(Mexp1)

Mexp2=mandelbrot(-2.2 + 2j, 1 - 2j, 240, 64, func=lambda z, c: exp(z**2) + c)

imshow(Mexp2)

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)

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)

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)

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)

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)

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)

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)

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)

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

Jj=julia(-1.2 + 1.6j, 1.2 - 1.6j, 240, 64, c=julia_c, func=julia_polynomial)

imshow(Jj)

Jjz = julia(0.2 + 0.1j, 0.4 - 0.05j, 320, 64, c=julia_c, func=julia_polynomial)

imshow(Jjz)

Jjzz = julia(0.366 + 0.09j, 0.386 + 0.075j, 320, 64, c=julia_c, func=julia_polynomial)

imshow(Jjzz)

Jjzzz = julia(0.366 + 0.09j, 0.386 + 0.075j, 320, 64, c=julia_c, func=julia_polynomial)

imshow(Jjzzz)