diff options
Diffstat (limited to 'src/solid_angle.c')
| -rw-r--r-- | src/solid_angle.c | 80 |
1 files changed, 79 insertions, 1 deletions
diff --git a/src/solid_angle.c b/src/solid_angle.c index 7201533..354547e 100644 --- a/src/solid_angle.c +++ b/src/solid_angle.c @@ -4,6 +4,7 @@ #include <gsl/gsl_interp2d.h> #include <gsl/gsl_spline2d.h> #include "vector.h" +#include "misc.h" static double hd[11] = {0.1,0.125,0.150,0.175,0.2,0.3,0.4,0.5,0.6,0.8,1.0}; static double Rd[13] = {0.1,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,2.0}; @@ -40,7 +41,8 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r) * regimes, the approximation is not very good and so it is modified by a * correction factor which comes from a lookup table. * - * This formula comes from "The Solid Angle Subtended at a Point by a * Circular Disk." Gardener et al. 1969. */ + * This formula comes from "The Solid Angle Subtended at a Point by a + * Circular Disk." Gardener et al. 1969. */ double dir[3]; double h, r0, D, a, u1, u2, F; static gsl_spline2d *spline; @@ -78,6 +80,82 @@ double get_solid_angle_approx(double *pos, double *pmt, double *n, double r) } } +double get_solid_angle2(double L2, double r2) +{ + double k = 2*sqrt(r2/(pow(L2,2)+pow(1+r2,2))); + double a2 = 4*r2/pow(1+r2,2); + double L_Rmax = L2/sqrt(pow(L2,2)+pow(1+r2,2)); + double r0_r = (r2-1)/(r2+1); + double K, P; + + if (fabs(r2-1) < 1e-5) { + /* 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); + return M_PI - 2*L_Rmax*K; + } else if (r2 <= 1) { + 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_Rmax*K + (2*L_Rmax)*r0_r*P; + } else { + K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE); + P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE); + return -2*L_Rmax*K + (2*L_Rmax)*r0_r*P; + } +} + +double get_solid_angle_lookup(double L2, double r2) +{ + size_t i, j; + static int initialized = 0; + static double xs[1000]; + static double ys[1000]; + static double zs[1000*1000]; + static double xlo = 0.0; + static double xhi = 1000; + static double ylo = 0.0; + static double yhi = 1000; + + if (!initialized) { + for (i = 0; i < LEN(xs); i++) { + xs[i] = xlo + (xhi-xlo)*i/(LEN(xs)-1); + } + for (i = 0; i < LEN(ys); i++) { + ys[i] = ylo + (yhi-ylo)*i/(LEN(ys)-1); + } + for (i = 0; i < LEN(xs); i++) { + for (j = 0; j < LEN(ys); j++) { + zs[i*LEN(ys) + j] = get_solid_angle2(xs[i],ys[j]); + } + } + initialized = 1; + } + + if (L2 < xlo || L2 > xhi || r2 < ylo || r2 > yhi) { + /* If the arguments are out of bounds of the lookup table, just call + * get_solid_angle2(). */ + return get_solid_angle2(L2,r2); + } + + return interp2d(L2,r2,xs,ys,zs,LEN(xs),LEN(ys)); +} + +double get_solid_angle_fast(double *pos, double *pmt, double *n, double r) +{ + double dir[3]; + double L, r0, R; + + dir[0] = pos[0] - pmt[0]; + dir[1] = pos[1] - pmt[1]; + dir[2] = pos[2] - pmt[2]; + + L = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]); + R = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); + r0 = sqrt(R*R - L*L); + + return get_solid_angle_lookup(L/r,r0/r); +} + double get_solid_angle(double *pos, double *pmt, double *n, double r) { /* Returns the solid angle subtended by a circular disk of radius r at a |