# Code from last time:
import numpy
# "Forward iteration method", or Boundary Scanning Method
def julia_color(c, max_iter, xmin=-2.0, xmax=2.0, ymin=-2.0,ymax=2.0, width=400, height=400):
    """Draw Julia set J(f) for f(z)=z^2+c, c a fixed complex number with |c|<2,
    in window where x ranges in (xmin,xmax), y ranges in (ymin,ymax), 
    grid of size width x height (in pixels)
    Test max_iter number iterates of f.
    Color based on how long it takes to escape"""

    xinc = (xmax-xmin)/width
    yinc = (ymax-ymin)/height

    escape = numpy.zeros((height,width))
    # array will eventually represent escaping

    for x in range(width):
        for y in range(height): 
            z = complex(xmin + x*xinc, ymin + y*yinc)
            # Want: set escape[y,x]=1 if z escapes i.e. |z|, |f(z)|, |f^{\circ 2}(z)| ->infty
            # Special thing about f(z) = z^2+c, when |c|<2:
            # if |z|>2, |f^{\circ n}(z)|-> infty as n ->infty
            # So: if for any z, |f^{\circ k}(z)|>2, z escapes 
            
            # compute iterates of f applied to z 
            n=0 # which iterate we're currently on
            while n < max_iter: 
                z = z^2 + c # take next iterate i.e. replace z by f(z)
                n+=1
                if abs(z)>2:
                    # then by above comments, know that original z escapes 
                    break # breaks out of containing loop 
            # at this point n is between 0 and max_iter,
            # bigger the longer it took escape
            escape[y,x]=n/max_iter
            # if didn't escape, n=max_iter
    
    return matrix_plot(escape, origin='lower', cmap='hot')

# Julia set for z^2-1
julia_color(complex(-1,0), 20)

# Mandlebrot set
import numpy
def mandlebrot(max_iter, xmin=-2.0, xmax=2.0, ymin=-2.0,ymax=2.0, width=400, height=400):
    """Draw Mandlebrot set in window specified by parameters. Represent
    complex plane by grid. Point c is in Mandlebrot if 0 is not in escaping set
    for f_c(z) = z^2 + c. Test max_iter number of iterations."""
    
    xinc = (xmax-xmin)/width
    yinc = (ymax-ymin)/height
    
    mand = numpy.zeros((height,width)) # array representing mandlebrot set
    
    for x in range(width):
        for y in range(height):
            c = complex(xmin + x*xinc, ymin + y*yinc)
            
            # determine whether 0 lies in escaping set for z^2 + c
            z = 0 # starting value; will update
            n = 0 # number of iterations done so far
            while n<max_iter:
                z = z*z + c # apply f to z, update
                n += 1
                if abs(z)>2:
                    break
            # at this point n = max_iter if 0 hasn't escaped yet
            # if 0 does escape, n is smaller the faster it has escaped 
            mand[y,x] = n/max_iter
    return matrix_plot(mand, origin='lower', cmap='hot')
                

mandlebrot(70)

# Feigenbaum point 
c=complex(-1.401155,0)
cr = real(c)
w=0.5 # width of window is 2w
mandlebrot(70, cr-w, cr+w, 0-w, 0+w)

# Feigenbaum point 
c=complex(-1.401155,0)
cr = real(c)
w=0.02 # width of window is 2w
mandlebrot(130, cr-w, cr+w, 0-w, 0+w)
# A sequence of shrinking "bulbs" starting from the main cardiod on the right
# and moving to the left converges to the Feigenbaum point. 

# Feigenbaum Julia set
julia_color(c,40)

See http://www.malinc.se/m/JuliaSets.php for a streamlined Julia/Mandlebrot interactive graphic

# code from last time for drawing orbits:
def draw_orbit(z, c, num_iter, col='blue'):
    """Return a plot of points z,f(z),f(f(z)),.. up to num_iter iterations,
    where f(z)=z^2+c. Can set color"""
    orbit = [z] # list that will store the orbit 
    for n in range(num_iter):
        new = orbit[-1]^2 + c 
        orbit.append(new)
    return list_plot(orbit, size=20, color=col)\
    + line([(real(z),imag(z)) for z in orbit], color=col)

draw_orbit(complex(-0.6,0.6), complex(0,0), 2)