# BakersMap.py: Plot out iterates of the Dissipative Baker's Map
#
# The Dissipative Baker's Map is given by
#     x_n+1 = f( x_n , y_n )
#     y_n+1 = g( x_n , y_n )
# with
#     f( x_n , y_n ) = 2 * x_n (mod 1)
#     g( x_n , y_n ) = a * y_n,       if x_n <= 1/2
#                    = 1/2 + a * y_n, if x_n > 1/2
#
# The state space is the unit square:
#       (x,y) in [0,1] x [0,1]
# and the control parameter ranges in
#       a in [0,1/2]
# 

def StandardMap(k,x,y):
	# We'll assume, for now, that (x,y) is in [0,1] x [0,1]
	# Calculate next y first, since it depends on current x
	# a '-' in the y equation shifts and is used
	# for the manifold plot for better visual
	y=(y-k/(2*pi)*sin(2*pi*x))
	x=(x+y)
	return  x%(1),y#%(1)



# Import plotting routines
from pylab import *
#from matplotlib import axes3d


def StandardMapIterate(k=.15, N=5000, n_init=20, box = [0.0,1.0,0.0,1.0]):
	# finds random initial points in the box
	xlen = box[1]-box[0]
	ylen = box[3]-box[2]
	inits=[]
	for i in xrange(n_init):
		xtemp = (random()*xlen+box[0])
		ytemp = (random()*ylen+box[2])
		inits.append([xtemp,ytemp])

	for nums in inits:
		# Set up arrays of iterates (x_n,y_n) and set the initital condition
		x=[nums[0]]
		y=[nums[1]]

		# The main loop that generates iterates and stores them
		for n in xrange(0,N):
		  # at each iteration calculate (x_n+1,y_n+1)
		  xtemp, ytemp = StandardMap(k,x[n],y[n])
		  # and append to lists x and y
		  x.append( xtemp )
		  y.append( ytemp )
		plot(x,y,',')
		  
	# Setup the plot
	xlabel('theta(n)') # set horizontal axis label
	ylabel('p(n)') # set vertical axis label
	title("Standard Map with k =" + str(k)) # set plot title
	# Plot the time series
		
	# Make sure the iterates appear in a unit square
	#axis('equal')
	axis(box)

	# Display the plot in a window
	show()
	
StandardMapIterate(k=1.0, N = 20000, n_init = 80, box =[.63,.64,.44,.455])
	
def StdMapLineIterate(k=0.8, nICs=5000, niterates=20, endpts = ((0.3,0.2),(0.3,0.8))):
	# create line segment
	xmin = endpts[0][0]
	ymin = endpts[0][1]
	xmax = endpts[1][0]
	ymax = endpts[1][1]
	
	xinc = (xmax - xmin) / nICs 
	yinc = (ymax - ymin) / nICs
	# An array of initial conditions on a straight line
	x = []
	y = []
	
	for i in xrange(0,nICs):
		x.append(xmin + xinc * i)
		y.append(ymin + yinc * i)
	# plot first line segment
	plot(x,y,',',label='0')
	
	# iterate the line segment
	for n in xrange(niterates):
	
		Tx=[]
		Ty=[]
		# now for the mapping
		for j in xrange(len(x)):
			xtemp,ytemp=StandardMap(k,x[j],y[j])
			Tx.append(xtemp)
			Ty.append(ytemp)
		# plot iterated line segment (different color)	
		x = Tx
		y = Ty
		
	# indent to plot each iteration 
	plot(x,y, ',', label=str(n+1))
	
	# Setup the plot
	xlabel('theta(n)') # set x-axis label
	ylabel('p(n)') # set y-axis label
	title("Standard Map with k =" + str(k)) # set plot title
	legend()
		
	# Make sure the iterates appear in a unit square
	#axis('equal')
	#axis([0.0, 1.0, 0.0, 1.0])

	# Use command below to save figure
	#savefig('BakersMapIterates', dpi=600)

	# Display the plot in a window
	show()
	
