chem/KPP/kpp/kpp-2.1/int/oldies/exqssa.c

   1 void INTEGRATE( double DT )
   2 {
   3 KPP_REAL P_VAR[NVAR], D_VAR[NVAR], V1[NVAR], V2[NVAR];
   4 int IsReject;
   5 KPP_REAL T, Tnext, STEP, STEPold, Told, SUP;
   6 KPP_REAL ERR, ERRold, ratio, factor, facmax, tmp;
   7 int i;
   8
   9   T = TIME;
  10   Tnext = TIME + DT;
  11   STEP = STEPMIN;
  12   Told = T;
  13   SUP  = 1e-14;
  14   IsReject = 0;
  15   ERR = .5;
  16
  17 /* -- BELOW THIS LIMIT USE TAYLOR INSTEAD OF EXP --- */
  18
  19   while ( T < Tnext ) {
  20
  21     T = Told + STEP;
  22     if ( T > Tnext ) {
  23       STEP = Tnext - Told;
  24       T = Tnext;
  25     }
  26
  27     FSPLIT_VAR ( VAR,  P_VAR, D_VAR );
  28
  29     for( i = 0; i < NVAR; i++ ) {
  30       if ( fabs(D_VAR[i]) > SUP ) {
  31         ratio = P_VAR[i] / D_VAR[i];
  32         tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) );
  33         V1[i] = tmp * tmp * (VAR[i] - ratio) + ratio;
  34         V2[i] = tmp * (VAR[i] - ratio) + ratio;
  35       } else {
  36         tmp = D_VAR[i] * STEP * 0.5;
  37         V1[i] = VAR[i] + P_VAR[i] * STEP * ( 1 - tmp *
  38                ( 1 - 2.0 / 3.0 * tmp ) );
  39         V2[i] = VAR[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp *
  40                ( 1 - 1.0 / 3.0 * tmp ) );
  41       }
  42     }
  43
  44     FSPLIT_VAR( V2,  P_VAR, D_VAR );
  45
  46     for( i = 0; i < NVAR; i++ ) {
  47       if ( fabs(D_VAR[i]) > SUP ) {
  48         ratio = P_VAR[i] / D_VAR[i];
  49         tmp = (KPP_REAL)exp( (double)(-D_VAR[i] * STEP * 0.5) );
  50         V2[i] = tmp * (V2[i] - ratio) + ratio;
  51       } else {
  52         tmp = D_VAR[i] * STEP * 0.5;
  53         V2[i] = V2[i] + P_VAR[i] * 0.5 * STEP * ( 1 - 0.5 * tmp *
  54                ( 1 - 1.0 / 3.0 * tmp ) );
  55       }
  56     }
  57 /* ==== Extrapolation and error estimation ======== */
  58
  59     ERRold=ERR;
  60     ERR=0.;
  61     for( i = 0; i < NVAR; i++ ) {
  62       V1[i] = 2.*V2[i] - V1[i];
  63       tmp = (V2[i] - V1[i]) / (ATOL[i] + RTOL[i]*V2[i]);
  64       ERR = ERR + tmp*tmp;
  65     }
  66     ERR = sqrt(ERR/NVAR);
  67     STEPold = STEP;
  68
  69 /* ===== choosing the stepsize ==================== */
  70
  71     factor = 0.9 / pow(ERR,0.35) * pow(ERRold,0.2);
  72     facmax = IsReject ? 1.0 : 5.0;
  73
  74     factor = max( 0.2, min(factor,facmax) );
  75     STEP = min( STEPMAX, max(STEPMIN,factor*STEP) );
  76
  77 /*================================================= */
  78
  79     if ( (ERR > 1) && (STEPold > STEPMIN) ) {
  80       T = Told;
  81       IsReject = 1;
  82     } else {
  83       IsReject = 0;
  84       Told = T;
  85       for( i = 0; i < NVAR; i++ )
  86         VAR[i] = max( V1[i], 0.0 );
  87       TIME = Tnext;
  88     }
  89   }
  90 }