5 #include "NumericalApproximation.h"
14 HapkeAtm1::HapkeAtm1(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
61 double xmunot, ymunot;
71 if(p_atmosTau == 0.0) {
82 p_wha2 = 0.5 * p_atmosWha;
92 p_delta = (1.0 - (p_x0 + p_y0) - (1.0 - p_atmosWha) / (1.0 - (p_x0 - p_y0))) / (p_atmosWha * (0.5 - p_e3));
95 p_alpha0 = 1.0 + p_delta * (0.5 - p_e3);
96 p_alpha1 = 0.5 + p_delta * ((1.0 / 3.0) - p_e4);
97 p_beta0 = p_e2 + p_delta * (0.5 - p_e3);
98 p_beta1 = p_e3 + p_delta * ((1.0 / 3.0) - p_e4);
101 if(p_atmosWha == 1.0) {
103 p_alpha2 = (1.0 / 3.0) + p_delta * (0.25 - p_e5);
104 p_beta2 = p_e4 + p_delta * (0.25 - p_e5);
105 p_fixcon = (p_beta0 * p_atmosTau - p_alpha1 + p_beta1) /
106 ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
111 p_gammax = p_wha2 * p_beta0;
112 p_gammay = 1.0 - p_wha2 * p_alpha0;
122 p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1) + p_atmosHahgsb;
124 SetOldTau(p_atmosTau);
125 SetOldWha(p_atmosWha);
131 if(incidence == 90.0) {
135 munot = cos((
PI / 180.0) * incidence);
138 maxval = max(1.0e-30, hpsq1 + munot * munot);
140 munotp = max(munotp, p_atmosTau / 69.0);
142 if(emission == 90.0) {
146 mu = cos((
PI / 180.0) * emission);
149 maxval = max(1.0e-30, hpsq1 + mu * mu);
151 mup = max(mup, p_atmosTau / 69.0);
153 maxval = max(1.0e-30, munotp);
154 xx = -p_atmosTau / maxval;
163 emunot = exp(-p_atmosTau / munotp);
166 maxval = max(1.0e-30, mup);
167 xx = -p_atmosTau / maxval;
176 emu = exp(-p_atmosTau / mup);
179 xmunot = 1.0 + p_delta * munotp * (1.0 - emunot);
180 ymunot = emunot + p_delta * munotp * (1.0 - emunot);
181 xmu = 1.0 + p_delta * mup * (1.0 - emu);
182 ymu = emu + p_delta * mup * (1.0 - emu);
185 if(p_atmosWha == 1.0) {
186 fix = p_fixcon * munotp * (xmunot + ymunot);
187 xmunot = xmunot + fix;
188 ymunot = ymunot + fix;
189 fix = p_fixcon * mup * (xmu + ymu);
200 p_atmosAtmSwitch = 1;
202 p_atmosInc = incidence;
203 p_atmosMunot = cos((
PI / 180.0) * incidence);
204 p_atmosSini = sin((
PI / 180.0) * incidence);
206 gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt *
AtmosWha() / 360.0;
207 p_atmosInc = emission;
208 p_atmosMunot = cos((
PI / 180.0) * emission);
209 p_atmosSini = sin((
PI / 180.0) * emission);
211 gmu = p_gammax * xmu + p_gammay * ymu + hahgt *
AtmosWha() / 360.0;
214 gmunot = p_gammax * xmunot + p_gammay * ymunot + hahgt;
216 gmu = p_gammax * xmu + p_gammay * ymu + hahgt;
222 phasefn = (1.0 - p_atmosHga * p_atmosHga) / pow(1.0 + 2.0 * p_atmosHga * 0.0 + p_atmosHga * p_atmosHga, 1.5);
225 phasefn = (1.0 - p_atmosHga * p_atmosHga) /
226 pow(1.0 + 2.0 * p_atmosHga * cos((
PI / 180.0) * phase) + p_atmosHga * p_atmosHga, 1.5);
229 p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) *
230 ((xmunot * xmu - ymunot * ymu) + (phasefn - 1.0) * (1.0 - emu * emunot));
242 p_atmosAtmSwitch = 3;
244 p_atmosInc = incidence;
245 p_atmosMunot = cos((
PI / 180.0) * incidence);
246 p_atmosSini = sin((
PI / 180.0) * incidence);
248 hahgt0 = hahgt0 *
AtmosWha() * p_atmosMunot / (360.0 * p_atmosSini);
261 p_atmosAtmSwitch = 1;
264 hahgt = .5 * (p_gammax * xmunot + p_gammay * ymunot - emunot) + hahgt *