physics/scattering.cpp
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // SPDX-License-Identifier: BSD-3-Clause | ||
| 2 | /// @file scattering.cpp | ||
| 3 | /// @brief Rayleigh and Mie scattering, refractivities and ozone. | ||
| 4 | #include "ernie/physics/scattering.h" | ||
| 5 | |||
| 6 | #include "ernie/util/math.h" | ||
| 7 | |||
| 8 | #include <algorithm> | ||
| 9 | #include <array> | ||
| 10 | #include <cmath> | ||
| 11 | #include <numbers> | ||
| 12 | #include <stdexcept> | ||
| 13 | |||
| 14 | namespace ernie { | ||
| 15 | |||
| 16 | namespace { | ||
| 17 | |||
| 18 | using std::numbers::pi; | ||
| 19 | |||
| 20 | /// Inverse square of a wavelength in nm, in um^-2, the variable of | ||
| 21 | /// the dispersion formulas. | ||
| 22 | 297 | double InverseSquareMicrons(double wavelengthNm) { | |
| 23 | 297 | const double um = wavelengthNm * 1e-3; | |
| 24 | 297 | return 1.0 / (um * um); | |
| 25 | } | ||
| 26 | |||
| 27 | /// Wiscombe's (1980) number of Mie series terms for size parameter | ||
| 28 | /// `x`. | ||
| 29 | 25068 | int MieTerms(double x) { | |
| 30 | 25068 | const double cubeRoot = std::cbrt(x); | |
| 31 | 25068 | const double terms = x <= 8.0 ? x + 4.0 * cubeRoot + 1.0 | |
| 32 | 9539 | : x < 4200.0 ? x + 4.05 * cubeRoot + 2.0 | |
| 33 | ✗ | : x + 4.0 * cubeRoot + 2.0; | |
| 34 | 25068 | return std::max(1, static_cast<int>(terms)); | |
| 35 | } | ||
| 36 | |||
| 37 | /// The logarithmic derivative D_n(z) = psi_n'(z) / psi_n(z) of the | ||
| 38 | /// Riccati-Bessel function, by Lentz's continued fraction for the | ||
| 39 | /// ratio of spherical Bessel functions: D_n = -n / z + a_1 a_2 ..., | ||
| 40 | /// with partial terms a_j = (-1)^(j+1) (2n + 2j - 1) / z and | ||
| 41 | /// convergents taken as ratios of two running products. | ||
| 42 | 50136 | std::complex<double> LentzLogDerivative(std::complex<double> z, int n) { | |
| 43 | 50136 | using Complex = std::complex<double>; | |
| 44 | 50136 | Complex inverse = 2.0 / z; | |
| 45 | 50136 | const Complex alpha = (n + 0.5) * inverse; | |
| 46 | 50136 | Complex term = -(n + 1.5) * inverse; | |
| 47 | 50136 | Complex numerator = term + 1.0 / alpha; | |
| 48 | 50136 | Complex denominator = term; | |
| 49 | 50136 | Complex ratio = numerator / denominator; | |
| 50 | 50136 | Complex product = alpha * ratio; | |
| 51 | 635142 | while (std::abs(std::abs(ratio) - 1.0) > 1e-12) { | |
| 52 | 585006 | term = inverse - term; | |
| 53 | 585006 | numerator = 1.0 / numerator + term; | |
| 54 | 585006 | denominator = 1.0 / denominator + term; | |
| 55 | 585006 | ratio = numerator / denominator; | |
| 56 | 585006 | inverse = -inverse; | |
| 57 | 585006 | product *= ratio; | |
| 58 | } | ||
| 59 | 50136 | return -static_cast<double>(n) / z + product; | |
| 60 | } | ||
| 61 | |||
| 62 | } // namespace | ||
| 63 | |||
| 64 | 170 | double AirRefractivity(double wavelengthNm, double co2Fraction) { | |
| 65 | 170 | const double s = InverseSquareMicrons(wavelengthNm); | |
| 66 | 170 | const double n300 = (8060.51 + 2480990.0 / (132.274 - s) + 17455.7 / (39.32957 - s)) * 1e-8; | |
| 67 | 170 | return n300 * (1.0 + 0.54 * (co2Fraction - 0.0003)); | |
| 68 | } | ||
| 69 | |||
| 70 | 8 | double CarbonDioxideRefractivity(double wavelengthNm) { | |
| 71 | 8 | const double s = InverseSquareMicrons(wavelengthNm); | |
| 72 | 8 | return (22822.1 + 117.8 * s + 2406030.0 / (130.0 - s) + 15997.0 / (38.9 - s)) * 1e-8; | |
| 73 | } | ||
| 74 | |||
| 75 | 119 | double AirKingFactor(double wavelengthNm, double co2Fraction) { | |
| 76 | 119 | const double s = InverseSquareMicrons(wavelengthNm); | |
| 77 | 119 | const double nitrogen = 1.034 + 3.17e-4 * s; | |
| 78 | 119 | const double oxygen = 1.096 + 1.385e-3 * s + 1.448e-4 * s * s; | |
| 79 | // Volume percentages of N2, O2 and Ar in dry air, and CO2's. | ||
| 80 | 119 | const double co2 = co2Fraction * 100.0; | |
| 81 | 119 | return (78.084 * nitrogen + 20.946 * oxygen + 0.934 * 1.0 + co2 * kCarbonDioxideKingFactor) / | |
| 82 | 119 | (78.084 + 20.946 + 0.934 + co2); | |
| 83 | } | ||
| 84 | |||
| 85 | 122 | double RayleighCrossSection(double wavelengthNm, double refractivity, double kingFactor) { | |
| 86 | 122 | const double lambda = wavelengthNm * 1e-9; | |
| 87 | 122 | const double n2 = (1.0 + refractivity) * (1.0 + refractivity); | |
| 88 | 122 | const double ratio = (n2 - 1.0) / (n2 + 2.0); | |
| 89 | 122 | const double lambda2 = lambda * lambda; | |
| 90 | 122 | return 24.0 * pi * pi * pi * ratio * ratio / | |
| 91 | 122 | (lambda2 * lambda2 * kStandardNumberDensity * kStandardNumberDensity) * kingFactor; | |
| 92 | } | ||
| 93 | |||
| 94 | 72398 | double RayleighPhase(double cosTheta, double kingFactor) { | |
| 95 | // F = (6 + 3 rho) / (6 - 7 rho) solved for the depolarisation | ||
| 96 | // ratio rho. | ||
| 97 | 72398 | const double rho = 6.0 * (kingFactor - 1.0) / (3.0 + 7.0 * kingFactor); | |
| 98 | 72398 | const double gamma = rho / (2.0 - rho); | |
| 99 | 72398 | return 3.0 / (16.0 * pi) * ((1.0 + 3.0 * gamma) + (1.0 - gamma) * cosTheta * cosTheta) / | |
| 100 | 72398 | (1.0 + 2.0 * gamma); | |
| 101 | } | ||
| 102 | |||
| 103 | 25068 | MieEfficiencies MieSphere(double sizeParameter, std::complex<double> refractiveIndex, | |
| 104 | std::span<const double> cosines, std::span<double> phase) { | ||
| 105 | 25068 | using Complex = std::complex<double>; | |
| 106 | 25068 | if (!phase.empty() && phase.size() != cosines.size()) { | |
| 107 | ✗ | throw std::invalid_argument("MieSphere needs one phase value per cosine"); | |
| 108 | } | ||
| 109 | 25068 | const double x = sizeParameter; | |
| 110 | 25068 | const Complex m = refractiveIndex; | |
| 111 | 25068 | const int terms = MieTerms(x); | |
| 112 | |||
| 113 | // The logarithmic derivative D_n(mx) by downward recurrence | ||
| 114 | // (Bohren and Huffman, eq. 4.89), stable for absorbing spheres | ||
| 115 | // where upward recurrence is not, started from its exact value | ||
| 116 | // by Lentz's continued fraction: Bohren and Huffman's start at | ||
| 117 | // 0 a few orders up loses digits for large |mx|. | ||
| 118 | 25068 | const Complex y = m * x; | |
| 119 | 25068 | const std::size_t count = static_cast<std::size_t>(terms); | |
| 120 | 25068 | std::vector<Complex> d(count + 1); | |
| 121 | 25068 | std::vector<double> dx(count + 1); | |
| 122 | 25068 | d[count] = LentzLogDerivative(y, terms); | |
| 123 | 25068 | dx[count] = LentzLogDerivative(Complex(x), terms).real(); | |
| 124 | 816204 | for (std::size_t n = count; n > 0; --n) { | |
| 125 | 791136 | const Complex nOverY = static_cast<double>(n) / y; | |
| 126 | 791136 | d[n - 1] = nOverY - 1.0 / (d[n] + nOverY); | |
| 127 | 791136 | const double nOverX = static_cast<double>(n) / x; | |
| 128 | 791136 | dx[n - 1] = nOverX - 1.0 / (dx[n] + nOverX); | |
| 129 | } | ||
| 130 | |||
| 131 | // Riccati-Bessel functions psi_n = x j_n(x) and chi_n = -x | ||
| 132 | // y_n(x) from n = -1, 0. chi grows with n and recurs upward | ||
| 133 | // stably, as does psi up to n = x; past that psi decays and | ||
| 134 | // upward recurrence loses digits, so it follows from D_n(x) | ||
| 135 | // instead: psi_n = psi_{n-1} / (D_n(x) + n / x). | ||
| 136 | 25068 | double psiPrev = std::cos(x); | |
| 137 | 25068 | double psi = std::sin(x); | |
| 138 | 25068 | double chiPrev = -std::sin(x); | |
| 139 | 25068 | double chi = std::cos(x); | |
| 140 | 25068 | Complex xiPrev(psi, -chi); | |
| 141 | |||
| 142 | // Angular functions pi_n and the amplitude functions S1, S2 per | ||
| 143 | // cosine. | ||
| 144 | 25068 | const std::size_t angles = phase.empty() ? 0 : cosines.size(); | |
| 145 | 25068 | std::vector<double> piPrev(angles, 0.0); | |
| 146 | 25068 | std::vector<double> piCurrent(angles, 1.0); | |
| 147 | 25068 | std::vector<Complex> s1(angles, Complex(0.0)); | |
| 148 | 25068 | std::vector<Complex> s2(angles, Complex(0.0)); | |
| 149 | |||
| 150 | 25068 | double extinction = 0.0; | |
| 151 | 25068 | double scattering = 0.0; | |
| 152 | 25068 | double asymmetry = 0.0; | |
| 153 | 25068 | Complex aPrev(0.0); | |
| 154 | 25068 | Complex bPrev(0.0); | |
| 155 | 816204 | for (std::size_t n = 1; n <= count; ++n) { | |
| 156 | 791136 | const double fn = static_cast<double>(n); | |
| 157 | 791136 | const double psiNext = | |
| 158 | 791136 | fn <= x ? (2.0 * fn - 1.0) / x * psi - psiPrev : psi / (dx[n] + fn / x); | |
| 159 | 791136 | const double chiNext = (2.0 * fn - 1.0) / x * chi - chiPrev; | |
| 160 | 791136 | psiPrev = psi; | |
| 161 | 791136 | psi = psiNext; | |
| 162 | 791136 | chiPrev = chi; | |
| 163 | 791136 | chi = chiNext; | |
| 164 | 791136 | const Complex xi(psi, -chi); | |
| 165 | |||
| 166 | // Bohren and Huffman, eq. 4.88. | ||
| 167 | 791136 | const Complex da = d[n] / m + fn / x; | |
| 168 | 791136 | const Complex db = d[n] * m + fn / x; | |
| 169 | 791136 | const Complex a = (da * psi - psiPrev) / (da * xi - xiPrev); | |
| 170 | 791136 | const Complex b = (db * psi - psiPrev) / (db * xi - xiPrev); | |
| 171 | |||
| 172 | 791136 | extinction += (2.0 * fn + 1.0) * (a.real() + b.real()); | |
| 173 | 791136 | scattering += (2.0 * fn + 1.0) * (std::norm(a) + std::norm(b)); | |
| 174 | 791136 | asymmetry += (2.0 * fn + 1.0) / (fn * (fn + 1.0)) * (a * std::conj(b)).real(); | |
| 175 | 791136 | if (n > 1) { | |
| 176 | 766068 | asymmetry += | |
| 177 | 766068 | (fn - 1.0) * (fn + 1.0) / fn * (aPrev * std::conj(a) + bPrev * std::conj(b)).real(); | |
| 178 | } | ||
| 179 | |||
| 180 | 791136 | const double weight = (2.0 * fn + 1.0) / (fn * (fn + 1.0)); | |
| 181 | 189188832 | for (std::size_t j = 0; j < angles; ++j) { | |
| 182 | 188397696 | const double mu = cosines[j]; | |
| 183 | 188397696 | const double tau = fn * mu * piCurrent[j] - (fn + 1.0) * piPrev[j]; | |
| 184 | 188397696 | s1[j] += weight * (a * piCurrent[j] + b * tau); | |
| 185 | 188397696 | s2[j] += weight * (a * tau + b * piCurrent[j]); | |
| 186 | 188397696 | const double piNext = | |
| 187 | 188397696 | ((2.0 * fn + 1.0) * mu * piCurrent[j] - (fn + 1.0) * piPrev[j]) / fn; | |
| 188 | 188397696 | piPrev[j] = piCurrent[j]; | |
| 189 | 188397696 | piCurrent[j] = piNext; | |
| 190 | } | ||
| 191 | |||
| 192 | 791136 | xiPrev = xi; | |
| 193 | 791136 | aPrev = a; | |
| 194 | 791136 | bPrev = b; | |
| 195 | } | ||
| 196 | |||
| 197 | 25068 | MieEfficiencies result; | |
| 198 | 25068 | result.extinction = 2.0 / (x * x) * extinction; | |
| 199 | 25068 | result.scattering = 2.0 / (x * x) * scattering; | |
| 200 | 50136 | result.asymmetry = | |
| 201 | 25068 | result.scattering > 0.0 ? 4.0 / (x * x * result.scattering) * asymmetry : 0.0; | |
| 202 | // The unpolarised intensity (|S1|^2 + |S2|^2) / 2 integrates to | ||
| 203 | // pi x^2 Q_sca over the sphere. | ||
| 204 | 25068 | const double normalisation = 2.0 * pi * x * x * result.scattering; | |
| 205 | 3370793 | for (std::size_t j = 0; j < angles; ++j) { | |
| 206 | 3345725 | phase[j] = normalisation > 0.0 ? (std::norm(s1[j]) + std::norm(s2[j])) / normalisation | |
| 207 | : kIsotropicPhase; | ||
| 208 | } | ||
| 209 | 25068 | return result; | |
| 210 | 25068 | } | |
| 211 | |||
| 212 | 8 | LogNormalMode LogNormalMode::FromEffective(double effectiveRadius, double effectiveVariance, | |
| 213 | double relativeNumber) { | ||
| 214 | 8 | const double variance = std::log1p(effectiveVariance); | |
| 215 | 8 | return {effectiveRadius * std::exp(-2.5 * variance), std::sqrt(variance), relativeNumber}; | |
| 216 | } | ||
| 217 | |||
| 218 | 1 | double LogNormalMode::EffectiveRadius() const { | |
| 219 | 1 | return medianRadius * std::exp(2.5 * sigma * sigma); | |
| 220 | } | ||
| 221 | |||
| 222 | 1 | double LogNormalMode::EffectiveVariance() const { | |
| 223 | 1 | return std::expm1(sigma * sigma); | |
| 224 | } | ||
| 225 | |||
| 226 | 3 | RefractiveIndex::RefractiveIndex(std::complex<double> value) : wavelengths_{1.0}, values_{value} {} | |
| 227 | |||
| 228 | 7 | RefractiveIndex::RefractiveIndex(std::vector<double> wavelengthsNm, | |
| 229 | 7 | std::vector<std::complex<double>> values) | |
| 230 | 7 | : wavelengths_(std::move(wavelengthsNm)), values_(std::move(values)) { | |
| 231 | 7 | if (wavelengths_.size() != values_.size()) { | |
| 232 | ✗ | throw std::invalid_argument("A refractive index needs one value per wavelength"); | |
| 233 | } | ||
| 234 | 398 | for (std::size_t i = 1; i < wavelengths_.size(); ++i) { | |
| 235 | 392 | if (!(wavelengths_[i] > wavelengths_[i - 1])) { | |
| 236 | 1 | throw std::invalid_argument("Refractive index wavelengths must be strictly " | |
| 237 | 1 | "increasing"); | |
| 238 | } | ||
| 239 | } | ||
| 240 | 7 | } | |
| 241 | |||
| 242 | 227 | std::complex<double> RefractiveIndex::operator()(double wavelengthNm) const { | |
| 243 | 227 | return Interpolate<std::complex<double>>(wavelengths_, values_, wavelengthNm); | |
| 244 | } | ||
| 245 | |||
| 246 | 161 | std::complex<double> MaxwellGarnett(std::complex<double> matrix, std::complex<double> inclusion, | |
| 247 | double fraction) { | ||
| 248 | 161 | const std::complex<double> em = matrix * matrix; | |
| 249 | 161 | const std::complex<double> ei = inclusion * inclusion; | |
| 250 | 161 | const std::complex<double> mixed = | |
| 251 | 161 | em * (ei + 2.0 * em + 2.0 * fraction * (ei - em)) / (ei + 2.0 * em - fraction * (ei - em)); | |
| 252 | 161 | return std::sqrt(mixed); | |
| 253 | } | ||
| 254 | |||
| 255 | 79 | ParticleOptics MieParticles(std::span<const LogNormalMode> modes, | |
| 256 | std::complex<double> refractiveIndex, double wavelengthNm, | ||
| 257 | std::span<const double> cosines, int radiusSamples) { | ||
| 258 | 79 | double totalNumber = 0.0; | |
| 259 | 222 | for (const LogNormalMode& mode : modes) { | |
| 260 | 143 | totalNumber += mode.relativeNumber; | |
| 261 | } | ||
| 262 | 79 | if (modes.empty() || !(totalNumber > 0.0) || radiusSamples < 2) { | |
| 263 | ✗ | throw std::invalid_argument("Mie particles need size modes with a positive number"); | |
| 264 | } | ||
| 265 | |||
| 266 | 79 | const double lambda = wavelengthNm * 1e-9; | |
| 267 | 79 | ParticleOptics result; | |
| 268 | 79 | result.phase.assign(cosines.size(), 0.0); | |
| 269 | 79 | std::vector<double> phase(cosines.size()); | |
| 270 | 222 | for (const LogNormalMode& mode : modes) { | |
| 271 | 143 | const double fraction = mode.relativeNumber / totalNumber; | |
| 272 | 143 | const double mean = std::log(mode.medianRadius); | |
| 273 | // Numbers fall off as exp(-t^2/2) in t = (ln r - mean) / | ||
| 274 | // sigma, and cross sections weight them by r^2, which moves | ||
| 275 | // their mass up by 2 sigma^2. | ||
| 276 | 143 | const double low = mean - 4.0 * mode.sigma; | |
| 277 | 143 | const double high = mean + 2.0 * mode.sigma * mode.sigma + 4.0 * mode.sigma; | |
| 278 | 143 | const int samples = mode.sigma > 0.0 ? radiusSamples : 1; | |
| 279 | 143 | double weightSum = 0.0; | |
| 280 | 143 | double modeScattering = 0.0; | |
| 281 | 143 | double modeAbsorption = 0.0; | |
| 282 | 143 | std::vector<double> modePhase(cosines.size(), 0.0); | |
| 283 | 13200 | for (int i = 0; i < samples; ++i) { | |
| 284 | 13057 | const double lnR = samples == 1 ? mean : low + (high - low) * i / (samples - 1); | |
| 285 | 13057 | const double t = mode.sigma > 0.0 ? (lnR - mean) / mode.sigma : 0.0; | |
| 286 | 13057 | const double trapezoid = (i == 0 || i == samples - 1) && samples > 1 ? 0.5 : 1.0; | |
| 287 | 13057 | const double weight = trapezoid * std::exp(-0.5 * t * t); | |
| 288 | 13057 | const double r = std::exp(lnR); | |
| 289 | 13057 | const MieEfficiencies q = | |
| 290 | 13057 | MieSphere(2.0 * pi * r / lambda, refractiveIndex, cosines, phase); | |
| 291 | 13057 | const double area = pi * r * r; | |
| 292 | 13057 | const double scattering = weight * q.scattering * area; | |
| 293 | 13057 | modeScattering += scattering; | |
| 294 | 13057 | modeAbsorption += weight * q.Absorption() * area; | |
| 295 | 3346756 | for (std::size_t j = 0; j < cosines.size(); ++j) { | |
| 296 | 3333699 | modePhase[j] += scattering * phase[j]; | |
| 297 | } | ||
| 298 | 13057 | weightSum += weight; | |
| 299 | } | ||
| 300 | 143 | result.scattering += fraction * modeScattering / weightSum; | |
| 301 | 143 | result.absorption += fraction * modeAbsorption / weightSum; | |
| 302 | 41016 | for (std::size_t j = 0; j < cosines.size(); ++j) { | |
| 303 | 40873 | result.phase[j] += fraction * modePhase[j] / weightSum; | |
| 304 | } | ||
| 305 | 143 | } | |
| 306 | 26661 | for (double& p : result.phase) { | |
| 307 | 26582 | p = result.scattering > 0.0 ? p / result.scattering : kIsotropicPhase; | |
| 308 | } | ||
| 309 | 79 | return result; | |
| 310 | 79 | } | |
| 311 | |||
| 312 | 60 | double OzoneCrossSection(double wavelengthNm) { | |
| 313 | // 10 nm averages starting at 360 nm (CosmoScout VR, | ||
| 314 | // ozoneMode.cpp, MIT). | ||
| 315 | 60 | static constexpr std::array<double, 48> kAverages = { | |
| 316 | 1.18006e-27, 2.18205e-28, 2.81764e-28, 6.63629e-28, 1.52685e-27, 2.76259e-27, 5.51975e-27, | ||
| 317 | 8.45102e-27, 1.58232e-26, 2.31555e-26, 3.66625e-26, 4.92413e-26, 7.76088e-26, 9.02900e-26, | ||
| 318 | 1.48333e-25, 1.60547e-25, 2.14349e-25, 2.76161e-25, 3.09823e-25, 3.50934e-25, 4.27703e-25, | ||
| 319 | 4.68477e-25, 4.40965e-25, 4.71385e-25, 5.03275e-25, 4.31623e-25, 3.74999e-25, 3.22324e-25, | ||
| 320 | 2.66819e-25, 2.24367e-25, 1.85651e-25, 1.47571e-25, 1.21105e-25, 9.43730e-26, 7.46292e-26, | ||
| 321 | 6.57117e-26, 5.10619e-26, 4.14823e-26, 4.22622e-26, 3.23257e-26, 2.44425e-26, 2.79549e-26, | ||
| 322 | 2.52744e-26, 1.61447e-26, 1.45506e-26, 2.07028e-26, 1.37295e-26, 6.98672e-27}; | ||
| 323 | 60 | constexpr double kFirstNm = 360.0; | |
| 324 | 60 | constexpr double kBinNm = 10.0; | |
| 325 | 60 | const double lastNm = kFirstNm + kBinNm * kAverages.size(); | |
| 326 | 60 | if (wavelengthNm < kFirstNm || wavelengthNm > lastNm) { | |
| 327 | return 0.0; | ||
| 328 | } | ||
| 329 | // Linear between bin centres, constant in the outer half bins. | ||
| 330 | 58 | const double position = (wavelengthNm - kFirstNm) / kBinNm - 0.5; | |
| 331 | 58 | const double clamped = std::clamp(position, 0.0, static_cast<double>(kAverages.size() - 1)); | |
| 332 | 58 | const std::size_t i = std::min(static_cast<std::size_t>(clamped), kAverages.size() - 2); | |
| 333 | 58 | const double t = clamped - static_cast<double>(i); | |
| 334 | 58 | return kAverages[i] + (kAverages[i + 1] - kAverages[i]) * t; | |
| 335 | } | ||
| 336 | |||
| 337 | 3 | const RefractiveIndex& HematiteRefractiveIndex() { | |
| 338 | // CosmoScout VR, | ||
| 339 | // scattering-table-generator/mie-settings/mars_realistic.json | ||
| 340 | // (CC0): every 10 nm from 220 to 1000 nm, n and k. | ||
| 341 | 3 | static constexpr std::array<std::array<double, 2>, 79> kTable = {{ | |
| 342 | {1.349, 1.127}, {1.517, 1.145}, {1.635, 1.145}, {1.733, 1.149}, {1.828, 1.148}, | ||
| 343 | {1.917, 1.133}, {1.992, 1.107}, {2.049, 1.076}, {2.090, 1.053}, {2.134, 1.041}, | ||
| 344 | {2.185, 1.018}, {2.218, 0.979}, {2.219, 0.947}, {2.206, 0.951}, {2.211, 0.987}, | ||
| 345 | {2.243, 1.034}, {2.305, 1.075}, {2.387, 1.097}, {2.475, 1.096}, {2.560, 1.075}, | ||
| 346 | {2.638, 1.040}, {2.708, 0.993}, {2.766, 0.936}, {2.809, 0.873}, {2.837, 0.813}, | ||
| 347 | {2.858, 0.758}, {2.876, 0.707}, {2.890, 0.654}, {2.893, 0.602}, {2.887, 0.560}, | ||
| 348 | {2.881, 0.532}, {2.885, 0.512}, {2.901, 0.492}, {2.927, 0.460}, {2.951, 0.413}, | ||
| 349 | {2.963, 0.351}, {2.956, 0.286}, {2.931, 0.229}, {2.898, 0.186}, {2.862, 0.156}, | ||
| 350 | {2.830, 0.134}, {2.800, 0.119}, {2.775, 0.110}, {2.753, 0.103}, {2.736, 0.096}, | ||
| 351 | {2.719, 0.089}, {2.703, 0.083}, {2.688, 0.078}, {2.675, 0.075}, {2.662, 0.072}, | ||
| 352 | {2.652, 0.068}, {2.641, 0.066}, {2.631, 0.063}, {2.621, 0.061}, {2.612, 0.060}, | ||
| 353 | {2.604, 0.059}, {2.596, 0.058}, {2.589, 0.057}, {2.582, 0.057}, {2.575, 0.057}, | ||
| 354 | {2.570, 0.058}, {2.566, 0.059}, {2.562, 0.059}, {2.559, 0.059}, {2.555, 0.058}, | ||
| 355 | {2.552, 0.058}, {2.549, 0.057}, {2.547, 0.056}, {2.544, 0.054}, {2.541, 0.053}, | ||
| 356 | {2.537, 0.053}, {2.535, 0.052}, {2.533, 0.051}, {2.531, 0.051}, {2.529, 0.049}, | ||
| 357 | {2.527, 0.048}, {2.525, 0.046}, {2.522, 0.045}, {2.520, 0.043}, | ||
| 358 | }}; | ||
| 359 | 3 | static const RefractiveIndex index = [] { | |
| 360 | 3 | std::vector<double> wavelengths; | |
| 361 | 3 | std::vector<std::complex<double>> values; | |
| 362 | 240 | for (std::size_t i = 0; i < kTable.size(); ++i) { | |
| 363 | 237 | wavelengths.push_back(220.0 + 10.0 * static_cast<double>(i)); | |
| 364 | 237 | values.emplace_back(kTable[i][0], kTable[i][1]); | |
| 365 | } | ||
| 366 | 6 | return RefractiveIndex(std::move(wavelengths), std::move(values)); | |
| 367 | 6 | }(); | |
| 368 | 3 | return index; | |
| 369 | } | ||
| 370 | |||
| 371 | } // namespace ernie | ||
| 372 |