aboutsummaryrefslogtreecommitdiff
path: root/scattering.c
diff options
context:
space:
mode:
authortlatorre <tlatorre@uchicago.edu>2018-08-14 10:08:27 -0500
committertlatorre <tlatorre@uchicago.edu>2018-08-14 10:08:27 -0500
commit24c8bcfe7f76b20124e2862ea050f815c0f768e7 (patch)
treee5bdbd638a2c7f38f1c094cc9e95cbdfe05b9481 /scattering.c
parent0b7f199c0d93074484ea580504485a32dc29f5e2 (diff)
downloadsddm-24c8bcfe7f76b20124e2862ea050f815c0f768e7.tar.gz
sddm-24c8bcfe7f76b20124e2862ea050f815c0f768e7.tar.bz2
sddm-24c8bcfe7f76b20124e2862ea050f815c0f768e7.zip
move everything to src directory
Diffstat (limited to 'scattering.c')
-rw-r--r--scattering.c121
1 files changed, 0 insertions, 121 deletions
diff --git a/scattering.c b/scattering.c
deleted file mode 100644
index a7fa717..0000000
--- a/scattering.c
+++ /dev/null
@@ -1,121 +0,0 @@
-#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);
-}