Isis 3.0 Programmer Reference
Back | Home
AlbedoAtm.cpp
1 #include <cmath>
2 #include "AlbedoAtm.h"
3 #include "NumericalApproximation.h"
4 #include "IException.h"
5 #include "IString.h"
6 
7 #define MIN(x,y) (((x) < (y)) ? (x) : (y))
8 #define MAX(x,y) (((x) > (y)) ? (x) : (y))
9 
10 namespace Isis {
23  AlbedoAtm::AlbedoAtm(Pvl &pvl, PhotoModel &pmodel, AtmosModel &amodel) :
24  NormModel(pvl, pmodel, amodel) {
25  PvlGroup &algo = pvl.findObject("NormalizationModel")
26  .findGroup("Algorithm", Pvl::Traverse);
27  // Set default value
28  SetNormPharef(0.0);
29  SetNormIncref(0.0);
30  SetNormEmaref(0.0);
31 
32  // Get value from user
33  if(algo.hasKeyword("Incref")) {
34  SetNormIncref(algo["Incref"]);
35  }
36 
37  if(algo.hasKeyword("Pharef")) {
38  SetNormPharef(algo["Pharef"]);
39  } else {
40  p_normPharef = p_normIncref;
41  }
42 
43  if(algo.hasKeyword("Emaref")) {
44  SetNormEmaref(algo["Emaref"]);
45  }
46 
47  // First-time setup:
48  // Calculate normalization at standard conditions
49  GetPhotoModel()->SetStandardConditions(true);
50  p_normPsurfref = GetPhotoModel()->CalcSurfAlbedo(p_normPharef, p_normIncref, p_normEmaref);
51  GetPhotoModel()->SetStandardConditions(false);
52 
53  // Get reference hemispheric albedo
54  GetAtmosModel()->GenerateAhTable();
55 
56  p_normAhref = (GetAtmosModel()->AtmosAhSpline()).Evaluate(p_normIncref, NumericalApproximation::Extrapolate);
57 
58  p_normMunotref = cos((PI / 180.0) * p_normIncref);
59 
60  // Now calculate atmosphere at standard conditions
61  GetAtmosModel()->SetStandardConditions(true);
62  GetAtmosModel()->CalcAtmEffect(p_normPharef, p_normIncref, p_normEmaref,
63  &p_normPstdref, &p_normTransref,
64  &p_normTrans0ref, &p_normSbar,
65  &p_normTranss);
66  GetAtmosModel()->SetStandardConditions(false);
67  }
68 
82  void AlbedoAtm::NormModelAlgorithm(double phase, double incidence, double emission,
83  double demincidence, double dememission, double dn,
84  double &albedo, double &mult, double &base) {
85  static double psurf;
86  static double ahInterp;
87  static double munot;
88  double pstd;
89  double trans;
90  double trans0;
91  double transs;
92  double rho;
93  double dpo;
94  double q;
95  double dpm;
96  double firsterm;
97  double secondterm;
98  double thirdterm;
99  double fourthterm;
100  double fifthterm;
101 
102  static double old_phase = -9999;
103  static double old_incidence = -9999;
104  static double old_emission = -9999;
105  static double old_demincidence = -9999;
106  static double old_dememission = -9999;
107 
108  if (old_phase != phase || old_incidence != incidence || old_emission != emission ||
109  old_demincidence != demincidence || old_dememission != dememission) {
110 
111  psurf = GetPhotoModel()->CalcSurfAlbedo(phase, demincidence,
112  dememission);
113 
114  ahInterp = (GetAtmosModel()->AtmosAhSpline()).Evaluate(incidence, NumericalApproximation::Extrapolate);
115 
116  munot = cos(incidence * (PI / 180.0));
117 
118  old_phase = phase;
119  old_incidence = incidence;
120  old_emission = emission;
121  old_demincidence = demincidence;
122  old_dememission = dememission;
123  }
124 
125  GetAtmosModel()->CalcAtmEffect(phase, incidence, emission, &pstd, &trans, &trans0, &p_normSbar,
126  &transs);
127 
128  // With model at actual geometry, calculate rho from dn
129  dpo = dn - pstd;
130  dpm = (psurf - ahInterp * munot) * trans0;
131  q = ahInterp * munot * trans + GetAtmosModel()->AtmosAb() * p_normSbar * dpo + dpm;
132 
133  if(dpo <= 0.0 && GetAtmosModel()->AtmosNulneg()) {
134  rho = 0.0;
135  }
136  else {
137  firsterm = GetAtmosModel()->AtmosAb() * p_normSbar;
138  secondterm = dpo * dpm;
139  thirdterm = firsterm * secondterm;
140  fourthterm = pow(q, 2.0) - 4.0 * thirdterm;
141 
142  if(fourthterm < 0.0) {
143  QString msg = "Square root of negative (math) encountered";
145  }
146  else {
147  fifthterm = q + sqrt(fourthterm);
148  }
149 
150  rho = 2 * dpo / fifthterm;
151  }
152 
153  // Now use rho and reference geometry to calculate output dnout
154  if((1.0 - rho * GetAtmosModel()->AtmosAb()*p_normSbar) <= 0.0) {
155  QString msg = "Divide by zero (math) encountered";
157  }
158  else {
159  albedo = p_normPstdref + rho * (p_normAhref * p_normMunotref *
160  p_normTransref / (1.0 - rho * GetAtmosModel()->AtmosAb() *
161  p_normSbar) + (p_normPsurfref - p_normAhref *
162  p_normMunotref) * p_normTrans0ref);
163  }
164  }
165 
175  void AlbedoAtm::SetNormPharef(const double pharef) {
176  if(pharef < 0.0 || pharef >= 180.0) {
177  QString msg = "Invalid value of normalization pharef [" +
178  toString(pharef) + "]";
180  }
181 
182  p_normPharef = pharef;
183  }
184 
194  void AlbedoAtm::SetNormIncref(const double incref) {
195  if(incref < 0.0 || incref >= 90.0) {
196  QString msg = "Invalid value of normalization incref [" +
197  toString(incref) + "]";
199  }
200 
201  p_normIncref = incref;
202  }
203 
213  void AlbedoAtm::SetNormEmaref(const double emaref) {
214  if(emaref < 0.0 || emaref >= 90.0) {
215  QString msg = "Invalid value of normalization emaref [" +
216  toString(emaref) + "]";
218  }
219 
220  p_normEmaref = emaref;
221  }
222 }
223 
224 extern "C" Isis::NormModel *AlbedoAtmPlugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel, Isis::AtmosModel &amodel) {
225  return new Isis::AlbedoAtm(pvl, pmodel, amodel);
226 }
double AtmosAb() const
Return atmospheric Ab value.
Definition: AtmosModel.h:155
double CalcSurfAlbedo(double pha, double inc, double ema)
Calculate the surface brightness using photometric angle information.
Definition: PhotoModel.cpp:171
void SetNormPharef(const double pharef)
Set parameters needed for albedo normalization.
Definition: AlbedoAtm.cpp:175
void SetNormEmaref(const double emaref)
Set the normalization function parameter.
Definition: AlbedoAtm.cpp:213
virtual void SetStandardConditions(bool standard)
Used to calculate atmosphere at standard conditions.
Definition: AtmosModel.cpp:492
void CalcAtmEffect(double pha, double inc, double ema, double *pstd, double *trans, double *trans0, double *sbar, double *transs)
Calculate the atmospheric scattering effect using photometric angle information.
Definition: AtmosModel.cpp:469
PvlObjectIterator findObject(const QString &name, PvlObjectIterator beg, PvlObjectIterator end)
Find the index of object with a specified name, between two indexes.
Definition: PvlObject.h:286
const double PI(3.14159265358979323846)
The mathematical constant PI.
Search child objects.
Definition: PvlObject.h:170
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition: IString.cpp:226
void SetNormIncref(const double incref)
Set the normalization function parameter.
Definition: AlbedoAtm.cpp:194
NumericalApproximation AtmosAhSpline()
If GenerateAhTable() has been called this returns a clamped cubic spline of the data set (p_atmosIncT...
Definition: AtmosModel.h:199
AlbedoAtm(Pvl &pvl, PhotoModel &pmodel, AtmosModel &amodel)
Constructs AlbedoAtm object using a Pvl, PhotoModel, and AtmosModel.
Definition: AlbedoAtm.cpp:23
Isotropic atmos scattering model.
Definition: AtmosModel.h:76
Contains multiple PvlContainers.
Definition: PvlGroup.h:57
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:38
A type of error that could only have occurred due to a mistake on the user&#39;s part (e...
Definition: IException.h:134
A type of error that cannot be classified as any of the other error types.
Definition: IException.h:126
Container for cube-like labels.
Definition: Pvl.h:135
Isis exception class.
Definition: IException.h:99
void GenerateAhTable()
This method computes the values of the atmospheric Ah table and sets the properties of the atmospheri...
Definition: AtmosModel.cpp:531
Albedo normalization with atmosphere.
Definition: AlbedoAtm.h:53
Evaluate() attempts to extrapolate if a is outside of the domain. &lt; This is only valid for NumericalA...
virtual void SetStandardConditions(bool standard)
Sets whether standard conditions will be used.
Definition: PhotoModel.cpp:44
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.

U.S. Department of the Interior | U.S. Geological Survey
ISIS | Privacy & Disclaimers | Astrogeology Research Program
To contact us, please post comments and questions on the ISIS Support Center
File Modified: 07/12/2023 23:14:07