|
11 | 11 |
|
12 | 12 | namespace |
13 | 13 | { |
14 | | - //Shared implementation of the Mie coefficient (an, bn) and qSca/qExt/qBack |
15 | | - //computation. |
16 | | - template <typename RefIndexType> |
17 | | - void ComputeCoefficientsImpl(MieCoefficients &coeff, double xPara, RefIndexType relRef) |
| 14 | +//Shared implementation of the Mie coefficient (an, bn) and qSca/qExt/qBack |
| 15 | +//computation. |
| 16 | +template <typename RefIndexType> |
| 17 | +void ComputeCoefficientsImpl(MieCoefficients &coeff, double xPara, RefIndexType relRef) |
| 18 | +{ |
| 19 | + Utilities util; |
| 20 | + |
| 21 | + //use conventional symbols |
| 22 | + double x = xPara; |
| 23 | + RefIndexType m = relRef; |
| 24 | + RefIndexType mx = m * x; |
| 25 | + double xStop = x + 4.05 * (pow(x,(1.0 / 3.0))) + 2.0; |
| 26 | + |
| 27 | + unsigned int nStop = static_cast<unsigned int>(ceil(xStop)); |
| 28 | + unsigned int yMod = static_cast<unsigned int>(ceil(std::abs(mx))); |
| 29 | + unsigned int nMx = static_cast<unsigned int>(MAX(xStop, yMod) + 15); |
| 30 | + unsigned int arraySize = nStop + 1; |
| 31 | + double x2 = x * x; |
| 32 | + |
| 33 | + std::vector<RefIndexType> dnMx(nMx); |
| 34 | + dnMx[nMx-1] = 0; |
| 35 | + |
| 36 | + for (unsigned int n = nMx - 1; n>0; n--) |
18 | 37 | { |
19 | | - Utilities util; |
20 | | - |
21 | | - //use conventional symbols |
22 | | - double x = xPara; |
23 | | - RefIndexType m = relRef; |
24 | | - RefIndexType mx = m * x; |
25 | | - double xStop = x + 4.05 * (pow(x,(1.0 / 3.0))) + 2.0; |
26 | | - |
27 | | - unsigned int nStop = static_cast<unsigned int>(ceil(xStop)); |
28 | | - unsigned int yMod = static_cast<unsigned int>(ceil(std::abs(mx))); |
29 | | - unsigned int nMx = static_cast<unsigned int>(MAX(xStop, yMod) + 15); |
30 | | - unsigned int arraySize = nStop + 1; |
31 | | - double x2 = x * x; |
32 | | - |
33 | | - std::vector<RefIndexType> dnMx(nMx); |
34 | | - dnMx[nMx-1] = 0; |
35 | | - |
36 | | - for (unsigned int n = nMx - 1; n>0; n--) |
37 | | - { |
38 | | - dnMx[n-1] = (double(n)/mx)-(1.0/(dnMx[n]+double(n)/mx)); |
39 | | - } |
40 | | - |
41 | | - // at the sphere boundary |
42 | | - double jX0 = cos(x); // phi(-1) |
43 | | - double yX0 = -sin(x); // kai(-1) |
44 | | - double jX1 = sin(x); // phi(0) |
45 | | - double yX1 = cos(x); // kai(0) |
46 | | - |
47 | | - std::vector<double> jX(arraySize); |
48 | | - std::vector<double> yX(arraySize); |
49 | | - std::vector<std::complex<double>> xi_x(arraySize); |
50 | | - jX[0] = jX1; |
51 | | - yX[0] = yX1; |
52 | | - xi_x[0]=std::complex<double> (jX1,-yX1); // xi(1) |
53 | | - |
54 | | - //Initialize temp holders |
55 | | - std::complex<double> tempQback = 0.0; |
56 | | - double tempQsca = 0.0; |
57 | | - double tempQext = 0.0; |
58 | | - double sign = -1.0; // tracks (-1)^(n-1), toggled each iteration instead of calling pow() |
59 | | - |
60 | | - coeff.an.assign(arraySize, std::complex<double>(0.0, 0.0)); |
61 | | - coeff.bn.assign(arraySize, std::complex<double>(0.0, 0.0)); |
62 | | - |
63 | | - for (unsigned int n = 1; n <= nStop; n++) |
64 | | - { |
65 | | - double fac2 = 2.0 * double(n) + 1.0; // 2n+1 |
66 | | - double fac3 = fac2 - 2.0; // 2n-1 |
67 | | - |
68 | | - //Update riccati Bessel functions for x (Array indices = n+1) |
69 | | - jX[n] = fac3 * jX1 / x - jX0; // phi recurrence |
70 | | - yX[n] = fac3 * yX1 / x - yX0; // kai recurrence |
71 | | - xi_x[n] = std::complex<double> (jX[n], -yX[n]); |
72 | | - |
73 | | - jX0 = jX1; |
74 | | - jX1 = jX[n]; |
75 | | - yX0 = yX1; |
76 | | - yX1 = yX[n]; |
77 | | - |
78 | | - // Calculate an and bn (According to Bohren and Huffman book) |
79 | | - // Remark: GouGouesbet uses size parameter as "ka" instead of "kx" |
80 | | - RefIndexType dervDn1 = (dnMx[n]/m) + (double(n) / x); |
81 | | - RefIndexType dervDn2 = (m*dnMx[n]) + (double(n) / x); |
82 | | - |
83 | | - std::complex<double> an = (dervDn1*jX[n] - jX[n-1])/ (dervDn1*xi_x[n] - xi_x[n-1]); |
84 | | - std::complex<double> bn = (dervDn2*jX[n] - jX[n-1])/ (dervDn2*xi_x[n] - xi_x[n-1]); |
85 | | - coeff.an[n-1] = an; |
86 | | - coeff.bn[n-1] = bn; |
87 | | - |
88 | | - sign = -sign; // (-1)^(n-1) |
89 | | - tempQback += fac2 * sign * (an-bn); |
90 | | - tempQsca += fac2 * (util.ComplexAbs(an) * util.ComplexAbs(an) + util.ComplexAbs(bn) * util.ComplexAbs(bn)); |
91 | | - tempQext += fac2 * (an + bn).real(); |
92 | | - } |
93 | | - coeff.nStop = nStop; |
94 | | - coeff.qBack = util.ComplexAbsSquared(tempQback)/x2; //back scattering efficiency |
95 | | - coeff.qSca = 2.0 * tempQsca / x2; //scattering efficiency |
96 | | - coeff.qExt = 2.0 * tempQext / x2; //extinction efficiency |
| 38 | + dnMx[n-1] = (double(n)/mx)-(1.0/(dnMx[n]+double(n)/mx)); |
97 | 39 | } |
| 40 | + |
| 41 | + // at the sphere boundary |
| 42 | + double jX0 = cos(x); // phi(-1) |
| 43 | + double yX0 = -sin(x); // kai(-1) |
| 44 | + double jX1 = sin(x); // phi(0) |
| 45 | + double yX1 = cos(x); // kai(0) |
| 46 | + |
| 47 | + std::vector<double> jX(arraySize); |
| 48 | + std::vector<double> yX(arraySize); |
| 49 | + std::vector<std::complex<double>> xi_x(arraySize); |
| 50 | + jX[0] = jX1; |
| 51 | + yX[0] = yX1; |
| 52 | + xi_x[0]=std::complex<double> (jX1,-yX1); // xi(1) |
| 53 | + |
| 54 | + //Initialize temp holders |
| 55 | + std::complex<double> tempQback = 0.0; |
| 56 | + double tempQsca = 0.0; |
| 57 | + double tempQext = 0.0; |
| 58 | + double sign = -1.0; // tracks (-1)^(n-1), toggled each iteration instead of calling pow() |
| 59 | + |
| 60 | + coeff.an.assign(arraySize, std::complex<double>(0.0, 0.0)); |
| 61 | + coeff.bn.assign(arraySize, std::complex<double>(0.0, 0.0)); |
| 62 | + |
| 63 | + for (unsigned int n = 1; n <= nStop; n++) |
| 64 | + { |
| 65 | + double fac2 = 2.0 * double(n) + 1.0; // 2n+1 |
| 66 | + double fac3 = fac2 - 2.0; // 2n-1 |
| 67 | + |
| 68 | + //Update riccati Bessel functions for x (Array indices = n+1) |
| 69 | + jX[n] = fac3 * jX1 / x - jX0; // phi recurrence |
| 70 | + yX[n] = fac3 * yX1 / x - yX0; // kai recurrence |
| 71 | + xi_x[n] = std::complex<double> (jX[n], -yX[n]); |
| 72 | + |
| 73 | + jX0 = jX1; |
| 74 | + jX1 = jX[n]; |
| 75 | + yX0 = yX1; |
| 76 | + yX1 = yX[n]; |
| 77 | + |
| 78 | + // Calculate an and bn (According to Bohren and Huffman book) |
| 79 | + // Remark: GouGouesbet uses size parameter as "ka" instead of "kx" |
| 80 | + RefIndexType dervDn1 = (dnMx[n]/m) + (double(n) / x); |
| 81 | + RefIndexType dervDn2 = (m*dnMx[n]) + (double(n) / x); |
| 82 | + |
| 83 | + std::complex<double> an = (dervDn1*jX[n] - jX[n-1])/ (dervDn1*xi_x[n] - xi_x[n-1]); |
| 84 | + std::complex<double> bn = (dervDn2*jX[n] - jX[n-1])/ (dervDn2*xi_x[n] - xi_x[n-1]); |
| 85 | + coeff.an[n-1] = an; |
| 86 | + coeff.bn[n-1] = bn; |
| 87 | + |
| 88 | + sign = -sign; // (-1)^(n-1) |
| 89 | + tempQback += fac2 * sign * (an-bn); |
| 90 | + tempQsca += fac2 * (util.ComplexAbs(an) * util.ComplexAbs(an) + util.ComplexAbs(bn) * util.ComplexAbs(bn)); |
| 91 | + tempQext += fac2 * (an + bn).real(); |
| 92 | + } |
| 93 | + coeff.nStop = nStop; |
| 94 | + coeff.qBack = util.ComplexAbsSquared(tempQback)/x2; //back scattering efficiency |
| 95 | + coeff.qSca = 2.0 * tempQsca / x2; //scattering efficiency |
| 96 | + coeff.qExt = 2.0 * tempQext / x2; //extinction efficiency |
| 97 | +} |
98 | 98 | } |
99 | 99 |
|
100 | 100 | //Computes the Mie coefficients (an, bn) for a real relative refractive index. |
|
0 commit comments