#Tasia Raymer
# 250 Project

#This program plots the time series x1 vs x2 for a fixed value of r1 and r2 for 6 different d values.
#It also calculates numerical estimates of the LCEs



#dynamic equations:
from CpledLogClass import *

import pylab as pl
from RandomArray import *
from math import sqrt
from math import log
from Scientific.Geometry import Vector

#set parameter values:
r1=3.8
r2=3.8
dvec=[0.0,0.05,0.09,0.1,0.11,0.15,0.25,0.4]


	
# intial conditions are choose randomly from uniform [0,1] distribution:
seed()
uniform(0,1)
x1_0=uniform(0,1)
x2_0=uniform(0,1)



print 'r1=%g r2=%g X_0=(%g,%g)'%(r1,r2,x1_0,x2_0)
nTrans=300
nIterations=1000

# Initial tangent vectors for LCE computation:
e1x1=1.0
e1x2=0.0
e2x1=0.0
e2x2=1.0

f=pl.figure(1,(15,10))
count=0
for i in range(len(dvec)):
	count +=1
	x1=x1_0
	x2=x2_0
	d=dvec[i]
	
	#set dynamics equations from the class with the above parameters:
	F=CpledLogEqs(r1,r2,d)
	
	for n in range(nTrans):
	#evolve flow:
		x1temp=F.x1dot(x1,x2)
		x2=F.x2dot(x1,x2)
		x1=x1temp
	#evolve tangent vector e1:
		tempe1x1=F.dx1dot(x1,x2,e1x1,e1x2)
		e1x2=F.dx2dot(x1,x2,e1x1,e1x2)
		e1x1=tempe1x1
	# normalize:
		mag=sqrt(e1x1**2+e1x2**2)
		e1x1=e1x1/mag
		e1x2=e1x2/mag
	# evolve tangent vector e2:
		tempe2x1=F.dx1dot(x1,x2,e2x1,e2x2)
		e2x2=F.dx2dot(x1,x2,e2x1,e2x2)
		e2x1=tempe2x1
	# remove e1 component from e2:
		e1dote2=e1x1*e2x1+e1x2*e2x2
		e2x1 -=e1dote2*e1x1
		e2x2 -=e1dote2*e1x2
	# normalize:
		mag=sqrt(e2x1**2+e2x2**2)
		e2x1=e2x1/mag
		e2x2=e2x2/mag
		
	# We will use the next nIterations to estimate the LCEs
	#initialize:
	minLCE=0.0
	maxLCE=0.0
	Cont=0.0
	# We will also use the iterations to plot the time series:
	x1vals=[]
	x2vals=[]
	
	for n in range(nIterations):
	#evolve flow:
		x1temp=F.x1dot(x1,x2)
		x2=F.x2dot(x1,x2)
		x1=x1temp
		#store values:
		x1vals.append(x1)
		x2vals.append(x2)
	# calculate the contraction:
		Cont +=log(abs(F.DetJac(x1,x2)))
	
	#evolve tangent vector e1:
		tempe1x1=F.dx1dot(x1,x2,e1x1,e1x2)
		e1x2=F.dx2dot(x1,x2,e1x1,e1x2)
		e1x1=tempe1x1
	# normalize:
		mag=sqrt(e1x1**2+e1x2**2)
		e1x1=e1x1/mag
		e1x2=e1x2/mag
	# accumulate stretch of e1:
		maxLCE +=log(mag)
		
	# evolve tangent vector e2:
		tempe2x1=F.dx1dot(x1,x2,e2x1,e2x2)
		e2x2=F.dx2dot(x1,x2,e2x1,e2x2)
		e2x1=tempe2x1
	# remove e1 component from e2:
		e1dote2=e1x1*e2x1+e1x2*e2x2
		e2x1 -=e1dote2*e1x1
		e2x2 -=e1dote2*e1x2
	# normalize:
		mag=sqrt(e2x1**2+e2x2**2)
		e2x1=e2x1/mag
		e2x2=e2x2/mag
	# accumulate shrink of e2:
		minLCE +=log(mag)
		
	# per-iterate:
	maxLCE=maxLCE/float(nIterations)
	minLCE=minLCE/float(nIterations)
	Cont=Cont/float(nIterations)

	print 'd=%g:  LCE= (%g,%g)   Cont=%g  diff=%g\n'%(d,minLCE,maxLCE,Cont,abs(minLCE+maxLCE-Cont))
	f.text(.5,.95,'r1=%g r2=%g X_0=(%g,%g)'%(r1,r2,x1_0,x2_0))
	pl.subplot(2,4,count)
	pl.title('d=%g'%d)
	if count>=4:
		pl.xlabel('x1(n)')
	pl.ylabel('x2(n)')
	pl.plot(x1vals,x2vals,'b.')
pl.show()