diff options
Diffstat (limited to 'solid_angle.c')
| -rw-r--r-- | solid_angle.c | 77 |
1 files changed, 76 insertions, 1 deletions
diff --git a/solid_angle.c b/solid_angle.c index d84ff22..4e44413 100644 --- a/solid_angle.c +++ b/solid_angle.c @@ -1,6 +1,81 @@ #include "solid_angle.h" #include <gsl/gsl_sf_ellint.h> #include <math.h> +#include <gsl/gsl_interp2d.h> +#include <gsl/gsl_spline2d.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}; + +static double lookupTable[13][11] = { +{0.9996,0.9996,0.9996,0.9996,0.9996,0.9996,0.9996,0.9996,0.9997,1.0000,1.0003}, +{1.0340,1.0334,1.0325,1.0316,1.0306,1.0257,1.0202,1.0150,1.0108,1.0059,1.0058}, +{1.0836,1.0814,1.0788,1.0758,1.0725,1.0578,1.0429,1.0300,1.0201,1.0091,1.0072}, +{1.1789,1.1718,1.1635,1.1546,1.1452,1.1067,1.0733,1.0480,1.0303,1.0122,1.0083}, +{1.3626,1.3375,1.3107,1.2837,1.2572,1.1669,1.1045,1.0642,1.0388,1.0147,1.0092}, +{1.7390,1.6319,1.5384,1.4590,1.3923,1.2167,1.1212,1.0737,1.0437,1.0163,1.0097}, +{2.2382,1.9036,1.6924,1.5489,1.4458,1.2241,1.1262,1.0742,1.0444,1.0170,1.0101}, +{1.7336,1.5812,1.4741,1.3945,1.3331,1.1850,1.1103,1.0676,1.0418,1.0170,1.0102}, +{1.2190,1.2218,1.2151,1.2032,1.1889,1.1314,1.0876,1.0575,1.0375,1.0165,1.0102}, +{1.0806,1.0920,1.0990,1.1023,1.1026,1.0880,1.0662,1.0471,1.0325,1.0157,1.0102}, +{1.0390,1.0464,1.0523,1.0566,1.0594,1.0591,1.0494,1.0380,1.0279,1.0148,1.0100}, +{1.0222,1.0270,1.0311,1.0345,1.0371,1.0408,1.0371,1.0305,1.0237,1.0138,1.0098}, +{1.0041,1.0051,1.0061,1.0070,1.0078,1.0106,1.0120,1.0122,1.0116,1.0097,1.0084} +}; + +static double A(double u, double a, double h) +{ + return atan(((a*a-h*h)*u - 2*a*a*h*h)/(2*a*h*sqrt((u-h*h)*(u+a*a)))); +} + +double get_solid_angle_approx(double *pos, double *pmt, double *n, double r) +{ + /* Returns the approximate solid angle subtended by a circular disk of + * radius r at a position `pmt` with a normal vector `n` from a position + * `pos`. + * + * This approximation comes from calculating the solid angle subtended by a + * square of equal area, which has a closed form solution. For certain + * 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. */ + double dir[3]; + double h, r0, D, a, u1, u2, F; + static gsl_spline2d *spline; + static gsl_interp_accel *xacc; + static gsl_interp_accel *yacc; + + dir[0] = pos[0] - pmt[0]; + dir[1] = pos[1] - pmt[1]; + dir[2] = pos[2] - pmt[2]; + + h = fabs(dir[0]*n[0] + dir[1]*n[1] + dir[2]*n[2]); + D = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]); + r0 = sqrt(D*D - h*h); + + a = sqrt(M_PI)*r/2; + + u1 = h*h + pow(r0-a,2); + u2 = h*h + pow(r0+a,2); + + F = 1.0; + if (h/D > 0.1 && h/D < 1.0 && r/D > 0.1 && r/D < 2.0) { + if (!spline) { + spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, sizeof(hd)/sizeof(double), sizeof(Rd)/sizeof(double)); + gsl_spline2d_init(spline, hd, Rd, (double *) lookupTable, sizeof(hd)/sizeof(double), sizeof(Rd)/sizeof(double)); + } + if (!xacc) xacc = gsl_interp_accel_alloc(); + if (!yacc) xacc = gsl_interp_accel_alloc(); + F = gsl_spline2d_eval(spline,h/D,r/D,xacc,yacc); + } + + if (r0 < a) { + return F*(A(u1,a,h) + A(u2,a,h) + M_PI); + } else { + return F*(A(u2,a,h) - A(u1,a,h)); + } +} double get_solid_angle(double *pos, double *pmt, double *n, double r) { @@ -23,7 +98,7 @@ 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 (fabs(r0-r) < 1e-9) { + if (fabs(r0-r) < 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); |