#StdMapLineIterate(k=.8,niterates=15,endpts = ((0.504,1.0023),(0.5,1.0)))
def StandardMapIteration(k=0, N=5000, n_init=((0,1./4.),(0,1./2.),(0,3./4.),(0,1./5.),(0,3./20.))):
	# allows user to designate initial conditions for the standard map
	
	for init in n_init:
		# Set up arrays of iterates (x_n,y_n) and set the initital condition
		x=[init[0]]
		y=[init[1]]
		
		# throw away transients
		xtemp, ytemp = StandardMap(k,x[0],y[0])
		
		# uncomment to throw away transients
		'''for n in xrange(10):
		  # at each iteration calculate (x_n+1,y_n+1)
		  xtemp, ytemp = StandardMap(k,xtemp,ytemp)		
		
		x[0],y[0] = StandardMap(k,xtemp,ytemp)
		'''
		# The main loop that generates iterates and stores them
		for n in xrange(0,N):
		  # at each iteration calculate (x_n+1,y_n+1)
		  xtemp, ytemp = StandardMap(k,x[n],y[n])
		  # and append to lists x and y
		  x.append( xtemp )
		  y.append( ytemp )
		  #plot each trjectory, top one includes label 
		  #plot(x,y,',',label = 'P_0: ' + str(y[0]))
		  plot(x,y,',')
		  
	# Setup the plot
	xlabel('theta(n)') # set horizontal axis label
	ylabel('p(n)') # set vertical axis label
	title("Standard Map with k =" + str(k)) # set plot title
	# uncomment with line 136
	#legend()

	# Then run the standard map iterate to fill in some quasiperiodic orbits
	# this function will also plot
	#StandardMapIterate()
	show()

#StandardMapIteration()
	
def UnstableManifold(k=1.0, fxpt = (0.5,0)):
	# computes and plots the unstable manifold
	# currently only works for the (0.5,0) fixed point because the 
	# eigvalues are easier to compute, but it can be easily modified
	# for larger N the manifolds look more continuous
	delta = 1e-4
	# initial fixed point
	x = fxpt[0]
	y = fxpt[1]
	# slope is along the unstable eigenvector  
	#eigvec = (1,-.5*k*cos(2*pi*x0)+.5*(k**2*cos(2*pi*x0)**2+4*k*cos(2*pi*x0))**.5)
	eigval = .5*(2+k-(k*(k+4))**.5)
	eigvec = (1, 1-eigval)
	# small displacement from fixed point along stable manifold
	p1 = x - delta*eigvec[0]
	q1 = y - delta*eigvec[1]
	# include both directions of manifold
	p2 = x + delta*eigvec[0]
	q2 = y + delta*eigvec[1]
	# iterate line segment
	StdMapLineIterate(k=k,niterates=15,endpts = ((p1,q1),(p2,q2)))
	
def StableManifold(k=.8, fxpt = (0.5,0)):
	# computes and plots the unstable manifold
	# currently only works for the (0.5,0) fixed point because the 
	# eigvalues are easier to compute, but it can be easily modified
	# for larger N the manifolds look more continuous
	delta = 1e-4
	# initial fixed point
	x = fxpt[0]
	y = fxpt[1]
	# slope is along the unstable eigenvector  
	#eigvec = (1,-.5*k*cos(2*pi*x0)+.5*(k**2*cos(2*pi*x0)**2+4*k*cos(2*pi*x0))**.5)
	eigval = .5*(2+k+(k*(k+4))**.5)
	eigvec = (1, 1-eigval)
	# small displacement from fixed point along stable manifold
	p1 = x - delta*eigvec[0]
	q1 = y - delta*eigvec[1]
	# include both directions of manifold
	p2 = x + delta*eigvec[0]
	q2 = y + delta*eigvec[1]
	# iterate line segment
	StdMapLineIterate(k=.8,niterates=15,endpts = ((p1,q1),(p2,q2)))