Fix subtle rare bug involving next timestep = infinity. Also prevent first rxn firing...

diff --git a/ss.cpp b/ss.cpp
index ed39d8c..e6f4138 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -269,11 +269,11 @@ void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long
     // stop if time runs out or too many steps or no more reactions will occur
     for (; (max_steps == 0 || step < max_steps) && (t < t_end) && (total_propensity > 0.0); ++step) {
 
-        double rand = (1.0-gsl_rng_uniform_pos(rnd)) * total_propensity;  // uniform on (0,R]
+        double rand = (1.0-gsl_rng_uniform(rnd)) * total_propensity;  // uniform on (0,R], because reasons
         double a_sum = 0.0;
 
         // time to reaction
-        dt = -log(gsl_rng_uniform(rnd)) / total_propensity;
+        dt = -log(gsl_rng_uniform_pos(rnd)) / total_propensity;     // disallow 0.0 when calculating log(r), r in (0,1)
         t += dt;
 
         // determine next reaction