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
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
#include "solid_angle.h"
#include <gsl/gsl_sf_ellint.h>
#include <math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include "vector.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}
};
double get_solid_angle_fast(double *pos, double *pmt, double *n, double r)
{
/* Returns a *very* fast approximation of the solid angle subtended by a
* circular disk of radius r at a position `pmt` with a normal vector `n` from
* a position `pos`. */
double dir[3];
double R, cos_theta_pmt;
SUB(dir,pos,pmt);
/* Distance to the PMT. */
R = NORM(dir);
cos_theta_pmt = fabs(DOT(dir,n)/R);
return pow(r,2)*cos_theta_pmt/pow(R,2);
}
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) yacc = 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);
}
}
|
