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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
#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"
#include "misc.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}
};
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_angle2(double L2, double r2)
{
double k = 2*sqrt(r2/(pow(L2,2)+pow(1+r2,2)));
double a2 = 4*r2/pow(1+r2,2);
double L_Rmax = L2/sqrt(pow(L2,2)+pow(1+r2,2));
double r0_r = (r2-1)/(r2+1);
double K, P;
if (fabs(r2-1) < 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_Rmax*K;
} else if (r2 <= 1) {
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_Rmax*K + (2*L_Rmax)*r0_r*P;
} else {
K = gsl_sf_ellint_Kcomp(k, GSL_PREC_DOUBLE);
P = gsl_sf_ellint_Pcomp(k, -a2, GSL_PREC_DOUBLE);
return -2*L_Rmax*K + (2*L_Rmax)*r0_r*P;
}
}
double get_solid_angle_lookup(double L2, double r2)
{
size_t i, j;
static int initialized = 0;
static double xs[1000];
static double ys[1000];
static double zs[1000*1000];
static double xlo = 0.0;
static double xhi = 1000;
static double ylo = 0.0;
static double yhi = 1000;
if (!initialized) {
for (i = 0; i < LEN(xs); i++) {
xs[i] = xlo + (xhi-xlo)*i/(LEN(xs)-1);
}
for (i = 0; i < LEN(ys); i++) {
ys[i] = ylo + (yhi-ylo)*i/(LEN(ys)-1);
}
for (i = 0; i < LEN(xs); i++) {
for (j = 0; j < LEN(ys); j++) {
zs[i*LEN(ys) + j] = get_solid_angle2(xs[i],ys[j]);
}
}
initialized = 1;
}
if (L2 < xlo || L2 > xhi || r2 < ylo || r2 > yhi) {
/* If the arguments are out of bounds of the lookup table, just call
* get_solid_angle2(). */
return get_solid_angle2(L2,r2);
}
return interp2d(L2,r2,xs,ys,zs,LEN(xs),LEN(ys));
}
double get_solid_angle_fast(double *pos, double *pmt, double *n, double r)
{
double dir[3];
double L, r0, R;
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);
return get_solid_angle_lookup(L/r,r0/r);
}
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);
}
}
|
