13 Hapke::Hapke(Pvl &pvl) : PhotoModel(pvl) {
18 p_photoThetaold = -999.0;
23 PvlGroup &algorithm = pvl.findObject(
"PhotometricModel").findGroup(
"Algorithm",
Pvl::Traverse);
25 p_algName = AlgorithmName().toUpper();
27 if(algorithm.hasKeyword(
"Hg1")) {
28 SetPhotoHg1(algorithm[
"Hg1"]);
31 if(algorithm.hasKeyword(
"Hg2")) {
32 SetPhotoHg2(algorithm[
"Hg2"]);
35 if(algorithm.hasKeyword(
"Bh")) {
36 SetPhotoBh(algorithm[
"Bh"]);
39 if(algorithm.hasKeyword(
"Ch")) {
40 SetPhotoCh(algorithm[
"Ch"]);
43 if(algorithm.hasKeyword(
"ZeroB0Standard")) {
44 SetPhoto0B0Standard(algorithm[
"ZeroB0Standard"][0]);
45 }
else if (algorithm.hasKeyword(
"ZeroB0St")) {
46 SetPhoto0B0Standard(algorithm[
"ZeroB0St"][0]);
48 SetPhoto0B0Standard(
"TRUE");
51 if(algorithm.hasKeyword(
"Wh")) {
52 SetPhotoWh(algorithm[
"Wh"]);
55 if(algorithm.hasKeyword(
"Hh")) {
56 SetPhotoHh(algorithm[
"Hh"]);
59 if(algorithm.hasKeyword(
"B0")) {
60 SetPhotoB0(algorithm[
"B0"]);
63 p_photoB0save = p_photoB0;
65 if(algorithm.hasKeyword(
"Theta")) {
66 SetPhotoTheta(algorithm[
"Theta"]);
103 double Hapke::PhotoModelAlgorithm(
double phase,
double incidence,
105 static double pht_hapke;
153 static double old_phase = -9999;
154 static double old_incidence = -9999;
155 static double old_emission= -9999;
157 if (old_phase == phase && old_incidence == incidence && old_emission == emission) {
162 old_incidence = incidence;
163 old_emission = emission;
165 pharad = phase *
PI / 180.0;
166 incrad = incidence *
PI / 180.0;
167 emarad = emission *
PI / 180.0;
171 if(p_photoTheta != p_photoThetaold) {
172 cost = cos(p_photoTheta *
PI / 180.0);
173 sint = sin(p_photoTheta *
PI / 180.0);
174 p_photoCott = cost / max(1.0e-10, sint);
175 p_photoCot2t = p_photoCott * p_photoCott;
176 p_photoTant = sint / cost;
177 tan2t = p_photoTant * p_photoTant;
178 p_photoSr = sqrt(1.0 +
PI * tan2t);
179 p_photoOsr = 1.0 / p_photoSr;
183 if(incidence >= 90.0) {
188 gamma = sqrt(1.0 - p_photoWh);
189 hgs = p_photoHg1 * p_photoHg1;
192 tang2 = tan(pharad/2.0);
194 if(p_photoHh == 0.0) {
198 bg = p_photoB0 / (1.0 + tang2 / p_photoHh);
201 if (p_algName ==
"HAPKEHEN") {
202 pg1 = (1.0 - p_photoHg2) * (1.0 - hgs) / pow((1.0 + hgs + 2.0 *
203 p_photoHg1 * cosg), 1.5);
204 pg2 = p_photoHg2 * (1.0 - hgs) / pow((1.0 + hgs - 2.0 *
205 p_photoHg1 * cosg), 1.5);
208 pg = 1.0 + p_photoBh * cosg + p_photoCh * (1.5 * pow(cosg, 2.0) - .5);
212 if(p_photoTheta <= 0.0) {
213 pht_hapke = p_photoWh / 4.0 * munot / (munot + mu) * ((1.0 + bg) *
214 pg - 1.0 +
Hfunc(munot, gamma) *
Hfunc(mu, gamma));
219 coti = munot / max(1.0e-10, sini);
221 ecoti = exp(min(-p_photoCot2t * cot2i /
PI , 23.0));
222 ecot2i = exp(min(-2.0 * p_photoCott * coti /
PI, 23.0));
223 u0p0 = p_photoOsr * (munot + sini * p_photoTant * ecoti / (2.0 - ecot2i));
226 cote = mu / max(1.0e-10, sine);
237 caz = (cosg - cosei) / sinei;
245 az = acos(caz) * 180.0 /
PI;
254 tanaz2 = tan(az2 *
PI / 180.0);
255 faz = exp(min(-2.0 * tanaz2, 23.0));
258 sin2a2 = pow(sin(az2 *
PI / 180.0), 2.0);
261 ecote = exp(min(-p_photoCot2t * cot2e /
PI, 23.0));
262 ecot2e = exp(min(-2.0 * p_photoCott * cote /
PI, 23.0));
263 up0 = p_photoOsr * (mu + sine * p_photoTant * ecote / (2.0 - ecot2e));
265 if(incidence <= emission) {
266 q = p_photoOsr * munot / u0p0;
269 q = p_photoOsr * mu / up0;
272 if(incidence <= emission) {
273 ecei = (2.0 - ecot2e - api * ecot2i);
274 s2ei = sin2a2 * ecoti;
275 u0p = p_photoOsr * (munot + sini * p_photoTant * (caz * ecote + s2ei) / ecei);
276 up = p_photoOsr * (mu + sine * p_photoTant * (ecote - s2ei) / ecei);
279 ecee = (2.0 - ecot2i - api * ecot2e);
280 s2ee = sin2a2 * ecote;
281 u0p = p_photoOsr * (munot + sini * p_photoTant * (ecoti - s2ee) / ecee);
282 up = p_photoOsr * (mu + sine * p_photoTant * (caz * ecoti + s2ee) / ecee);
285 rr1 = p_photoWh / 4.0 * u0p / (u0p + up) * ((1.0 + bg) * pg -
287 rr2 = up * munot / (up0 * u0p0 * p_photoSr * (1.0 - faz + faz * q));
288 pht_hapke = rr1 * rr2;
302 if(hg1 <= -1.0 || hg1 >= 1.0) {
303 string msg =
"Invalid value of Hapke Henyey Greenstein hg1 [" +
319 if(hg2 < 0.0 || hg2 > 1.0) {
320 string msg =
"Invalid value of Hapke Henyey Greenstein hg2 [" +
336 if(bh < -1.0 || bh > 1.0) {
337 string msg =
"Invalid value of Hapke Legendre bh [" +
353 if(ch < -1.0 || ch > 1.0) {
354 string msg =
"Invalid value of Hapke Legendre ch [" +
369 if(wh <= 0.0 || wh > 1.0) {
370 string msg =
"Invalid value of Hapke wh [" +
386 string msg =
"Invalid value of Hapke hh [" +
402 string msg =
"Invalid value of Hapke b0 [" +
421 if(temp !=
"NO" && temp !=
"YES" && temp !=
"FALSE" && temp !=
"TRUE") {
422 string msg =
"Invalid value of Hapke ZeroB0Standard [" +
427 if (temp ==
"NO" || temp ==
"FALSE") p_photo0B0Standard =
"FALSE";
428 if (temp ==
"YES" || temp ==
"TRUE") p_photo0B0Standard =
"TRUE";
439 if(theta < 0.0 || theta > 90.0) {
440 string msg =
"Invalid value of Hapke theta [" +
444 p_photoTheta = theta;
452 p_photoB0save = p_photoB0;
453 if (p_photo0B0Standard ==
"TRUE") p_photoB0 = 0.0;
456 p_photoB0 = p_photoB0save;
const double PI
The mathematical constant PI.
void SetPhotoHg2(const double hg2)
Return photometric Hg1 value.
void SetOldTheta(double theta)
Return photometric Theta value.
Namespace for the standard library.
void SetPhotoBh(const double bh)
Return photometric Hg2 value.
void SetPhotoHg1(const double hg1)
Set the Hapke Henyey Greenstein coefficient for the single particle phase function.
void SetPhoto0B0Standard(const QString &b0standard)
Determine if the Hapke opposition surge component is initialized to zero during the SetStandardCondit...
void SetPhotoCh(const double ch)
Return photometric Bh value.
#define _FILEINFO_
Macro for the filename and line number.
A type of error that could only have occurred due to a mistake on the user's part (e...
void SetPhotoHh(const double hh)
Return photometric Wh value.
Container for cube-like labels.
void SetStandardConditions(bool standard)
Sets whether standard conditions will be used.
Adds specific functionality to C++ strings.
Namespace for ISIS/Bullet specific routines.
double Hfunc(double u, double gamma)
Hapke's approximation to Chandra's H function.
void SetPhotoTheta(const double theta)
Return photometric B0 value.
IString UpCase()
Converst any lower case characters in the object IString with uppercase characters.
virtual void SetStandardConditions(bool standard)
Sets whether standard conditions will be used.
void SetPhotoWh(const double wh)
Return photometric Ch value.
Hapke-Henyey-Greenstein photometric model.
void SetPhotoB0(const double b0)
Return photometric Hh value.