# Benchmark 1 Code starts here
# 3 body problem near L2

M=1.0 # set the mass of the star
m=1.0/100.0 # set the mass of the planet

def f(Y,t):
# x1,y1 are the coordinates of the planet      
       x1 = Y[0]
       y1 = Y[1]
       vx1 = Y[4]
       vy1 = Y[5]
# x2, y2 are the coordinates of the sattellite      
       x2 = Y[2]
       y2 = Y[3]
       vx2 = Y[6]
       vy2 = Y[7]
       r13 = (x1**2+y1**2)**1.5 #star-planet distance
       r23 = (x2**2+y2**2)**1.5 #star-spacecraft distance
       r213 = ((x2-x1)**2+(y2-y1)**2)**1.5 #planet-spacecraft distance
       return array([vx1,vy1,vx2,vy2,-M*x1/r13,-M*y1/r13,-M*x2/r23-m*(x2-x1)/r213,-M*y2/r23-m*(y2-y1)/r213])

# initial conditions
Y = array([0.25,0.0,0.289118496229226,0.0,0.0,2.0,0.0,2.31244137],float)

#Solve for these tpoints
ti = 0.0
tf = 0.83
N=1000
H=(tf-ti)/N
tpoints = arange(ti,tf,H)
# for fixed step methods use 40 substeps
# for adapative use param=1e6 or tol=1e-6

#Benchmark 2: Forced Damped Pendulum
#Damped pendulum that receives a period impulse kick.

g=1.0
l=1.0

impulsetime=0.01
def force(t):
   return exp(-(1+cos(t*0.25))**2/(2*impulsetime)**2)

def f(Y,t):
       theta = Y[0]
       omega = Y[1]
       return array([omega,-g/l*sin(theta)-0.1*omega+force(t)])

# init conditions
theta0 = 3.14159
omega0 = -0.01
Yi= array([theta0,omega0],float)

# solve for these times
a = 0.0
b = 200.0
N = 400         # Number of Big Steps
H = (b-a)/N     # Size of Big Steps
tpoints = arange(a,b,H)
# use 8 substeps for fixed point methods
# param = 1e4 or tol=1d-4 for adaptive