Isis 3.0 Object Programmers' Reference

Home

---

AtmosModel.cpp

#include <cmath>
#include <string>
#include "Pvl.h"
#include "iString.h"
#include "AtmosModel.h"
#include "NumericalApproximation.h"
#include "NumericalAtmosApprox.h"
#include "PhotoModel.h"
#include "Minnaert.h"
#include "LunarLambert.h"
#include "Plugin.h"
#include "iException.h"
#include "Filename.h"

using namespace std;

#define MAX(x,y) (((x) > (y)) ? (x) : (y))
namespace Isis {
  AtmosModel::AtmosModel (Pvl &pvl, PhotoModel &pmodel) {
    p_atmosAlgorithmName = "Unknown";
    p_atmosPM = &pmodel;

    p_atmosIncTable.resize(91);
    p_atmosAhTable.resize(91);
    p_atmosHahgtTable.resize(91);
    p_atmosHahgt0Table.resize(91);
    p_atmosAb = 0.0;
    p_atmosCosphi = 0.0;
    p_atmosAtmSwitch = 0;
    p_atmosEulgam = 0.5772156;
    p_atmosHahgsb = 0.0;
    p_atmosHga = 0.68;
    p_atmosInc = 0.0;
    p_atmosMunot = 0.0;
    p_atmosNinc = 91;
    p_atmosPhi = 0.0;
    p_atmosSini = 0.0;
    p_atmosTau = 0.28;
    p_atmosTauref = 0.0;
    p_atmosTauold = -1.0;
    p_atmosWha = 0.95;
    p_atmosWhaold = 1.0e30;
    p_atmosBha = 0.85;
    p_pstd = 0.0;
    p_sbar = 0.0;
    p_trans = 0.0;
    p_trans0 = 0.0;
    p_standardConditions = false;

    PvlGroup &algorithm = pvl.FindObject("AtmosphericModel").FindGroup("Algorithm",Pvl::Traverse);

    if (algorithm.HasKeyword("Nulneg")) {
      SetAtmosNulneg( ((string)algorithm["Nulneg"]).compare("YES") == 0 );
    }
    else {
      p_atmosNulneg = false;
    }

    if (algorithm.HasKeyword("Tau")) {
      SetAtmosTau( algorithm["Tau"] );
    }
    p_atmosTausave = p_atmosTau;

    if (algorithm.HasKeyword("Tauref")) {
      SetAtmosTauref( algorithm["Tauref"] );
    }

    if (algorithm.HasKeyword("Wha")) {
      SetAtmosWha( algorithm["Wha"] );
    }
    p_atmosWhasave = p_atmosWha;

    if (algorithm.HasKeyword("Wharef")) {
      SetAtmosWharef( algorithm["Wharef"] );
    } else {
      p_atmosWharef = p_atmosWha;
    }

    if (algorithm.HasKeyword("Hga")) {
      SetAtmosHga( algorithm["Hga"] );
    }
    p_atmosHgasave = p_atmosHga;

    if (algorithm.HasKeyword("Hgaref")) {
      SetAtmosHgaref( algorithm["Hgaref"] );
    } else {
      p_atmosHgaref = p_atmosHga;
    }

    if (algorithm.HasKeyword("Bha")) {
      SetAtmosBha( algorithm["Bha"] );
    }
    p_atmosBhasave = p_atmosBha;

    if (algorithm.HasKeyword("Bharef")) {
      SetAtmosBharef( algorithm["Bharef"] );
    } else {
      p_atmosBharef = p_atmosBha;
    }

    if (algorithm.HasKeyword("Inc")) {
      SetAtmosInc( algorithm["Inc"] );
    }

    if (algorithm.HasKeyword("Phi")) {
      SetAtmosPhi( algorithm["Phi"] );
    }
  }


