#
# pendulum_demo.py
#
#
# Python ODE solver for the vector field named: pendulum
#
# This script uses the scipy odeint function to solve the differential equations.
#
# This file was generated by the program VFGEN (Version:2.4.0)
# Generated on 16-Jul-2008 at 16:43
#

import sys
from math import pi
from scipy.integrate import odeint
import pendulum


def use():
    print 'use: ',sys.argv[0],' [options]'
    print 'options:'
    print '    -h    Print this message.'
    for i in range(N_):
        print '    '+varnames_[i]+'=<initial_condition>  Default value is ', def_y_[i]
    for i in range(P_):
        print '    '+parnames_[i]+'=<parameter_value>    Default value is ', def_p_[i]
    print '    abserr=<absolute_error_tolerance>   Default value is ', abserr
    print '    relerr=<relative_error_tolerance>   Default value is ', relerr
    print '    stoptime=<stop_time>                Default value is ', stoptime
    print '    numpoints=<number_of_output_points> Default value is ', numpoints

#
# Main script begins here...
#
Pi = pi
N_ = 2
P_ = 4
# Default values for the initial conditions, parameters and solver parameters
def_y_ = [ Pi-1.0000000000000000e-02, 0.0000000000000000e+00]
def_p_ = [9.8100000000000005e+00, 0.0000000000000000e+00, 1.0000000000000000e+00, 1.0000000000000000e+00]
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 250
varnames_ = ["theta", "v"]
parnames_ = ["g", "b", "L", "m"]

# Create a dict of all the options that can be given on the command line.
# Set the values to the default value of option.
options = {}
for i in range(N_):
    options[varnames_[i]] = def_y_[i]
for i in range(P_):
    options[parnames_[i]] = def_p_[i]
options['abserr'] = abserr
options['relerr'] = relerr
options['stoptime'] = stoptime
options['numpoints'] = numpoints

# Process the command line arguments.
for a in sys.argv[1:]:
    if a == '-h' or a == '-help' or a == '--help' or a == 'help':
        use()
        sys.exit()
    eqloc = a.find('=')
    if (eqloc == -1):
        print 'Invalid argument (missing =): ', a
        use()
        sys.exit()
    else:
        var = a[0:eqloc]
        val = a[eqloc+1:]
        if var in options:
            options[var] = float(val)
        else:
            print 'Unknown argument: ', a
            use()
            sys.exit()

# Get the values for the initial conditions and parameters from the options dict.
y_ = []
for i in range(N_):
    y_ = y_ + [options[varnames_[i]]]
p_ = []
for i in range(P_):
    p_ = p_ + [options[parnames_[i]]]

# Create the time samples for the output of the ODE solver.
tfinal = options['stoptime']
N = int(options['numpoints'])
if N < 2:
    print 'The number of points must be at least 2.'
    sys.exit()
t = [tfinal*float(i)/(N-1) for i in range(N)]

# Call the ODE solver.
ysol = odeint(pendulum.vectorfield,y_,t,args=(p_,),Dfun=pendulum.jacobian,atol=options['abserr'],rtol=options['relerr'])

# Print the solution.
for t1,y1 in zip(t,ysol):
    print t1,y1[0],y1[1],
    print pendulum.energy(y1,t1,p_)