Subway ride (Muscod)
The differential equations for the Lotka Volterra fishing problem in C code read as follows
#define V4 21.0/GAMMA
#define S4 1000.0
#define V5 22.0/GAMMA
// LOCAL
#define S 2112.0
#define TEND 65.0
//#define TEND 70.0
// EXPRESS
//#define S 6336.0
//#define TEND 151.0
/* User changes end */
// Weight
#define W 78000.0
#define WEFF (W + 72000.0 * 0.1)
/* 3600.0 / 5280.0 */
#define GAMMA (3600.0 / 5280.0)
#define A 100.0
#define N 10.0
#define B 0.045
/* [ 0.24 + 0.034 * (N-1) ] / (100*N) */
#define C 0.000546
#define CC (0.367 * GAMMA)
#define G (32.2 * GAMMA)
#define E 1.0
#define V1 (1.03 - (W - 72000) / ( 110000 - 72000) * (1.03 - 0.71) ) / GAMMA
#define V2 (6.86 - (W - 72000) / ( 110000 - 72000) * (6.86 - 6.05) ) / GAMMA
#define V3 (14.49 - (W - 72000) / ( 110000 - 72000) * (14.49 - 13.07) ) / GAMMA
#define A1 (5998.6162 + (W - 72000) / ( 110000 - 72000) * (6118.9179 - 5998.6162) )
#define A2 (11440.7968 + (W - 72000) / ( 110000 - 72000) * (17188.6252 - 11440.7968) )
#define A3 (10280.0514 + (W - 72000) / ( 110000 - 72000) * (15629.0954 - 10280.0514) )
#define P1 (105.880645 + (W - 72000) / ( 110000 - 72000) * (107.872258 - 105.880645) )
#define P2 (168.931957 + (W - 72000) / ( 110000 - 72000) * (245.209888 - 168.931957) )
#define P3 (334.626716 + (W - 72000) / ( 110000 - 72000) * (458.188550 - 334.626716) )
#define B01 -0.1983670410E02
#define B11 0.1952738055E03
#define B21 0.2061789974E04
#define B31 -0.7684409308E03
#define B41 0.2677869201E03
#define B51 -0.3159629687E02
#define B02 -0.1577169936E03
#define B12 0.3389010339E04
#define B22 0.6202054610E04
#define B32 -0.4608734450E04
#define B42 0.2207757061E04
#define B52 -0.3673344160E03
#define C01 0.3629738340E02
#define C11 -0.2115281047E03
#define C21 0.7488955419E03
#define C31 -0.9511076467E03
#define C41 0.5710015123E03
#define C51 -0.1221306465E03
#define C02 0.4120568887E02
#define C12 0.3408049202E03
#define C22 -0.1436283271E03
#define C32 0.8108316584E02
#define C42 -0.5689703073E01
#define C52 -0.2191905731E01
/* -------------------------------------------- */
static void mfcn(double *ts, double *xd, double *xa, double *p, double *pr,
double *mval, long *dpnd, long *info)
{
if (*dpnd) { *dpnd = MFCN_DPND(0, 0, *xd, 0, 0); return; }
*mval = xd[2];
}
/* -------------------------------------------- */
/* Seriell */
static void ffcn1a( double *xd, double *rhs )
{
rhs[0] = xd[1];
rhs[1] = G * E * A1 / WEFF / GAMMA;
rhs[2] = P1;
}
static void ffcn1b( double *xd, double *rhs )
{
rhs[0] = xd[1];
rhs[1] = G * E * A2 / WEFF / GAMMA;
rhs[2] = P2;
}
static void ffcn1c( double *xd, double *rhs )
{
double T, R, temp, temp2, tempb, tempb2;
rhs[0] = xd[1];
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 0.3);
temp2 = temp * temp;
T = B01
+ B11 * temp
+ B21 * temp2
+ B31 * temp2 * temp
+ B41 * temp2 * temp2
+ B51 * temp2 * temp2 * temp
;
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;
tempb = 1.0 / ( 0.1 * GAMMA * xd[1]);
tempb2 = tempb * tempb;
rhs[2] = ( C01
+ C11 * tempb
+ C21 * tempb2
+ C31 * tempb2 * tempb
+ C41 * tempb2 * tempb2
+ C51 * tempb2 * tempb2 * tempb )
;
if (xd[1] < V2 + 1e-2) {
rhs[1] = G * E * A2 / WEFF / GAMMA;
rhs[2] = P2;
}
}
/* -------------------------------------------- */
/* Parallel */
static void ffcn2b( double *xd, double *rhs )
{
rhs[0] = xd[1];
rhs[1] = G * E * A3 / WEFF / GAMMA;
rhs[2] = P3;
}
static void ffcn2c( double *xd, double *rhs )
{
double T, R, temp, temp2, tempb, tempb2;
rhs[0] = xd[1];
temp = 1.0 / ( 0.1 * GAMMA * xd[1] - 1.0);
temp2 = temp * temp;
T = B02
+ B12 * temp
+ B22 * temp2
+ B32 * temp2 * temp
+ B42 * temp2 * temp2
+ B52 * temp2 * temp2 * temp
;
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;
rhs[1] = G * ( E * T - R ) / WEFF / GAMMA;
tempb = 1.0 / ( 0.1 * GAMMA * xd[1] - 1);
tempb2 = tempb * tempb;
rhs[2] = ( C02
+ C12 * tempb
+ C22 * tempb2
+ C32 * tempb2 * tempb
+ C42 * tempb2 * tempb2
+ C52 * tempb2 * tempb2 * tempb )
;
// Achtung
if (xd[1] < V3 + 1e-2) {
rhs[1] = G * E * A3 / WEFF / GAMMA;
rhs[2] = P3;
}
}
/* -------------------------------------------- */
/* Rollen */
static void ffcn3( double *xd, double *rhs )
{
double R;
rhs[0] = xd[1];
R = C*A*GAMMA*GAMMA*xd[1]*xd[1] + B*W*GAMMA*xd[1] / 2000 + 1.3 / 2000 * W + 116;
rhs[1] = ( - G * R / WEFF - CC / GAMMA ) ;
rhs[2] = 0;
}
/* -------------------------------------------- */
/* Bremsen */
static void ffcn4( double *xd, double *rhs )
{
rhs[0] = xd[1];
rhs[1] = - 3.0 / GAMMA;
rhs[2] = 0;
}
/* -------------------------------------------- */
/* Overall right hand side for v < V1*/
static void ffcn_A(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, long *info)
{
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];
long i;
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }
ffcn1a( xd, rhs0 );
ffcn3( xd, rhs2 );
ffcn4( xd, rhs3 );
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];
}
/* -------------------------------------------- */
/* Overall right hand side for V1 < v < V2*/
static void ffcn_B(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, long *info)
{
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];
long i;
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }
ffcn1b( xd, rhs0 );
ffcn3( xd, rhs2 );
ffcn4( xd, rhs3 );
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];
}
/* -------------------------------------------- */
/* Overall right hand side for V2 < v < V3*/
static void ffcn_C(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, long *info)
{
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];
long i;
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }
ffcn1c( xd, rhs0 );
ffcn2b( xd, rhs1 );
ffcn3( xd, rhs2 );
ffcn4( xd, rhs3 );
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];
}
/* -------------------------------------------- */
/* Overall right hand side for V3 < v */
static void ffcn_D(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, long *info)
{
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];
long i;
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }
ffcn1c( xd, rhs0 );
ffcn2c( xd, rhs1 );
ffcn3( xd, rhs2 );
ffcn4( xd, rhs3 );
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];
}
/* -------------------------------------------- */
/* Pure coasting resp. braking on last stage to avoid nondifferentiabilities*/
static void ffcn_E(double *t, double *xd, double *xa, double *u,
double *p, double *rhs, double *rwh, long *iwh, long *info)
{
double rhs0[NXD], rhs1[NXD], rhs2[NXD], rhs3[NXD];
long i;
for(i=0;i<NXD;i++) { rhs0[i] = rhs1[i] = rhs2[i] = rhs3[i] = 0; }
ffcn3( xd, rhs2 );
ffcn4( xd, rhs3 );
rhs[0] = u[0]*rhs0[0] + u[1]*rhs1[0] + u[2]*rhs2[0] + u[3]*rhs3[0];
rhs[1] = u[0]*rhs0[1] + u[1]*rhs1[1] + u[2]*rhs2[1] + u[3]*rhs3[1];
rhs[2] = u[0]*rhs0[2] + u[1]*rhs1[2] + u[2]*rhs2[2] + u[3]*rhs3[2];
}