File failed to load: https://isis.astrogeology.usgs.gov/9.0.0/Object/assets/jax/output/NativeMML/config.js
Isis 3 Programmer Reference
HapkeAtm1.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include <cmath>
8#include "AtmosModel.h"
9#include "Constants.h"
10#include "HapkeAtm1.h"
11#include "NumericalApproximation.h"
12#include "Pvl.h"
13#include "PvlGroup.h"
14#include "IException.h"
15#include "IString.h"
16
17using std::max;
18
19namespace Isis {
20 HapkeAtm1::HapkeAtm1(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
21 }
22
63 void HapkeAtm1::AtmosModelAlgorithm(double phase, double incidence, double emission) {
64 double munot, mu;
65 double xx;
66 double emunot, emu;
67 double xmunot, ymunot;
68 double xmu, ymu;
69 double gmunot, gmu;
70 double hpsq1;
71 double munotp, mup;
72 double fix;
73 double hahgt, hahgt0;
74 double phasefn;
75 double maxval;
76
77 if(p_atmosTau == 0.0) {
78 p_pstd = 0.0;
79 p_trans = 1.0;
80 p_trans0 = 1.0;
81 p_sbar = 0.0;
82 p_transs = 1.0;
83 return;
84 }
85
86 if(TauOrWhaChanged()) {
87 // preparation includes exponential integrals e sub 2 through 4
88 p_wha2 = 0.5 * p_atmosWha;
89 p_e2 = AtmosModel::En(2, p_atmosTau);
90 p_e3 = AtmosModel::En(3, p_atmosTau);
91 p_e4 = AtmosModel::En(4, p_atmosTau);
92
93 // zeroth moments of (uncorrected) x and y times characteristic fn
94 p_x0 = p_wha2;
95 p_y0 = p_wha2 * p_e2;
96
97 // higher-order correction term for x and y
98 p_delta = (1.0 - (p_x0 + p_y0) - (1.0 - p_atmosWha) / (1.0 - (p_x0 - p_y0))) / (p_atmosWha * (0.5 - p_e3));
99
100 // moments of (corrected) x and y
101 p_alpha0 = 1.0 + p_delta * (0.5 - p_e3);
102 p_alpha1 = 0.5 + p_delta * ((1.0 / 3.0) - p_e4);
103 p_beta0 = p_e2 + p_delta * (0.5 - p_e3);
104 p_beta1 = p_e3 + p_delta * ((1.0 / 3.0) - p_e4);
105
106 // prepare to find correct mixture of x and y in conservative case
107 if(p_atmosWha == 1.0) {
108 p_e5 = AtmosModel::En(5, p_atmosTau);
109 p_alpha2 = (1.0 / 3.0) + p_delta * (0.25 - p_e5);
110 p_beta2 = p_e4 + p_delta * (0.25 - p_e5);
111 p_fixcon = (p_beta0 * p_atmosTau - p_alpha1 + p_beta1) /
112 ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
113 }
114
115
116 // gamma will be a weighted sum of x and y functions
117 p_gammax = p_wha2 * p_beta0;
118 p_gammay = 1.0 - p_wha2 * p_alpha0;
119
120
121 // sbar is total diffuse illumination
122 // isotropic part comes from moments, correction is numerical integral
123 if (p_atmosEstTau) {
125 } else {
127 }
128 p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1) + p_atmosHahgsb;
129
130 SetOldTau(p_atmosTau);
131 SetOldWha(p_atmosWha);
132 }
133
134 // correct the path lengths for planetary curvature
135 hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
136
137 if(incidence == 90.0) {
138 munot = 0.0;
139 }
140 else {
141 munot = cos((PI / 180.0) * incidence);
142 }
143
144 maxval = max(1.0e-30, hpsq1 + munot * munot);
145 munotp = p_atmosHnorm / (sqrt(maxval) - munot);
146 munotp = max(munotp, p_atmosTau / 69.0);
147
148 if(emission == 90.0) {
149 mu = 0.0;
150 }
151 else {
152 mu = cos((PI / 180.0) * emission);
153 }
154
155 maxval = max(1.0e-30, hpsq1 + mu * mu);
156 mup = p_atmosHnorm / (sqrt(maxval) - mu);
157 mup = max(mup, p_atmosTau / 69.0);
158// build the x and y functions of mu0 and mu
159 maxval = max(1.0e-30, munotp);
160 xx = -p_atmosTau / maxval;
161
162 if(xx < -69.0) {
163 emunot = 0.0;
164 }
165 else if(xx > 69.0) {
166 emunot = 1.0e30;
167 }
168 else {
169 emunot = exp(-p_atmosTau / munotp);
170 }
171
172 maxval = max(1.0e-30, mup);
173 xx = -p_atmosTau / maxval;
174
175 if(xx < -69.0) {
176 emu = 0.0;
177 }
178 else if(xx > 69.0) {
179 emu = 1.0e30;
180 }
181 else {
182 emu = exp(-p_atmosTau / mup);
183 }
184
185 xmunot = 1.0 + p_delta * munotp * (1.0 - emunot);
186 ymunot = emunot + p_delta * munotp * (1.0 - emunot);
187 xmu = 1.0 + p_delta * mup * (1.0 - emu);
188 ymu = emu + p_delta * mup * (1.0 - emu);
189
190 // mix the x and y as required in the conservative case
191 if(p_atmosWha == 1.0) {
192 fix = p_fixcon * munotp * (xmunot + ymunot);
193 xmunot = xmunot + fix;
194 ymunot = ymunot + fix;
195 fix = p_fixcon * mup * (xmu + ymu);
196 xmu = xmu + fix;
197 ymu = ymu + fix;
198 }
199
200 // gamma1 functions come from x and y, with a correction for
201 // highly forward-scattered light as tabulated in hahgtTable
202 if (p_atmosEstTau) {
205 NumericalAtmosApprox qromb;
206 p_atmosAtmSwitch = 1;
207 qromb.Reset();
208 p_atmosInc = incidence;
209 p_atmosMunot = cos((PI / 180.0) * incidence);
210 p_atmosSini = sin((PI / 180.0) * incidence);
211 hahgt = qromb.RombergsMethod(this, sub, 0, 180);
212 gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt * AtmosWha() / 360.0;
213 p_atmosInc = emission;
214 p_atmosMunot = cos((PI / 180.0) * emission);
215 p_atmosSini = sin((PI / 180.0) * emission);
216 hahgt = qromb.RombergsMethod(this, sub, 0, 180);
217 gmu = p_gammax * xmu + p_gammay * ymu + hahgt * AtmosWha() / 360.0;
218 } else {
219 hahgt = p_atmosHahgtSpline.Evaluate(incidence, NumericalApproximation::Extrapolate);
220 gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt;
222 gmu = p_gammax * xmu + p_gammay * ymu + hahgt;
223 }
224
225 // purely atmos term uses x and y (plus single-particle phase
226 // function correction)
227 if(phase == 90.0) {
228 phasefn = (1.0 - p_atmosHga * p_atmosHga) / pow(1.0 + 2.0 * p_atmosHga * 0.0 + p_atmosHga * p_atmosHga, 1.5);
229 }
230 else {
231 phasefn = (1.0 - p_atmosHga * p_atmosHga) /
232 pow(1.0 + 2.0 * p_atmosHga * cos((PI / 180.0) * phase) + p_atmosHga * p_atmosHga, 1.5);
233 }
234
235 p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) *
236 ((xmunot * xmu - ymunot * ymu) + (phasefn - 1.0) * (1.0 - emu * emunot));
237
238 // xmitted surface term uses gammas
239 p_trans = gmunot * gmu;
240
241 // finally, never-scattered term is given by pure attenuation, with
242 // a correction for highly forward-scattered light (on the way down
243 // but not on the way up) as tabulated in hahgt0Table
244 if (p_atmosEstTau) {
247 NumericalAtmosApprox qromb;
248 p_atmosAtmSwitch = 3;
249 qromb.Reset();
250 p_atmosInc = incidence;
251 p_atmosMunot = cos((PI / 180.0) * incidence);
252 p_atmosSini = sin((PI / 180.0) * incidence);
253 hahgt0 = qromb.RombergsMethod(this, sub, 0, 180);
254 hahgt0 = hahgt0 * AtmosWha() * p_atmosMunot / (360.0 * p_atmosSini);
255 } else {
256 hahgt0 = p_atmosHahgt0Spline.Evaluate(incidence, NumericalApproximation::Extrapolate);
257 }
258 p_trans0 = (emunot + hahgt0) * emu;
259
260 // Calculate the transmission of light that must be subtracted from a shadow. This
261 // includes direct flux and the scattered flux in the upsun half of the sky
262 // downwelling onto the surface, and the usual transmission upward.
263 if (p_atmosEstTau) {
266 NumericalAtmosApprox qromb;
267 p_atmosAtmSwitch = 1;
268 qromb.Reset();
269 hahgt = qromb.RombergsMethod(this, sub, 90, 180);
270 hahgt = .5 * (p_gammax * xmunot + p_gammay * ymunot - emunot) + hahgt *
271 AtmosWha() / 360.0;
272 } else {
273 hahgt = p_atmosHahgtSpline.Evaluate(incidence, NumericalApproximation::Extrapolate);
274 }
275 p_transs = (emunot + hahgt) * emu;
276 }
277}
278
279extern "C" Isis::AtmosModel *HapkeAtm1Plugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel) {
280 return new Isis::HapkeAtm1(pvl, pmodel);
281}
Isotropic atmos scattering model.
Definition AtmosModel.h:60
double p_trans
Transmission of surface reflected light through the atmosphere overall.
Definition AtmosModel.h:259
static double En(unsigned int n, double x)
This routine evaluates the generalized exponential integral, En(x).
double p_atmosHnorm
Atmospheric shell thickness normalized to planet radius.
Definition AtmosModel.h:267
double p_sbar
Illumination of the ground by the sky.
Definition AtmosModel.h:262
bool p_atmosEstTau
Estimate optical depth tau using shadows.
Definition AtmosModel.h:269
double p_transs
Transmission of light that must be subtracted from the flat surface model to get the shadow model.
Definition AtmosModel.h:261
NumericalApproximation p_atmosHahgt0Spline
Spline object for the atmospheric Hahg0 Table. Properties are set in GenerateHahgTables().
Definition AtmosModel.h:287
void GenerateHahgTablesShadow()
This method is a modified version of the GenerateHahgTables method and is used solely for shadow mode...
double p_pstd
Pure atmospheric-scattering term.
Definition AtmosModel.h:258
double p_trans0
Transmission of surface reflected light through the atmosphere with no scatterings in the atmosphere.
Definition AtmosModel.h:260
NumericalApproximation p_atmosHahgtSpline
Spline object for the atmospheric Hahg Table. Properties are set in GenerateHahgTables().
Definition AtmosModel.h:285
double AtmosWha() const
Return atmospheric Wha value.
Definition AtmosModel.h:123
void GenerateHahgTables()
This method computes the values of the atmospheric Hahg and Hahg0 tables and sets the properties of t...
bool TauOrWhaChanged() const
Checks whether tau or wha have changed.
Implements the Hapke Atmospheric Model.
Definition HapkeAtm1.h:40
virtual void AtmosModelAlgorithm(double phase, double incidence, double emission)
Henyey-Greenstein atmos scattering in the 1st approximation.
Definition HapkeAtm1.cpp:63
void Reset()
Resets the state of the object.
@ Extrapolate
Evaluate() attempts to extrapolate if a is outside of the domain. This is only valid for NumericalApp...
double RombergsMethod(AtmosModel *am, IntegFunc sub, double a, double b)
This variation on the NumericalApproximation method integrates a specified AtmosModel function rather...
IntegFunc
This enum defines function to be integrated by Romberg's method.
@ OuterFunction
Indicates that Romberg's method will integrate the function OutrFunc2Bint()
Container for cube-like labels.
Definition Pvl.h:119
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
const double PI
The mathematical constant PI.
Definition Constants.h:40