Isis 3 Programmer Reference
MoonAlbedo.cpp
1
6/* SPDX-License-Identifier: CC0-1.0 */
7#include <cmath>
8#include "MoonAlbedo.h"
9#include "SpecialPixel.h"
10#include "IException.h"
11
12namespace Isis {
13 MoonAlbedo::MoonAlbedo(Pvl &pvl, PhotoModel &pmodel) :
14 NormModel(pvl, pmodel) {
15 PvlGroup &algo = pvl.findObject("NormalizationModel")
16 .findGroup("Algorithm", Pvl::Traverse);
17 // Set default values
18 SetNormD(0.14);
19 SetNormE(-0.4179);
20 SetNormF(0.55);
21 SetNormG2(0.02);
22 SetNormXmul(1.0);
23 SetNormWl(p_normWavelength);
24 SetNormH(0.048);
25 SetNormBsh1(0.08);
26 SetNormXb1(-0.0817);
27 SetNormXb2(0.0081);
28
29 // Get values from user
30 if(algo.hasKeyword("D")) {
31 SetNormD(algo["D"]);
32 }
33
34 if(algo.hasKeyword("Wl")) {
35 SetNormWl(algo["Wl"]);
36 }
37
38 if(algo.hasKeyword("E")) {
39 SetNormE(algo["E"]);
40 }
41 else {
42 if(p_normWl < 1.0) {
43 p_normE = -0.3575 * p_normWl - 0.0607;
44 }
45 else {
46 p_normE = -0.4179;
47 }
48 }
49
50 if(algo.hasKeyword("F")) {
51 SetNormF(algo["F"]);
52 }
53
54 if(algo.hasKeyword("G2")) {
55 SetNormG2(algo["G2"]);
56 }
57 else {
58 if(p_normWl < 1.0) {
59 p_normG2 = -0.9585 * p_normWl + 0.98;
60 }
61 else {
62 p_normG2 = 0.02;
63 }
64 }
65
66 if(algo.hasKeyword("Xmul")) {
67 SetNormXmul(algo["Xmul"]);
68 }
69
70 if(algo.hasKeyword("H")) {
71 SetNormH(algo["H"]);
72 }
73
74 if(algo.hasKeyword("Bsh1")) {
75 SetNormBsh1(algo["Bsh1"]);
76 }
77 else {
78 p_normBsh1 = 19.89 - 59.58 * p_normWl + 59.86 * pow(p_normWl, 2) -
79 20.09 * pow(p_normWl, 3);
80 if(p_normBsh1 < 0.0) {
81 p_normBsh1 = 0.0;
82 }
83 }
84
85 if(algo.hasKeyword("Xb1")) {
86 SetNormXb1(algo["Xb1"]);
87 }
88
89 if(algo.hasKeyword("Xb2")) {
90 SetNormXb2(algo["Xb2"]);
91 }
92
93 // Initialize values that will be needed to normalize to a
94 // Buratti function at phase = 2.0 degrees
95 p_normF1 = 1.0 - p_normF;
96 double g1 = p_normD * 0.1 + p_normE;
97 double g1sq = g1 * g1;
98 p_normG2sq = p_normG2 * p_normG2;
99 double c30 = cos(30.0 * Isis::PI / 180.0);
100 if(1.0 + g1sq + 2.0 * g1 *c30 <= 0.0) {
101 std::string msg = "Error while initializing Buratti function";
102 throw IException(IException::Unknown, msg, _FILEINFO_);
103 }
104 double pg130 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * c30), 1.5));
105 if(1.0 + p_normG2sq + 2.0 * p_normG2 *c30 <= 0.0) {
106 std::string msg = "Error while initializing Buratti function";
107 throw IException(IException::Unknown, msg, _FILEINFO_);
108 }
109 double pg230 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * c30), 1.5));
110 if(p_normBsh1 < 0.0) p_normBsh1 = 0.0;
111 if(1.0 + tan(15.0 * Isis::PI / 180.0) / p_normH == 0.0) {
112 std::string msg = "Error while initializing Buratti function";
113 throw IException(IException::Unknown, msg, _FILEINFO_);
114 }
115 double bshad30 = 1.0 + p_normBsh1 / (1.0 + tan(15.0 * Isis::PI / 180.0) / p_normH);
116 p_normPg30 = (pg130 + pg230) * bshad30;
117 p_normBc1 = p_normXb1 + p_normXb2 * p_normWl;
118 p_normFbc3 = 1.0 + p_normBc1 * 2.0;
119 if(p_normFbc3 == 0.0) {
120 std::string msg = "Error while initializing Buratti function";
121 throw IException(IException::Unknown, msg, _FILEINFO_);
122 }
123 p_normC3 = cos(2.0 * Isis::PI / 180.0);
124 if(1.0 + p_normG2sq + 2.0 * p_normG2 *p_normC3 <= 0.0) {
125 std::string msg = "Error while initializing Buratti function";
126 throw IException(IException::Unknown, msg, _FILEINFO_);
127 }
128 p_normPg32 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * p_normC3), 1.5));
129 if(1.0 + tan(Isis::PI / 180.0) / p_normH == 0.0) {
130 std::string msg = "Error while initializing Buratti function";
131 throw IException(IException::Unknown, msg, _FILEINFO_);
132 }
133 p_normBshad3 = 1.0 + p_normBsh1 / (1.0 + tan(Isis::PI / 180.0) / p_normH);
134 }
135
136 void MoonAlbedo::NormModelAlgorithm(double phase, double incidence, double emission,
137 double demincidence, double dememission, double dn,
138 double &albedo, double &mult, double &base) {
139 double a;
140 double cosa;
141 double pg2;
142 double bshad;
143 double r;
144 double g1;
145 double g1sq;
146 double pg1;
147 double pg;
148 double fbc;
149 double pg31;
150 double pg3;
151
152 a = GetPhotoModel()->CalcSurfAlbedo(phase, demincidence, dememission);
153
154 if(a != 0.0) {
155 cosa = cos(phase * Isis::PI / 180.0);
156 if(1.0 + p_normG2sq + 2.0 * p_normG2 *cosa <= 0.0) {
157 albedo = NULL8;
158 return;
159 }
160 pg2 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * cosa), 1.5));
161 if(1.0 + tan(phase * .5 * Isis::PI / 180.0) / p_normH == 0.0) {
162 albedo = NULL8;
163 return;
164 }
165 bshad = 1.0 + p_normBsh1 / (1.0 + tan(phase * .5 * Isis::PI / 180.0) / p_normH);
166 r = dn * p_normXmul;
167 // Estimate the albedo at 0.1, then iterate
168 albedo = 0.1;
169 for(int i = 0; i < 6; i++) {
170 g1 = p_normD * albedo + p_normE;
171 g1sq = g1 * g1;
172 if(1.0 + g1sq + 2.0 * g1 *cosa <= 0.0) {
173 albedo = NULL8;
174 return;
175 }
176 pg1 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * cosa), 1.5));
177 pg = (pg1 + pg2) * bshad;
178 if(phase <= 2.0) {
179 fbc = 1.0 + p_normBc1 * phase;
180 if(1.0 + g1sq + 2.0 * g1 *p_normC3 <= 0.0) {
181 albedo = NULL8;
182 return;
183 }
184 pg31 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * p_normC3), 1.5));
185 pg3 = (pg31 + p_normPg32) * p_normBshad3;
186 pg = fbc * (pg3 / p_normFbc3);
187 }
188 if(pg == 0.0) {
189 albedo = NULL8;
190 return;
191 }
192 albedo = r * a * p_normPg30 / pg;
193 }
194 }
195 else {
196 albedo = NULL8;
197 return;
198 }
199 }
200
208 void MoonAlbedo::SetNormD(const double d) {
209 p_normD = d;
210 }
211
220 void MoonAlbedo::SetNormWl(const double wl) {
221 p_normWl = wl;
222 }
223
231 void MoonAlbedo::SetNormE(const double e) {
232 p_normE = e;
233 }
234
242 void MoonAlbedo::SetNormF(const double f) {
243 p_normF = f;
244 }
245
253 void MoonAlbedo::SetNormG2(const double g2) {
254 p_normG2 = g2;
255 }
256
264 void MoonAlbedo::SetNormXmul(const double xmul) {
265 p_normXmul = xmul;
266 }
267
275 void MoonAlbedo::SetNormH(const double h) {
276 if(h == 0.0) {
277 std::string msg = "Invalid value of normalization h [" +
278 IString(h) + "]";
279 throw IException(IException::User, msg, _FILEINFO_);
280 }
281
282 p_normH = h;
283 }
284
292 void MoonAlbedo::SetNormBsh1(const double bsh1) {
293 if(bsh1 < 0.0) {
294 std::string msg = "Invalid value of normalization bsh1 [" +
295 IString(bsh1) + "]";
296 throw IException(IException::User, msg, _FILEINFO_);
297 }
298
299 p_normBsh1 = bsh1;
300 }
301
309 void MoonAlbedo::SetNormXb1(const double xb1) {
310 p_normXb1 = xb1;
311 }
312
320 void MoonAlbedo::SetNormXb2(const double xb2) {
321 p_normXb2 = xb2;
322 }
323}
324
325extern "C" Isis::NormModel *MoonAlbedoPlugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel) {
326 return new Isis::MoonAlbedo(pvl, pmodel);
327}
Isis::IException
Isis exception class.
Definition IException.h:91
Isis::IString
Adds specific functionality to C++ strings.
Definition IString.h:165
Isis::MoonAlbedo
Albedo dependent phase function normalization for the Moon.
Definition MoonAlbedo.h:25
Isis::Pvl
Container for cube-like labels.
Definition Pvl.h:119
Isis
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
Isis::PI
const double PI
The mathematical constant PI.
Definition Constants.h:40