#
# pendulum_demo.py
#
#
# Python ODE solver for the vector field named: pendulum
#
# This script uses the PyGSL odeiv functions to solve the differential equations.
#
# This file was generated by the program VFGEN (Version:2.4.0)
# Generated on 10-Jul-2008 at 17:45
#

import sys
from math import pi
from pygsl import odeiv
import pendulum

def use():
    print 'use: ',sys.argv[0],' [options]'
    print 'options:'
    print '    -h    Print this message.'
    for vname in varnames_:
        print '    '+vname+'=<initial_condition>  Default value is ',def_options[vname]
    for pname in parnames_:
        print '    '+pname+'=<parameter_value>    Default value is ',def_options[pname]
    print '    abserr=<absolute_error_tolerance>  Default value is ',def_options['abserr']
    print '    relerr=<relative_error_tolerance>  Default value is ',def_options['relerr']
    print '    stoptime=<stop_time>               Default value is ',def_options['stoptime']

#
# Main script begins here...
#
Pi = Numeric.pi
N_ = 2
varnames_ = ["theta", "v"]
parnames_ = ["g", "b", "L", "m"]
solver_param_names_ = ["abserr","relerr","stoptime"]

#
# Set the default values of everything in options dict.
#
options = {}
# Default initial conditions
options['theta'] =  Pi-1.0000000000000000e-02
options['v'] = 0.0000000000000000e+00
# Default vector field parameter values
options['g'] = 9.8100000000000005e+00
options['b'] = 0.0000000000000000e+00
options['L'] = 1.0000000000000000e+00
options['m'] = 1.0000000000000000e+00
# Default ODE solver parameters:
options['abserr'] = 1.0e-8
options['relerr'] = 1.0e-6
options['stoptime'] = 10.0

# Make a copy of the default options so use() can show the default values
def_options = options.copy()

#
# 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()

y_ = (options['theta'],options['v'])
p_ = (options['g'],options['b'],options['L'],options['m'])

# Create the GSL ODE solver
step    = odeiv.step_rk8pd(N_,pendulum.vectorfield,pendulum.jacobian,args=p_)
control = odeiv.control_y_new(step,options['abserr'],options['relerr'])
evolve  = odeiv.evolve(step,control,N_)

stoptime = options['stoptime']
# Initial step size is stoptime/500
h = stoptime/500.0
t = 0
# Call evolve.apply(...) until the solution reaches stoptime
while t < stoptime:
    t, h, y_ = evolve.apply(t,stoptime,h,y_)
    y_ = y_[0]
    print t,y_[0],y_[1],pendulum.energy(t,y_,args=p_)