13 Isotropic2::Isotropic2(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
65 double f1munot, f1mmunot, f1mmu, f1mu;
66 double xmunot, ymunot;
71 if(p_atmosTau == 0.0) {
82 p_wha2 = 0.5 * p_atmosWha;
100 p_f1m = log(2.0) - p_em * p_e1 + p_e1_2;
101 p_f2m = -1.0 * (p_f1m + p_em * p_e2 - 1.0);
102 p_f3m = -1.0 * (p_f2m + p_em * p_e3 - 0.5);
103 p_g12 = (p_atmosTau * p_e1 * p_e2 + p_f1m + p_f2m) * 0.5;
104 p_g13 = (p_atmosTau * p_e1 * p_e3 + p_f1m + p_f3m) * (1.0 / 3.0);
118 p_f1 = Eulgam() + log(p_atmosTau) + p_e * p_e1;
119 p_f2 = p_f1 + p_e * p_e2 - 1.0;
120 p_f3 = p_f2 + p_e * p_e3 - 0.5;
122 p_g12p = (p_atmosTau * (p_e1 - p_g11p) + p_em * (p_f1 + p_f2)) * 0.25;
123 p_g13p = (p_atmosTau * (0.5 * p_e1 - p_g12p) + p_em * (p_f1 + p_f3)) * 0.2;
126 p_x0 = p_wha2 * (1.0 + p_wha2 * p_g12);
127 p_y0 = p_wha2 * (p_e2 + p_wha2 * p_g12p);
130 p_delta = (1.0 - (p_x0 + p_y0) - (1.0 - p_atmosWha) / (1.0 - (p_x0 - p_y0))) / (p_atmosWha * (0.5 - p_e3));
133 p_alpha0 = 1.0 + p_wha2 * p_g12 + p_delta * (0.5 - p_e3);
134 p_alpha1 = 0.5 + p_wha2 * p_g13 + p_delta * ((1.0 / 3.0) - p_e4);
135 p_beta0 = p_e2 + p_wha2 * p_g12p + p_delta * (0.5 - p_e3);
136 p_beta1 = p_e3 + p_wha2 * p_g13p + p_delta * ((1.0 / 3.0) - p_e4);
139 if(p_atmosWha == 1.0) {
141 p_f4m = -1.0 * (p_f3m + p_em * p_e4 - (1.0 / 3.0));
142 p_g14 = (p_atmosTau * p_e1 * p_e4 + p_f1m + p_f4m) * 0.25;
143 p_f4 = p_f3 + p_e * p_e4 - (1.0 / 3.0);
144 p_g14p = (p_atmosTau * (0.5 * p_e1 - p_g13p) + p_em * (p_f1 + p_f4)) * (1.0 / 6.0);
145 p_alpha2 = (1.0 / 3.0) + p_wha2 * p_g14 + p_delta * (0.25 - p_e5);
146 p_beta2 = p_e4 + p_wha2 * p_g14p + p_delta * (0.25 - p_e5);
147 p_fixcon = (p_beta0 * p_atmosTau - p_alpha1 + p_beta1) / ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
151 p_gammax = p_wha2 * p_beta0;
152 p_gammay = 1.0 - p_wha2 * p_alpha0;
155 p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1);
157 SetOldTau(p_atmosTau);
158 SetOldWha(p_atmosWha);
163 munot = cos((
PI / 180.0) * incidence);
164 maxval = max(1.0e-30, hpsq1 + munot * munot);
166 munotp = max(munotp, p_atmosTau / 69.0);
167 mu = cos((
PI / 180.0) * emission);
168 maxval = max(1.0e-30, hpsq1 + mu * mu);
170 mup = max(mup, p_atmosTau / 69.0);
173 xx = -p_atmosTau / max(munotp, 1.0e-30);
181 emunot = exp(-p_atmosTau / munotp);
184 xx = -p_atmosTau / max(mup, 1.0e-30);
192 emu = exp(-p_atmosTau / mup);
197 if(fabs(xx - 1.0) < 1.0e-10) {
199 f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * emunot +
203 f1munot = xx * (log(xx / (1.0 - xx)) + p_e1 / emunot +
205 f1mmunot = xx * (log(1.0 + 1.0 / xx) - p_e1 * emunot +
209 std::string msg =
"Negative length of planetary curvature ";
210 msg +=
"encountered";
215 if(fabs(xx - 1.0) < 1.0e-10) {
217 f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu +
AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
220 f1mu = xx * (log(xx / (1.0 - xx)) + p_e1 / emu +
AtmosModel::Ei(p_atmosTau * (1.0 / xx - 1.0)));
221 f1mmu = xx * (log(1.0 + 1.0 / xx) - p_e1 * emu +
AtmosModel::En(1, p_atmosTau * (1.0 + 1.0 / xx)));
224 std::string msg =
"Negative length of planetary curvature encountered";
228 xmunot = 1.0 + p_wha2 * f1mmunot + p_delta * munotp * (1.0 - emunot);
229 ymunot = emunot * (1.0 + p_wha2 * f1munot) + p_delta * munotp * (1.0 - emunot);
230 xmu = 1.0 + p_wha2 * f1mmu + p_delta * mup * (1.0 - emu);
231 ymu = emu * (1.0 + p_wha2 * f1mu) + p_delta * mup * (1.0 - emu);
234 if(p_atmosWha == 1.0) {
235 fix = p_fixcon * munotp * (xmunot + ymunot);
236 xmunot = xmunot + fix;
237 ymunot = ymunot + fix;
238 fix = p_fixcon * mup * (xmu + ymu);
244 gmunot = p_gammax * xmunot + p_gammay * ymunot;
245 gmu = p_gammax * xmu + p_gammay * ymu;
248 p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * (xmunot * xmu - ymunot * ymu);
257 p_transs = (emunot + 0.5 * (p_gammax * xmunot + p_gammay * ymunot - emunot)) * emu;
static double G11Prime(double tau)
Perform Chandra and Van de Hulst's series approximation for the g'11 function needed in second order ...
double p_atmosHnorm
Atmospheric shell thickness normalized to planet radius.
virtual void AtmosModelAlgorithm(double phase, double incidence, double emission)
Isotropic atmospheric scattering in the first approximation The model for scattering for a general...
double p_trans0
Transmission of surface reflected light through the atmosphere with no scatterings in the atmosphere...
const double PI
The mathematical constant PI.
bool TauOrWhaChanged() const
Checks whether tau or wha have changed.
double p_trans
Transmission of surface reflected light through the atmosphere overall.
Isotropic atmos scattering model.
double p_sbar
Illumination of the ground by the sky.
static double En(unsigned int n, double x)
This routine evaluates the generalized exponential integral, En(x).
#define _FILEINFO_
Macro for the filename and line number.
A type of error that cannot be classified as any of the other error types.
Container for cube-like labels.
static double Ei(double x)
This routine computes the exponential integral, Ei(x).
double p_transs
Transmission of light that must be subtracted from the flat surface model to get the shadow model...
double p_pstd
Pure atmospheric-scattering term.
Namespace for ISIS/Bullet specific routines.