1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83
| #define S_FUNCTION_NAME adm1_ODE
#include "simstruc.h"
#include <math.h>
#define XINIT ssGetSFcnParam(S,0)
#define PAR ssGetSFcnParam(S,1)
#define V ssGetSFcnParam(S,2)
/*
* mdlInitializeSizes - initialize the sizes array
*/
static void mdlInitializeSizes(SimStruct * S)
{
ssSetNumContStates( S, 42); /* number of continuous states */
ssSetNumDiscStates( S, 0); /* number of discrete states */
ssSetNumInputPorts( S, 33); /* number of inputs */
ssSetNumOutputPorts( S, 51); /* number of outputs */
ssSetDirectFeedThrough(S, 1); /* direct feedthrough flag */
ssSetNumSampleTimes( S, 1); /* number of sample times */
ssSetNumSFcnParams( S, 3); /* number of input arguments */
ssSetNumRWork( S, 0); /* number of real work vector elements */
ssSetNumIWork( S, 0); /* number of integer work vector elements*/
ssSetNumPWork( S, 0); /* number of pointer work vector elements*/
}
/*
* mdlInitializeSampleTimes - initialize the sample times array
*/
static void mdlInitializeSampleTimes(SimStruct *S)
{
ssSetSampleTime(S, 0, CONTINUOUS_SAMPLE_TIME);
ssSetOffsetTime(S, 0, 0.0);
}
/*
* mdlInitializeConditions - initialize the states
*/
static void mdlInitializeConditions(double *x0, SimStruct *S)
{
int i;
for (i = 0; i < 42; i++) {
x0[i] = mxGetPr(XINIT)[i];
}
}
/*
* mdlOutputs - compute the outputs
*/
static void mdlOutputs(double *y, double *x, double *u, SimStruct *S, int tid)
{
double R, T_base, T_op, P_atm, p_gas_h2o, P_gas, k_P, q_gas, V_liq, procT8, procT9, procT10, p_gas_h2, p_gas_ch4, p_gas_co2;
double kLa, K_H_h2_base, K_H_ch4_base, K_H_co2_base, phi, S_H_ion, pK_w_base, K_H_h2o_base, K_H_h2, K_H_ch4, K_H_co2, K_w, factor;
int i;
R = mxGetPr(PAR)[77];
T_base = mxGetPr(PAR)[78];
T_op = mxGetPr(PAR)[79];
P_atm = mxGetPr(PAR)[93];
V_liq = mxGetPr(V)[0];
kLa = mxGetPr(PAR)[94];
K_H_h2_base = mxGetPr(PAR)[98];
K_H_ch4_base = mxGetPr(PAR)[97];
K_H_co2_base = mxGetPr(PAR)[96];
K_H_h2o_base = mxGetPr(PAR)[95];
pK_w_base = mxGetPr(PAR)[80];
k_P = mxGetPr(PAR)[99];
factor = (1.0/T_base - 1.0/T_op)/(100.0*R);
K_H_h2 = K_H_h2_base*exp(-4180.0*factor); /* T adjustment for K_H_h2 */
K_H_ch4 = K_H_ch4_base*exp(-14240.0*factor); /* T adjustment for K_H_ch4 */
K_H_co2 = K_H_co2_base*exp(-19410.0*factor); /* T adjustment for K_H_co2 */
K_w = pow(10,-pK_w_base)*exp(55900.0*factor); /* T adjustment for K_w */
p_gas_h2o = K_H_h2o_base*exp(5290.0*(1.0/T_base - 1.0/T_op)); /* T adjustement for water vapour saturation pressure */
for (i = 0; i < 26; i++) {
y[i] = x[i];
} |
Partager