| 1 | /* #include <math.h> */ |
|---|
| 2 | #include "mex.h" |
|---|
| 3 | #include <stdarg.h> |
|---|
| 4 | |
|---|
| 5 | #define MAX_BUF 200 |
|---|
| 6 | |
|---|
| 7 | void Usage() |
|---|
| 8 | { |
|---|
| 9 | mexPrintf(" \n" |
|---|
| 10 | "To get this help message use: KPP_ROOT ? \n" |
|---|
| 11 | " \n" |
|---|
| 12 | "To initialize default values use: KPP_ROOT \n" |
|---|
| 13 | " (type who to see the variables created) \n" |
|---|
| 14 | " \n" |
|---|
| 15 | "To integrate the model use: \n" |
|---|
| 16 | " [ c, m, f ] = KPP_ROOT( t, c0, k, p, fn, tfn ); \n" |
|---|
| 17 | " \n" |
|---|
| 18 | " input : \n" |
|---|
| 19 | " t - Time vector, contains the time at which results \n" |
|---|
| 20 | " should be reported; \n" |
|---|
| 21 | " c0 - Vector with the initial concentrations for all \n" |
|---|
| 22 | " species; \n" |
|---|
| 23 | " k - Vector with all rate constants; \n" |
|---|
| 24 | " p - Vector of parameters for the integration; \n" |
|---|
| 25 | " p(1) holds the relative tolerance \n" |
|---|
| 26 | " p(2) holds the absolute tolerance \n" |
|---|
| 27 | " p(3) holds the minimum step size allowed \n" |
|---|
| 28 | " p(4) holds the maximum step size allowed \n" |
|---|
| 29 | " If any of the above is zero the default value is \n" |
|---|
| 30 | " used; \n" |
|---|
| 31 | " fn - (optional) Name of a matlab function to be called \n" |
|---|
| 32 | " to update the values of k's and concentrations \n" |
|---|
| 33 | " If not present no update is performed. \n" |
|---|
| 34 | " \n" |
|---|
| 35 | " tfn - (optional) Time at which the fn function should \n" |
|---|
| 36 | " be called. If missing <t> is assumed. \n" |
|---|
| 37 | " \n" |
|---|
| 38 | " output: \n" |
|---|
| 39 | " c - Matrix of concentrations of all species vs. time; \n" |
|---|
| 40 | " m - (optional) Mass conservation of all atoms vs. time; \n" |
|---|
| 41 | " f - (optional) Matrix of fluxes of all species vs. time; \n" |
|---|
| 42 | " \n" |
|---|
| 43 | ); |
|---|
| 44 | } |
|---|
| 45 | |
|---|
| 46 | int giveusage; |
|---|
| 47 | |
|---|
| 48 | void F9Error( char *fmt, ... ) |
|---|
| 49 | { |
|---|
| 50 | va_list args; |
|---|
| 51 | char buf[ MAX_BUF ]; |
|---|
| 52 | char errmsg[ MAX_BUF ]; |
|---|
| 53 | |
|---|
| 54 | va_start( args, fmt ); |
|---|
| 55 | vsprintf( buf, fmt, args ); |
|---|
| 56 | va_end( args ); |
|---|
| 57 | |
|---|
| 58 | if( giveusage ) Usage(); |
|---|
| 59 | |
|---|
| 60 | mexPrintf("Error: %s\n", buf); |
|---|
| 61 | mexErrMsgTxt( 0 ); |
|---|
| 62 | } |
|---|
| 63 | |
|---|
| 64 | char allvars[1000]; |
|---|
| 65 | |
|---|
| 66 | int CreateVar(char *name, KPP_REAL val) |
|---|
| 67 | { |
|---|
| 68 | mxArray *GA; |
|---|
| 69 | double *pga; |
|---|
| 70 | |
|---|
| 71 | sprintf(allvars, "%s %s",allvars, name); |
|---|
| 72 | GA = mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 73 | pga = mxGetPr(GA); |
|---|
| 74 | *pga = (double)val; |
|---|
| 75 | mxSetName(GA,name); |
|---|
| 76 | return mexPutArray(GA,"global"); |
|---|
| 77 | } |
|---|
| 78 | |
|---|
| 79 | int CreateVec(char *name, int len, KPP_REAL *val) |
|---|
| 80 | { |
|---|
| 81 | mxArray *GA; |
|---|
| 82 | double *pga; |
|---|
| 83 | int i; |
|---|
| 84 | |
|---|
| 85 | sprintf(allvars, "%s %s",allvars, name); |
|---|
| 86 | GA = mxCreateDoubleMatrix(1,len,mxREAL); |
|---|
| 87 | pga = mxGetPr(GA); |
|---|
| 88 | if( sizeof(KPP_REAL) == sizeof(double) ) { |
|---|
| 89 | memmove( pga, val, len*sizeof(double) ); |
|---|
| 90 | } else { |
|---|
| 91 | for( i = 0; i < len; i++ ) pga[i] = (double)val[i]; |
|---|
| 92 | } |
|---|
| 93 | mxSetName(GA,name); |
|---|
| 94 | return mexPutArray(GA,"global"); |
|---|
| 95 | } |
|---|
| 96 | |
|---|
| 97 | int CreateStrVec(char *name, int len, char **val) |
|---|
| 98 | { |
|---|
| 99 | mxArray *GA; |
|---|
| 100 | |
|---|
| 101 | sprintf(allvars, "%s %s",allvars, name); |
|---|
| 102 | GA = mxCreateCharMatrixFromStrings( len, (const char **)val ); |
|---|
| 103 | mxSetName(GA,name); |
|---|
| 104 | return mexPutArray(GA,"global"); |
|---|
| 105 | } |
|---|
| 106 | |
|---|
| 107 | int CreateStr(char *name, char *val) |
|---|
| 108 | { |
|---|
| 109 | mxArray *GA; |
|---|
| 110 | |
|---|
| 111 | sprintf(allvars, "%s %s",allvars, name); |
|---|
| 112 | GA = mxCreateString( val ); |
|---|
| 113 | mxSetName(GA,name); |
|---|
| 114 | return mexPutArray(GA,"global"); |
|---|
| 115 | } |
|---|
| 116 | |
|---|
| 117 | #define T_PRM prhs[0] |
|---|
| 118 | #define C0_PRM prhs[1] |
|---|
| 119 | #define K_PRM prhs[2] |
|---|
| 120 | #define P_PRM prhs[3] |
|---|
| 121 | #define FN_PRM prhs[4] |
|---|
| 122 | #define TFN_PRM prhs[5] |
|---|
| 123 | |
|---|
| 124 | #define C_PRM plhs[0] |
|---|
| 125 | #define M_PRM plhs[1] |
|---|
| 126 | #define F_PRM plhs[2] |
|---|
| 127 | |
|---|
| 128 | #define HAS_FN (nrhs >= 5) |
|---|
| 129 | #define HAS_TFN (nrhs >= 6) |
|---|
| 130 | |
|---|
| 131 | #define HAS_M (nlhs >= 2) |
|---|
| 132 | #define HAS_F (nlhs >= 3) |
|---|
| 133 | |
|---|
| 134 | #define DBL (sizeof(KPP_REAL) == sizeof(double)) |
|---|
| 135 | |
|---|
| 136 | void mexFunction( |
|---|
| 137 | int nlhs, mxArray *plhs[], |
|---|
| 138 | int nrhs, const mxArray *prhs[] |
|---|
| 139 | ) |
|---|
| 140 | { |
|---|
| 141 | double * tp; |
|---|
| 142 | double * c0p; |
|---|
| 143 | double * kp; |
|---|
| 144 | double * pp; |
|---|
| 145 | char fnp[ MAX_BUF ]; |
|---|
| 146 | double *tfnp; |
|---|
| 147 | |
|---|
| 148 | double * cp; |
|---|
| 149 | double * mp; |
|---|
| 150 | double * fp; |
|---|
| 151 | double ATOLS; |
|---|
| 152 | double dval[ NMASS+NSPEC ]; |
|---|
| 153 | |
|---|
| 154 | mxArray *Carr, *Karr, *Tarr; |
|---|
| 155 | double *fcp; |
|---|
| 156 | double *fkp; |
|---|
| 157 | double *ftp, *ftp1; |
|---|
| 158 | |
|---|
| 159 | int i,j,m,n,nd,t; |
|---|
| 160 | int nsteps, nspc, nreact, ncb; |
|---|
| 161 | int tcb, CallBack; |
|---|
| 162 | KPP_REAL prm[4]; |
|---|
| 163 | |
|---|
| 164 | giveusage = 1; |
|---|
| 165 | |
|---|
| 166 | if(nrhs == 0) { |
|---|
| 167 | |
|---|
| 168 | InitVal(); |
|---|
| 169 | Update_RCONST(); |
|---|
| 170 | |
|---|
| 171 | prm[0] = 1e-4; |
|---|
| 172 | prm[1] = 1.0E-18; |
|---|
| 173 | prm[2] = 0.01; |
|---|
| 174 | prm[3] = 900; |
|---|
| 175 | |
|---|
| 176 | sprintf(allvars,"global "); |
|---|
| 177 | |
|---|
| 178 | CreateVec("PRM",4, prm); |
|---|
| 179 | |
|---|
| 180 | CreateVar("NSPEC",NSPEC); |
|---|
| 181 | CreateVar("NREACT",NREACT); |
|---|
| 182 | CreateVar("NMASS",NMASS); |
|---|
| 183 | |
|---|
| 184 | CreateVec("C0", NSPEC, C); |
|---|
| 185 | CreateVec("K0", NREACT, RCONST); |
|---|
| 186 | |
|---|
| 187 | for( i = 0; i < NLOOKAT; i++ ) |
|---|
| 188 | CreateVar( SLOOKAT[i], (double)(i+1) ); |
|---|
| 189 | |
|---|
| 190 | for( i = 0; i < NMASS; i++ ) |
|---|
| 191 | CreateVar( SMASS[i], (double)(i+1) ); |
|---|
| 192 | |
|---|
| 193 | CreateStrVec("SSPEC", NSPEC, SLOOKAT); |
|---|
| 194 | CreateStrVec("SMASS", NMASS, SMASS); |
|---|
| 195 | CreateStrVec("SEQN", NREACT, SEQN); |
|---|
| 196 | |
|---|
| 197 | CreateStr("GLOBALCMD", allvars); |
|---|
| 198 | |
|---|
| 199 | mexEvalString(allvars); |
|---|
| 200 | /* |
|---|
| 201 | mexPrintf("The KPP_ROOT model parameters were sucessfully initialized.\n"); |
|---|
| 202 | */ |
|---|
| 203 | return; |
|---|
| 204 | } |
|---|
| 205 | |
|---|
| 206 | if( nrhs < 4 ) |
|---|
| 207 | F9Error("First 4 parameters are REQUIRED only %d received.", nrhs); |
|---|
| 208 | if( nlhs < 1 ) |
|---|
| 209 | F9Error("At least one output parameter REQUIRED."); |
|---|
| 210 | |
|---|
| 211 | if(! mxIsDouble(T_PRM)) F9Error("<t> must be of type double."); |
|---|
| 212 | if(! mxIsDouble(C0_PRM)) F9Error("<c0> must be of type double."); |
|---|
| 213 | if(! mxIsDouble(K_PRM)) F9Error("<k> must be of type double."); |
|---|
| 214 | if(! mxIsDouble(P_PRM)) F9Error("<p> must be of type double."); |
|---|
| 215 | if((nrhs > 4) && (! mxIsChar(FN_PRM))) F9Error("<fn> must be of type char."); |
|---|
| 216 | if((nrhs > 5) && (! mxIsDouble(TFN_PRM))) F9Error("<tfn> must be of type double."); |
|---|
| 217 | |
|---|
| 218 | nd = mxGetNumberOfDimensions( T_PRM ); |
|---|
| 219 | m = mxGetM( T_PRM ); |
|---|
| 220 | n = mxGetN( T_PRM ); |
|---|
| 221 | if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<t> must be a column vector."); |
|---|
| 222 | nsteps = (m == 1) ? n : m; |
|---|
| 223 | tp = mxGetPr( T_PRM ); |
|---|
| 224 | |
|---|
| 225 | nd = mxGetNumberOfDimensions( C0_PRM ); |
|---|
| 226 | m = mxGetM( C0_PRM ); |
|---|
| 227 | n = mxGetN( C0_PRM ); |
|---|
| 228 | if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<c0> must be a column vector."); |
|---|
| 229 | nspc = (m == 1) ? n : m; |
|---|
| 230 | c0p = mxGetPr( C0_PRM ); |
|---|
| 231 | |
|---|
| 232 | nd = mxGetNumberOfDimensions( K_PRM ); |
|---|
| 233 | m = mxGetM( K_PRM ); |
|---|
| 234 | n = mxGetN( K_PRM ); |
|---|
| 235 | if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<k> must be a column vector."); |
|---|
| 236 | nreact = (m == 1) ? n : m; |
|---|
| 237 | kp = mxGetPr( K_PRM ); |
|---|
| 238 | |
|---|
| 239 | nd = mxGetNumberOfDimensions( P_PRM ); |
|---|
| 240 | m = mxGetM( P_PRM ); |
|---|
| 241 | n = mxGetN( P_PRM ); |
|---|
| 242 | if( !( (nd == 2) && ((m == 1) || (n == 1)) && (n*m == 4) ) ) |
|---|
| 243 | F9Error("<p> must be a column vectorof length 4."); |
|---|
| 244 | pp = mxGetPr( P_PRM ); |
|---|
| 245 | |
|---|
| 246 | *fnp = 0; |
|---|
| 247 | if( HAS_FN ) { |
|---|
| 248 | nd = mxGetNumberOfDimensions( FN_PRM ); |
|---|
| 249 | m = mxGetM( FN_PRM ); |
|---|
| 250 | n = mxGetN( FN_PRM ); |
|---|
| 251 | if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<fn> must be a character string."); |
|---|
| 252 | if( mxGetString( FN_PRM, fnp, MAX_BUF ) ) |
|---|
| 253 | F9Error("Can not read function mane (too long?)"); |
|---|
| 254 | |
|---|
| 255 | Carr = mxCreateDoubleMatrix(1,NSPEC,mxREAL); |
|---|
| 256 | fcp = mxGetPr(Carr); |
|---|
| 257 | mxSetName(Carr,"C"); |
|---|
| 258 | mexPutArray(Carr,"base"); |
|---|
| 259 | |
|---|
| 260 | Karr = mxCreateDoubleMatrix(1,NREACT,mxREAL); |
|---|
| 261 | fkp = mxGetPr(Karr); |
|---|
| 262 | mxSetName(Karr,"K"); |
|---|
| 263 | mexPutArray(Karr,"base"); |
|---|
| 264 | |
|---|
| 265 | Tarr = mxCreateDoubleMatrix(1,1,mxREAL); |
|---|
| 266 | ftp = mxGetPr(Tarr); |
|---|
| 267 | mxSetName(Tarr,"T"); |
|---|
| 268 | mexPutArray(Tarr,"base"); |
|---|
| 269 | } |
|---|
| 270 | |
|---|
| 271 | tfnp = 0; ncb = 0; |
|---|
| 272 | if( HAS_TFN ) { |
|---|
| 273 | nd = mxGetNumberOfDimensions( TFN_PRM ); |
|---|
| 274 | m = mxGetM( TFN_PRM ); |
|---|
| 275 | n = mxGetN( TFN_PRM ); |
|---|
| 276 | if( !( (nd == 2) && ((m == 1) || (n == 1)) ) ) F9Error("<tfn> must be a column vector."); |
|---|
| 277 | ncb = (m == 1) ? n : m; |
|---|
| 278 | tfnp = mxGetPr( TFN_PRM ); |
|---|
| 279 | } |
|---|
| 280 | |
|---|
| 281 | giveusage = 0; |
|---|
| 282 | |
|---|
| 283 | if( !((nspc == NSPEC) && (nreact == NREACT)) ) { |
|---|
| 284 | F9Error("Size of parameters do not match the model:\n\n" |
|---|
| 285 | " Number of species was %d and should be %d;\n" |
|---|
| 286 | " Number of rections (rate constants) was %d and should be %d;\n", |
|---|
| 287 | nspc, NSPEC, nreact, NREACT); |
|---|
| 288 | } |
|---|
| 289 | |
|---|
| 290 | if( DBL ) { memmove( C, c0p, sizeof(double)*NSPEC ); |
|---|
| 291 | memmove( RCONST, kp, sizeof(double)*NREACT ); } |
|---|
| 292 | else { for( i = 0; i < NSPEC; i++ ) C[i] = (KPP_REAL)c0p[i]; |
|---|
| 293 | for( i = 0; i < NREACT; i++ ) RCONST[i] = (KPP_REAL)kp[i]; } |
|---|
| 294 | |
|---|
| 295 | RTOLS = 1e-4; |
|---|
| 296 | ATOLS = 1e-18; |
|---|
| 297 | STEPMIN = 0.01; |
|---|
| 298 | STEPMAX = 900.0; |
|---|
| 299 | |
|---|
| 300 | if( pp[0] ) RTOLS = pp[0]; |
|---|
| 301 | if( pp[1] ) ATOLS = pp[1]; |
|---|
| 302 | if( pp[2] ) STEPMIN = pp[2]; |
|---|
| 303 | if( pp[3] ) STEPMAX = pp[3]; |
|---|
| 304 | |
|---|
| 305 | for( i = 0; i < NVAR; i++ ) { |
|---|
| 306 | RTOL[i] = RTOLS; |
|---|
| 307 | ATOL[i] = ATOLS; |
|---|
| 308 | } |
|---|
| 309 | |
|---|
| 310 | C_PRM = mxCreateDoubleMatrix(NSPEC,nsteps,mxREAL); |
|---|
| 311 | cp = mxGetPr(C_PRM); |
|---|
| 312 | |
|---|
| 313 | if( HAS_M ) { |
|---|
| 314 | M_PRM = mxCreateDoubleMatrix(NMASS,nsteps,mxREAL); |
|---|
| 315 | mp = mxGetPr(M_PRM); |
|---|
| 316 | } |
|---|
| 317 | |
|---|
| 318 | if( HAS_F ) { |
|---|
| 319 | F_PRM = mxCreateDoubleMatrix(NSPEC,nsteps,mxREAL); |
|---|
| 320 | fp = mxGetPr(F_PRM); |
|---|
| 321 | } |
|---|
| 322 | |
|---|
| 323 | tcb = 0; |
|---|
| 324 | |
|---|
| 325 | for( t = 0; t < nsteps; t++ ) { |
|---|
| 326 | if( t ) { |
|---|
| 327 | TIME = tp[t-1]; |
|---|
| 328 | |
|---|
| 329 | CallBack = 0; |
|---|
| 330 | if( HAS_TFN ) { |
|---|
| 331 | if( tcb < ncb ) |
|---|
| 332 | if( tfnp[tcb] <= TIME ) { CallBack = 1; tcb++; } |
|---|
| 333 | } else { |
|---|
| 334 | CallBack = HAS_FN; |
|---|
| 335 | } |
|---|
| 336 | |
|---|
| 337 | if( CallBack ) { |
|---|
| 338 | if( DBL ) { memmove( fcp, C, sizeof(double)*NSPEC ); |
|---|
| 339 | memmove( fkp, RCONST, sizeof(double)*NREACT ); } |
|---|
| 340 | else { for( i = 0; i < NSPEC; i++ ) fcp[i] = (double)C[i]; |
|---|
| 341 | for( i = 0; i < NREACT; i++ ) fkp[i] = (double)RCONST[i]; } |
|---|
| 342 | *ftp = TIME; |
|---|
| 343 | |
|---|
| 344 | mexPutArray(Carr,"base"); |
|---|
| 345 | mexPutArray(Karr,"base"); |
|---|
| 346 | mexPutArray(Tarr,"base"); |
|---|
| 347 | |
|---|
| 348 | mexCallMATLAB( 0, 0, 0, 0, fnp ); |
|---|
| 349 | |
|---|
| 350 | mxDestroyArray(Carr); Carr = mexGetArray("C","base"); fcp = mxGetPr(Carr); |
|---|
| 351 | mxDestroyArray(Karr); Karr = mexGetArray("K","base"); fkp = mxGetPr(Karr); |
|---|
| 352 | mxDestroyArray(Tarr); Tarr = mexGetArray("T","base"); ftp = mxGetPr(Tarr); |
|---|
| 353 | |
|---|
| 354 | if( DBL ) { memmove( C, fcp, sizeof(double)*NSPEC ); |
|---|
| 355 | memmove( RCONST, fkp, sizeof(double)*NREACT ); } |
|---|
| 356 | else { for( i = 0; i < NSPEC; i++ ) C[i] = (KPP_REAL)fcp[i]; |
|---|
| 357 | for( i = 0; i < NREACT; i++ ) RCONST[i] = (KPP_REAL)fkp[i]; } |
|---|
| 358 | |
|---|
| 359 | } |
|---|
| 360 | |
|---|
| 361 | INTEGRATE( tp[t-1], tp[t] ); |
|---|
| 362 | } |
|---|
| 363 | if( DBL ) { memmove( cp, C, sizeof(double)*NSPEC ); cp += NSPEC; } |
|---|
| 364 | else { for( i = 0; i < NSPEC; i++ ) *cp++ = (double)C[i]; } |
|---|
| 365 | if( HAS_M ) { |
|---|
| 366 | if( DBL ) { GetMass( mp ); mp += NMASS; } |
|---|
| 367 | else { GetMass( dval ); |
|---|
| 368 | for( i = 0; i < NMASS; i++ ) *mp++ = (double)dval[i]; } |
|---|
| 369 | } |
|---|
| 370 | if( HAS_F ) { |
|---|
| 371 | if( DBL ) { FLUX( fp ); fp += NSPEC; } |
|---|
| 372 | else { FLUX( dval ); |
|---|
| 373 | for( i = 0; i < NSPEC; i++ ) *fp++ = (double)dval[i]; } |
|---|
| 374 | } |
|---|
| 375 | } |
|---|
| 376 | } |
|---|
| 377 | |
|---|