  double AtmosModel::G11Prime(double tau) {
    double sum;
    int icnt;
    double fac;
    double term;
    double tol;
    double fi;
    double elog;
    double eulgam;
    double e1_2;
    double result;

    tol = 1.0e-6;
    eulgam = 0.5772156;

    sum = 0.0;
    icnt = 1;
    fac = -tau;
    term = fac;
    while (fabs(term) > fabs(sum)*tol) {
      sum = sum + term;
      icnt = icnt + 1;
      fi = (double) icnt;
      fac = fac * (-tau) / fi;
      term = fac / (fi * fi);
    }
    elog = log(MAX(1.0e-30,tau)) + eulgam;
    e1_2 = sum + PI * PI / 12.0 + 0.5 *
        pow(elog,2.0);
    result = 2.0 * (AtmosModel::En(1,tau) 
                    + elog * AtmosModel::En(2,tau) 
                    - tau * e1_2);
    return(result);
  } 

  double AtmosModel::Ei (double x) throw (iException &){
    //static double r8ei(double x) in original NumericalMethods class
    // This method was derived from an algorithm in the text 
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 6.3 by Flannery, Press, Teukolsky, and Vetterling
    double result;
    double fpmin; // close to smallest representable floating-point number
    double euler; // Euler's constant, lower-case gamma
    double epsilon;   // desired relative error (tolerance)
    int maxit;    // maximum number of iterations allowed
    double sum,fact,term,prev;

    fpmin = 1.0e-30;
    maxit = 100;
    epsilon = 6.0e-8;
    euler = 0.57721566;

    if (x <= 0.0) {
      throw iException::Message(iException::Programmer,
                                "AtmosModel::Ei() - Invalid argument. Definition requires x > 0.0. Entered x = "
                                + iString(x),
                                _FILEINFO_);
    }
    if (x < fpmin) { //special case: avoid failure of convergence test because underflow
      result = log(x) + euler;
    } 
    else if (x <= -log(epsilon)) { //Use power series
      sum = 0.0;
      fact = 1.0;
      for (int k=1; k<=maxit; k++) {
        fact = fact * x / k;
        term = fact / k;
        sum = sum + term;
        if (term < epsilon*sum) {
          result = sum + log(x) + euler;
          return(result);
        }
      }
      throw iException::Message(iException::Math,
                                "AtmosModel::Ei() - Power series failed to converge in " 
                                + iString(maxit) + " iterations. Unable to calculate exponential integral.",
                                _FILEINFO_);
    } 
    else { // Use asymptotic series
      sum = 0.0;
      term = 1.0;
      for (int k=1; k<=maxit; k++) {
        prev = term;
        term = term * k / x;
        if (term < epsilon) {
          result = exp(x) * (1.0 + sum) / x;
          return(result);
        } 
        else {
          if (term < prev) { // still converging: add new term
            sum = sum + term;
          } 
          else { // diverging: subtract previous term and exit
            sum = sum - prev;
            result = exp(x) * (1.0 + sum) / x;
            return(result);
          }
        }
      }
      result = exp(x) * (1.0 + sum) / x;
    }
    return(result);
  }

  double AtmosModel::En (unsigned int n, double x) throw (iException &){
    //static double r8expint(int n, double x) in original NumericalMethods class
    // This method was derived from an algorithm in the text 
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 6.3 by Flannery, Press, Teukolsky, and Vetterling
    int nm1;
    double result;
    double a,b,c,d,h;
    double delta;
    double fact;
    double psi;
    double fpmin;  // close to smallest representable floating-point number 
    double maxit;  // maximum number of iterations allowed                  
    double epsilon;    // desired relative error (tolerance)                    
    double euler;  // Euler's constant, gamma                               

    fpmin = 1.0e-30;
    maxit = 100;
    epsilon = 1.0e-7;
    euler = 0.5772156649;

    nm1 = n - 1;
    if (n < 0 || x < 0.0 || (x == 0.0 && (n == 0 || n == 1))) {
      iString msg = "AtmosModel::En() - Invalid arguments. ";
        msg+= "Definition requires (x > 0.0 and n >=0 ) or (x >= 0.0 and n > 1). ";
        msg+= "Entered x = " + iString(x) + " and n = " + iString((int) n);
      throw iException::Message(iException::Programmer,msg,
                                _FILEINFO_);
    } 
    else if (n == 0) { // special case, this implies x > 0 by logic above
      result = exp(-x) / x;
    } 
    else if (x == 0.0) { // special case, this implies n > 1
      result = 1.0 / nm1;
    } 
    else if (x > 1.0) { // Lentz's algorithm, this implies n > 0
      b = x + n;
      c = 1.0 / fpmin;
      d = 1.0 / b;
      h = d;
      for (int i=1; i<=maxit; i++) {
        a = -i * (nm1 + i);
        b = b + 2.0;
        d = 1.0 / (a * d + b);
        c = b + a / c;
        delta = c * d;
        h = h * delta;
        if (fabs(delta-1.0) < epsilon) {
          result = h * exp(-x);
          return(result);
        }
      }
      throw iException::Message(iException::Math,
                                "AtmosModel::En() - Continued fraction failed to converge in "
                                + iString(maxit) + " iterations. Unable to calculate exponential integral.", 
                                _FILEINFO_);
    } 
    else { // evaluate series
      if (nm1 != 0) {
        result = 1.0 / nm1;
      } 
      else {
        result = -log(x) - euler;
      }
      fact = 1.0;
      for (int i=1; i<=maxit; i++) {
        fact = -fact * x / i;
        if (i != nm1) {
          delta = -fact / (i - nm1);
        } 
        else {
          psi = -euler;
          for (int ii=1; ii<=nm1; ii++) {
            psi = psi + 1.0 / ii;
          }
          delta = fact * (-log(x) + psi);
        }
        result = result + delta;
        if (fabs(delta) < fabs(result)*epsilon) {
          return(result);
        }
      }
      throw iException::Message(iException::Math,
                                "AtmosModel::En() - Series representation failed to converge in "
                                + iString(maxit) + " iterations. Unable to calculate exponential integral.",
                                _FILEINFO_);
    }
    return(result);
  }

