#include <math.h>
#include "auto_f2c.h"

#define PI 3.14159265358979

#define EPSILON par[0]
#define A       par[1]
#define OMEGA   par[2]
#define Y_0      par[3]
#define PHI_0    par[4]
#define T       par[10]


#define X     u[0]
#define Y     u[1]
#define PHI   u[2]
#define dx1   u[3]   /* dx1 = dx/dy0 */  
#define dy1   u[4]   /* dy1 = dy/dy0 */
#define dphi1 u[5]   /* dphi1 = dphi/dy0 */
#define dx2   u[6]   /* dx2 = dx/dtheta0 */



/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int func (integer ndim, const doublereal *u, const integer *icp,
          const doublereal *par, integer ijac,
          doublereal *f, doublereal *dfdu, doublereal *dfdp)
{
  integer dfdu_dim1, dfdp_dim1;

  dfdu_dim1 = ndim;
  dfdp_dim1 = ndim;

  /* f[0] - f[2] are the original FVDP */
  f[0] = T * (Y + X - X*X*X/3)/EPSILON;
  f[1] = T * (-X + A*sin(PHI));
  f[2] = T * (2.0*PI*OMEGA) ;

  /* d[dx/dy0]/dt */
  f[3] = T * (1/EPSILON) * ((1-X*X)*u[3] + u[4]);
   
  /* d[dx/dtheta0]/dt */
  f[4] = T * (1/EPSILON) * ((1-X*X)*u[6] + u[7]);

  /* d[dy/dy0]/dt */
  f[5] = -T*u[3]; 

  /* d[dy/dtheta0]/dt */
  f[6] = T* (-1*u[6] + A*cos(PHI));  

  
  /*  if (ijac == 0)
     {
     return 0;
     }
  */
  /* Jacobian */

  /*  ARRAY2D(dfdu,0,0) = (1.0-X*X)/EPSILON;
  ARRAY2D(dfdu,0,1) = 1.0/EPSILON;
  ARRAY2D(dfdu,0,2) = (doublereal) 0.0;
  ARRAY2D(dfdu,1,0) = (doublereal) -1.0;
  ARRAY2D(dfdu,1,1) = (doublereal) 0.0;
  ARRAY2D(dfdu,1,2) = A*cos(PHI);
  ARRAY2D(dfdu,2,0) = (doublereal) 0.0;
  ARRAY2D(dfdu,2,1) = (doublereal) 0.0;
  ARRAY2D(dfdu,2,2) = (doublereal) 0.0;
  */

  /* if (ijac == 1)
     {
     return 0;
     }
  */
  /* Parameter derivatives */

  /* ARRAY2D(dfdp,0,0) = -(Y+X-X*X*X/3)/(EPSILON*EPSILON);
  ARRAY2D(dfdp,0,1) = (doublereal) 0.0;
  ARRAY2D(dfdp,0,2) = (doublereal) 0.0;
  ARRAY2D(dfdp,1,0) = (doublereal) 0.0;
  ARRAY2D(dfdp,1,1) = sin(PHI);
  ARRAY2D(dfdp,1,2) = (doublereal) 0.0;
  ARRAY2D(dfdp,2,0) = (doublereal) 0.0;
  ARRAY2D(dfdp,2,1) = (doublereal) 0.0;
  ARRAY2D(dfdp,2,2) = (doublereal) 2.0*PI;
  */
  return 0;
}

/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int stpnt (integer ndim, doublereal t,
           doublereal *u, doublereal *par)
{
  /* Parameter values for the starting orbit in fvdp.dat : */

  EPSILON = (doublereal) 0.1;
  A       = (doublereal) 1.1;
  OMEGA   = (doublereal) 1.55;
  Y_0      = (doublereal) -0.63;
  PHI_0    = (doublereal) 2.26194671058465;
  
  return 0;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int bcnd (integer ndim, const doublereal *par, const integer *icp,
          integer nbc, const doublereal *u0, const doublereal *u1, integer ijac,
          doublereal *fb, doublereal *dbc)
{
  
  fb[0] = u0[0]+1;  /* x(0) = -1 */
  fb[1] = u1[0]-1;  /* x(1) = 1 */
  fb[2] = u0[1]-Y_0; /* Y(0) = Y_0 */
  fb[3] = u0[2]-PHI_0; /* phi(0) = phi_0 */
  fb[4] = u0[3];    
  fb[5] = u0[4];
  fb[6] = u0[5]-1;
  fb[7] = u0[6];
  
  
  return 0;


}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int icnd (integer ndim, const doublereal *par, const integer *icp,
          integer nint, const doublereal *u, const doublereal *uold,
          const doublereal *udot, const doublereal *upold, integer ijac,
          doublereal *fi, doublereal *dint)
{
  return 0;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int fopt (integer ndim, const doublereal *u, const integer *icp,
          const doublereal *par, integer ijac,
          doublereal *fs, doublereal *dfdu, doublereal *dfdp)
{
  return 0;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */
int pvls (integer ndim, const doublereal *u,
          doublereal *par)
{
  return 0;
}
/* ---------------------------------------------------------------------- */
/* ---------------------------------------------------------------------- */