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