/*
 *  pendulum_gvf.c
 *
 *  GSL 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 13:45
 */

#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>

/*
 *  The vector field.
 */

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

    p_ = (double *) params;

    theta      = y_[0];
    v          = y_[1];

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

    f_[0] = v;
    f_[1] = -1.0/m*v/(L*L)*b-sin(theta)*g/L;

    return GSL_SUCCESS;
    }

/*
 *  The Jacobian.
 */

int pendulum_jac(double t, const double y_[], double *jac_, double *dfdt_, void *params)
    {
    const double Pi = M_PI;
    double theta, v;
    double g, b, L, m;
    double *p_;

    p_ = (double *) params;

    theta      = y_[0];
    v          = y_[1];

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

    gsl_matrix_view dfdy_mat = gsl_matrix_view_array(jac_,2,2);
    gsl_matrix *m_ = &dfdy_mat.matrix;

    gsl_matrix_set(m_, 0, 0, 0.0);
    gsl_matrix_set(m_, 0, 1, 1.0);
    gsl_matrix_set(m_, 1, 0, -g/L*cos(theta));
    gsl_matrix_set(m_, 1, 1, -1.0/m/(L*L)*b);

    dfdt_[0] = 0.0;
    dfdt_[1] = 0.0;

    return GSL_SUCCESS;
    }

/*
 *  The Jacobian with respect to the parameters.
 */

int pendulum_jacp(double t, const double y_[], double *jacp_, void *params)
    {
    const double Pi = M_PI;
    double theta, v;
    double g, b, L, m;
    double *p_;

    p_ = (double *) params;

    theta      = y_[0];
    v          = y_[1];

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

    gsl_matrix_view dfdp_mat = gsl_matrix_view_array(jacp_,2,4);
    gsl_matrix *m_ = &dfdp_mat.matrix;

    gsl_matrix_set(m_, 0, 0, 0.0);
    gsl_matrix_set(m_, 0, 1, 0.0);
    gsl_matrix_set(m_, 0, 2, 0.0);
    gsl_matrix_set(m_, 0, 3, 0.0);
    gsl_matrix_set(m_, 1, 0, -sin(theta)/L);
    gsl_matrix_set(m_, 1, 1, -1.0/m*v/(L*L));
    gsl_matrix_set(m_, 1, 2,  sin(theta)*g/(L*L)+2.0*1.0/m*v/(L*L*L)*b);
    gsl_matrix_set(m_, 1, 3, 1.0/(m*m)*v/(L*L)*b);

    return GSL_SUCCESS;
    }

/*
 *  User function: energy
 */

double pendulum_energy(double t, const double y_[], void *params)
    {
    const double Pi = M_PI;
    double theta, v;
    double g, b, L, m;
    double *p_;

    p_ = (double *) params;

    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;
    }