GCC Code Coverage Report

ernie/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
3 / 3
100.0%
Branches:
2 / 4
50.0%
Functions:
-
-%
physics/atmosphere.h
Line Branch Exec Source
1 // SPDX-License-Identifier: BSD-3-Clause
2 /// @file atmosphere.h
3 /// @brief Planetary atmospheres as physical constituents, and their
4 /// reduction to a camera's spectral bands.
5 /// Altitudes and radii are in m, number densities m^-3, cross
6 /// sections m^2 and coefficients m^-1.
7 #pragma once
8
9 #include "ernie/physics/scattering.h"
10 #include "ernie/physics/spectrum.h"
11
12 #include <functional>
13 #include <span>
14 #include <string>
15 #include <vector>
16
17 namespace ernie {
18
19 /// @brief One term of a vertical density profile, relative to its
20 /// reference density.
21 struct DensityLayer {
22 /// @brief The layer's shape over altitude h.
23 enum class Shape {
24 /// exp(-h / height): an isothermal gas.
25 Exponential,
26 /// exp(1 - exp(h / height)), 1 at the ground: dust falling
27 /// off faster than exponentially (Conrath; Schneegans et
28 /// al.).
29 DoubleExponential,
30 /// 1 - |h - height| / width, clamped at 0: peaks at `height`.
31 Tent,
32 };
33
34 /// The shape.
35 Shape shape = Shape::Exponential;
36 /// Scale height or peak altitude, in metres.
37 double height = 8000.0;
38 /// A tent's half width, in metres.
39 double width = 0.0;
40 /// The layer's weight in the profile's mean.
41 double weight = 1.0;
42 };
43
44 /// @brief A vertical density profile: the weighted mean of its
45 /// layers, or a table.
46 /// Exponential layers give 1 at the ground.
47 ///
48 /// References:
49 /// - B. J. Conrath, "Thermal structure of the Martian atmosphere
50 /// during the dissipation of the dust storm of 1971", Icarus
51 /// 24(1):36-46, 1975. https://doi.org/10.1016/0019-1035(75)90156-6
52 /// - S. Schneegans, T. Meyran, I. Ginkel, G. Zachmann and A. Gerndt,
53 /// "Physically Based Real-Time Rendering of Atmospheres using Mie
54 /// Theory", Computer Graphics Forum 43(2), 2024, eq. (7).
55 /// https://doi.org/10.1111/cgf.15010
56
2/4
✓ Branch 3 taken 1 time.
✗ Branch 4 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 time.
 19 struct DensityProfile {
57 /// The layers, when `table` is empty.
58 std::vector<DensityLayer> layers;
59 /// Relative densities every `tableStep` m from the ground, used
60 /// in place of the layers when not empty; linear between, last
61 /// value beyond.
62 std::vector<double> table;
63 /// The table's altitude step, in metres.
64 double tableStep = 0.0;
65
66 /// @param scaleHeight The scale height, in metres.
67 /// @return One exponential layer.
68 static DensityProfile Exponential(double scaleHeight);
69 /// @param peakAltitude The peak's altitude, in metres.
70 /// @param halfWidth The half width, in metres.
71 /// @return One tent layer.
72 static DensityProfile Tent(double peakAltitude, double halfWidth);
73
74 /// @brief A hydrostatic gas under the U.S. Standard Atmosphere
75 /// 1976's lapse rates, shifted to a surface temperature.
76 /// Tabulated every 100 m to 120 km; heights are treated as
77 /// geometric (the standard's are geopotential, 1.6% apart at 100
78 /// km).
79 /// @param surfaceTemperature Ground temperature, in kelvin.
80 /// @param molarMass Molar mass, in kg/mol.
81 /// @param gravity Surface gravity, in m/s^2.
82 /// @return Density relative to the ground's.
83 ///
84 /// References:
85 /// - NOAA, NASA and USAF, "U.S. Standard Atmosphere, 1976",
86 /// NOAA-S/T 76-1562, 1976, table 4.
87 /// https://ntrs.nasa.gov/citations/19770009539
88 static DensityProfile StandardAtmosphere(double surfaceTemperature, double molarMass,
89 double gravity);
90
91 /// @param altitude The altitude, in metres.
92 /// @return The relative density; 0 without layers or a table.
93 double operator()(double altitude) const;
94
95 /// @param top The upper altitude, in metres.
96 /// @return The relative density integrated from the ground.
97 double Column(double top) const;
98 };
99
100 /// Optical properties of one particle or molecule at a wavelength
101 /// in nm, with the phase function at the given cosines.
102 using OpticsModel = std::function<ParticleOptics(double wavelengthNm, std::span<const double>)>;
103
104 /// @brief Dry air molecules, Rayleigh after Bodhaine et al. (1999).
105 /// @param co2Fraction The CO2 volume fraction.
106 /// @return The optics model.
107 OpticsModel AirMolecules(double co2Fraction = 360e-6);
108
109 /// @brief Molecules of a CO2 atmosphere such as Mars's.
110 /// CO2 scatters per Owens and Bates; the remainder (mostly N2 and Ar
111 /// on Mars, within 7% of air) scatters as air.
112 /// @param co2Fraction The CO2 volume fraction.
113 /// @return The optics model.
114 OpticsModel CarbonDioxideMolecules(double co2Fraction);
115
116 /// @return Ozone molecules, absorbing only (`OzoneCrossSection`).
117 OpticsModel OzoneMolecules();
118
119 /// @brief Mie spheres with a log-normal size distribution
120 /// (`MieParticles`).
121 /// @param modes The size distribution.
122 /// @param refractiveIndex The index, possibly varying by wavelength.
123 /// @param radiusSamples Integration samples per mode.
124 /// @return The optics model.
125 OpticsModel MieAerosol(std::vector<LogNormalMode> modes, RefractiveIndex refractiveIndex,
126 int radiusSamples = 96);
127
128 /// @brief One constituent of an atmosphere.
129 struct Constituent {
130 /// Its name.
131 std::string name;
132 /// Particles per m^3 where the profile is 1.
133 double numberDensity = 0.0;
134 /// Its vertical profile.
135 DensityProfile density;
136 /// Its optics per particle.
137 OpticsModel optics;
138 /// Refractivity n - 1 by wavelength in nm where the profile is 1,
139 /// scaling with density (Gladstone-Dale); empty for constituents
140 /// too sparse to bend light, such as aerosols and ozone.
141 std::function<double(double)> refractivity;
142 };
143
144 /// @brief A spherical shell of constituents around a planet.
145 7 struct AtmosphereModel {
146 /// The planet's radius, in metres.
147 double bottomRadius = 6371e3;
148 /// The shell's thickness, in metres.
149 double topAltitude = 100e3;
150 /// The constituents.
151 std::vector<Constituent> constituents;
152
153 /// @param wavelengthNm The wavelength, in nm.
154 /// @return Extinction optical depth from the ground straight up.
155 double VerticalOpticalDepth(double wavelengthNm) const;
156 };
157
158 /// @brief Earth's atmosphere: an isothermal dry-air column, ozone in
159 /// a tent profile, and haze.
160 struct EarthAtmosphereOptions {
161 /// Sea-level pressure, in Pa; with the temperature, gives the
162 /// molecules' number density p / (k T).
163 double surfacePressure = 101325.0;
164 /// Sea-level temperature, in kelvin.
165 double surfaceTemperature = 288.15;
166 /// Surface gravity, in m/s^2.
167 double gravity = 9.80665;
168 /// The CO2 volume fraction.
169 double co2Fraction = 360e-6;
170 /// Total column ozone in Dobson units, spread over a 10-40 km
171 /// tent peaking at 25 km as in Bruneton's model.
172 double ozoneColumn = 300.0;
173 /// Vertical aerosol optical depth at 550 nm.
174 double aerosolOpticalDepth = 0.1;
175 /// The aerosol scale height, in metres.
176 double aerosolScaleHeight = 1200.0;
177 };
178
179 /// @brief Builds Earth's atmosphere.
180 /// The haze is Schneegans et al.'s bimodal fit to AERONET, index
181 /// 1.4 + 0.015i.
182 /// @param radius The planet's radius, in metres.
183 /// @param options The atmosphere's settings.
184 /// @return The model.
185 ///
186 /// References:
187 /// - NOAA, NASA and USAF, "U.S. Standard Atmosphere, 1976", NOAA-S/T
188 /// 76-1562, 1976. https://ntrs.nasa.gov/citations/19770009539
189 /// - E. Bruneton, "Precomputed Atmospheric Scattering: a New
190 /// Implementation", 2017 (ozone profile and column).
191 /// https://ebruneton.github.io/precomputed_atmospheric_scattering/
192 /// - Schneegans et al. 2024 (see `DensityProfile`), section 3 (haze).
193 /// - B. N. Holben et al., "AERONET - A Federated Instrument Network
194 /// and Data Archive for Aerosol Characterization", Remote Sensing
195 /// of Environment 66(1):1-16, 1998.
196 /// https://doi.org/10.1016/S0034-4257(98)00031-5
197 AtmosphereModel EarthAtmosphere(double radius, const EarthAtmosphereOptions& options = {});
198
199 /// @brief Mars's atmosphere: an isothermal CO2 column and dust.
200 struct MarsAtmosphereOptions {
201 /// Mean surface pressure, in Pa.
202 double surfacePressure = 610.0;
203 /// Mean surface temperature, in kelvin.
204 double surfaceTemperature = 210.0;
205 /// Surface gravity, in m/s^2.
206 double gravity = 3.71;
207 /// The CO2 volume fraction, after Curiosity's SAM.
208 double co2Fraction = 0.951;
209 /// Mean molar mass, in kg/mol.
210 double molarMass = 43.42e-3;
211 /// Vertical dust optical depth at 880 nm, the rovers' reference
212 /// (clear seasons 0.3-0.9, dust storms above 1).
213 double dustOpticalDepth = 0.5;
214 };
215
216 /// @brief Builds Mars's atmosphere.
217 /// Dust is feldspar (n = 1.5) with 3% hematite by Maxwell Garnett
218 /// (Ehlers et al.), Schneegans et al.'s bimodal sizes (0.04 and 1.4
219 /// um, 1000:1) in two double-exponential layers.
220 /// @param radius The planet's radius, in metres.
221 /// @param options The atmosphere's settings.
222 /// @return The model.
223 ///
224 /// References:
225 /// - H. B. Franz et al., "Initial SAM calibration gas experiments on
226 /// Mars: Quadrupole mass spectrometer results and implications",
227 /// Planetary and Space Science 138:44-54, 2017.
228 /// https://doi.org/10.1016/j.pss.2017.01.014
229 /// - M. T. Lemmon et al., "Dust aerosol, clouds, and the atmospheric
230 /// optical depth record over 5 Mars years of the Mars Exploration
231 /// Rover mission", Icarus 251:96-111, 2015.
232 /// https://doi.org/10.1016/j.icarus.2014.03.029
233 /// - K. Ehlers, R. Chakrabarty and H. Moosmueller, "Blue moons and
234 /// Martian sunsets", Applied Optics 53(9):1808-1819, 2014.
235 /// https://doi.org/10.1364/AO.53.001808
236 /// - Schneegans et al. 2024 (see `DensityProfile`), section 3 and
237 /// table 1.
238 AtmosphereModel MarsAtmosphere(double radius, const MarsAtmosphereOptions& options = {});
239
240 /// @brief Scattering angles for tabulated phase functions.
241 /// theta_i = pi (i / (count - 1))^2, dense near forward where
242 /// diffraction peaks are narrow.
243 /// @param count The number of angles.
244 /// @return The angles' cosines.
245 std::vector<double> PhaseCosines(int count);
246
247 /// @brief A constituent averaged over a camera's bands.
248 struct BandedConstituent {
249 /// Its name.
250 std::string name;
251 /// Its vertical profile.
252 DensityProfile density;
253 /// Scattering coefficient per band where the profile is 1, in
254 /// m^-1.
255 std::vector<float> scattering;
256 /// Absorption coefficient per band where the profile is 1, in
257 /// m^-1.
258 std::vector<float> absorption;
259 /// Phase function per band at `PhaseCosines`, band-major:
260 /// phase[band * angles + i].
261 std::vector<float> phase;
262 /// Refractivity n - 1 per band where the profile is 1; zeros
263 /// without one.
264 std::vector<float> refractivity;
265 };
266
267 /// @brief An atmosphere averaged over a camera's bands, ready to
268 /// render.
269 2 struct BandedAtmosphere {
270 /// The planet's radius, in metres.
271 double bottomRadius = 0.0;
272 /// The shell's thickness, in metres.
273 double topAltitude = 0.0;
274 /// The bands averaged over.
275 SpectralBands bands;
276 /// Phase function samples (`PhaseCosines`).
277 int phaseAngles = 0;
278 /// The constituents.
279 std::vector<BandedConstituent> constituents;
280
281 /// @param altitude The altitude, in metres.
282 /// @param band The band.
283 /// @return The extinction coefficient, in m^-1.
284 double Extinction(double altitude, int band) const;
285
286 /// @brief Transmittance to the top along a ray above the
287 /// horizon, by the midpoint rule.
288 /// The reference the renderer's transmittance table is tested
289 /// against.
290 /// @param r The start radius, in metres.
291 /// @param mu The ray's cosine from the zenith.
292 /// @param band The band.
293 /// @param steps Integration steps.
294 /// @return The transmittance.
295 double TransmittanceToTop(double r, double mu, int band, int steps = 4000) const;
296
297 /// @param altitude The altitude, in metres.
298 /// @param band The band.
299 /// @return 1 plus each constituent's refractivity times density.
300 double RefractiveIndex(double altitude, int band) const;
301
302 /// @brief Where a refracted ray ends (`Refract`).
303 struct Refraction {
304 /// Angle between the first and last directions, in radians:
305 /// its bending towards the ground.
306 double deviation = 0.0;
307 /// Whether it reached the ground rather than the top.
308 bool ground = false;
309 };
310 /// @brief Traces a refracted ray out of the atmosphere by the ray
311 /// equation d(n t)/ds = grad n.
312 /// @param r The start radius, in metres.
313 /// @param mu The apparent zenith cosine.
314 /// @param band The band.
315 /// @param step The largest step, in metres.
316 /// @return The ray's deviation and where it ended.
317 Refraction Refract(double r, double mu, int band, double step = 100.0) const;
318 };
319
320 /// @brief Averages an atmosphere over a camera's bands.
321 /// Samples each band at 10 nm or finer; phase functions are weighted
322 /// by scattering to preserve scattered radiance.
323 /// @param atmosphere The model.
324 /// @param bands The bands.
325 /// @param phaseAngles Phase function samples (`PhaseCosines`).
326 /// @return The banded atmosphere.
327 BandedAtmosphere BandAverage(const AtmosphereModel& atmosphere, const SpectralBands& bands,
328 int phaseAngles = 256);
329
330 } // namespace ernie
331