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
|
#include "solid_angle.h"
#include <math.h>
#include <stdio.h>
/* Table of the values of solid angle for various values of r0/r and L/r.
*
* See Table 1 in http://www.umich.edu/~ners312/CourseLibrary/SolidAngleOfADiskOffAxis.pdf. */
struct solid_angle_results {
double L;
double r0;
double omega;
} solid_angle_results[] = {
{0.5,0.0,3.4732594},
{0.5,0.2,3.4184435},
{0.5,0.4,3.2435434},
{0.5,0.6,2.9185178},
{0.5,0.8,2.4122535},
{0.5,1.0,1.7687239},
{0.5,1.2,1.1661307},
{0.5,1.4,0.7428889},
{0.5,1.6,0.4841273},
{0.5,1.8,0.3287007},
{0.5,2.0,0.2324189},
{1.0,0.0,1.8403024},
{1.0,0.2,1.8070933},
{1.0,0.4,1.7089486},
{1.0,0.6,1.5517370},
{1.0,0.8,1.3488367},
{1.0,1.0,1.1226876},
{1.0,1.2,0.9003572},
{1.0,1.4,0.7039130},
{1.0,1.6,0.5436956},
{1.0,1.8,0.4195415},
{1.0,2.0,0.3257993},
{1.5,0.0,1.0552591},
{1.5,0.2,1.0405177},
{1.5,0.4,0.9975504},
{1.5,0.6,0.9301028},
{1.5,0.8,0.8441578},
{1.5,1.0,0.7472299},
{1.5,1.2,0.6472056},
{1.5,1.4,0.5509617},
{1.5,1.6,0.4632819},
{1.5,1.8,0.3866757},
{1.5,2.0,0.3217142},
{2.0,0.0,0.6633335},
{2.0,0.2,0.6566352},
{2.0,0.4,0.6370508},
{2.0,0.6,0.6060694},
{2.0,0.8,0.5659755},
{2.0,1.0,0.5195359},
{2.0,1.2,0.4696858},
{2.0,1.4,0.4191714},
{2.0,1.6,0.3702014},
{2.0,1.8,0.3243908},
{2.0,2.0,0.282707}
};
int isclose(double a, double b, double rel_tol, double abs_tol)
{
/* Returns 1 if a and b are "close". This algorithm is taken from Python's
* math.isclose() function.
*
* See https://www.python.org/dev/peps/pep-0485/. */
return fabs(a-b) <= fmax(rel_tol*fmax(fabs(a),fabs(b)),abs_tol);
}
int test_solid_angle(char *err)
{
/* Tests the get_solid_angle() function. */
int i;
double pmt[3] = {0,0,0};
double pos[3] = {0,0,1};
double n[3] = {0,0,1};
double r = 1.0;
double solid_angle;
solid_angle = get_solid_angle(pos,pmt,n,r);
if (!isclose(solid_angle, 2*M_PI*(1-1/sqrt(2)), 1e-9, 0)) {
sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, 2*M_PI*(1-1/sqrt(2)));
return 1;
}
for (i = 0; i < sizeof(solid_angle_results)/sizeof(struct solid_angle_results); i++) {
pos[0] = solid_angle_results[i].r0*r;
pos[2] = solid_angle_results[i].L*r;
solid_angle = get_solid_angle(pos,pmt,n,r);
if (!isclose(solid_angle, solid_angle_results[i].omega, 1e-4, 0)) {
sprintf(err, "solid angle = %.5f, but expected %.5f", solid_angle, solid_angle_results[i].omega);
return 1;
}
}
return 0;
}
int main(int argc, char **argv)
{
char err[256];
if (!test_solid_angle(err)) {
printf("[\033[92mok\033[0m] test_solid_angle\n");
return 0;
} else {
printf("[\033[91mfail\033[0m] test_solid_angle: %s\n", err);
return 1;
}
}
|
