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