move everything to src directory

1 files changed, 121 insertions, 0 deletions

diff --git a/src/scattering.c b/src/scattering.c
new file mode 100644
index 0000000..a7fa717
--- /dev/null
+++ b/src/scattering.c
@@ -0,0 +1,121 @@
+#include "scattering.h"
+#include "quantum_efficiency.h"
+#include <gsl/gsl_sf_bessel.h>
+#include <gsl/gsl_integration.h>
+#include <gsl/gsl_interp2d.h>
+#include <gsl/gsl_spline2d.h>
+#include "optics.h"
+#include "sno.h"
+#include <stddef.h> /* for size_t */
+#include <stdlib.h> /* for exit() */
+#include <stdio.h> /* for fprintf() */
+#include <gsl/gsl_errno.h> /* for gsl_strerror() */
+
+static double xlo = -1.0;
+static double xhi = 1.0;
+static size_t nx = 200;
+static double ylo = 0.0;
+static double yhi = 1.0;
+static size_t ny = 1000;
+
+static double *x;
+static double *y;
+static double *z;
+
+static gsl_spline2d *spline;
+static gsl_interp_accel *xacc;
+static gsl_interp_accel *yacc;
+
+static double prob_scatter(double wavelength, void *params)
+{
+    /* Calculate the number of photons emitted per unit solid angle per cm at
+     * an angle `theta` for a particle travelling with velocity `beta` with an
+     * angular distribution of width `theta0`. */
+    double qe, delta, index;
+
+    double beta_cos_theta = ((double *) params)[0];
+    double beta_sin_theta_theta0 = ((double *) params)[1];
+
+    qe = get_quantum_efficiency(wavelength);
+
+    index = get_index(HEAVY_WATER_DENSITY, wavelength, 10.0);
+
+    delta = (1.0/index - beta_cos_theta)/(2*beta_sin_theta_theta0);
+
+    /* FIXME: ignore GSL error for underflow here. */
+    if (delta*delta > 500) return 0.0;
+    else if (delta*delta == 0.0) return INFINITY;
+    return qe*exp(-delta*delta)*gsl_sf_bessel_K0(delta*delta)/pow(wavelength,2)*1e7/(4*M_PI*M_PI);
+}
+
+void init_interpolation(void)
+{
+    size_t i, j;
+    double params[2];
+    double result, error;
+    size_t nevals;
+    int status;
+
+    x = malloc(nx*sizeof(double));
+    y = malloc(ny*sizeof(double));
+    z = malloc(nx*ny*sizeof(double));
+
+    spline = gsl_spline2d_alloc(gsl_interp2d_bilinear, nx, ny);
+    xacc = gsl_interp_accel_alloc();
+    yacc = gsl_interp_accel_alloc();
+
+    for (i = 0; i < nx; i++) {
+        x[i] = xlo + (xhi-xlo)*i/(nx-1);
+    }
+
+    for (i = 0; i < ny; i++) {
+        y[i] = ylo + (yhi-ylo)*i/(ny-1);
+    }
+
+    gsl_integration_cquad_workspace *w = gsl_integration_cquad_workspace_alloc(100);
+
+    gsl_function F;
+    F.function = &prob_scatter;
+    F.params = params;
+
+    for (i = 0; i < nx; i++) {
+        for (j = 0; j < ny; j++) {
+            params[0] = x[i];
+            params[1] = y[j];
+
+            status = gsl_integration_cquad(&F, 200, 800, 0, 1e-2, w, &result, &error, &nevals);
+
+            if (status) {
+                fprintf(stderr, "error integrating photon angular distribution: %s\n", gsl_strerror(status));
+                exit(1);
+            }
+            gsl_spline2d_set(spline, z, i, j, result);
+        }
+    }
+
+    gsl_spline2d_init(spline, x, y, z, nx, ny);
+
+    gsl_integration_cquad_workspace_free(w);
+}
+
+double get_probability(double beta, double cos_theta, double theta0)
+{
+    double sin_theta;
+
+    /* Technically this isn't defined up to a sign, but it doesn't matter since
+     * we are going to square it everywhere. */
+    sin_theta = fabs(sin(acos(cos_theta)));
+
+    return gsl_spline2d_eval(spline, beta*cos_theta, beta*sin_theta*theta0, xacc, yacc)/sin_theta;
+}
+
+void free_interpolation(void)
+{
+    free(x);
+    free(y);
+    free(z);
+
+    gsl_spline2d_free(spline);
+    gsl_interp_accel_free(xacc);
+    gsl_interp_accel_free(yacc);
+}