GCC Code Coverage Report

ernie/
Coverage:
low: ≥ 0%
medium: ≥ 75.0%
high: ≥ 90.0%
Lines:
276 / 300
92.0%
Branches:
144 / 234
61.5%
Functions:
27 / 29
93.1%
physics/atmosphere.cpp
Line Branch Exec Source
1 // SPDX-License-Identifier: BSD-3-Clause
2 /// @file atmosphere.cpp
3 /// @brief Atmosphere models, density profiles and band averaging.
4 #include "ernie/physics/atmosphere.h"
5
6 #include <algorithm>
7 #include <cmath>
8 #include <iterator>
9 #include <numbers>
10 #include <stdexcept>
11 #include <utility>
12
13 namespace ernie {
14
15 namespace {
16
17 /// One Dobson unit, in molecules m^-2.
18 constexpr double kDobsonUnit = 2.6867e20;
19
20 /// The molar mass of dry air (kg/mol) in the U.S. Standard Atmosphere
21 /// 1976.
22 constexpr double kAirMolarMass = 28.9644e-3;
23
24 /// Intervals of the Simpson rule integrating density profiles.
25 constexpr int kColumnIntervals = 4000;
26
27 /// The longest spacing between samples when averaging over a band
28 /// (nm).
29 constexpr double kBandSampleNm = 10.0;
30
31 335190 double LayerDensity(const DensityLayer& layer, double altitude) {
32
3/4
✓ Branch 0 taken 134242 times.
✓ Branch 1 taken 24010 times.
✓ Branch 2 taken 176938 times.
✗ Branch 3 not taken.
 335190 switch (layer.shape) {
33 134242 case DensityLayer::Shape::Exponential:
34 134242 return std::exp(-altitude / layer.height);
35 24010 case DensityLayer::Shape::DoubleExponential:
36 24010 return std::exp(1.0 - std::exp(altitude / layer.height));
37 176938 case DensityLayer::Shape::Tent:
38
2/2
✓ Branch 0 taken 55048 times.
✓ Branch 1 taken 121890 times.
 231986 return std::max(0.0, 1.0 - std::abs(altitude - layer.height) / layer.width);
39 }
40 return 0.0;
41 }
42
43 /// An isothermal gas column: number density p / (k T) at the ground
44 /// and scale height R T / (M g).
45 3 std::pair<double, DensityProfile> IsothermalGas(double pressure, double temperature,
46 double molarMass, double gravity) {
47 3 const double numberDensity = pressure / (kBoltzmann * temperature);
48 3 const double scaleHeight = kGasConstant * temperature / (molarMass * gravity);
49 6 return {numberDensity, DensityProfile::Exponential(scaleHeight)};
50 }
51
52 /// The number density that gives `profile` the vertical optical depth
53 /// `opticalDepth` at `wavelengthNm` up to `top`.
54 5 double NumberDensityForDepth(const OpticsModel& optics, const DensityProfile& profile, double top,
55 double wavelengthNm, double opticalDepth) {
56 5 const ParticleOptics sample = optics(wavelengthNm, {});
57 5 const double extinction = sample.scattering + sample.absorption;
58 5 const double column = profile.Column(top);
59
2/4
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 5 if (!(extinction > 0.0) || !(column > 0.0)) {
60 ✗ throw std::invalid_argument("An aerosol needs extinction and "
61
0/2
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 ✗ "a column to set its depth");
62 }
63 5 return opticalDepth / (extinction * column);
64 5 }
65
66 } // namespace
67
68 8 DensityProfile DensityProfile::Exponential(double scaleHeight) {
69 8 DensityProfile profile;
70
1/2
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
 8 profile.layers = {{DensityLayer::Shape::Exponential, scaleHeight, 0.0, 1.0}};
71 8 return profile;
72 ✗ }
73
74 6 DensityProfile DensityProfile::Tent(double peakAltitude, double halfWidth) {
75 6 DensityProfile profile;
76
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 profile.layers = {{DensityLayer::Shape::Tent, peakAltitude, halfWidth, 1.0}};
77 6 return profile;
78 ✗ }
79
80 5 DensityProfile DensityProfile::StandardAtmosphere(double surfaceTemperature, double molarMass,
81 double gravity) {
82 // U.S. Standard Atmosphere 1976, table 4: the layers' base
83 // heights (m) and lapse rates (K/m) to 84.852 km.
84 5 struct Layer {
85 double base;
86 double lapse;
87 };
88 5 static constexpr Layer kLayers[] = {{0.0, -6.5e-3}, {11000.0, 0.0}, {20000.0, 1.0e-3},
89 {32000.0, 2.8e-3}, {47000.0, 0.0}, {51000.0, -2.8e-3},
90 {71000.0, -2.0e-3}, {84852.0, 0.0}};
91 66005 const auto temperature = [&](double h) {
92 66000 double t = surfaceTemperature;
93
2/2
✓ Branch 0 taken 528000 times.
✓ Branch 1 taken 66000 times.
 594000 for (std::size_t i = 0; i < std::size(kLayers); ++i) {
94
2/2
✓ Branch 0 taken 462000 times.
✓ Branch 1 taken 66000 times.
 528000 const double end = i + 1 < std::size(kLayers) ? kLayers[i + 1].base : 1e30;
95
2/2
✓ Branch 0 taken 174235 times.
✓ Branch 1 taken 353765 times.
 1056000 t += kLayers[i].lapse * (std::clamp(h, kLayers[i].base, end) - kLayers[i].base);
96 }
97 66000 return t;
98 5 };
99 5 constexpr double kStep = 100.0;
100 5 constexpr int kSamples = 1201; // to 120 km
101 5 constexpr int kSubSteps = 10;
102 5 const double scale = molarMass * gravity / kGasConstant;
103 5 DensityProfile profile;
104 5 profile.tableStep = kStep;
105
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 profile.table.reserve(kSamples);
106 5 double logPressure = 0.0; // relative to the ground's
107
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 profile.table.push_back(1.0);
108
2/2
✓ Branch 0 taken 6000 times.
✓ Branch 1 taken 5 times.
 6005 for (int i = 1; i < kSamples; ++i) {
109 // d ln p / dh = -M g / (R T), by the midpoint rule.
110
2/2
✓ Branch 0 taken 60000 times.
✓ Branch 1 taken 6000 times.
 66000 for (int s = 0; s < kSubSteps; ++s) {
111 60000 const double h = (i - 1) * kStep + (s + 0.5) * kStep / kSubSteps;
112 60000 logPressure -= scale / temperature(h) * kStep / kSubSteps;
113 }
114 6000 const double h = i * kStep;
115
1/2
✓ Branch 1 taken 6000 times.
✗ Branch 2 not taken.
 6000 profile.table.push_back(std::exp(logPressure) * surfaceTemperature / temperature(h));
116 }
117 5 return profile;
118 ✗ }
119
120 414263 double DensityProfile::operator()(double altitude) const {
121
2/2
✓ Branch 0 taken 91077 times.
✓ Branch 1 taken 323186 times.
 414263 if (!table.empty()) {
122 91077 const double x =
123
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 91077 times.
 91077 std::clamp(altitude / tableStep, 0.0, static_cast<double>(table.size() - 1));
124
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 91077 times.
 91077 const auto i = std::min(static_cast<std::size_t>(x), table.size() - 2);
125 91077 const double t = x - static_cast<double>(i);
126 91077 return table[i] + (table[i + 1] - table[i]) * t;
127 }
128 323186 double sum = 0.0;
129 323186 double weights = 0.0;
130
2/2
✓ Branch 0 taken 335190 times.
✓ Branch 1 taken 323186 times.
 658376 for (const DensityLayer& layer : layers) {
131 335190 sum += layer.weight * LayerDensity(layer, altitude);
132 335190 weights += layer.weight;
133 }
134
2/2
✓ Branch 0 taken 323185 times.
✓ Branch 1 taken 1 time.
 323186 return weights > 0.0 ? sum / weights : 0.0;
135 }
136
137 23 double DensityProfile::Column(double top) const {
138 23 const double step = top / kColumnIntervals;
139 23 double sum = (*this)(0.0) + (*this)(top);
140
2/2
✓ Branch 0 taken 91977 times.
✓ Branch 1 taken 23 times.
 92000 for (int i = 1; i < kColumnIntervals; ++i) {
141
2/2
✓ Branch 0 taken 45977 times.
✓ Branch 1 taken 46000 times.
 137954 sum += (i % 2 == 1 ? 4.0 : 2.0) * (*this)(i * step);
142 }
143 23 return sum * step / 3.0;
144 }
145
146 5 OpticsModel AirMolecules(double co2Fraction) {
147 63 return [co2Fraction](double wavelengthNm, std::span<const double> cosines) {
148 58 const double king = AirKingFactor(wavelengthNm, co2Fraction);
149 58 ParticleOptics optics;
150 116 optics.scattering =
151
2/4
✓ Branch 1 taken 58 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 58 times.
✗ Branch 5 not taken.
 58 RayleighCrossSection(wavelengthNm, AirRefractivity(wavelengthNm, co2Fraction), king);
152
2/2
✓ Branch 0 taken 12320 times.
✓ Branch 1 taken 58 times.
 12378 for (const double mu : cosines) {
153
2/4
✓ Branch 1 taken 12320 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 12320 times.
✗ Branch 5 not taken.
 12320 optics.phase.push_back(RayleighPhase(mu, king));
154 }
155 58 return optics;
156 5 };
157 }
158
159 3 OpticsModel CarbonDioxideMolecules(double co2Fraction) {
160 7 return [co2Fraction](double wavelengthNm, std::span<const double> cosines) {
161 4 const double co2 = co2Fraction * RayleighCrossSection(
162 wavelengthNm, CarbonDioxideRefractivity(wavelengthNm),
163 4 kCarbonDioxideKingFactor);
164 4 const double airKing = AirKingFactor(wavelengthNm);
165 4 const double rest =
166 4 (1.0 - co2Fraction) *
167 4 RayleighCrossSection(wavelengthNm, AirRefractivity(wavelengthNm), airKing);
168 4 ParticleOptics optics;
169 4 optics.scattering = co2 + rest;
170 // Each gas scatters with its own depolarisation; the
171 // mixture's phase function is their mean weighted by
172 // scattering.
173
2/2
✓ Branch 0 taken 32 times.
✓ Branch 1 taken 4 times.
 36 for (const double mu : cosines) {
174
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 optics.phase.push_back((co2 * RayleighPhase(mu, kCarbonDioxideKingFactor) +
175
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 64 rest * RayleighPhase(mu, airKing)) /
176
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 optics.scattering);
177 }
178 4 return optics;
179 3 };
180 }
181
182 5 OpticsModel OzoneMolecules() {
183 58 return [](double wavelengthNm, std::span<const double> cosines) {
184 53 ParticleOptics optics;
185
1/2
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
 53 optics.absorption = OzoneCrossSection(wavelengthNm);
186
1/2
✓ Branch 1 taken 53 times.
✗ Branch 2 not taken.
 53 optics.phase.assign(cosines.size(), kIsotropicPhase);
187 53 return optics;
188 5 };
189 }
190
191 5 OpticsModel MieAerosol(std::vector<LogNormalMode> modes, RefractiveIndex refractiveIndex,
192 int radiusSamples) {
193
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 return [modes = std::move(modes), refractiveIndex = std::move(refractiveIndex),
194 radiusSamples](double wavelengthNm, std::span<const double> cosines) {
195 63 return MieParticles(modes, refractiveIndex(wavelengthNm), wavelengthNm, cosines,
196 63 radiusSamples);
197
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 };
198 }
199
200 1 double AtmosphereModel::VerticalOpticalDepth(double wavelengthNm) const {
201 1 double depth = 0.0;
202
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
 4 for (const Constituent& constituent : constituents) {
203 3 const ParticleOptics optics = constituent.optics(wavelengthNm, {});
204 3 depth += constituent.numberDensity * (optics.scattering + optics.absorption) *
205 3 constituent.density.Column(topAltitude);
206 3 }
207 1 return depth;
208 }
209
210 5 AtmosphereModel EarthAtmosphere(double radius, const EarthAtmosphereOptions& options) {
211 5 AtmosphereModel atmosphere;
212 5 atmosphere.bottomRadius = radius;
213 5 atmosphere.topAltitude = 100e3;
214
215 // Air with the U.S. Standard Atmosphere's temperatures: their
216 // fall through the troposphere shapes both the scattering and,
217 // through the density's gradient at the ground, the refraction
218 // near the horizon (some 3 arcminutes less than an isothermal
219 // column's).
220 5 const double airDensity = options.surfacePressure / (kBoltzmann * options.surfaceTemperature);
221 5 const DensityProfile airProfile = DensityProfile::StandardAtmosphere(
222
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 options.surfaceTemperature, kAirMolarMass, options.gravity);
223 // Air's refractivity at the reference state scaled to its density
224 // (refractivity is proportional to it, Gladstone-Dale).
225 5 const double airRatio = airDensity / kStandardNumberDensity;
226
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 10 atmosphere.constituents.push_back({"molecules", airDensity, airProfile,
227 5 AirMolecules(options.co2Fraction),
228
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
 5 [airRatio, co2 = options.co2Fraction](double nm) {
229 50 return AirRefractivity(nm, co2) * airRatio;
230 }});
231
232 // A tent from 10 to 40 km holds the column in a triangle of base
233 // 30 km: peak * 15 km.
234 5 const DensityProfile ozoneProfile = DensityProfile::Tent(25e3, 15e3);
235 5 const double ozonePeak = options.ozoneColumn * kDobsonUnit / 15e3;
236
3/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 2 times.
 10 atmosphere.constituents.push_back({"ozone", ozonePeak, ozoneProfile, OzoneMolecules(), {}});
237
238 5 if (options.aerosolOpticalDepth > 0.0) {
239 // Schneegans et al. 2024, CosmoScout VR
240 // mie-settings/earth_haze.json (CC0): median radii
241 // exp(-16.316) and exp(-13.616) m, sigma 0.6, 10000 small
242 // particles per large one.
243
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 OpticsModel haze =
244 MieAerosol({{std::exp(-16.316), 0.6, 10000.0}, {std::exp(-13.616), 0.6, 1.0}},
245
3/6
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
 3 RefractiveIndex({1.4, 0.015}));
246
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 const DensityProfile hazeProfile = DensityProfile::Exponential(options.aerosolScaleHeight);
247 6 const double hazeDensity = NumberDensityForDepth(haze, hazeProfile, atmosphere.topAltitude,
248
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 550.0, options.aerosolOpticalDepth);
249
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 atmosphere.constituents.push_back(
250 3 {"aerosols", hazeDensity, hazeProfile, std::move(haze), {}});
251 3 }
252 5 return atmosphere;
253
12/22
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 5 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 5 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 5 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 5 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 5 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 3 times.
✓ Branch 22 taken 2 times.
✓ Branch 24 taken 3 times.
✗ Branch 25 not taken.
✓ Branch 27 taken 3 times.
✗ Branch 28 not taken.
✓ Branch 31 taken 3 times.
✗ Branch 32 not taken.
 31 }
254
255 3 AtmosphereModel MarsAtmosphere(double radius, const MarsAtmosphereOptions& options) {
256 3 AtmosphereModel atmosphere;
257 3 atmosphere.bottomRadius = radius;
258 3 atmosphere.topAltitude = 100e3;
259
260
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 const auto [gasDensity, gasProfile] = IsothermalGas(
261
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 3 options.surfacePressure, options.surfaceTemperature, options.molarMass, options.gravity);
262 // CO2's refractivity and the remainder's (as air's) at the
263 // reference state, scaled to the gas's density (Gladstone-Dale).
264 3 const double gasRatio = gasDensity / kStandardNumberDensity;
265
1/2
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
 6 atmosphere.constituents.push_back(
266 9 {"molecules", gasDensity, gasProfile, CarbonDioxideMolecules(options.co2Fraction),
267
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
 3 [gasRatio, co2 = options.co2Fraction](double nm) {
268 2 return (co2 * CarbonDioxideRefractivity(nm) + (1.0 - co2) * AirRefractivity(nm)) *
269 2 gasRatio;
270 }});
271
272 3 if (options.dustOpticalDepth > 0.0) {
273 // Feldspar with 3% hematite (Ehlers et al. 2014), indexed
274 // where the hematite table is.
275
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const RefractiveIndex& hematite = HematiteRefractiveIndex();
276 2 std::vector<double> wavelengths;
277 2 std::vector<std::complex<double>> values;
278
2/2
✓ Branch 0 taken 158 times.
✓ Branch 1 taken 2 times.
 160 for (double nm = 220.0; nm <= 1000.0; nm += 10.0) {
279
1/2
✓ Branch 1 taken 158 times.
✗ Branch 2 not taken.
 158 wavelengths.push_back(nm);
280
2/4
✓ Branch 1 taken 158 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 158 times.
✗ Branch 5 not taken.
 316 values.push_back(MaxwellGarnett({1.5, 0.0}, hematite(nm), 0.03));
281 }
282
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 OpticsModel dust = MieAerosol({LogNormalMode::FromEffective(0.04e-6, 0.1, 1000.0),
283
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 LogNormalMode::FromEffective(1.4e-6, 0.35, 1.0)},
284
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 6 RefractiveIndex(std::move(wavelengths), std::move(values)));
285 // Dense low dust and thin high dust (Schneegans et al. 2024,
286 // table 1).
287 2 DensityProfile dustProfile;
288 4 dustProfile.layers = {{DensityLayer::Shape::DoubleExponential, 4000.0, 0.0, 3.0},
289
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 {DensityLayer::Shape::DoubleExponential, 20000.0, 0.0, 1.0}};
290 4 const double dustDensity = NumberDensityForDepth(dust, dustProfile, atmosphere.topAltitude,
291
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 880.0, options.dustOpticalDepth);
292
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 atmosphere.constituents.push_back({"dust", dustDensity, dustProfile, std::move(dust), {}});
293 2 }
294 3 return atmosphere;
295
8/14
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 3 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✓ Branch 10 taken 1 time.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
 13 }
296
297 5 std::vector<double> PhaseCosines(int count) {
298
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 4 times.
 5 if (count < 2) {
299
1/2
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 1 throw std::invalid_argument("A phase table needs at least two angles");
300 }
301 4 std::vector<double> cosines(static_cast<std::size_t>(count));
302
2/2
✓ Branch 0 taken 293 times.
✓ Branch 1 taken 4 times.
 297 for (int i = 0; i < count; ++i) {
303 293 const double u = static_cast<double>(i) / (count - 1);
304 293 cosines[static_cast<std::size_t>(i)] = std::cos(std::numbers::pi * u * u);
305 }
306 4 return cosines;
307 }
308
309 ✗ double BandedAtmosphere::Extinction(double altitude, int band) const {
310 ✗ const auto b = static_cast<std::size_t>(band);
311 ✗ double extinction = 0.0;
312
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 ✗ for (const BandedConstituent& constituent : constituents) {
313 ✗ extinction += (static_cast<double>(constituent.scattering[b]) +
314 ✗ static_cast<double>(constituent.absorption[b])) *
315 ✗ constituent.density(altitude);
316 }
317 ✗ return extinction;
318 }
319
320 ✗ double BandedAtmosphere::TransmittanceToTop(double r, double mu, int band, int steps) const {
321 ✗ const double top = bottomRadius + topAltitude;
322 // Where the ray leaves the sphere of radius `top`.
323 ✗ const double discriminant = r * r * (mu * mu - 1.0) + top * top;
324
0/4
✗ Branch 0 not taken.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 ✗ const double length = std::max(0.0, -r * mu + std::sqrt(std::max(discriminant, 0.0)));
325 ✗ const double step = length / steps;
326 ✗ double depth = 0.0;
327
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
 ✗ for (int i = 0; i < steps; ++i) {
328 ✗ const double t = (i + 0.5) * step;
329 ✗ const double radius = std::sqrt(r * r + t * t + 2.0 * r * mu * t);
330 ✗ depth += Extinction(radius - bottomRadius, band) * step;
331 }
332 ✗ return std::exp(-depth);
333 }
334
335 165302 double BandedAtmosphere::RefractiveIndex(double altitude, int band) const {
336 165302 const auto b = static_cast<std::size_t>(band);
337 165302 double n = 1.0;
338
2/2
✓ Branch 0 taken 322233 times.
✓ Branch 1 taken 165302 times.
 487535 for (const BandedConstituent& constituent : constituents) {
339 322233 n += static_cast<double>(constituent.refractivity[b]) * constituent.density(altitude);
340 }
341 165302 return n;
342 }
343
344 6 BandedAtmosphere::Refraction BandedAtmosphere::Refract(double r, double mu, int band,
345 double step) const {
346 // In the plane of the ray and the planet's centre: x across, z
347 // up through the start.
348 6 const double top = bottomRadius + topAltitude;
349 165300 const auto index = [&](double radius) { return RefractiveIndex(radius - bottomRadius, band); };
350 // The gradient of n, radial: dn/dr times the unit position.
351 55106 const auto bend = [&](double px, double pz, double tx, double tz, double& ax, double& az) {
352 55100 const double radius = std::sqrt(px * px + pz * pz);
353 55100 const double h = 1.0;
354 55100 const double dndr = (index(radius + h) - index(radius - h)) / (2.0 * h);
355 55100 const double gx = dndr * px / radius;
356 55100 const double gz = dndr * pz / radius;
357 55100 const double along = gx * tx + gz * tz;
358 55100 const double n = index(radius);
359 55100 ax = (gx - along * tx) / n;
360 55100 az = (gz - along * tz) / n;
361 55106 };
362 6 double px = 0.0;
363 6 double pz = r;
364
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 6 const double sine = std::sqrt(std::max(1.0 - mu * mu, 0.0));
365 6 double tx = sine;
366 6 double tz = mu;
367 6 Refraction result;
368
1/2
✓ Branch 0 taken 27556 times.
✗ Branch 1 not taken.
 27556 for (int i = 0; i < 1000000; ++i) {
369 27556 const double radius = std::sqrt(px * px + pz * pz);
370 27556 const bool outwards = px * tx + pz * tz > 0.0;
371
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 27551 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 27556 if (radius > top && outwards) {
372 break;
373 }
374
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 27550 times.
 27551 if (radius < bottomRadius) {
375 1 result.ground = true;
376 1 break;
377 }
378 // The midpoint rule on the ray equation: t' = (grad n -
379 // (grad n . t) t) / n.
380 27550 double ax = 0.0;
381 27550 double az = 0.0;
382 27550 bend(px, pz, tx, tz, ax, az);
383 27550 const double mx = px + 0.5 * step * tx;
384 27550 const double mz = pz + 0.5 * step * tz;
385 27550 double mtx = tx + 0.5 * step * ax;
386 27550 double mtz = tz + 0.5 * step * az;
387 27550 bend(mx, mz, mtx, mtz, ax, az);
388 27550 px += step * mtx;
389 27550 pz += step * mtz;
390 27550 tx += step * ax;
391 27550 tz += step * az;
392 27550 const double length = std::sqrt(tx * tx + tz * tz);
393 27550 tx /= length;
394 27550 tz /= length;
395 }
396 // atan2 of the cross and dot products keeps small angles (Mars's
397 // microradians) that acos rounds away.
398 6 result.deviation = std::atan2(std::abs(sine * tz - mu * tx), tx * sine + tz * mu);
399 6 return result;
400 }
401
402 3 BandedAtmosphere BandAverage(const AtmosphereModel& atmosphere, const SpectralBands& bands,
403 int phaseAngles) {
404 3 bands.Check();
405 3 const std::vector<double> cosines = PhaseCosines(phaseAngles);
406
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 const int samples = std::max(1, static_cast<int>(std::ceil(bands.Width() / kBandSampleNm)));
407
408 3 BandedAtmosphere result;
409 3 result.bottomRadius = atmosphere.bottomRadius;
410 3 result.topAltitude = atmosphere.topAltitude;
411 3 result.bands = bands;
412 3 result.phaseAngles = phaseAngles;
413
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 3 times.
 9 for (const Constituent& constituent : atmosphere.constituents) {
414 6 BandedConstituent banded;
415
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 banded.name = constituent.name;
416
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 banded.density = constituent.density;
417
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 banded.phase.resize(static_cast<std::size_t>(bands.count) * cosines.size());
418
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 6 times.
 57 for (int band = 0; band < bands.count; ++band) {
419 51 double scattering = 0.0;
420 51 double absorption = 0.0;
421 51 double refractivity = 0.0;
422
1/2
✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
 51 std::vector<double> phase(cosines.size(), 0.0);
423
2/2
✓ Branch 0 taken 150 times.
✓ Branch 1 taken 51 times.
 201 for (int s = 0; s < samples; ++s) {
424
1/2
✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
 150 const double nm = bands.Lower(band) + (s + 0.5) * bands.Width() / samples;
425
1/2
✓ Branch 1 taken 150 times.
✗ Branch 2 not taken.
 150 const ParticleOptics optics = constituent.optics(nm, cosines);
426
2/2
✓ Branch 0 taken 52 times.
✓ Branch 1 taken 98 times.
 150 if (constituent.refractivity) {
427
1/2
✓ Branch 1 taken 52 times.
✗ Branch 2 not taken.
 52 refractivity += constituent.refractivity(nm);
428 }
429 150 scattering += optics.scattering;
430 150 absorption += optics.absorption;
431
2/2
✓ Branch 0 taken 36960 times.
✓ Branch 1 taken 150 times.
 37110 for (std::size_t i = 0; i < cosines.size(); ++i) {
432 36960 phase[i] += optics.scattering * optics.phase[i];
433 }
434 150 }
435 102 banded.scattering.push_back(
436
1/2
✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
 51 static_cast<float>(constituent.numberDensity * scattering / samples));
437 102 banded.absorption.push_back(
438
1/2
✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
 51 static_cast<float>(constituent.numberDensity * absorption / samples));
439
1/2
✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
 51 banded.refractivity.push_back(static_cast<float>(refractivity / samples));
440 51 float* out = banded.phase.data() + static_cast<std::size_t>(band) * cosines.size();
441
2/2
✓ Branch 0 taken 12336 times.
✓ Branch 1 taken 51 times.
 12387 for (std::size_t i = 0; i < cosines.size(); ++i) {
442 12336 out[i] =
443
2/2
✓ Branch 0 taken 8224 times.
✓ Branch 1 taken 4112 times.
 12336 static_cast<float>(scattering > 0.0 ? phase[i] / scattering : kIsotropicPhase);
444 }
445 51 }
446
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
 6 result.constituents.push_back(std::move(banded));
447 6 }
448 3 return result;
449 3 }
450
451 } // namespace ernie
452