initial commit of a function to calculate the solid angle subtended by a circular disk

4 files changed, 41 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
index d17415b..9e63049 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,3 @@
 *.swp
 calculate_limits
+*.o
diff --git a/Makefile b/Makefile
index 12de948..6148848 100644
--- a/Makefile
+++ b/Makefile
@@ -3,6 +3,8 @@ LDLIBS=-lm -lgsl -lgslcblas
 
 calculate_limits: calculate_limits.c
 
+solid_angle.o: solid_angle.c
+
 clean:
 	rm calculate_limits
 
diff --git a/solid_angle.c b/solid_angle.c
new file mode 100644
index 0000000..5fc6033
--- /dev/null
+++ b/solid_angle.c
@@ -0,0 +1,35 @@
+#include "solid_angle.h"
+#include <gsl/gsl_sf_ellint.h>
+#include <math.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, a, k, Rmax, K, P;
+
+    dir[0] = pos[0] - pmt[0];
+    dir[1] = pos[1] - pmt[1];
+    dir[2] = pos[2] - pmt[2];
+
+    L = 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);
+
+    a = 4*r0*r/pow(r0+r,2);
+    Rmax = L*L + (r0+r)*(r0+r);
+    k = sqrt(4*r0*r/(L*L + pow(r0+r,2)));
+
+    if (r0 <= r) {
+        K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
+        P = gsl_sf_ellint_Pcomp(k, a, 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, a, GSL_PREC_DOUBLE);
+        return -2*L*K/Rmax + (2*L/Rmax)*(r0-r)*P/(r0+r);
+    }
+}
diff --git a/solid_angle.h b/solid_angle.h
new file mode 100644
index 0000000..99fda37
--- /dev/null
+++ b/solid_angle.h
@@ -0,0 +1,3 @@
+#include <math.h>
+
+double get_solid_angle(double *pos, double *pmt, double *n, double r);