#
# pendulum.py
#
# Python file for the vector field named: pendulum
#

"""
This module implements the vector field name "pendulum" as
the function vectorfield().  The Jacobian of the vector field
is computed by jacobian().  These functions can be used with
the SciPy odeint function.
This module also defines the function energy().

For example:

    from scipy.integrate import odeint
    import pendulum

    params = [g,b,L,m]   # Assume the parameters have been set elsewhere
    t = [i/10.0 for i in range(0,101)]
    ic = [1.0,0.0]
    sol = odeint(pendulum.vectorfield,ic,t,args=(params,),Dfun=pendulum.jacobian)

This file was generated by the program VFGEN (Version:2.4.0)
Generated on 16-Jul-2008 at 16:43

"""

from math import *
import numpy

#
# The vector field.
#

def vectorfield(y_,t_,p_):
    """
    The vector field function for the vector field "pendulum"
    Arguments:
        y_ :  vector of the state variables:
                  y_[0] is theta
                  y_[1] is v
        t_ :  time
        p_ :  vector of the parameters
                  p_[0] is g
                  p_[1] is b
                  p_[2] is L
                  p_[3] is m
    """
    Pi = pi
    theta      = y_[0]
    v          = y_[1]

    g          = p_[0]
    b          = p_[1]
    L          = p_[2]
    m          = p_[3]

    f_ = numpy.zeros((2,))
    f_[0] = v
    f_[1] = -1.0/m*v/(L*L)*b-sin(theta)*g/L

    return f_

#
#  The Jacobian.
#

def jacobian(y_, t_, p_):
    """
    The Jacobian of the vector field "pendulum"
    Arguments:
        y_ :  vector of the state variables:
                  y_[0] is theta
                  y_[1] is v
        t_ :  time
        p_ :  vector of the parameters
                  p_[0] is g
                  p_[1] is b
                  p_[2] is L
                  p_[3] is m
    """
    Pi = pi
    theta      = y_[0]
    v          = y_[1]
    g          = p_[0]
    b          = p_[1]
    L          = p_[2]
    m          = p_[3]

    # Create the Jacobian matrix:
    jac_ = numpy.zeros((2,2))
    jac_[0,1] = 1.0
    jac_[1,0] = -g/L*cos(theta)
    jac_[1,1] = -1.0/m/(L*L)*b
    return jac_


#
# User function: energy
#

def energy(y_, t_, p_):
    """
    The user-defined function "energy" for the vector field "pendulum"
    Arguments:
        y_ :  vector of the state variables:
                  y_[0] is theta
                  y_[1] is v
        t_ :  time
        p_ :  vector of the parameters
                  p_[0] is g
                  p_[1] is b
                  p_[2] is L
                  p_[3] is m
    """
    Pi = pi
    theta      = y_[0]
    v          = y_[1]

    g          = p_[0]
    b          = p_[1]
    L          = p_[2]
    m          = p_[3]

    return -m*g*L*cos(theta)+m*(v*v)*(L*L)/2.0