move everything to src directory

1 files changed, 115 insertions, 0 deletions

diff --git a/src/solid_angle.c b/src/solid_angle.c
new file mode 100644
index 0000000..4e44413
--- /dev/null
+++ b/src/solid_angle.c
@@ -0,0 +1,115 @@
+#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)
+{
+    /* Returns the solid angle subtended by a circular disk of radius r at a
+     * position `pmt` with a normal vector `n` from a position `pos`.
+     *
+     * See http://www.umich.edu/~ners312/CourseLibrary/SolidAngleOfADiskOffAxis.pdf. */
+    double dir[3];
+    double L, r0, R, a2, k, Rmax, K, P;
+
+    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);
+
+    a2 = 4*r0*r/pow(r0+r,2);
+    Rmax = sqrt(L*L + (r0+r)*(r0+r));
+    k = sqrt(4*r0*r/(L*L + pow(r0+r,2)));
+
+    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);
+        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);
+        return -2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r);
+    }
+}