USGS

Isis 3.0 Object Programmers' Reference

 Home
Hapke.cpp
1 #include <cmath>
2 #include "Constants.h"
3 #include "Hapke.h"
4 #include "Pvl.h"
5 #include "PvlGroup.h"
6 #include "IException.h"
7 #include "IString.h"
8 
9 using namespace std;
10 
11 namespace Isis {
12 
13  Hapke::Hapke(Pvl &pvl) : PhotoModel(pvl) {
14  p_photoHh = 0.0;
15  p_photoB0 = 0.0;
16  p_photoTheta = 0.0;
17  p_photoWh = 0.5;
18  p_photoThetaold = -999.0;
19 
20  p_photoHg1 = 0.0;
21  p_photoHg2 = 0.0;
22 
23  PvlGroup &algorithm = pvl.findObject("PhotometricModel").findGroup("Algorithm", Pvl::Traverse);
24 
25  p_algName = AlgorithmName().toUpper();
26 
27  if(algorithm.hasKeyword("Hg1")) {
28  SetPhotoHg1(algorithm["Hg1"]);
29  }
30 
31  if(algorithm.hasKeyword("Hg2")) {
32  SetPhotoHg2(algorithm["Hg2"]);
33  }
34 
35  if(algorithm.hasKeyword("Bh")) {
36  SetPhotoBh(algorithm["Bh"]);
37  }
38 
39  if(algorithm.hasKeyword("Ch")) {
40  SetPhotoCh(algorithm["Ch"]);
41  }
42 
43  if(algorithm.hasKeyword("ZeroB0Standard")) {
44  SetPhoto0B0Standard(algorithm["ZeroB0Standard"][0]);
45  } else if (algorithm.hasKeyword("ZeroB0St")) {
46  SetPhoto0B0Standard(algorithm["ZeroB0St"][0]);
47  } else {
48  SetPhoto0B0Standard("TRUE");
49  }
50 
51  if(algorithm.hasKeyword("Wh")) {
52  SetPhotoWh(algorithm["Wh"]);
53  }
54 
55  if(algorithm.hasKeyword("Hh")) {
56  SetPhotoHh(algorithm["Hh"]);
57  }
58 
59  if(algorithm.hasKeyword("B0")) {
60  SetPhotoB0(algorithm["B0"]);
61  }
62 
63  p_photoB0save = p_photoB0;
64 
65  if(algorithm.hasKeyword("Theta")) {
66  SetPhotoTheta(algorithm["Theta"]);
67  }
68  }
69 
70  /*
71  * Computes normal albedo mult factor w/o opp surge from
72  * Hapke input parameters: W,H,BO,HG,THETA.
73  * Full-up Hapke's Law with macroscopic roughness. The photometric
74  * function multiplied back in will be modified to take out oppos-
75  * ition effect. This requires saving the actual value of B0 while
76  * temporarily setting it to zero in the common block to compute
77  * the overall normalization.
78  *
79  * @param phase Value of phase angle, in degrees.
80  * @param incidence Value of incidence angle, in degrees.
81  * @param emission Value of emission angle, in degrees.
82  * @returns <b>double</b>
83  *
84  *
85  * @history 1989-08-02 Unknown author in Isis2 under name pht_hapke
86  * @history 1991-08-07 Tammy Becker relinked hapke to new photompr
87  * @history 1997-02-16 James M Anderson - changed nonstandard degree trig
88  * to use radians
89  * @history 1999-01-11 KTT - Remove mu,munot,and alpha from the argument
90  * list and pass in only ema,inc, and phase. Remove
91  * lat and lon from argument list because they aren't
92  * used.
93  * @history 1999-03-01 K Teal Thompson Implement 1999-01-08 Randy Kirk
94  * Original Specs & Code. Declare vars, add implicit none.
95  * @history 1999-11-18 Randy Kirk - fixed minor typos, implemented return with
96  * smooth Hapke (Theta=0) result before doing rough Hapke
97  * calculations, allow single-particle-phase params = 0
98  * @history 2008-01-14 Janet Barrett - Imported into Isis3. Changed name from pht_hapke to PhotoModelAlgorithm()
99  * @history 2008-11-05 Jeannie Walldren - Added documentation
100  * from Isis2 files. Replaced Isis::PI with PI since this is in Isis namespace.
101  *
102  */
103  double Hapke::PhotoModelAlgorithm(double phase, double incidence,
104  double emission) {
105  static double pht_hapke;
106  double pharad; //phase angle in radians
107  double incrad; // incidence angle in radians
108  double emarad; // emission angle in radians
109  double munot;
110  double mu;
111  double cost;
112  double sint;
113  double tan2t;
114  double gamma;
115  double hgs;
116  double cosg;
117  double tang2;
118  double bg;
119  double pg;
120  double pg1;
121  double pg2;
122  double sini;
123  double coti;
124  double cot2i;
125  double ecoti;
126  double ecot2i;
127  double u0p0;
128  double sine;
129  double cote;
130  double cot2e;
131  double cosei;
132  double sinei;
133  double caz;
134  double az;
135  double az2;
136  double faz;
137  double tanaz2;
138  double sin2a2;
139  double api;
140  double ecote;
141  double ecot2e;
142  double up0;
143  double q;
144  double ecei;
145  double s2ei;
146  double u0p;
147  double up;
148  double ecee;
149  double s2ee;
150  double rr1;
151  double rr2;
152 
153  static double old_phase = -9999;
154  static double old_incidence = -9999;
155  static double old_emission= -9999;
156 
157  if (old_phase == phase && old_incidence == incidence && old_emission == emission) {
158  return pht_hapke;
159  }
160 
161  old_phase = phase;
162  old_incidence = incidence;
163  old_emission = emission;
164 
165  pharad = phase * PI / 180.0;
166  incrad = incidence * PI / 180.0;
167  emarad = emission * PI / 180.0;
168  munot = cos(incrad);
169  mu = cos(emarad);
170 
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;
180  SetOldTheta(p_photoTheta);
181  }
182 
183  if(incidence >= 90.0) {
184  pht_hapke = 0.0;
185  return pht_hapke;
186  }
187 
188  gamma = sqrt(1.0 - p_photoWh);
189  hgs = p_photoHg1 * p_photoHg1;
190 
191  cosg = cos(pharad);
192  tang2 = tan(pharad/2.0);
193 
194  if(p_photoHh == 0.0) {
195  bg = 0.0;
196  }
197  else {
198  bg = p_photoB0 / (1.0 + tang2 / p_photoHh);
199  }
200 
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);
206  pg = pg1 + pg2;
207  } else { // Hapke Legendre
208  pg = 1.0 + p_photoBh * cosg + p_photoCh * (1.5 * pow(cosg, 2.0) - .5);
209  }
210 
211  // If smooth Hapke is wanted then set Theta<=0.0
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));
215  return pht_hapke;
216  }
217 
218  sini = sin(incrad);
219  coti = munot / max(1.0e-10, sini);
220  cot2i = coti * coti;
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));
224 
225  sine = sin(emarad);
226  cote = mu / max(1.0e-10, sine);
227  cot2e = cote * cote;
228 
229  cosei = mu * munot;
230  sinei = sine * sini;
231 
232  if(sinei == 0.0) {
233  caz = 1.0;
234  az = 0.0;
235  }
236  else {
237  caz = (cosg - cosei) / sinei;
238  if(caz <= -1.0) {
239  az = 180.0;
240  }
241  else if(caz > 1.0) {
242  az = 0.0;
243  }
244  else {
245  az = acos(caz) * 180.0 / PI;
246  }
247  }
248 
249  az2 = az / 2.0;
250  if(az2 >= 90.0) {
251  faz = 0.0;
252  }
253  else {
254  tanaz2 = tan(az2 * PI / 180.0);
255  faz = exp(min(-2.0 * tanaz2, 23.0));
256  }
257 
258  sin2a2 = pow(sin(az2 * PI / 180.0), 2.0);
259  api = az / 180.0;
260 
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));
264 
265  if(incidence <= emission) {
266  q = p_photoOsr * munot / u0p0;
267  }
268  else {
269  q = p_photoOsr * mu / up0;
270  }
271 
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);
277  }
278  else {
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);
283  }
284 
285  rr1 = p_photoWh / 4.0 * u0p / (u0p + up) * ((1.0 + bg) * pg -
286  1.0 + Hfunc(u0p, gamma) * Hfunc(up, gamma));
287  rr2 = up * munot / (up0 * u0p0 * p_photoSr * (1.0 - faz + faz * q));
288  pht_hapke = rr1 * rr2;
289 
290  return pht_hapke;
291  }
292 
301  void Hapke::SetPhotoHg1(const double hg1) {
302  if(hg1 <= -1.0 || hg1 >= 1.0) {
303  string msg = "Invalid value of Hapke Henyey Greenstein hg1 [" +
304  IString(hg1) + "]";
305  throw IException(IException::User, msg, _FILEINFO_);
306  }
307  p_photoHg1 = hg1;
308  }
309 
318  void Hapke::SetPhotoHg2(const double hg2) {
319  if(hg2 < 0.0 || hg2 > 1.0) {
320  string msg = "Invalid value of Hapke Henyey Greenstein hg2 [" +
321  IString(hg2) + "]";
322  throw IException(IException::User, msg, _FILEINFO_);
323  }
324  p_photoHg2 = hg2;
325  }
326 
335  void Hapke::SetPhotoBh(const double bh) {
336  if(bh < -1.0 || bh > 1.0) {
337  string msg = "Invalid value of Hapke Legendre bh [" +
338  IString(bh) + "]";
339  throw IException(IException::User, msg, _FILEINFO_);
340  }
341  p_photoBh = bh;
342  }
343 
352  void Hapke::SetPhotoCh(const double ch) {
353  if(ch < -1.0 || ch > 1.0) {
354  string msg = "Invalid value of Hapke Legendre ch [" +
355  IString(ch) + "]";
356  throw IException(IException::User, msg, _FILEINFO_);
357  }
358  p_photoCh = ch;
359  }
360 
368  void Hapke::SetPhotoWh(const double wh) {
369  if(wh <= 0.0 || wh > 1.0) {
370  string msg = "Invalid value of Hapke wh [" +
371  IString(wh) + "]";
372  throw IException(IException::User, msg, _FILEINFO_);
373  }
374  p_photoWh = wh;
375  }
376 
384  void Hapke::SetPhotoHh(const double hh) {
385  if(hh < 0.0) {
386  string msg = "Invalid value of Hapke hh [" +
387  IString(hh) + "]";
388  throw IException(IException::User, msg, _FILEINFO_);
389  }
390  p_photoHh = hh;
391  }
392 
400  void Hapke::SetPhotoB0(const double b0) {
401  if(b0 < 0.0) {
402  string msg = "Invalid value of Hapke b0 [" +
403  IString(b0) + "]";
404  throw IException(IException::User, msg, _FILEINFO_);
405  }
406  p_photoB0 = b0;
407  }
408 
416  void Hapke::SetPhoto0B0Standard(const QString &b0standard) {
417  IString temp(b0standard);
418  temp = temp.UpCase();
419 
420  if(temp != "NO" && temp != "YES" && temp != "FALSE" && temp != "TRUE") {
421  string msg = "Invalid value of Hapke ZeroB0Standard [" +
422  temp + "]";
423  throw IException(IException::User, msg, _FILEINFO_);
424  }
425  if (temp == "NO" || temp == "FALSE") p_photo0B0Standard = "FALSE";
426  if (temp == "YES" || temp == "TRUE") p_photo0B0Standard = "TRUE";
427  }
428 
435  void Hapke::SetPhotoTheta(const double theta) {
436  if(theta < 0.0 || theta > 90.0) {
437  string msg = "Invalid value of Hapke theta [" +
438  IString(theta) + "]";
439  throw IException(IException::User, msg, _FILEINFO_);
440  }
441  p_photoTheta = theta;
442  }
443 
444 
445  void Hapke::SetStandardConditions(bool standard) {
446  PhotoModel::SetStandardConditions(standard);
447 
448  if(standard) {
449  p_photoB0save = p_photoB0;
450  if (p_photo0B0Standard == "TRUE") p_photoB0 = 0.0;
451  }
452  else {
453  p_photoB0 = p_photoB0save;
454  }
455  }
456 }
457 
458 extern "C" Isis::PhotoModel *HapkePlugin(Isis::Pvl &pvl) {
459  return new Isis::Hapke(pvl);
460 }