1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
#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_snoman_d2o(wavelength);
delta = (1.0/index - beta_cos_theta)/beta_sin_theta_theta0;
return qe*exp(-pow(delta,2)/2.0)/pow(wavelength,2)*1e7/sqrt(2*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)/(theta0*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);
}
|
