move everything to src directory

1 files changed, 0 insertions, 115 deletions

diff --git a/solid_angle.c b/solid_angle.c
deleted file mode 100644
index 4e44413..0000000
--- a/solid_angle.c
+++ /dev/null
@@ -1,115 +0,0 @@
-#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);
-    }
-}