#include #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; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */