from pylab import *
from scipy.integrate import odeint
from scipy.interpolate import interp1d
from numpy import *
def rhs(xvector,t):
global prevx
global prevy
global prevvx
global prevvy
global prevt
global poincaresection
#****************
# Paste the contents of your rhs function from 3c here
#****************
x1dot=?????
x2dot=?????
x3dot=?????
x4dot=?????
#Here I'm going to rename the xvector variables in terms of their
# original meanings to make the remainder of this function easier
# to read.
#****************
x=xvector[0]
y=xvector[2]
vx=xvector[1]
vy=xvector[3]
# This first checks to make sure vy>0, we only want to plot
# cases when vy>0
if(vy>0.0):
# If y>0 now and was <0 on the previous timestep then we know
# it crossed the y=0 surface between then and now, we don't know
# exactly when or where yet though
if y>0.0 and (prevy<0.0):
interpx=[prevx, x]
interpy=[prevy, y]
interpvx=[prevvx, vx]
interpvy=[prevvy, vy]
interpt=[prevt,t]
# interpolate all our variables between the previous timestep
# and this one. This will let us get the value of any variable
# at any time between the last timestep and now
fy=interp1d(interpy,interpt,kind='linear')
fx=interp1d(interpt,interpx,kind='linear')
fvx=interp1d(interpt,interpvx,kind='linear')
fvy=interp1d(interpt,interpvy,kind='linear')
# Get the time when y was equal to zero exactly
crossingtime=fy(0.0)
# Find the values of all the other variables at this time
poincarepointvx=fvx(crossingtime)
poincarepointx=fx(crossingtime)
poincarepointvy=fvy(crossingtime)
# Add the coordinates of the point we want (in this case [x,vx]) to the poincaresection
# list
poincaresection=vstack((poincaresection,array([poincarepointx,poincarepointvx])))
# Make our current timestep the previous timestep
prevx=x
prevy=y
prevvy=vy
prevvx=vx
prevt=t
return [x1dot, x2dot, x3dot, x4dot]
#Define how long we will run for and how many steps we want to take
start=.0
end=17000.0
num=8000
t=linspace(start,end,num)
sigma=0.25
# Define the list we will store our poincare section in
poincaresection=array([0,0])
# Our initial conditions. You will need to fill these in.
E=?????
x0=?????
vx0=?????
y0=?????
vy0=?????
# Initially, we will say our previous timestep is just our current timestep
# Don't worry too much about this
prevx=x0
prevy=y0
prevvx=vx0
prevvy=vy0
# Call odeint
solution=odeint(rhs,[x0,y0,vx0,vy0],t)
# Plot our poincare section
figure()
plot(poincaresection[:,0],poincaresection[:,1],'o',markersize=0.6)
xlabel('x')
ylabel('vx')
title('E='+str(E)+', sigma='+str(sigma))
#savefig('poincare.pdf')
show()