diff options
| author | tlatorre <tlatorre@uchicago.edu> | 2018-07-04 14:31:33 -0400 |
|---|---|---|
| committer | tlatorre <tlatorre@uchicago.edu> | 2018-07-04 14:31:33 -0400 |
| commit | ac9bc989fea328db6db1820d2676cbbf9335584f (patch) | |
| tree | 2b11eec051b9ca4f1e0922400735d54cb0d793fe /solid_angle.c | |
| parent | e61f3f60d73b2c0618f1c87533bfc756ff217ee5 (diff) | |
| download | sddm-ac9bc989fea328db6db1820d2676cbbf9335584f.tar.gz sddm-ac9bc989fea328db6db1820d2676cbbf9335584f.tar.bz2 sddm-ac9bc989fea328db6db1820d2676cbbf9335584f.zip | |
update solid angle test and fix a bug in the solid angle calculation
Diffstat (limited to 'solid_angle.c')
| -rw-r--r-- | solid_angle.c | 11 |
1 files changed, 8 insertions, 3 deletions
diff --git a/solid_angle.c b/solid_angle.c index f89d65e..d84ff22 100644 --- a/solid_angle.c +++ b/solid_angle.c @@ -23,13 +23,18 @@ double get_solid_angle(double *pos, double *pmt, double *n, double r) Rmax = sqrt(L*L + (r0+r)*(r0+r)); k = sqrt(4*r0*r/(L*L + pow(r0+r,2))); - if (r0 <= r) { + if (fabs(r0-r) < 1e-9) { + /* If r0 is very close to r, the GSL routines below will occasionally + * throw a domain error. */ K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); - P = gsl_sf_ellint_Pcomp(k, a2, GSL_PREC_DOUBLE); + return M_PI - 2*L*K/Rmax; + } else if (r0 <= r) { + K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); + P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE); return 2*M_PI - 2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r); } else { K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); - P = gsl_sf_ellint_Pcomp(k, a2, GSL_PREC_DOUBLE); + P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE); return -2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r); } } |