USGS

Isis 3.0 Object Programmers' Reference

 Home
HapkeAtm2.cpp
1 #include <cmath>
2 #include "AtmosModel.h"
3 #include "Constants.h"
4 #include "HapkeAtm2.h"
5 #include "Pvl.h"
6 #include "PvlGroup.h"
7 #include "NumericalApproximation.h"
8 #include "IException.h"
9 #include "IString.h"
10 
11 using std::max;
12 
13 namespace Isis {
14  HapkeAtm2::HapkeAtm2(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
15  }
56  void HapkeAtm2::AtmosModelAlgorithm(double phase, double incidence,
57  double p_emission) {
58  double xx;
59  double maxval;
60  double munot, mu;
61  double p_emunot, emu;
62  double munotp, mup;
63  double f1munot, f1mmunot;
64  double xmunot, ymunot;
65  double xmu, ymu;
66  double gmunot, gmu;
67  double hpsq1;
68  double fix;
69  double f1mu, f1mmu;
70  double hahgt, hahgt0;
71  double phasefn;
72 
73  if(p_atmosTau == 0.0) {
74  p_pstd = 0.0;
75  p_trans = 1.0;
76  p_trans0 = 1.0;
77  p_sbar = 0.0;
78  p_transs = 1.0;
79  return;
80  }
81 
82  if(TauOrWhaChanged()) {
83 
84  // preparation includes exponential integrals p_e sub 2 through 4
85  p_wha2 = 0.5 * p_atmosWha;
86  p_e1 = AtmosModel::En(1, p_atmosTau);
87  p_e1_2 = AtmosModel::En(1, 2.0 * p_atmosTau);
88  p_e2 = AtmosModel::En(2, p_atmosTau);
89  p_e3 = AtmosModel::En(3, p_atmosTau);
90  p_e4 = AtmosModel::En(4, p_atmosTau);
91 
92  // chandra's gmn functions require fm and fn at mu=-1
93  xx = -p_atmosTau;
94  if(xx < -69.0) {
95  p_em = 0.0;
96  }
97  else if(xx > 69.0) {
98  p_em = 1.0e30;
99  }
100  else {
101  p_em = exp(xx);
102  }
103 
104  p_f1m = log(2.0) - p_em * p_e1 + p_e1_2;
105  p_f2m = -1.0 * (p_f1m + p_em * p_e2 - 1.0);
106  p_f3m = -1.0 * (p_f2m + p_em * p_e3 - 0.5);
107  p_g12 = (p_atmosTau * p_e1 * p_e2 + p_f1m + p_f2m) * 0.5;
108  p_g13 = (p_atmosTau * p_e1 * p_e3 + p_f1m + p_f3m) * (1.0 / 3.0);
109 
110  // chandra's g'mn functions require g'11 and f at mu=+1
111  xx = p_atmosTau;
112  if(xx < -69.0) {
113  p_e = 0.0;
114  }
115  else if(xx > 69.0) {
116  p_e = 1.0e30;
117  }
118  else {
119  p_e = exp(xx);
120  }
121 
122  p_f1 = Eulgam() + log(p_atmosTau) + p_e * p_e1;
123  p_f2 = p_f1 + p_e * p_e2 - 1.0;
124  p_f3 = p_f2 + p_e * p_e3 - 0.5;
125  p_g11p = AtmosModel::G11Prime(p_atmosTau);
126  p_g12p = (p_atmosTau * (p_e1 - p_g11p) + p_em * (p_f1 + p_f2)) * 0.25;
127  p_g13p = (p_atmosTau * (0.5 * p_e1 - p_g12p) + p_em * (p_f1 + p_f3)) * 0.2;
128 
129  // zeroth moments of (uncorrected) x and y times characteristic fn
130  p_x0 = p_wha2 * (1.0 + p_wha2 * p_g12);
131  p_y0 = p_wha2 * (p_e2 + p_wha2 * p_g12p);
132 
133  // higher-order correction term for x and y
134  p_delta = (1.0 - (p_x0 + p_y0) - (1.0 - p_atmosWha) / (1.0 - (p_x0 - p_y0))) / (p_atmosWha * (0.5 - p_e3));
135 
136  // moments of (corrected) x and y
137  p_alpha0 = 1.0 + p_wha2 * p_g12 + p_delta * (0.5 - p_e3);
138  p_alpha1 = 0.5 + p_wha2 * p_g13 + p_delta * ((1.0 / 3.0) - p_e4);
139  p_beta0 = p_e2 + p_wha2 * p_g12p + p_delta * (0.5 - p_e3);
140  p_beta1 = p_e3 + p_wha2 * p_g13p + p_delta * ((1.0 / 3.0) - p_e4);
141 
142  // prepare to find correct mixture of x and y in conservative case
143  if(p_atmosWha == 1.0) {
144  p_e5 = AtmosModel::En(5, p_atmosTau);
145  p_f4m = -1.0 * (p_f3m + p_em * p_e4 - (1.0 / 3.0));
146  p_g14 = (p_atmosTau * p_e1 * p_e4 + p_f1m + p_f4m) * 0.25;
147  p_f4 = p_f3 + p_e * p_e4 - (1.0 / 3.0);
148  p_g14p = (p_atmosTau * (0.5 * p_e1 - p_g13p) + p_em * (p_f1 + p_f4)) * (1.0 / 6.0);
149  p_alpha2 = (1.0 / 3.0) + p_wha2 * p_g14 + p_delta * (0.25 - p_e5);
150  p_beta2 = p_e4 + p_wha2 * p_g14p + p_delta * (0.25 - p_e5);
151  p_fixcon = (p_beta0 * p_atmosTau - p_alpha1 + p_beta1) /
152  ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
153  }
154 
155  // gamma will be a weighted sum of x and y functions
156  p_gammax = p_wha2 * p_beta0;
157  p_gammay = 1.0 - p_wha2 * p_alpha0;
158 
159  // sbar is total diffuse illumination
160  // isotropic part comes from moments, correction is numerical integral
161  if (p_atmosEstTau) {
162  GenerateHahgTablesShadow();
163  } else {
164  GenerateHahgTables();
165  }
166  p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1) + p_atmosHahgsb;
167 
168  SetOldTau(p_atmosTau);
169  SetOldWha(p_atmosWha);
170  }
171 
172  // correct the path lengths for planetary curvature
173  hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
174  munot = cos((PI / 180.0) * incidence);
175  maxval = max(1.0e-30, hpsq1 + munot * munot);
176  munotp = p_atmosHnorm / (sqrt(maxval) - munot);
177  munotp = max(munotp, p_atmosTau / 69.0);
178  mu = cos((PI / 180.0) * p_emission);
179  maxval = max(1.0e-30, hpsq1 + mu * mu);
180  mup = p_atmosHnorm / (sqrt(maxval) - mu);
181  mup = max(mup, p_atmosTau / 69.0);
182 
183  // build the x and y functions of mu0 and mu
184  maxval = max(1.0e-30, munotp);
185  xx = -p_atmosTau / maxval;
186 
187  if(xx < -69.0) {
188  p_emunot = 0.0;
189  }
190  else if(xx > 69.0) {
191  p_emunot = 1.0e30;
192  }
193  else {
194  p_emunot = exp(-p_atmosTau / munotp);
195  }
196 
197  maxval = max(1.0e-30, mup);
198  xx = -p_atmosTau / maxval;
199  if(xx < -69.0) {
200  emu = 0.0;
201  }
202  else if(xx > 69.0) {
203  emu = 1.0e30;
204  }
205  else {
206  emu = exp(-p_atmosTau / mup);
207  }
208 
209  // in the second approximation the x and y include the p_f1 function
210  xx = munotp;
211  if(fabs(xx - 1.0) < 1.0e-10) {
212  f1munot = p_f1;
213  f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * p_emunot + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
214  }
215  else if(xx > 0.0) {
216  f1munot = xx * (log(xx / (1.0 - xx)) + p_e1 / p_emunot + AtmosModel::Ei(p_atmosTau * (1.0 / xx - 1.0)));
217  f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * p_emunot + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
218  }
219  else {
220  std::string msg = "Negative length of planetary curvature ";
221  msg += "encountered";
222  throw IException(IException::Unknown, msg, _FILEINFO_);
223  }
224 
225  xx = mup;
226  if(fabs(xx - 1.0) < 1.0e-10) {
227  f1mu = p_f1;
228  f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
229  }
230  else if(xx > 0.0) {
231  f1mu = xx * (log(xx / (1.0 - xx)) + p_e1 / emu + AtmosModel::Ei(p_atmosTau * (1.0 / xx - 1.0)));
232  f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
233  }
234  else {
235  std::string msg = "Negative length of planetary curvature ";
236  msg += "encountered";
237  throw IException(IException::Unknown, msg, _FILEINFO_);
238  }
239 
240  xmunot = 1.0 + p_wha2 * f1mmunot + p_delta * munotp * (1.0 - p_emunot);
241  ymunot = p_emunot * (1.0 + p_wha2 * f1munot) + p_delta * munotp * (1.0 - p_emunot);
242  xmu = 1.0 + p_wha2 * f1mmu + p_delta * mup * (1.0 - emu);
243  ymu = emu * (1.0 + p_wha2 * f1mu) + p_delta * mup * (1.0 - emu);
244 
245  // mix the x and y as required in the conservative case
246  if(p_atmosWha == 1.0) {
247  fix = p_fixcon * munotp * (xmunot + ymunot);
248  xmunot = xmunot + fix;
249  ymunot = ymunot + fix;
250  fix = p_fixcon * mup * (xmu + ymu);
251  xmu = xmu + fix;
252  ymu = ymu + fix;
253  }
254 
255  // gamma1 functions come from x and y, with a correction for
256  // highly forward-scattered light as tabulated in hahgtTable
257  if (p_atmosEstTau) {
258  NumericalAtmosApprox::IntegFunc sub;
259  sub = NumericalAtmosApprox::OuterFunction;
260  NumericalAtmosApprox qromb;
261  p_atmosAtmSwitch = 1;
262  qromb.Reset();
263  p_atmosInc = incidence;
264  p_atmosMunot = cos((PI / 180.0) * incidence);
265  p_atmosSini = sin((PI / 180.0) * incidence);
266  hahgt = qromb.RombergsMethod(this, sub, 0, 180);
267  gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt * AtmosWha() / 360.0;
268  p_atmosInc = p_emission;
269  p_atmosMunot = cos((PI / 180.0) * p_emission);
270  p_atmosSini = sin((PI / 180.0) * p_emission);
271  hahgt = qromb.RombergsMethod(this, sub, 0, 180);
272  gmu = p_gammax * xmu + p_gammay * ymu + hahgt * AtmosWha() / 360.0;
273  } else {
274  hahgt = p_atmosHahgtSpline.Evaluate(incidence, NumericalApproximation::Extrapolate);
275  gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt;
276  hahgt = p_atmosHahgtSpline.Evaluate(p_emission, NumericalApproximation::Extrapolate);
277  gmu = p_gammax * xmu + p_gammay * ymu + hahgt;
278  }
279 
280  // purely atmos term uses x and y (plus single-particle phase
281  // function correction)
282  phasefn = (1.0 - p_atmosHga * p_atmosHga) / pow(1.0 + 2.0 * p_atmosHga *
283  cos((PI / 180.0) * phase) + p_atmosHga * p_atmosHga, 1.5);
284  p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * ((xmunot * xmu - ymunot * ymu) +
285  (phasefn - 1.0) * (1.0 - emu * p_emunot));
286 
287  // xmitted surface term uses gammas
288  p_trans = gmunot * gmu;
289 
290  // finally, never-scattered term is given by pure attenuation, with
291  // a correction for highly forward-scattered light (on the way down
292  // but not on the way up) as tabulated in hahgt0Table
293  if (p_atmosEstTau) {
294  NumericalAtmosApprox::IntegFunc sub;
295  sub = NumericalAtmosApprox::OuterFunction;
296  NumericalAtmosApprox qromb;
297  p_atmosAtmSwitch = 3;
298  qromb.Reset();
299  p_atmosInc = incidence;
300  p_atmosMunot = cos((PI / 180.0) * incidence);
301  p_atmosSini = sin((PI / 180.0) * incidence);
302  hahgt0 = qromb.RombergsMethod(this, sub, 0, 180);
303  hahgt0 = hahgt0 * AtmosWha() * p_atmosMunot / (360.0 * p_atmosSini);
304  } else {
305  hahgt0 = p_atmosHahgt0Spline.Evaluate(incidence, NumericalApproximation::Extrapolate);
306  }
307  p_trans0 = (p_emunot + hahgt0) * emu;
308 
309  // Calculate the transmission of light that must be subtracted from a shadow. This
310  // includes direct flux and the scattered flux in the upsun half of the sky
311  // downwelling onto the surface, and the usual transmission upward.
312  if (p_atmosEstTau) {
313  NumericalAtmosApprox::IntegFunc sub;
314  sub = NumericalAtmosApprox::OuterFunction;
315  NumericalAtmosApprox qromb;
316  p_atmosAtmSwitch = 1;
317  qromb.Reset();
318  hahgt = qromb.RombergsMethod(this, sub, 90, 180);
319  hahgt = .5 * (p_gammax * xmunot + p_gammay * ymunot - p_emunot) + hahgt *
320  AtmosWha() / 360.0;
321  } else {
322  hahgt = p_atmosHahgtSpline.Evaluate(incidence, NumericalApproximation::Extrapolate);
323  }
324  p_transs = (p_emunot + hahgt) * emu;
325 
326  }
327 }
328 
329 extern "C" Isis::AtmosModel *HapkeAtm2Plugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel) {
330  return new Isis::HapkeAtm2(pvl, pmodel);
331 }