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
12#define MIN(x,y) (((x) < (y)) ? (x) : (y))
13#define MAX(x,y) (((x) > (y)) ? (x) : (y))
14
15namespace Isis {
16 MoonAlbedo::MoonAlbedo(Pvl &pvl, PhotoModel &pmodel) :
17 NormModel(pvl, pmodel) {
18 PvlGroup &algo = pvl.findObject("NormalizationModel")
19 .findGroup("Algorithm", Pvl::Traverse);
20 // Set default values
21 SetNormD(0.14);
22 SetNormE(-0.4179);
23 SetNormF(0.55);
24 SetNormG2(0.02);
25 SetNormXmul(1.0);
26 SetNormWl(p_normWavelength);
27 SetNormH(0.048);
28 SetNormBsh1(0.08);
29 SetNormXb1(-0.0817);
30 SetNormXb2(0.0081);
31
32 // Get values from user
33 if(algo.hasKeyword("D")) {
34 SetNormD(algo["D"]);
35 }
36
37 if(algo.hasKeyword("Wl")) {
38 SetNormWl(algo["Wl"]);
39 }
40
41 if(algo.hasKeyword("E")) {
42 SetNormE(algo["E"]);
43 }
44 else {
45 if(p_normWl < 1.0) {
46 p_normE = -0.3575 * p_normWl - 0.0607;
47 }
48 else {
49 p_normE = -0.4179;
50 }
51 }
52
53 if(algo.hasKeyword("F")) {
54 SetNormF(algo["F"]);
55 }
56
57 if(algo.hasKeyword("G2")) {
58 SetNormG2(algo["G2"]);
59 }
60 else {
61 if(p_normWl < 1.0) {
62 p_normG2 = -0.9585 * p_normWl + 0.98;
63 }
64 else {
65 p_normG2 = 0.02;
66 }
67 }
68
69 if(algo.hasKeyword("Xmul")) {
70 SetNormXmul(algo["Xmul"]);
71 }
72
73 if(algo.hasKeyword("H")) {
74 SetNormH(algo["H"]);
75 }
76
77 if(algo.hasKeyword("Bsh1")) {
78 SetNormBsh1(algo["Bsh1"]);
79 }
80 else {
81 p_normBsh1 = 19.89 - 59.58 * p_normWl + 59.86 * pow(p_normWl, 2) -
82 20.09 * pow(p_normWl, 3);
83 if(p_normBsh1 < 0.0) {
84 p_normBsh1 = 0.0;
85 }
86 }
87
88 if(algo.hasKeyword("Xb1")) {
89 SetNormXb1(algo["Xb1"]);
90 }
91
92 if(algo.hasKeyword("Xb2")) {
93 SetNormXb2(algo["Xb2"]);
94 }
95
96 // Initialize values that will be needed to normalize to a
97 // Buratti function at phase = 2.0 degrees
98 p_normF1 = 1.0 - p_normF;
99 double g1 = p_normD * 0.1 + p_normE;
100 double g1sq = g1 * g1;
101 p_normG2sq = p_normG2 * p_normG2;
102 double c30 = cos(30.0 * Isis::PI / 180.0);
103 if(1.0 + g1sq + 2.0 * g1 *c30 <= 0.0) {
104 std::string msg = "Error while initializing Buratti function";
105 throw IException(IException::Unknown, msg, _FILEINFO_);
106 }
107 double pg130 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * c30), 1.5));
108 if(1.0 + p_normG2sq + 2.0 * p_normG2 *c30 <= 0.0) {
109 std::string msg = "Error while initializing Buratti function";
110 throw IException(IException::Unknown, msg, _FILEINFO_);
111 }
112 double pg230 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * c30), 1.5));
113 if(p_normBsh1 < 0.0) p_normBsh1 = 0.0;
114 if(1.0 + tan(15.0 * Isis::PI / 180.0) / p_normH == 0.0) {
115 std::string msg = "Error while initializing Buratti function";
116 throw IException(IException::Unknown, msg, _FILEINFO_);
117 }
118 double bshad30 = 1.0 + p_normBsh1 / (1.0 + tan(15.0 * Isis::PI / 180.0) / p_normH);
119 p_normPg30 = (pg130 + pg230) * bshad30;
120 p_normBc1 = p_normXb1 + p_normXb2 * p_normWl;
121 p_normFbc3 = 1.0 + p_normBc1 * 2.0;
122 if(p_normFbc3 == 0.0) {
123 std::string msg = "Error while initializing Buratti function";
124 throw IException(IException::Unknown, msg, _FILEINFO_);
125 }
126 p_normC3 = cos(2.0 * Isis::PI / 180.0);
127 if(1.0 + p_normG2sq + 2.0 * p_normG2 *p_normC3 <= 0.0) {
128 std::string msg = "Error while initializing Buratti function";
129 throw IException(IException::Unknown, msg, _FILEINFO_);
130 }
131 p_normPg32 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * p_normC3), 1.5));
132 if(1.0 + tan(Isis::PI / 180.0) / p_normH == 0.0) {
133 std::string msg = "Error while initializing Buratti function";
134 throw IException(IException::Unknown, msg, _FILEINFO_);
135 }
136 p_normBshad3 = 1.0 + p_normBsh1 / (1.0 + tan(Isis::PI / 180.0) / p_normH);
137 }
138
139 void MoonAlbedo::NormModelAlgorithm(double phase, double incidence, double emission,
140 double demincidence, double dememission, double dn,
141 double &albedo, double &mult, double &base) {
142 double a;
143 double cosa;
144 double pg2;
145 double bshad;
146 double r;
147 double g1;
148 double g1sq;
149 double pg1;
150 double pg;
151 double fbc;
152 double pg31;
153 double pg3;
154
155 a = GetPhotoModel()->CalcSurfAlbedo(phase, demincidence, dememission);
156
157 if(a != 0.0) {
158 cosa = cos(phase * Isis::PI / 180.0);
159 if(1.0 + p_normG2sq + 2.0 * p_normG2 *cosa <= 0.0) {
160 albedo = NULL8;
161 return;
162 }
163 pg2 = p_normF * (1.0 - p_normG2sq) / (pow((1.0 + p_normG2sq + 2.0 * p_normG2 * cosa), 1.5));
164 if(1.0 + tan(phase * .5 * Isis::PI / 180.0) / p_normH == 0.0) {
165 albedo = NULL8;
166 return;
167 }
168 bshad = 1.0 + p_normBsh1 / (1.0 + tan(phase * .5 * Isis::PI / 180.0) / p_normH);
169 r = dn * p_normXmul;
170 // Estimate the albedo at 0.1, then iterate
171 albedo = 0.1;
172 for(int i = 0; i < 6; i++) {
173 g1 = p_normD * albedo + p_normE;
174 g1sq = g1 * g1;
175 if(1.0 + g1sq + 2.0 * g1 *cosa <= 0.0) {
176 albedo = NULL8;
177 return;
178 }
179 pg1 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * cosa), 1.5));
180 pg = (pg1 + pg2) * bshad;
181 if(phase <= 2.0) {
182 fbc = 1.0 + p_normBc1 * phase;
183 if(1.0 + g1sq + 2.0 * g1 *p_normC3 <= 0.0) {
184 albedo = NULL8;
185 return;
186 }
187 pg31 = p_normF1 * (1.0 - g1sq) / (pow((1.0 + g1sq + 2.0 * g1 * p_normC3), 1.5));
188 pg3 = (pg31 + p_normPg32) * p_normBshad3;
189 pg = fbc * (pg3 / p_normFbc3);
190 }
191 if(pg == 0.0) {
192 albedo = NULL8;
193 return;
194 }
195 albedo = r * a * p_normPg30 / pg;
196 }
197 }
198 else {
199 albedo = NULL8;
200 return;
201 }
202 }
203
211 void MoonAlbedo::SetNormD(const double d) {
212 p_normD = d;
213 }
214
223 void MoonAlbedo::SetNormWl(const double wl) {
224 p_normWl = wl;
225 }
226
234 void MoonAlbedo::SetNormE(const double e) {
235 p_normE = e;
236 }
237
245 void MoonAlbedo::SetNormF(const double f) {
246 p_normF = f;
247 }
248
256 void MoonAlbedo::SetNormG2(const double g2) {
257 p_normG2 = g2;
258 }
259
267 void MoonAlbedo::SetNormXmul(const double xmul) {
268 p_normXmul = xmul;
269 }
270
278 void MoonAlbedo::SetNormH(const double h) {
279 if(h == 0.0) {
280 std::string msg = "Invalid value of normalization h [" +
281 IString(h) + "]";
282 throw IException(IException::User, msg, _FILEINFO_);
283 }
284
285 p_normH = h;
286 }
287
295 void MoonAlbedo::SetNormBsh1(const double bsh1) {
296 if(bsh1 < 0.0) {
297 std::string msg = "Invalid value of normalization bsh1 [" +
298 IString(bsh1) + "]";
299 throw IException(IException::User, msg, _FILEINFO_);
300 }
301
302 p_normBsh1 = bsh1;
303 }
304
312 void MoonAlbedo::SetNormXb1(const double xb1) {
313 p_normXb1 = xb1;
314 }
315
323 void MoonAlbedo::SetNormXb2(const double xb2) {
324 p_normXb2 = xb2;
325 }
326}
327
328extern "C" Isis::NormModel *MoonAlbedoPlugin(Isis::Pvl &pvl, Isis::PhotoModel &pmodel) {
329 return new Isis::MoonAlbedo(pvl, pmodel);
330}
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