chem/KPP/kpp/kpp-2.1/util/mex.c

   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