/* Copyright (c) 2019, Anthony Latorre * * This program is free software: you can redistribute it and/or modify it * under the terms of the GNU General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) * any later version. * This program is distributed in the hope that it will be useful, but WITHOUT * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for * more details. * You should have received a copy of the GNU General Public License along with * this program. If not, see . */ #include "solid_angle.h" #include #include #include #include #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, LEN(hd), LEN(Rd)); gsl_spline2d_init(spline, hd, Rd, (double *) lookupTable, LEN(hd), LEN(Rd)); } 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, R2; dir[0] = pos[0] - pmt[0]; dir[1] = pos[1] - pmt[1]; dir[2] = pos[2] - pmt[2]; L = fabs(DOT(dir,n)); R2 = DOT(dir,dir); r0 = sqrt(R2 - 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); } }