Isis 3 Programmer Reference
Isotropic1.cpp
1 #include <cmath>
2 #include "AtmosModel.h"
3 #include "Constants.h"
4 #include "Isotropic1.h"
5 #include "Pvl.h"
6 #include "PvlGroup.h"
7 #include "IException.h"
8 #include "IString.h"
9 
10 using std::max;
11 
12 namespace Isis {
13  Isotropic1::Isotropic1(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
14  }
15 
53  void Isotropic1::AtmosModelAlgorithm(double phase, double incidence, double emission) {
54  double hpsq1;
55  double mu, munot;
56  double mup, munotp;
57  double xx;
58  double maxval;
59  double emu, emunot;
60  double xmu, xmunot;
61  double ymu, ymunot;
62  double fix;
63  double gmunot, gmu;
64 
65  if(p_atmosTau == 0.0) {
66  p_pstd = 0.0;
67  p_trans = 1.0;
68  p_trans0 = 1.0;
69  p_sbar = 0.0;
70  p_transs = 1.0;
71  return;
72  }
73 
74  if(TauOrWhaChanged()) {
75  // preparation includes exponential integrals e sub 2 through 4
76  p_wha2 = 0.5 * p_atmosWha;
77  p_e2 = AtmosModel::En(2, p_atmosTau);
78  p_e3 = AtmosModel::En(3, p_atmosTau);
79  p_e4 = AtmosModel::En(4, p_atmosTau);
80 
81  // zeroth moments of (uncorrected) x and y times characteristic fn
82  p_x0 = p_wha2;
83  p_y0 = p_wha2 * p_e2;
84 
85  // higher-order correction term for x and y
86  p_delta = (1.0 - (p_x0 + p_y0) - (1.0 - p_atmosWha) / (1.0 - (p_x0 - p_y0))) / (p_atmosWha * (0.5 - p_e3));
87 
88  // moments of (corrected) x and y
89  p_alpha0 = 1.0 + p_delta * (0.5 - p_e3);
90  p_alpha1 = 0.5 + p_delta * ((1.0 / 3.0) - p_e4);
91  p_beta0 = p_e2 + p_delta * (0.5 - p_e3);
92  p_beta1 = p_e3 + p_delta * ((1.0 / 3.0) - p_e4);
93 
94  // prepare to find correct mixture of x and y in conservative case
95  if(p_atmosWha == 1.0) {
96  p_e5 = AtmosModel::En(5, p_atmosTau);
97  p_alpha2 = (1.0 / 3.0) + p_delta * (0.25 - p_e5);
98  p_beta2 = p_e4 + p_delta * (0.25 - p_e5);
99  p_fixcon = (p_beta0 * p_atmosTau - p_alpha1 + p_beta1) / ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
100  }
101 
102  // gamma will be a weighted sum of x and y functions
103  p_gammax = p_wha2 * p_beta0;
104  p_gammay = 1.0 - p_wha2 * p_alpha0;
105 
106  // sbar is total diffuse illumination and comes from moments
107  p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1);
108 
109  SetOldTau(p_atmosTau);
110  SetOldWha(p_atmosWha);
111  }
112 
113  // correct the path lengths for planetary curvature
114  hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
115  munot = cos((PI / 180.0) * incidence);
116  maxval = max(1.0e-30, hpsq1 + munot * munot);
117  munotp = p_atmosHnorm / (sqrt(maxval) - munot);
118  munotp = max(munotp, p_atmosTau / 69.0);
119  mu = cos((PI / 180.0) * emission);
120  maxval = max(1.0e-30, hpsq1 + mu * mu);
121  mup = p_atmosHnorm / (sqrt(maxval) - mu);
122  mup = max(mup, p_atmosTau / 69.0);
123 
124  // build the x and y functions of mu0 and mu
125  maxval = max(1.0e-30, munotp);
126  xx = -p_atmosTau / maxval;
127  if(xx < -69.0) {
128  emunot = 0.0;
129  }
130  else if(xx > 69.0) {
131  emunot = 1.0e30;
132  }
133  else {
134  emunot = exp(-p_atmosTau / munotp);
135  }
136 
137  maxval = max(1.0e-30, mup);
138  xx = -p_atmosTau / maxval;
139 
140  if(xx < -69.0) {
141  emu = 0.0;
142  }
143  else if(xx > 69.0) {
144  emu = 1.0e30;
145  }
146  else {
147  emu = exp(-p_atmosTau / mup);
148  }
149 
150  xmunot = 1.0 + p_delta * munotp * (1.0 - emunot);
151  ymunot = emunot + p_delta * munotp * (1.0 - emunot);
152  xmu = 1.0 + p_delta * mup * (1.0 - emu);
153  ymu = emu + p_delta * mup * (1.0 - emu);
154 
155  // mix the x and y as required in the conservative case
156  if(p_atmosWha == 1.0) {
157  fix = p_fixcon * munotp * (xmunot + ymunot);
158  xmunot = xmunot + fix;
159  ymunot = ymunot + fix;
160  fix = p_fixcon * mup * (xmu + ymu);
161  xmu = xmu + fix;
162  ymu = ymu + fix;
163  }
164 
165  // gamma1 functions come from x and y
166  gmunot = p_gammax * xmunot + p_gammay * ymunot;
167  gmu = p_gammax * xmu + p_gammay * ymu;
168 
169  // purely atmos term uses x and y, xmitted surface term uses gammas
170  p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * (xmunot * xmu - ymunot * ymu);
171  p_trans = gmunot * gmu;
172 
173  // finally, never-scattered term is given by pure attenuation
174  p_trans0 = emunot * emu;
175 
176  // Calculate the transmission of light that must be subtracted from a shadow.
177  // This includes direct flux and the scattered flux in the upsun half of the
178  // sky downwelling onto the surface, and the usual transmission upward.
179  p_transs = (emunot + 0.5 * (p_gammax * xmunot + p_gammay * ymunot - emunot)) * emu;
180  }
181 }
182 
183 extern "C" Isis::AtmosModel *Isotropic1Plugin(Isis::Pvl &pvl,
184  Isis::PhotoModel &pmodel) {
185  return new Isis::Isotropic1(pvl, pmodel);
186 }
const double PI
The mathematical constant PI.
Definition: Constants.h:56
Isotropic atmos scattering model.
Definition: AtmosModel.h:76
Container for cube-like labels.
Definition: Pvl.h:135
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31