File failed to load: https://isis.astrogeology.usgs.gov/9.0.0/Object/assets/jax/output/NativeMML/config.js
Isis 3 Programmer Reference
NumericalAtmosApprox.cpp
1
5
6/* SPDX-License-Identifier: CC0-1.0 */
7#include <cmath>
8#include "AtmosModel.h"
9#include "NumericalAtmosApprox.h"
10#include "NumericalApproximation.h"
11#include "IException.h"
12#include "IString.h"
13
14using namespace std;
15namespace Isis {
52 double NumericalAtmosApprox::RombergsMethod(AtmosModel *am, IntegFunc sub, double a, double b) {
53 // This method was derived from an algorithm in the text
54 // Numerical Recipes in C: The Art of Scientific Computing
55 // Section 4.3 by Flannery, Press, Teukolsky, and Vetterling
56 int maxits = 20; // maximium number of iterations allowed to converge
57 double dss = 0; // error estimate for
58 double h[maxits+1]; // relative stepsizes for trap
59 double trap[maxits+1]; // successive trapeziodal approximations
60 double epsilon; // desired fractional accuracy
61 double epsilon2; // desired fractional accuracy
62 double ss; // result
63
64 epsilon = 1.0e-4;
65 epsilon2 = 1.0e-6;
66
67 h[0] = 1.0;
68 try {
70 for(int i = 0; i < maxits; i++) {
71 // i will determine number of trapezoidal partitions of area
72 // under curve for "integration" using refined trapezoidal rule
73 trap[i] = RefineExtendedTrap(am, sub, a, b, trap[i], i + 1); // validates data here
74 if(i >= 4) {
75 interp.AddData(5, &h[i-4], &trap[i-4]);
77 dss = interp.PolynomialNevilleErrorEstimate()[0];
78 interp.Reset();
79 // we work only until our necessary accuracy is achieved
80 if(fabs(dss) <= epsilon * fabs(ss)) return ss;
81 if(fabs(dss) <= epsilon2) return ss;
82 }
83 trap[i+1] = trap[i];
84 h[i+1] = 0.25 * h[i];
85 // This is a key step: the factor is 0.25d0 even though
86 // the stepsize is decreased by 0.5d0. This makes the
87 // extraplolation a polynomial in h-squared as allowed
88 // by the equation from Numerical Recipes 4.2.1 pg.132,
89 // not just a polynomial in h.
90 }
91 }
92 catch(IException &e) { // catch error from RefineExtendedTrap, Constructor, Evaluate, PolynomialNevilleErrorEstimate
93 throw IException(e,
94 e.errorType(),
95 "NumericalAtmosApprox::RombergsMethod() - Caught the following error: ",
96 _FILEINFO_);
97 }
99 "NumericalAtmosApprox::RombergsMethod() - Failed to converge in "
100 + IString(maxits) + " iterations.",
101 _FILEINFO_);
102 }
103
136 double NumericalAtmosApprox::RefineExtendedTrap(AtmosModel *am, IntegFunc sub, double a, double b, double s, unsigned int n) {
137 // This method was derived from an algorithm in the text
138 // Numerical Recipes in C: The Art of Scientific Computing
139 // Section 4.2 by Flannery, Press, Teukolsky, and Vetterling
140 try {
141 if(n == 1) {
142 double begin;
143 double end;
144 if(sub == InnerFunction) {
145 begin = InrFunc2Bint(am, a);
146 end = InrFunc2Bint(am, b);
147 }
148 else {
149 begin = OutrFunc2Bint(am, a);
150 end = OutrFunc2Bint(am, b);
151 }
152 return (0.5 * (b - a) * (begin + end));
153 }
154 else {
155 int it;
156 double delta, tnm, x, sum;
157 it = (int)(pow(2.0, (double)(n - 2)));
158 tnm = it;
159 delta = (b - a) / tnm; // spacing of the points to be added
160 x = a + 0.5 * delta;
161 sum = 0.0;
162 for(int i = 0; i < it; i++) {
163 if(sub == InnerFunction) {
164 sum = sum + InrFunc2Bint(am, x);
165 }
166 else {
167 sum = sum + OutrFunc2Bint(am, x);
168 }
169 x = x + delta;
170 }
171 return (0.5 * (s + (b - a) * sum / tnm));// replace s with refined value
172 }
173 }
174 catch(IException &e) { // catch exception from Evaluate()
175 throw IException(e,
176 e.errorType(),
177 "NumericalAtmosApprox::RefineExtendedTrap() - Caught the following error: ",
178 _FILEINFO_);
179 }
180 }
181
202 double result;
205
206 am->p_atmosPhi = phi;
207 am->p_atmosCosphi = cos((PI / 180.0) * phi);
208
210 try {
211 result = qromb.RombergsMethod(am, sub, 1.0e-6, 1.0);
212 return result;
213 }
214 catch(IException &e) { // catch exception from RombergsMethod()
215 throw IException(e,
216 e.errorType(),
217 "NumericalAtmosApprox::OutrFunc2Bint() - Caught the following error: ",
218 _FILEINFO_);
219 }
220 }
221
251 double ema; //pass in mu, get emission angle
252 double sine; //sin(ema)
253 double alpha;
254 double phase; //angle between sun and camera
255 double result;
256 double phasefn;//Henyey-Greenstein phase fn
257 double xx; //temp var to calculate emunot, emu
258 double emunot; //exp(-tau/munot)
259 double emu; //exp(-tau/mu)
260 double tfac; //factor that occurs in the integrals for transmission
261
262 // calculate the phase angle
263 // also calculate any of the other redundant parameters
264 ema = acos(mu) * (180.0 / PI);
265 sine = sin(ema * (PI / 180.0));
266 if((am->p_atmosAtmSwitch == 0) || (am->p_atmosAtmSwitch == 2)) { // Reflection phase <= 90
267 alpha = am->p_atmosSini * sine * am->p_atmosCosphi +
268 am->p_atmosMunot * mu;
269 }
270 else { // Transmission phase <= 90
271 alpha = am->p_atmosSini * sine * am->p_atmosCosphi -
272 am->p_atmosMunot * mu;
273 }
274 phase = acos(alpha) * (180.0 / PI);
275 // now evaluate the integrand; all needed parameters
276 // have been hidden separately and passed to it in COMMON.
277 if(am->p_atmosAtmSwitch == 0) {
278 // integrand for hemispheric albedo
279 result = mu * am->p_atmosPM->CalcSurfAlbedo(phase, am->p_atmosInc, ema);
280 }
281 else {
282 phasefn = (1.0 - am->AtmosHga() * am->AtmosHga()) / pow((1.0 + 2.0 * am->AtmosHga() * alpha + am->AtmosHga() * am->AtmosHga()), 1.5);
283 xx = -am->AtmosTau() / max(am->p_atmosMunot, 1.0e-30);
284 if(xx < -69.0) {
285 emunot = 0.0;
286 }
287 else if(xx > 69.0) {
288 emunot = 1.0e30;
289 }
290 else {
291 emunot = exp(xx);
292 }
293 xx = -am->AtmosTau() / max(mu, 1.0e-30);
294
295 if(xx < -69.0) {
296 emu = 0.0;
297 }
298 else if(xx > 69.0) {
299 emu = 1.0e30;
300 }
301 else {
302 emu = exp(xx);
303 }
304 if(mu == am->p_atmosMunot) {
305 tfac = am->AtmosTau() * emunot / (am->p_atmosMunot * am->p_atmosMunot);
306 }
307 else {
308 tfac = (emunot - emu) / (am->p_atmosMunot - mu);
309 }
310 if(am->p_atmosAtmSwitch == 1) {
311 result = mu * (phasefn - 1.0) * tfac;
312 }
313 else if(am->p_atmosAtmSwitch == 2) {
314 result = am->p_atmosMunot * mu * (phasefn - 1.0) * (1.0 - emunot * emu) / (am->p_atmosMunot + mu);
315 }
316 else if(am->p_atmosAtmSwitch == 3) {
317 result = -sine * am->p_atmosCosphi * (phasefn - 1.0) * tfac;
318 }
319 else {
320 string msg = "NumericalAtmosApprox::InrFunc2Bint() - Invalid value of atmospheric ";
321 msg += "switch used as argument to this function";
322 throw IException(IException::Programmer, msg, _FILEINFO_);
323 }
324 }
325 return result;
326 }
327}//end NumericalAtmosApprox.cpp
Isotropic atmos scattering model.
Definition AtmosModel.h:60
double AtmosTau() const
Return atmospheric Tau value.
Definition AtmosModel.h:119
double AtmosHga() const
Return atmospheric Hga value.
Definition AtmosModel.h:127
Isis exception class.
Definition IException.h:91
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
ErrorType errorType() const
Returns the source of the error for this exception.
Adds specific functionality to C++ strings.
Definition IString.h:165
double Evaluate(const double a, const ExtrapType &etype=ThrowError)
Calculates interpolated or extrapolated value of tabulated data set for given domain value.
@ PolynomialNeville
Polynomial interpolation using Neville's algorithm.
void Reset()
Resets the state of the object.
NumericalApproximation(const NumericalApproximation::InterpType &itype=CubicNatural)
Default constructor creates NumericalApproximation object.
@ Extrapolate
Evaluate() attempts to extrapolate if a is outside of the domain. This is only valid for NumericalApp...
vector< double > PolynomialNevilleErrorEstimate()
Retrieves the error estimate for the Neville's polynomial interpolation type.
void AddData(const double x, const double y)
Add a datapoint to the set.
static double InrFunc2Bint(AtmosModel *am, double mu)
Inner function to be integrated.
double RombergsMethod(AtmosModel *am, IntegFunc sub, double a, double b)
This variation on the NumericalApproximation method integrates a specified AtmosModel function rather...
NumericalAtmosApprox(const NumericalApproximation::InterpType &itype=CubicNatural)
Uses Isis::NumericalApproximation constructor.
static double OutrFunc2Bint(AtmosModel *am, double phi)
This function is the outer integrand over mu at specified phi.
double RefineExtendedTrap(AtmosModel *am, IntegFunc sub, double a, double b, double s, unsigned int n)
This variation on the NumericalApproximation method integrates a specified AtmosModel function rather...
IntegFunc
This enum defines function to be integrated by Romberg's method.
@ InnerFunction
Indicates that Romberg's method will integrate the function InrFunc2Bint()
double CalcSurfAlbedo(double pha, double inc, double ema)
Calculate the surface brightness using photometric angle information.
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
Namespace for the standard library.