/*
 *  pendulum_cv.c
 *
 *  CVODE C file for the vector field named: pendulum
 *
 *  This file was generated by the program VFGEN (Version:2.4.0)
 *  Generated on 10-Jul-2008 at 11:48
 */

#include <math.h>
/* Include headers for SUNDIALS v2.3.0, CVODE v2.5.0 */
#include <sundials/sundials_types.h>
#include <sundials/sundials_dense.h>
#include <sundials/sundials_nvector.h>
#include <nvector/nvector_serial.h>
#include <cvode/cvode.h>
#include <cvode/cvode_dense.h>
/*
 *  The vector field.
 */

int pendulum_vf(realtype t, N_Vector y_, N_Vector f_, void *params)
    {
    const realtype Pi = RCONST(M_PI);
    realtype theta, v;
    realtype g, b, L, m;
    realtype *p_;

    p_ = (realtype *) params;

    theta      = NV_Ith_S(y_,0);
    v          = NV_Ith_S(y_,1);

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

    NV_Ith_S(f_,0) = v;
    NV_Ith_S(f_,1) = -1.0/m*v/(L*L)*b-sin(theta)*g/L;
    return 0;
    }

/*
 *  The Jacobian.
 */

int pendulum_jac(long int N_, DenseMat jac_, realtype t,
                N_Vector y_, N_Vector fy_, void *params,
                N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
    {
    const realtype Pi = RCONST(M_PI);
    realtype theta, v;
    realtype g, b, L, m;
    realtype *p_;

    p_ = (realtype *) params;

    theta      = NV_Ith_S(y_,0);
    v          = NV_Ith_S(y_,1);

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

    DENSE_ELEM(jac_, 0, 1) = 1.0;
    DENSE_ELEM(jac_, 1, 0) = -g/L*cos(theta);
    DENSE_ELEM(jac_, 1, 1) = -1.0/m/(L*L)*b;
    return 0;
    }

/*
 *  User-defined functions. 
 */

int pendulum_func(realtype t, N_Vector y_, realtype *func_, void *params)
    {
    const realtype Pi = RCONST(M_PI);
    realtype theta, v;
    realtype g, b, L, m;
    realtype *p_;

    p_ = (realtype *) params;

    theta      = NV_Ith_S(y_,0);
    v          = NV_Ith_S(y_,1);

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

    /* energy:  */
    func_[0] = -m*g*L*cos(theta)+m*(v*v)*(L*L)/2.0;
    return 0;
    }