GCC Code Coverage Report

ernie/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
7 / 7
100.0%
Branches:
5 / 8
62.5%
Functions:
-
-%
physics/scattering.h
Line Branch Exec Source
1 // SPDX-License-Identifier: BSD-3-Clause
2 /// @file scattering.h
3 /// @brief Rayleigh and Mie scattering and absorption cross sections.
4 /// Wavelengths are vacuum nm, radii m, cross sections m^2. Phase
5 /// functions are in sr^-1, integrate to 1 and take the scattering
6 /// angle's cosine (1 is forward).
7 #pragma once
8
9 #include <complex>
10 #include <numbers>
11 #include <span>
12 #include <utility>
13 #include <vector>
14
15 namespace ernie {
16
17 /// The isotropic phase function, 1 / (4 pi) sr^-1: light scattered equally in every direction.
18 inline constexpr double kIsotropicPhase = 1.0 / (4.0 * std::numbers::pi);
19
20 /// Ideal gas number density at 288.15 K and 101325 Pa, in m^-3: the
21 /// reference state of the refractivities below.
22 ///
23 /// References:
24 /// - B. A. Bodhaine, N. B. Wood, E. G. Dutton and J. R. Slusser, "On
25 /// Rayleigh Optical Depth Calculations", Journal of Atmospheric and
26 /// Oceanic Technology 16(11):1854-1861, 1999, eq. (24).
27 /// https://doi.org/10.1175/1520-0426(1999)016<1854:ORODC>2.0.CO;2
28 inline constexpr double kStandardNumberDensity = 2.546899e25;
29
30 /// @brief Dry air's refractivity at 288.15 K and 101325 Pa.
31 /// Peck and Reeder's formula with Bodhaine et al.'s CO2 correction.
32 /// @param wavelengthNm The wavelength, in nm.
33 /// @param co2Fraction The CO2 volume fraction.
34 /// @return n - 1.
35 ///
36 /// References:
37 /// - E. R. Peck and K. Reeder, "Dispersion of Air", Journal of the
38 /// Optical Society of America 62(8):958-962, 1972.
39 /// https://doi.org/10.1364/JOSA.62.000958
40 /// - Bodhaine et al. 1999 (see `kStandardNumberDensity`), eqs. (18)
41 /// and (19).
42 double AirRefractivity(double wavelengthNm, double co2Fraction = 360e-6);
43
44 /// @brief CO2's refractivity at 288.15 K and 101325 Pa.
45 /// @param wavelengthNm The wavelength, in nm.
46 /// @return n - 1.
47 ///
48 /// References:
49 /// - J. C. Owens, "Optical Refractive Index of Air: Dependence on
50 /// Pressure, Temperature and Composition", Applied Optics
51 /// 6(1):51-59, 1967. https://doi.org/10.1364/AO.6.000051
52 /// - Bodhaine et al. 1999 (see `kStandardNumberDensity`), eq. (27).
53 double CarbonDioxideRefractivity(double wavelengthNm);
54
55 /// @brief Dry air's King factor for molecular anisotropy.
56 /// Bates's N2 and O2 factors, 1 for Ar and 1.15 for CO2, weighted
57 /// by volume fraction.
58 /// @param wavelengthNm The wavelength, in nm.
59 /// @param co2Fraction The CO2 volume fraction.
60 /// @return (6 + 3 rho) / (6 - 7 rho).
61 ///
62 /// References:
63 /// - D. R. Bates, "Rayleigh scattering by air", Planetary and Space
64 /// Science 32(6):785-790, 1984.
65 /// https://doi.org/10.1016/0032-0633(84)90102-8
66 /// - Bodhaine et al. 1999 (see `kStandardNumberDensity`), eqs. (5),
67 /// (6) and (23).
68 double AirKingFactor(double wavelengthNm, double co2Fraction = 360e-6);
69
70 /// CO2's King factor (Bates 1984, as used by Bodhaine et al. 1999).
71 inline constexpr double kCarbonDioxideKingFactor = 1.15;
72
73 /// @brief Rayleigh cross section of one molecule.
74 /// 24 pi^3 (n^2 - 1)^2 / (lambda^4 Ns^2 (n^2 + 2)^2) F; independent
75 /// of density since n - 1 scales with it.
76 /// @param wavelengthNm The wavelength, in nm.
77 /// @param refractivity n - 1 at `kStandardNumberDensity`.
78 /// @param kingFactor The King factor F.
79 /// @return The cross section, in m^2.
80 ///
81 /// References:
82 /// - Bodhaine et al. 1999 (see `kStandardNumberDensity`), eq. (22).
83 double RayleighCrossSection(double wavelengthNm, double refractivity, double kingFactor);
84
85 /// @brief The Rayleigh phase function with depolarisation.
86 /// 3 / (16 pi) ((1 + 3 g) + (1 - g) mu^2) / (1 + 2 g), with
87 /// g = rho / (2 - rho) and rho the depolarisation ratio.
88 /// @param cosTheta The scattering angle's cosine.
89 /// @param kingFactor The King factor; 1 for isotropic molecules.
90 /// @return The phase function, in sr^-1.
91 ///
92 /// References:
93 /// - S. Chandrasekhar, "Radiative Transfer", Oxford University Press,
94 /// 1950, chapter I.
95 /// - A. Bucholtz, "Rayleigh-scattering calculations for the
96 /// terrestrial atmosphere", Applied Optics 34(15):2765-2773, 1995,
97 /// eqs. (12)-(13). https://doi.org/10.1364/AO.34.002765
98 double RayleighPhase(double cosTheta, double kingFactor);
99
100 /// @brief Mie efficiencies of a homogeneous sphere.
101 /// Efficiencies are cross sections over pi r^2.
102 struct MieEfficiencies {
103 /// Extinction efficiency.
104 double extinction = 0.0;
105 /// Scattering efficiency.
106 double scattering = 0.0;
107 /// Asymmetry parameter g, the mean scattering cosine.
108 double asymmetry = 0.0;
109
110 /// @return The absorption efficiency, at least 0.
111 13058 double Absorption() const {
112 // Rounding can leave a non-absorbing sphere slightly
113 // negative.
114
3/4
✓ Branch 0 taken 1 time.
✓ Branch 0 taken 12193 times.
✗ Branch 1 not taken.
✓ Branch 1 taken 864 times.
 13058 return extinction > scattering ? extinction - scattering : 0.0;
115 }
116 };
117
118 /// @brief Mie scattering by a homogeneous sphere.
119 /// Sums the series to Wiscombe's number of terms.
120 /// @param sizeParameter 2 pi r / lambda.
121 /// @param refractiveIndex Relative index n + ik, k >= 0 absorbing.
122 /// @param cosines Cosines to evaluate the phase function at.
123 /// @param phase Receives the unpolarised phase function at
124 /// `cosines`; empty to skip.
125 /// @return The efficiencies.
126 ///
127 /// References:
128 /// - G. Mie, "Beitraege zur Optik trueber Medien, speziell
129 /// kolloidaler Metalloesungen", Annalen der Physik 330(3):377-445,
130 /// 1908. https://doi.org/10.1002/andp.19083300302
131 /// - C. F. Bohren and D. R. Huffman, "Absorption and Scattering of
132 /// Light by Small Particles", Wiley, 1983, chapter 4 and appendix A
133 /// (BHMIE). https://doi.org/10.1002/9783527618156
134 /// - W. J. Lentz, "Generating Bessel functions in Mie scattering
135 /// calculations using continued fractions", Applied Optics
136 /// 15(3):668-671, 1976. https://doi.org/10.1364/AO.15.000668
137 /// - W. J. Wiscombe, "Improved Mie scattering algorithms", Applied
138 /// Optics 19(9):1505-1509, 1980.
139 /// https://doi.org/10.1364/AO.19.001505
140 MieEfficiencies MieSphere(double sizeParameter, std::complex<double> refractiveIndex,
141 std::span<const double> cosines, std::span<double> phase);
142
143 /// @brief One mode of a log-normal particle size distribution.
144 /// Particles per unit ln(r) are Gaussian.
145 ///
146 /// References:
147 /// - J. E. Hansen and L. D. Travis, "Light scattering in planetary
148 /// atmospheres", Space Science Reviews 16(4):527-610, 1974,
149 /// section 2.4 (effective radius and variance).
150 /// https://doi.org/10.1007/BF00168069
151 struct LogNormalMode {
152 /// Median radius, in metres.
153 double medianRadius = 1e-6;
154 /// Standard deviation of ln(r).
155 double sigma = 0.5;
156 /// Particle count relative to the other modes.
157 double relativeNumber = 1.0;
158
159 /// @brief A mode from its effective radius and variance.
160 /// sigma^2 = ln(1 + v_eff); r_median = r_eff exp(-2.5 sigma^2).
161 /// @param effectiveRadius Area-weighted mean radius, in metres.
162 /// @param effectiveVariance The effective variance.
163 /// @param relativeNumber Count relative to other modes.
164 /// @return The mode.
165 static LogNormalMode FromEffective(double effectiveRadius, double effectiveVariance,
166 double relativeNumber = 1.0);
167
168 /// @return The effective radius, in metres.
169 double EffectiveRadius() const;
170 /// @return The effective variance.
171 double EffectiveVariance() const;
172 };
173
174 /// @brief A complex refractive index tabulated by wavelength.
175 /// Linear between samples, constant beyond the ends.
176
2/4
✓ Branch 8 taken 3 times.
✗ Branch 9 not taken.
✓ Branch 13 taken 5 times.
✗ Branch 14 not taken.
 11 class RefractiveIndex {
177 public:
178 RefractiveIndex() = default;
179
180 /// @brief A constant index.
181 /// @param value The index at every wavelength.
182 explicit RefractiveIndex(std::complex<double> value);
183
184 /// @brief Tabulates an index.
185 /// @param wavelengthsNm Strictly increasing wavelengths, in nm.
186 /// @param values The index at each.
187 /// @throws std::invalid_argument If the wavelengths are not
188 /// strictly increasing or the sizes differ.
189 RefractiveIndex(std::vector<double> wavelengthsNm, std::vector<std::complex<double>> values);
190
191 /// @param wavelengthNm The wavelength, in nm.
192 /// @return The interpolated index.
193 std::complex<double> operator()(double wavelengthNm) const;
194
195 private:
196 std::vector<double> wavelengths_;
197 std::vector<std::complex<double>> values_;
198 };
199
200 /// @brief Effective index of a matrix with spherical inclusions, by
201 /// the Maxwell Garnett rule on permittivities m^2.
202 /// @param matrix The matrix's index.
203 /// @param inclusion The inclusions' index.
204 /// @param fraction The inclusions' volume fraction.
205 /// @return The effective index.
206 ///
207 /// References:
208 /// - J. C. Maxwell Garnett, "Colours in Metal Glasses and in Metallic
209 /// Films", Philosophical Transactions of the Royal Society A
210 /// 203:385-420, 1904. https://doi.org/10.1098/rsta.1904.0024
211 /// - Bohren and Huffman 1983 (see `MieSphere`), section 8.5.
212 std::complex<double> MaxwellGarnett(std::complex<double> matrix, std::complex<double> inclusion,
213 double fraction);
214
215 /// @brief A particle population's mean optics per particle at one
216 /// wavelength.
217 192 struct ParticleOptics {
218 /// Scattering cross section, in m^2.
219 double scattering = 0.0;
220 /// Absorption cross section, in m^2.
221 double absorption = 0.0;
222 /// Phase function at the requested cosines, in sr^-1.
223 std::vector<double> phase;
224 };
225
226 /// @brief Mean Mie optics of a multi-mode log-normal population.
227 /// Each mode is integrated over ln(r) by trapezoids, which unlike
228 /// Monte Carlo keeps tables smooth across wavelengths. The phase
229 /// function is weighted by scattering cross section.
230 /// @param modes The size distribution's modes.
231 /// @param refractiveIndex The particles' index.
232 /// @param wavelengthNm The wavelength, in nm.
233 /// @param cosines Cosines to evaluate the phase function at.
234 /// @param radiusSamples Trapezoid samples per mode.
235 /// @return The mean optics per particle.
236 /// @throws std::invalid_argument If there are no modes or their
237 /// numbers do not sum above zero.
238 ///
239 /// References:
240 /// - Hansen and Travis 1974 (see `LogNormalMode`), section 2.
241 ParticleOptics MieParticles(std::span<const LogNormalMode> modes,
242 std::complex<double> refractiveIndex, double wavelengthNm,
243 std::span<const double> cosines, int radiusSamples = 96);
244
245 /// @brief Ozone's absorption cross section at 233 K.
246 /// 10 nm averages of the IUP Bremen data over 360-840 nm (Chappuis
247 /// band included), as tabulated by CosmoScout VR.
248 /// @param wavelengthNm The wavelength, in nm.
249 /// @return The cross section in m^2; 0 outside 360-840 nm.
250 ///
251 /// References:
252 /// - V. Gorshelev, A. Serdyuchenko, M. Weber, W. Chehade and J. P.
253 /// Burrows, "High spectral resolution ozone absorption
254 /// cross-sections - Part 1: Measurements, data analysis and
255 /// comparison with previous measurements around 293 K", Atmospheric
256 /// Measurement Techniques 7(2):609-624, 2014.
257 /// https://doi.org/10.5194/amt-7-609-2014
258 /// - A. Serdyuchenko, V. Gorshelev, M. Weber, W. Chehade and J. P.
259 /// Burrows, "High spectral resolution ozone absorption
260 /// cross-sections - Part 2: Temperature dependence", Atmospheric
261 /// Measurement Techniques 7(2):625-636, 2014.
262 /// https://doi.org/10.5194/amt-7-625-2014
263 /// - CosmoScout VR,
264 /// plugins/csp-atmospheres/scattering-table-generator/ozoneMode.cpp
265 /// (MIT), German Aerospace Center (DLR).
266 /// https://github.com/cosmoscout/cosmoscout-vr
267 double OzoneCrossSection(double wavelengthNm);
268
269 /// @brief Hematite's (alpha-Fe2O3) refractive index, 220-1000 nm.
270 /// As tabulated by CosmoScout VR (CC0) after Ehlers et al.; mixed
271 /// with feldspar for Martian dust.
272 /// @return The tabulated index.
273 ///
274 /// References:
275 /// - K. Ehlers, R. Chakrabarty and H. Moosmueller, "Blue moons and
276 /// Martian sunsets", Applied Optics 53(9):1808-1819, 2014.
277 /// https://doi.org/10.1364/AO.53.001808
278 /// - S. Schneegans, T. Meyran, I. Ginkel, G. Zachmann and A. Gerndt,
279 /// "Physically Based Real-Time Rendering of Atmospheres using Mie
280 /// Theory", Computer Graphics Forum 43(2), 2024.
281 /// https://doi.org/10.1111/cgf.15010
282 const RefractiveIndex& HematiteRefractiveIndex();
283
284 } // namespace ernie
285