  void AtmosModel::CalcAtmEffect(double pha, double inc, double ema, 
        double *pstd, double *trans, double *trans0, double *sbar) {

    // Check validity of photometric angles
    //if (pha > 180.0 || inc > 90.0 || ema > 90.0 || pha < 0.0 ||
    //    inc < 0.0 || ema < 0.0) {
    //  string msg = "Invalid photometric angles";
    //  throw iException::Message(iException::Programmer,msg,_FILEINFO_);
    //}

    // Apply atmospheric function
    AtmosModelAlgorithm(pha,inc,ema);
    *pstd = p_pstd;
    *trans = p_trans;
    *trans0 = p_trans0;
    *sbar = p_sbar;
  }

  void AtmosModel::SetStandardConditions(bool standard) {
    p_standardConditions = standard;
    if (p_standardConditions) {
      p_atmosBhasave = p_atmosBha;
      p_atmosBha = p_atmosBharef;
      p_atmosHgasave = p_atmosHga;
      p_atmosHga = p_atmosHgaref;
      p_atmosTausave = p_atmosTau;
      p_atmosTau = p_atmosTauref;
      p_atmosWhasave = p_atmosWha;
      p_atmosWha = p_atmosWharef;
    } else {
      p_atmosBha = p_atmosBhasave;
      p_atmosHga = p_atmosHgasave;
      p_atmosTau = p_atmosTausave;
      p_atmosWha = p_atmosWhasave;
    }
  }

