Official websites use .gov
A .gov website belongs to an official government organization in the United States.

Secure .gov websites use HTTPS
A lock ( ) or https:// means you’ve safely connected to the .gov website. Share sensitive information only on official, secure websites.

Isis 3 Programmer Reference
Isotropic1.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include <cmath>
8#include "AtmosModel.h"
9#include "Constants.h"
10#include "Isotropic1.h"
11#include "Pvl.h"
12#include "PvlGroup.h"
13#include "IException.h"
14#include "IString.h"
15
16using std::max;
17
18namespace Isis {
19 Isotropic1::Isotropic1(Pvl &pvl, PhotoModel &pmodel) : AtmosModel(pvl, pmodel) {
20 }
21
59 void Isotropic1::AtmosModelAlgorithm(double phase, double incidence, double emission) {
60 double hpsq1;
61 double mu, munot;
62 double mup, munotp;
63 double xx;
64 double maxval;
65 double emu, emunot;
66 double xmu, xmunot;
67 double ymu, ymunot;
68 double fix;
69 double gmunot, gmu;
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(TauOrWhaChanged()) {
81 // preparation includes exponential integrals e sub 2 through 4
82 p_wha2 = 0.5 * p_atmosWha;
83 p_e2 = AtmosModel::En(2, p_atmosTau);
84 p_e3 = AtmosModel::En(3, p_atmosTau);
85 p_e4 = AtmosModel::En(4, p_atmosTau);
86
87 // zeroth moments of (uncorrected) x and y times characteristic fn
88 p_x0 = p_wha2;
89 p_y0 = p_wha2 * p_e2;
90
91 // higher-order correction term for x and y
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));
93
94 // moments of (corrected) x and y
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);
99
100 // prepare to find correct mixture of x and y in conservative case
101 if(p_atmosWha == 1.0) {
102 p_e5 = AtmosModel::En(5, p_atmosTau);
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) / ((p_alpha1 + p_beta1) * p_atmosTau + 2.0 * (p_alpha2 + p_beta2));
106 }
107
108 // gamma will be a weighted sum of x and y functions
109 p_gammax = p_wha2 * p_beta0;
110 p_gammay = 1.0 - p_wha2 * p_alpha0;
111
112 // sbar is total diffuse illumination and comes from moments
113 p_sbar = 1.0 - ((2.0 - p_atmosWha * p_alpha0) * p_alpha1 + p_atmosWha * p_beta0 * p_beta1);
114
115 SetOldTau(p_atmosTau);
116 SetOldWha(p_atmosWha);
117 }
118
119 // correct the path lengths for planetary curvature
120 hpsq1 = pow((1.0 + p_atmosHnorm), 2.0) - 1.0;
121 munot = cos((PI / 180.0) * incidence);
122 maxval = max(1.0e-30, hpsq1 + munot * munot);
123 munotp = p_atmosHnorm / (sqrt(maxval) - munot);
124 munotp = max(munotp, p_atmosTau / 69.0);
125 mu = cos((PI / 180.0) * emission);
126 maxval = max(1.0e-30, hpsq1 + mu * mu);
127 mup = p_atmosHnorm / (sqrt(maxval) - mu);
128 mup = max(mup, p_atmosTau / 69.0);
129
130 // build the x and y functions of mu0 and mu
131 maxval = max(1.0e-30, munotp);
132 xx = -p_atmosTau / maxval;
133 if(xx < -69.0) {
134 emunot = 0.0;
135 }
136 else if(xx > 69.0) {
137 emunot = 1.0e30;
138 }
139 else {
140 emunot = exp(-p_atmosTau / munotp);
141 }
142
143 maxval = max(1.0e-30, mup);
144 xx = -p_atmosTau / maxval;
145
146 if(xx < -69.0) {
147 emu = 0.0;
148 }
149 else if(xx > 69.0) {
150 emu = 1.0e30;
151 }
152 else {
153 emu = exp(-p_atmosTau / mup);
154 }
155
156 xmunot = 1.0 + p_delta * munotp * (1.0 - emunot);
157 ymunot = emunot + p_delta * munotp * (1.0 - emunot);
158 xmu = 1.0 + p_delta * mup * (1.0 - emu);
159 ymu = emu + p_delta * mup * (1.0 - emu);
160
161 // mix the x and y as required in the conservative case
162 if(p_atmosWha == 1.0) {
163 fix = p_fixcon * munotp * (xmunot + ymunot);
164 xmunot = xmunot + fix;
165 ymunot = ymunot + fix;
166 fix = p_fixcon * mup * (xmu + ymu);
167 xmu = xmu + fix;
168 ymu = ymu + fix;
169 }
170
171 // gamma1 functions come from x and y
172 gmunot = p_gammax * xmunot + p_gammay * ymunot;
173 gmu = p_gammax * xmu + p_gammay * ymu;
174
175 // purely atmos term uses x and y, xmitted surface term uses gammas
176 p_pstd = 0.25 * p_atmosWha * munotp / (munotp + mup) * (xmunot * xmu - ymunot * ymu);
177 p_trans = gmunot * gmu;
178
179 // finally, never-scattered term is given by pure attenuation
180 p_trans0 = emunot * emu;
181
182 // Calculate the transmission of light that must be subtracted from a shadow.
183 // This includes direct flux and the scattered flux in the upsun half of the
184 // sky downwelling onto the surface, and the usual transmission upward.
185 p_transs = (emunot + 0.5 * (p_gammax * xmunot + p_gammay * ymunot - emunot)) * emu;
186 }
187}
188
189extern "C" Isis::AtmosModel *Isotropic1Plugin(Isis::Pvl &pvl,
190 Isis::PhotoModel &pmodel) {
191 return new Isis::Isotropic1(pvl, pmodel);
192}
Isotropic atmos scattering model.
Definition AtmosModel.h:60
double p_trans
Transmission of surface reflected light through the atmosphere overall.
Definition AtmosModel.h:259
static double En(unsigned int n, double x)
This routine evaluates the generalized exponential integral, En(x).
double p_atmosHnorm
Atmospheric shell thickness normalized to planet radius.
Definition AtmosModel.h:267
double p_sbar
Illumination of the ground by the sky.
Definition AtmosModel.h:262
double p_transs
Transmission of light that must be subtracted from the flat surface model to get the shadow model.
Definition AtmosModel.h:261
double p_pstd
Pure atmospheric-scattering term.
Definition AtmosModel.h:258
double p_trans0
Transmission of surface reflected light through the atmosphere with no scatterings in the atmosphere.
Definition AtmosModel.h:260
bool TauOrWhaChanged() const
Checks whether tau or wha have changed.
virtual void AtmosModelAlgorithm(double phase, double incidence, double emission)
Isotropic atmospheric scattering in the first approximation The model for scattering for a general,...
Container for cube-like labels.
Definition Pvl.h:119
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
const double PI
The mathematical constant PI.
Definition Constants.h:40