Subway ride (Muscod)
From mintOC
Revision as of 19:00, 21 November 2010 by SebastianSager (Talk | contribs)
The differential equations for the Subway ride 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]; }
