View source for Subway ride (Muscod)

The differential equations for the [[Subway ride]] in [[:Category:C | C code]] read as follows

<source lang="cpp">

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

</source>

[[Category:C]]