USGS

Isis 3.0 Object Programmers' Reference

 Home
Anisotropic2.cpp
1 #include <cmath>
2 #include "Anisotropic2.h"
3 #include "AtmosModel.h"
4 #include "Constants.h"
5 #include "Pvl.h"
6 #include "PvlGroup.h"
7 #include "IException.h"
8 #include "IString.h"
9 
10 using std::min;
11 using std::max;
12 
13 namespace Isis {
14  Anisotropic2::Anisotropic2(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
15  }
16 
45  void Anisotropic2::AtmosModelAlgorithm(double phase, double incidence, double emission) {
46  double xx;
47  double munot, mu;
48  double emunot, emu;
49  double munotp, mup;
50  double gmunot, gmu;
51  double hpsq1;
52  double f1munot, f2munot, f3munot;
53  double f1mmunot, f2mmunot, f3mmunot;
54  double maxval;
55  double f1mu, f2mu, f3mu;
56  double f1mmu, f2mmu, f3mmu;
57  double sum;
58  double prod;
59  double cosazss;
60  double xystuff;
61  double xmunot_0, ymunot_0;
62  double xmu_0, ymu_0;
63  double cxx, cyy;
64  double xmunot_1, ymunot_1;
65  double xmu_1, ymu_1;
66 
67  if(p_atmosBha == 0.0) {
68  p_atmosBha = 1.0e-6;
69  }
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(p_atmosWha == 1.0) {
81  QString msg = "Anisotropic conservative case not implemented yet - WHA parameter cannot be set to 1.0.";
82  msg += "This will cause negative planetary curvature to occur.";
83  throw IException(IException::User, msg, _FILEINFO_);
84  }
85 
86  if(TauOrWhaChanged()) {
87  // preparation includes exponential integrals e sub 2 through 5
88  p_wha2 = 0.5 * p_atmosWha;
89  p_wham = 1.0 - p_atmosWha;
90  p_e1 = AtmosModel::En(1, p_atmosTau);
91  p_e1_2 = AtmosModel::En(1, 2.0 * p_atmosTau);
92  p_e2 = AtmosModel::En(2, p_atmosTau);
93  p_e3 = AtmosModel::En(3, p_atmosTau);
94  p_e4 = AtmosModel::En(4, p_atmosTau);
95  p_e5 = AtmosModel::En(5, p_atmosTau);
96 
97  // chandra's gmn functions require fm and fn at mu=-1
98  xx = -p_atmosTau;
99  if(xx < -69.0) {
100  p_em = 0.0;
101  }
102  else if(xx > 69.0) {
103  p_em = 1.0e30;
104  }
105  else {
106  p_em = exp(xx);
107  }
108 
109  p_f1m = log(2.0) - p_em * p_e1 + p_e1_2;
110  p_f2m = -1.0 * (p_f1m + p_em * p_e2 - 1.0);
111  p_f3m = -1.0 * (p_f2m + p_em * p_e3 - 0.5);
112  p_f4m = -1.0 * (p_f3m + p_em * p_e4 - (1.0 / 3.0));
113  p_g12 = (p_atmosTau * p_e1 * p_e2 + p_f1m + p_f2m) * 0.5;
114  p_g13 = (p_atmosTau * p_e1 * p_e3 + p_f1m + p_f3m) * (1.0 / 3.0);
115  p_g14 = (p_atmosTau * p_e1 * p_e4 + p_f1m + p_f4m) * 0.25;
116  p_g32 = (p_atmosTau * p_e3 * p_e2 + p_f3m + p_f2m) * 0.25;
117  p_g33 = (p_atmosTau * p_e3 * p_e3 + p_f3m + p_f3m) * 0.2;
118  p_g34 = (p_atmosTau * p_e3 * p_e4 + p_f3m * p_f4m) * (1.0 / 6.0);
119 
120  // chandra's g'mn functions require g'11 and f at mu=+1
121  xx = p_atmosTau;
122  if(xx < -69.0) {
123  p_e = 0.0;
124  }
125  else if(xx > 69.0) {
126  p_e = 1.0e30;
127  }
128  else {
129  p_e = exp(xx);
130  }
131 
132  p_f1 = Eulgam() + log(p_atmosTau) + p_e * p_e1;
133  p_f2 = p_f1 + p_e * p_e2 - 1.0;
134  p_f3 = p_f2 + p_e * p_e3 - 0.5;
135  p_f4 = p_f3 + p_e * p_e4 - (1.0 / 3.0);
136  p_g11p = AtmosModel::G11Prime(p_atmosTau);
137  p_g12p = (p_atmosTau * (p_e1 - p_g11p) + p_em * (p_f1 + p_f2)) * 0.25;
138  p_g13p = (p_atmosTau * (0.5 * p_e1 - p_g12p) + p_em * (p_f1 + p_f3)) * 0.2;
139  p_g14p = (p_atmosTau * ((1.0 / 3.0) * p_e1 - p_g13p) + p_em * (p_f1 + p_f4)) * (1.0 / 6.0);
140  p_g32p = (p_atmosTau * (p_e1 - p_g13p) + p_em * (p_f3 + p_f2)) * (1.0 / 6.0);
141  p_g33p = (p_atmosTau * (0.5 * p_e1 - p_g32p) + p_em * (p_f3 + p_f3)) * 0.142857;
142  p_g34p = (p_atmosTau * ((1.0 / 3.0) * p_e1 - p_g33p) + p_em * (p_f3 + p_f4)) * 0.125;
143 
144  // first, get the required quantities for the axisymmetric m=0 part
145  // zeroth moments of (uncorrected) x and y times characteristic fn
146  p_x0_0 = p_wha2 * (1.0 + (1.0 / 3.0) * p_atmosBha * p_wham + p_wha2 * (p_g12 + p_atmosBha *
147  p_wham * (p_g14 + p_g32) + pow(p_atmosBha, 2.0) * pow(p_wham, 2.0) * p_g34));
148  p_y0_0 = p_wha2 * (p_e2 + p_atmosBha * p_wham * p_e4 + p_wha2 * (p_g12p
149  + p_atmosBha * p_wham * (p_g14p + p_g32p) +
150  pow(p_atmosBha, 2.0) * pow(p_wham, 2.0) * p_g34p));
151 
152  // higher-order correction term for x and y
153  p_delta_0 = (1.0 - (p_x0_0 + p_y0_0) - (1.0 - p_atmosWha * (1.0 + (1.0 / 3.0) * p_atmosBha *
154  p_wham)) / (1.0 - (p_x0_0 - p_y0_0))) / (p_atmosWha * (0.5 - p_e3 + p_atmosBha * p_wham * (0.25 - p_e5)));
155 
156  // moments of (corrected) x and y
157  p_alpha0_0 = 1.0 + p_wha2 * (p_g12 + p_atmosBha * p_wham * p_g32) + p_delta_0 * (0.5 - p_e3);
158  p_alpha1_0 = 0.5 + p_wha2 * (p_g13 + p_atmosBha * p_wham * p_g33) + p_delta_0 * ((1.0 / 3.0) - p_e4);
159  p_beta0_0 = p_e2 + p_wha2 * (p_g12p + p_atmosBha * p_wham * p_g32p) + p_delta_0 * (0.5 - p_e3);
160  p_beta1_0 = p_e3 + p_wha2 * (p_g13p + p_atmosBha * p_wham * p_g33p) + p_delta_0 * ((1.0 / 3.0) - p_e4);
161 
162  // gamma will be a weighted sum of m=0 x and y functions
163  // but unlike before, call the weights q1 and p1, and we also
164  // need additional weights q0 and p0
165  p_fac = 2.0 - p_atmosWha * p_alpha0_0;
166  p_den = pow(p_fac, 2.0) - pow((p_atmosWha * p_beta0_0), 2.0);
167  p_q0 = p_atmosBha * p_atmosWha * p_wham * (p_fac * p_alpha1_0 - p_atmosWha * p_beta0_0 * p_beta1_0) / p_den;
168  p_p0 = p_atmosBha * p_atmosWha * p_wham * (-p_fac * p_beta1_0 - p_atmosWha * p_beta0_0 * p_alpha1_0) / p_den;
169  p_q02p02 = p_q0 * p_q0 - p_p0 * p_p0;
170  p_q1 = (2.0 * p_wham * p_fac) / p_den;
171  p_p1 = (2.0 * p_wham * p_atmosWha * p_beta0_0) / p_den;
172  p_q12p12 = p_q1 * p_q1 - p_p1 * p_p1;
173 
174  // sbar is total diffuse illumination and comes from moments
175  p_sbar = 1.0 - 2.0 * (p_q1 * p_alpha1_0 + p_p1 * p_beta1_0);
176 
177  // we're not done yet! still have to calculate the m=1 portion
178  // zeroth moments of (uncorrected) x and y times characteristic fn
179  p_x0_1 = 0.5 * p_wha2 * p_atmosBha * (1.0 - (1.0 / 3.0) + 0.5 * p_wha2 * p_atmosBha *
180  (p_g12 - (p_g14 + p_g32) + p_g34));
181  p_y0_1 = 0.5 * p_wha2 * p_atmosBha * (p_e2 - p_e4 + 0.5 * p_wha2 * p_atmosBha *
182  (p_g12p - (p_g14p + p_g32p) + p_g34p));
183 
184  // higher-order correction term for x and y
185  p_delta_1 = (1.0 - (p_x0_1 + p_y0_1) - (1.0 - (1.0 / 3.0) * p_atmosWha * p_atmosBha) /
186  (1.0 - (p_x0_1 - p_y0_1))) / (p_wha2 * p_atmosBha *
187  ((0.5 - 0.25) - (p_e3 - p_e5)));
188 
189  // moments of (corrected) x and y are not needed for m=1, so we're done
190  SetOldTau(p_atmosTau);
191  SetOldWha(p_atmosWha);
192  }
193 
194  // correct the path lengths for planetary curvature
195  hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
196 
197  if(incidence == 90.0) {
198  munot = 0.0;
199  }
200  else {
201  munot = cos((PI / 180.0) * incidence);
202  }
203 
204  maxval = max(1.0e-30, hpsq1 + munot * munot);
205  munotp = p_atmosHnorm / (sqrt(maxval) - munot);
206  munotp = max(munotp, p_atmosTau / 69.0);
207 
208  if(emission == 90.0) {
209  mu = 0.0;
210  }
211  else {
212  mu = cos((PI / 180.0) * emission);
213  }
214 
215  maxval = max(1.0e-30, hpsq1 + mu * mu);
216  mup = p_atmosHnorm / (sqrt(maxval) - mu);
217  mup = max(mup, p_atmosTau / 69.0);
218 
219  // build the x and y functions of mu0 and mu
220  maxval = max(1.0e-30, munotp);
221  xx = -p_atmosTau / maxval;
222  if(xx < -69.0) {
223  emunot = 0.0;
224  }
225  else if(xx > 69.0) {
226  emunot = 1.0e30;
227  }
228  else {
229  emunot = exp(-p_atmosTau / munotp);
230  }
231 
232  maxval = max(1.0e-30, mup);
233  xx = -p_atmosTau / maxval;
234  if(xx < -69.0) {
235  emu = 0.0;
236  }
237  else if(xx > 69.0) {
238  emu = 1.0e30;
239  }
240  else {
241  emu = exp(-p_atmosTau / mup);
242  }
243 
244  // in the second approximation the x and y include the f1, f3 functions
245  xx = munotp;
246  if(fabs(xx - 1.0) < 1.0e-10) {
247  f1munot = p_f1;
248  f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * emunot +
249  AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
250  }
251  else if(xx > 0.0) {
252  f1munot = xx * (log(xx / (1.0 - xx)) + p_e1 / emunot +
253  AtmosModel::Ei(p_atmosTau * (1.0 / xx - 1.0)));
254  f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * emunot +
255  AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
256  }
257  else {
258  QString msg = "Negative length of planetary curvature encountered";
259  throw IException(IException::Unknown, msg, _FILEINFO_);
260  }
261 
262  f2munot = munotp * (f1munot + p_e2 / emunot - 1);
263  f2mmunot = -munotp * (f1mmunot + p_e2 * emunot - 1);
264  f3munot = munotp * (f2munot + p_e3 / emunot - 0.5);
265  f3mmunot = -munotp * (f2mmunot + p_e3 * emunot - 0.5);
266  xx = mup;
267 
268  if(fabs(xx - 1.0) < 1.0e-10) {
269  f1mu = p_f1;
270  f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
271  }
272  else if(xx > 0.0) {
273  f1mu = xx * (log(xx / (1.0 - xx)) + p_e1 / emu + AtmosModel::Ei(p_atmosTau * (1.0 / xx - 1.0)));
274  f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu + AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
275  }
276  else {
277  QString msg = "Negative length of planetary curvature encountered";
278  throw IException(IException::Unknown, msg, _FILEINFO_);
279  }
280 
281  f2mu = mup * (f1mu + p_e2 / emu - 1);
282  f2mmu = -mup * (f1mmu + p_e2 * emu - 1);
283  f3mu = mup * (f2mu + p_e3 / emu - 0.5);
284  f3mmu = -mup * (f2mmu + p_e3 * emu - 0.5);
285 
286  // first for m=0
287  xmunot_0 = 1.0 + p_wha2 * (f1mmunot + p_atmosBha * p_wham * f3mmunot) + p_delta_0 * munotp * (1.0 - emunot);
288  ymunot_0 = emunot * (1.0 + p_wha2 * (f1munot + p_atmosBha * p_wham * f3munot)) +
289  p_delta_0 * munotp * (1.0 - emunot);
290  xmu_0 = 1.0 + p_wha2 * (f1mmu + p_atmosBha * p_wham * f3mmu) + p_delta_0 * mup * (1.0 - emu);
291  ymu_0 = emu * (1.0 + p_wha2 * (f1mu + p_atmosBha * p_wham * f3mu)) + p_delta_0 * mup * (1.0 - emu);
292 
293  // then for m=1
294  xmunot_1 = 1.0 + 0.5 * p_wha2 * p_atmosBha * (f1mmunot - f3mmunot) + p_delta_1 * munotp * (1.0 - emunot);
295  ymunot_1 = emunot * (1.0 + 0.5 * p_wha2 * p_atmosBha *
296  (f1munot - f3munot)) + p_delta_1 * munotp * (1.0 - emunot);
297  xmu_1 = 1.0 + 0.5 * p_wha2 * p_atmosBha * (f1mmu - f3mmu) + p_delta_1 * mup * (1.0 - emu);
298  ymu_1 = emu * (1.0 + 0.5 * p_wha2 * p_atmosBha * (f1mu - f3mu)) + p_delta_1 * mup * (1.0 - emu);
299 
300  // gamma1 functions come from x and y with m=0
301  gmunot = p_p1 * xmunot_0 + p_q1 * ymunot_0;
302  gmu = p_p1 * xmu_0 + p_q1 * ymu_0;
303 
304  // purely atmos term uses x and y of both orders and is complex
305  sum = munot + mu;
306  prod = munot * mu;
307  cxx = 1.0 - p_q0 * sum + (p_q02p02 - p_atmosBha * p_q12p12) * prod;
308  cyy = 1.0 + p_q0 * sum + (p_q02p02 - p_atmosBha * p_q12p12) * prod;
309 
310  if(phase == 90.0) {
311  cosazss = 0.0 - munot * mu;
312  }
313  else {
314  cosazss = cos((PI / 180.0) * phase) - munot * mu;
315  }
316 
317  xystuff = cxx * xmunot_0 * xmu_0 - cyy * ymunot_0 *
318  ymu_0 - p_p0 * sum * (xmu_0 * ymunot_0 + ymu_0 * xmunot_0) + cosazss * p_atmosBha * (xmu_1 *
319  xmunot_1 - ymu_1 * ymunot_1);
320  p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * xystuff;
321 
322  // xmitted surface term uses gammas from m=0
323  p_trans = gmunot * gmu;
324 
325  // finally, never-scattered term is given by pure attenuation
326  p_trans0 = emunot * emu;
327 
328  // Calculate the transmission of light that must be subtracted from a
329  // shadow. This includes direct flux and the scattered flux in the
330  // upsun half of the sky downwelling onto the surface, and the usual
331  // transmission upward. NOTE: We need to derive the analytic expression
332  // for the light from half the sky in the Legendre scattering model. Until
333  // we do so, we are setting the shadow transmission to the purely
334  // unscattered part (same as trans0). This will give a result but is
335  // not fully consistent with how the other scattering models are
336  // implemented.
337  p_transs = p_trans0;
338  }
339 }
340 
341 extern "C" Isis::AtmosModel *Anisotropic2Plugin(Isis::Pvl &pvl,
342  Isis::PhotoModel &pmodel) {
343 
344  return new Isis::Anisotropic2(pvl, pmodel);
345 }