diff options
Diffstat (limited to 'src/likelihood.c')
| -rw-r--r-- | src/likelihood.c | 14 |
1 files changed, 14 insertions, 0 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 4a0e755..5ca9c8d 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -358,6 +358,13 @@ int get_delta_ray_weights(particle *p, double range, double distance_to_psup, do cdf = cdf1 + i*delta; cos_theta = gsl_spline_eval(p->spline_delta,cdf,p->acc_delta); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + + /* Occasionally due to floating point rounding error, it's possible + * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, + * we just make sure that xcdf is bounded by x1 and x2. */ + if (xcdf < x1) xcdf = x1; + if (xcdf > x2) xcdf = x2; + if (i > 0) { /* For each interval we approximate the integral as: * @@ -442,6 +449,13 @@ int get_shower_weights(particle *p, double distance_to_psup, double *x, double * cdf = cdf1 + i*delta; cos_theta = gsl_spline_eval(p->spline_shower,cdf,p->acc_shower); xcdf = cos_theta_to_x(cos_theta,R,cos_gamma,sin_gamma); + + /* Occasionally due to floating point rounding error, it's possible + * that cdf > cdf2 and so we end up with xcdf = inf. To prevent this, + * we just make sure that xcdf is bounded by x1 and x2. */ + if (xcdf < x1) xcdf = x1; + if (xcdf > x2) xcdf = x2; + if (i > 0) { /* For each interval we approximate the integral as: * |