  void AtmosModel::GenerateAhTable() { 
    int inccnt;     //for loop incidence angle count, 1:ninc
    double fun_temp;//temporary variable for integral
    double factor;  //factor for integration: 1 except at ends of interval where it's 1/2
    double yp1,yp2; //first derivatives of first and last x values of @a p_atmosIncTable

    NumericalAtmosApprox::IntegFunc sub;
    sub = NumericalAtmosApprox::OuterFunction;

    p_atmosAb = 0.0;

    //Create NumericalAtmosApprox object here for RombergsMethod used in for loop
    NumericalAtmosApprox qromb;

    for (inccnt=0; inccnt<p_atmosNinc; inccnt++) {
      p_atmosInc = (double) inccnt;
      p_atmosIncTable[inccnt] = p_atmosInc;
      p_atmosMunot = cos((PI/180.0)*p_atmosInc);
      p_atmosSini = sin((PI/180.0)*p_atmosInc);

      iString phtName = p_atmosPM->AlgorithmName();
      phtName.UpCase();
      if (p_atmosInc == 90.0) {
        p_atmosAhTable[inccnt] = 0.0;
      }
      else if (phtName == "LAMBERT") {
        p_atmosAhTable[inccnt] = 1.0;
      }
      else if (phtName == "LOMMELSEELIGER") {
        p_atmosAhTable[inccnt] = 2.0 * log((1.0/p_atmosMunot)/p_atmosMunot);
      }
      else if (phtName == "MINNAERT") {
        p_atmosAhTable[inccnt] = (pow(p_atmosMunot,((Minnaert*)p_atmosPM)->PhotoK()))/((Minnaert*)p_atmosPM)->PhotoK();
      }
      else if (phtName == "LUNARLAMBERT") {
        p_atmosAhTable[inccnt] = 2.0 * ((LunarLambert*)p_atmosPM)->PhotoL() * 
                                 log((1.0+p_atmosMunot)/p_atmosMunot) + 1.0 - 
                                 ((LunarLambert*)p_atmosPM)->PhotoL();
      }
      else {
        // numerically integrate the other photometric models    
        // outer integral is over phi from 0 to pi rad = 180 deg 
        p_atmosAtmSwitch = 0;
        // integrate AtmosModel function from 0 to 180
        fun_temp = qromb.RombergsMethod(this,sub,0,180);
        // the correct normalization with phi in degrees is:
        p_atmosAhTable[inccnt] = fun_temp / (90.0*p_atmosMunot);
      }
      // Let's get our estimate of Ab by integrating (summing)
      // A (i)sinicosi over our A  table                      
      if ((inccnt == 0) || (inccnt == p_atmosNinc-1)) {
        factor = 0.5;
      }
      else {
        factor = 1.0;
      }

      p_atmosAb = p_atmosAb + p_atmosAhTable[inccnt] * p_atmosMunot * p_atmosSini * factor;
    }

    factor = 2.0 * PI/180.0;
    p_atmosAb = p_atmosAb * factor;

    // set up clamped cubic spline
    yp1 = 1.0e30;
    yp2 = 1.0e30;
    p_atmosAhSpline.Reset();
    p_atmosAhSpline.SetInterpType(NumericalApproximation::CubicClamped);
    p_atmosAhSpline.AddData(p_atmosIncTable,p_atmosAhTable);
    p_atmosAhSpline.SetCubicClampedEndptDeriv(yp1,yp2);
  }

  void AtmosModel::GenerateHahgTables() { 
    int inccnt;         // for loop incidence angle count,1:ninc
    double fun_temp;    // temporary variable for integral
    double hahgsb_temp; // save increment to hahgsb
    double factor;      // factor for integration: 1 except at ends of interval where it's 1/2
    double yp1,yp2;     // derivatives of endpoints of data set

    NumericalAtmosApprox::IntegFunc sub;
    sub = NumericalAtmosApprox::OuterFunction;

    p_atmosHahgsb = 0.0;
    NumericalAtmosApprox qromb;

    for (inccnt = 0; inccnt < p_atmosNinc; inccnt++) {
      p_atmosInc = (double) inccnt;
      p_atmosIncTable[inccnt] = p_atmosInc;
      p_atmosMunot = cos((PI/180.0)*p_atmosInc);
      p_atmosSini = sin((PI/180.0)*p_atmosInc);

      p_atmosAtmSwitch = 1;

      qromb.Reset();
      fun_temp = qromb.RombergsMethod(this,sub,0,180);

      p_atmosHahgtTable[inccnt] = fun_temp * AtmosWha() / 360.0;
      p_atmosAtmSwitch = 2;

      fun_temp = qromb.RombergsMethod(this,sub,0,180);

      hahgsb_temp = fun_temp * AtmosWha() / 360.0;

      if ((inccnt == 0) || (inccnt == p_atmosNinc-1)) {
        factor = 0.5;
      }
      else {
        factor = 1.0;
      }

      p_atmosHahgsb = p_atmosHahgsb + p_atmosSini * factor * hahgsb_temp;
      if (p_atmosInc == 0.0) {
        p_atmosHahgt0Table[inccnt] = 0.0;
      }
      else {
        p_atmosAtmSwitch = 3;
        fun_temp = qromb.RombergsMethod(this,sub,0,180);
        p_atmosHahgt0Table[inccnt] = fun_temp * AtmosWha() * p_atmosMunot / (360.0 * p_atmosSini);
      }
    }

    factor = 2.0 * PI/180.0;
    p_atmosHahgsb = p_atmosHahgsb * factor;

    // set up clamped cubic splines
    yp1 = 1.0e30;
    yp2 = 1.0e30;
    p_atmosHahgtSpline.Reset();
    p_atmosHahgtSpline.SetInterpType(NumericalApproximation::CubicClamped);
    p_atmosHahgtSpline.AddData(p_atmosIncTable,p_atmosHahgtTable); 
    p_atmosHahgtSpline.SetCubicClampedEndptDeriv(yp1,yp2);

    p_atmosHahgt0Spline.Reset();
    p_atmosHahgt0Spline.SetInterpType(NumericalApproximation::CubicClamped);
    p_atmosHahgt0Spline.AddData(p_atmosIncTable,p_atmosHahgt0Table);
    p_atmosHahgt0Spline.SetCubicClampedEndptDeriv(yp1,yp2);
  }

