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

Last change on this file since 2898 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 1.5 KB
Line 
1void StochasticRates( double RCT[], double Volume, double SCT[] );   
2void Propensity ( int V[], int F[], double SCT[], double A[] );
3void MoleculeChange ( int j, int NmlcV[] );
4double CellMass(double T); 
5void Update_RCONST();
6
7/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
8void Gillespie(int Nevents, double Volume, double* T, int NmlcV[], int NmlcF[]) 
9/* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
10{
11   
12      int i, m=0, event;
13      double r1, r2;
14      double A[NREACT], SCT[NREACT], x;
15
16     /* Compute the stochastic reaction rates */
17      Update_RCONST();
18      StochasticRates( RCONST, Volume, SCT );   
19   
20      for (event = 1; event <= Nevents; event++) {
21
22          /* Uniformly distributed ranfor (m numbers */
23          r1 = (double)rand()/(double)RAND_MAX;
24          r2 = (double)rand()/(double)RAND_MAX;
25
26          /* Avoid log of zero */
27          r2 = (r2-1.0e-14) ? r2 : 1.0e-14;
28         
29          /* Propensity vector */
30          TIME = *T;
31          Propensity ( NmlcV, NmlcF, SCT, A );
32         
33          /* Cumulative sum of propensities */
34          for (i=1; i<NREACT; i++)
35            A[i] = A[i-1]+A[i];
36         
37          /* Index of next reaction */
38          x = r1*A[NREACT-1];
39          for ( i = 0; i<NREACT; i++)
40            if (A[i] >= x) {
41              m = i+1;
42              break;
43            }
44         
45          /* Update T with time to next reaction */
46          *T = *T - log(r2)/A[NREACT-1];
47
48          /* Update state vector after reaction m */
49          MoleculeChange( m, NmlcV );
50         
51     } /* for event */
52   
53} /* Gillespie */
Note: See TracBrowser for help on using the repository browser.