USGS

Isis 3.0 Object Programmers' Reference

 Home
Anisotropic1.cpp
1 #include <cmath>
2 #include "Anisotropic1.h"
3 #include "AtmosModel.h"
4 #include "IException.h"
5 #include "IString.h"
6 
7 using std::max;
8 
9 namespace Isis {
10  Anisotropic1::Anisotropic1(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
11 
12  // Set default values
13  p_atmosAlpha0_0 = 0.0;
14  p_atmosAlpha1_0 = 0.0;
15  p_atmosBeta0_0 = 0.0;
16  p_atmosBeta1_0 = 0.0;
17  p_atmosDelta_0 = 0.0;
18  p_atmosDelta_1 = 0.0;
19  p_atmosDen = 0.0;
20  p_atmosE2 = 0.0;
21  p_atmosE3 = 0.0;
22  p_atmosE4 = 0.0;
23  p_atmosE5 = 0.0;
24  p_atmosFac = 0.0;
25  p_atmosP0 = 0.0;
26  p_atmosP1 = 0.0;
27  p_atmosQ0 = 0.0;
28  p_atmosQ02p02 = 0.0;
29  p_atmosQ1 = 0.0;
30  p_atmosQ12p12 = 0.0;
31  p_atmosWha2 = 0.0;
32  p_atmosWham = 0.0;
33  p_atmosX0_0 = 0.0;
34  p_atmosX0_1 = 0.0;
35  p_atmosY0_0 = 0.0;
36  p_atmosY0_1 = 0.0;
37 
38  }
39 
66  void Anisotropic1::AtmosModelAlgorithm(double phase, double incidence, double emission) {
67  double hpsq1;
68  double munot, munotp;
69  double maxval;
70  double mu, mup;
71  double xx;
72  double emunot, emu;
73  double gmunot, gmu;
74  double sum, prod;
75  double cosazss;
76  double xmunot_0, ymunot_0;
77  double xmu_0, ymu_0;
78  double xmunot_1, ymunot_1;
79  double xmu_1, ymu_1;
80  double cxx, cyy;
81  double xystuff;
82 
83  if(p_atmosBha == 0.0) {
84  p_atmosBha = 1.0e-6;
85  }
86 
87  if(p_atmosTau == 0.0) {
88  p_pstd = 0.0;
89  p_trans = 1.0;
90  p_trans0 = 1.0;
91  p_sbar = 0.0;
92  p_transs = 1.0;
93  return;
94  }
95 
96  if(p_atmosWha == 1.0) {
97  QString msg = "Anisotropic conservative case not implemented yet - WHA parameter cannot be set to 1.0";
98  msg += "This will cause negative planetary curvature to occur.";
99  throw IException(IException::User, msg, _FILEINFO_);
100  }
101 
102  if(TauOrWhaChanged()) {
103  // preparation includes exponential integrals e sub 2 through 5
104  p_atmosWha2 = 0.5 * p_atmosWha;
105  p_atmosWham = 1.0 - p_atmosWha;
106  p_atmosE2 = AtmosModel::En(2, p_atmosTau);
107  p_atmosE3 = AtmosModel::En(3, p_atmosTau);
108  p_atmosE4 = AtmosModel::En(4, p_atmosTau);
109  p_atmosE5 = AtmosModel::En(5, p_atmosTau);
110  // first, get the required quantities for the axisymmetric m=0 part
111  // zeroth moments of (uncorrected) x and y times characteristic fn
112  p_atmosX0_0 = p_atmosWha2 * (1.0 + (1.0 / 3.0) * p_atmosBha *
113  p_atmosWham);
114  p_atmosY0_0 = p_atmosWha2 * (p_atmosE2 + p_atmosBha * p_atmosWham *
115  p_atmosE4);
116  // higher-order correction term for x and y
117  p_atmosDelta_0 = (1.0 - (p_atmosX0_0 + p_atmosY0_0) - (1.0 - p_atmosWha *
118  (1.0 + (1.0 / 3.0) * p_atmosBha * p_atmosWham)) / (1.0 -
119  (p_atmosX0_0 - p_atmosY0_0))) / (p_atmosWha * (0.5 - p_atmosE3 +
120  p_atmosBha * p_atmosWham * (0.25 - p_atmosE5)));
121 
122  // moments of (corrected) x and y
123  p_atmosAlpha0_0 = 1.0 + p_atmosDelta_0 * (0.5 - p_atmosE3);
124  p_atmosAlpha1_0 = 0.5 + p_atmosDelta_0 * ((1.0 / 3.0) - p_atmosE4);
125  p_atmosBeta0_0 = p_atmosE2 + p_atmosDelta_0 * (0.5 - p_atmosE3);
126  p_atmosBeta1_0 = p_atmosE3 + p_atmosDelta_0 * ((1.0 / 3.0) - p_atmosE4);
127  // gamma will be a weighted sum of m=0 x and y functions
128  // but unlike before, call the weights q1 and p1, and we also
129  // need additional weights q0 and p0
130  p_atmosFac = 2.0 - p_atmosWha * p_atmosAlpha0_0;
131  p_atmosDen = pow(p_atmosFac, 2.0) - pow(p_atmosWha * p_atmosBeta0_0, 2.0);
132  p_atmosQ0 = p_atmosBha * p_atmosWha * p_atmosWham * (p_atmosFac *
133  p_atmosAlpha1_0 - p_atmosWha * p_atmosBeta0_0 * p_atmosBeta1_0) /
134  p_atmosDen;
135  p_atmosP0 = p_atmosBha * p_atmosWha * p_atmosWham * (-p_atmosFac *
136  p_atmosBeta1_0 - p_atmosWha * p_atmosBeta0_0 * p_atmosAlpha1_0) /
137  p_atmosDen;
138  p_atmosQ02p02 = p_atmosQ0 * p_atmosQ0 - p_atmosP0 * p_atmosP0;
139  p_atmosQ1 = (2.0 * p_atmosWham * p_atmosFac) / p_atmosDen;
140  p_atmosP1 = (2.0 * p_atmosWham * p_atmosWha * p_atmosBeta0_0) /
141  p_atmosDen;
142  p_atmosQ12p12 = p_atmosQ1 * p_atmosQ1 - p_atmosP1 * p_atmosP1;
143  // sbar is total diffuse illumination and comes from moments
144  p_sbar = 1.0 - 2.0 * (p_atmosQ1 * p_atmosAlpha1_0 + p_atmosP1 *
145  p_atmosBeta1_0);
146  // we're not done yet! still have to calculate the m=1 portion
147  // zeroth moments of (uncorrected) x and y times characteristic fn
148  p_atmosX0_1 = 0.5 * p_atmosWha2 * p_atmosBha * (1.0 - (1.0 / 3.0));
149  p_atmosY0_1 = 0.5 * p_atmosWha2 * p_atmosBha * (p_atmosE2 - p_atmosE4);
150  // higher-order correction term for x and y
151  p_atmosDelta_1 = (1.0 - (p_atmosX0_1 + p_atmosY0_1) - (1.0 -
152  (1.0 / 3.0) * p_atmosWha * p_atmosBha) / (1.0 - (p_atmosX0_1 -
153  p_atmosY0_1))) / (p_atmosWha2 * p_atmosBha * ((0.5 - 0.25) -
154  (p_atmosE3 - p_atmosE5)));
155  // moments of (corrected) x and y are not needed for m=1, so we're done
156  SetOldTau(p_atmosTau);
157  SetOldWha(p_atmosWha);
158  }
159 
160  // correct the path lengths for planetary curvature
161  hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
162 
163  if(incidence == 90.0) {
164  munot = 0.0;
165  }
166  else {
167  munot = cos((PI / 180.0) * incidence);
168  }
169 
170  maxval = max(1.0e-30, hpsq1 + munot * munot);
171  munotp = p_atmosHnorm / (sqrt(maxval) - munot);
172  munotp = max(munotp, p_atmosTau / 69.0);
173  if(emission == 90.0) {
174  mu = 0.0;
175  }
176  else {
177  mu = cos((PI / 180.0) * emission);
178  }
179 
180  maxval = max(1.0e-30, hpsq1 + mu * mu);
181  mup = p_atmosHnorm / (sqrt(maxval) - mu);
182  mup = max(mup, p_atmosTau / 69.0);
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  emunot = 0.0;
189  }
190  else if(xx > 69.0) {
191  emunot = 1.0e30;
192  }
193  else {
194  emunot = exp(-p_atmosTau / munotp);
195  }
196 
197  maxval = max(1.0e-30, mup);
198  xx = -p_atmosTau / maxval;
199 
200  if(xx < -69.0) {
201  emu = 0.0;
202  }
203  else if(xx > 69.0) {
204  emu = 1.0e30;
205  }
206  else {
207  emu = exp(-p_atmosTau / mup);
208  }
209 
210  // first for m=0
211  xmunot_0 = 1.0 + p_atmosDelta_0 * munotp * (1.0 - emunot);
212  ymunot_0 = emunot + p_atmosDelta_0 * munotp * (1.0 - emunot);
213  xmu_0 = 1.0 + p_atmosDelta_0 * mup * (1.0 - emu);
214  ymu_0 = emu + p_atmosDelta_0 * mup * (1.0 - emu);
215 
216  // then for m=1
217  xmunot_1 = 1.0 + p_atmosDelta_1 * munotp * (1.0 - emunot);
218  ymunot_1 = emunot + p_atmosDelta_1 * munotp * (1.0 - emunot);
219  xmu_1 = 1.0 + p_atmosDelta_1 * mup * (1.0 - emu);
220  ymu_1 = emu + p_atmosDelta_1 * mup * (1.0 - emu);
221 
222  // gamma1 functions come from x and y with m=0
223  gmunot = p_atmosP1 * xmunot_0 + p_atmosQ1 * ymunot_0;
224  gmu = p_atmosP1 * xmu_0 + p_atmosQ1 * ymu_0;
225 
226  // purely atmos term uses x and y of both orders and is complex
227  sum = munot + mu;
228  prod = munot * mu;
229  cxx = 1.0 - p_atmosQ0 * sum + (p_atmosQ02p02 - p_atmosBha *
230  p_atmosQ12p12) * prod;
231  cyy = 1.0 + p_atmosQ0 * sum + (p_atmosQ02p02 - p_atmosBha *
232  p_atmosQ12p12) * prod;
233 
234  if(phase == 90.0) {
235  cosazss = 0.0 - munot * mu;
236  }
237  else {
238  cosazss = cos((PI / 180.0) * phase) - munot * mu;
239  }
240 
241  xystuff = cxx * xmunot_0 * xmu_0 - cyy * ymunot_0 *
242  ymu_0 - p_atmosP0 * sum * (xmu_0 * ymunot_0 + ymu_0 *
243  xmunot_0) + cosazss * p_atmosBha * (xmu_1 *
244  xmunot_1 - ymu_1 * ymunot_1);
245  p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * xystuff;
246 
247  // xmitted surface term uses gammas from m=0
248  p_trans = gmunot * gmu;
249 
250  // finally, never-scattered term is given by pure attenuation
251  p_trans0 = emunot * emu;
252 
253  // Calculate the transmission of light that must be subtracted from a
254  // shadow. This includes direct flux and the scattered flux in the
255  // upsun half of the sky downwelling onto the surface, and the usual
256  // transmission upward. NOTE: We need to derive the analytic expression
257  // for the light from half the sky in the Legendre scattering model. Until
258  // we do so, we are setting the shadow transmission to the purely
259  // unscattered part (same as trans0). This will give a result but is
260  // not fully consistent with how the other scattering models are
261  // implemented.
262  p_transs = p_trans0;
263  }
264 }
265 
266 extern "C" Isis::AtmosModel *Anisotropic1Plugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel) {
267  return new Isis::Anisotropic1(pvl, pmodel);
268 }