[2696] | 1 | void StochasticRates( double RCT[], double Volume, double SCT[] ); |
---|
| 2 | void Propensity ( int V[], int F[], double SCT[], double A[] ); |
---|
| 3 | void MoleculeChange ( int j, int NmlcV[] ); |
---|
| 4 | double CellMass(double T); |
---|
| 5 | void Update_RCONST(); |
---|
| 6 | |
---|
| 7 | /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ |
---|
| 8 | void 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 */ |
---|