  void AtmosModel::SetAtmosAtmSwitch (const int atmswitch) {
    if (atmswitch < 0 || atmswitch > 3) {
      string msg = "Invalid value of atmospheric atmswitch [" + iString(atmswitch) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosAtmSwitch = atmswitch;
  }

  void AtmosModel::SetAtmosBha (const double bha) {
    if (bha < -1.0 || bha > 1.0) {
      string msg = "Invalid value of Anisotropic atmospheric bha [" +
                        iString(bha) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosBha = bha;
  }

  void AtmosModel::SetAtmosBharef (const double bharef) {
    if (bharef < -1.0 || bharef > 1.0) {
      string msg = "Invalid value of Atmospheric bharef [" + iString(bharef) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosBharef = bharef;
  }

  void AtmosModel::SetAtmosHga (const double hga) {
    if (hga <= -1.0 || hga >= 1.0) {
      string msg = "Invalid value of Hapke atmospheric hga [" + iString(hga) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosHga = hga;
  }

  void AtmosModel::SetAtmosHgaref (const double hgaref) {
    if (hgaref <= -1.0 || hgaref >= 1.0) {
      string msg = "Invalid value of Atmospheric hgaref [" + iString(hgaref) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosHgaref = hgaref;
  }

  void AtmosModel::SetAtmosInc (const double inc) {
    if (inc < 0.0 || inc > 90.0) {
      string msg = "Invalid value of atmospheric inc [" + iString(inc) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosInc = inc;
    p_atmosMunot = cos(inc*PI/180.0);
    p_atmosSini = sin(inc*PI/180.0);
  }

  void AtmosModel::SetAtmosNulneg (const string nulneg) {
    iString temp(nulneg);
    temp = temp.UpCase();

    if (temp != "NO" && temp != "YES") {
      string msg = "Invalid value of Atmospheric nulneg [" + temp + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    SetAtmosNulneg( temp.compare("YES") == 0 );
  }

  void AtmosModel::SetAtmosPhi (const double phi) {
    if (phi < 0.0 || phi > 360.0) {
      string msg = "Invalid value of atmospheric phi [" + iString(phi) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosPhi = phi;
    p_atmosCosphi = cos(phi*PI/180.0);
  }

  void AtmosModel::SetAtmosTau (const double tau) {
    if (tau < 0.0) {
      string msg = "Invalid value of Atmospheric tau [" + iString(tau) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosTau = tau;
  }

  void AtmosModel::SetAtmosTauref (const double tauref) {
    if (tauref < 0.0) {
      string msg = "Invalid value of Atmospheric tauref [" + iString(tauref) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosTauref = tauref;
  }

  void AtmosModel::SetAtmosWha (const double wha) {
    if (wha <= 0.0 || wha > 1.0) {
      string msg = "Invalid value of Atmospheric wha [" + iString(wha) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosWha = wha;
  }

  void AtmosModel::SetAtmosWharef (const double wharef) {
    if (wharef <= 0.0 || wharef > 1.0) {
      string msg = "Invalid value of Atmospheric wharef [" + iString(wharef) + "]";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    p_atmosWharef = wharef;
  }

  bool AtmosModel::TauOrWhaChanged() const {
    return (AtmosTau() != p_atmosTauold) || (AtmosWha() != p_atmosWhaold);
  }
}