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
|
#include "solid_angle.h"
#include <gsl/gsl_sf_ellint.h>
#include <math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.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_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);
}
}
|
