source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/int/gillespie.c @ 3821

void StochasticRates( double RCT[], double Volume, double SCT[] );   
void Propensity ( int V[], int F[], double SCT[], double A[] );
void MoleculeChange ( int j, int NmlcV[] );
double CellMass(double T); 
void Update_RCONST();

/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
void Gillespie(int Nevents, double Volume, double* T, int NmlcV[], int NmlcF[]) 
/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
{
    
      int i, m=0, event;
      double r1, r2;
      double A[NREACT], SCT[NREACT], x;

     /* Compute the stochastic reaction rates */
      Update_RCONST();
      StochasticRates( RCONST, Volume, SCT );   
   
      for (event = 1; event <= Nevents; event++) {

          /* Uniformly distributed ranfor (m numbers */
          r1 = (double)rand()/(double)RAND_MAX;
          r2 = (double)rand()/(double)RAND_MAX;

          /* Avoid log of zero */
          r2 = (r2-1.0e-14) ? r2 : 1.0e-14;
          
          /* Propensity vector */
          TIME = *T;
          Propensity ( NmlcV, NmlcF, SCT, A );
          
          /* Cumulative sum of propensities */
          for (i=1; i<NREACT; i++)
            A[i] = A[i-1]+A[i];
          
          /* Index of next reaction */
          x = r1*A[NREACT-1];
          for ( i = 0; i<NREACT; i++)
            if (A[i] >= x) {
              m = i+1;
              break;
            }
          
          /* Update T with time to next reaction */
          *T = *T - log(r2)/A[NREACT-1];

          /* Update state vector after reaction m */
          MoleculeChange( m, NmlcV );
          
     } /* for event */
    
} /* Gillespie */