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_atmosHnorm = 0.003;
    p_pstd = 0.0;
    p_sbar = 0.0;
    p_trans = 0.0;
    p_trans0 = 0.0;
    p_transs = 0.0;
    p_standardConditions = false;

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

    if(algorithm.hasKeyword("Nulneg")) {
      SetAtmosNulneg(algorithm["Nulneg"][0] == "YES");
    }
    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("Hga")) {
      SetAtmosHga(algorithm["Hga"]);
    }
    p_atmosHgasave = p_atmosHga;

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

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

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

    if(algorithm.hasKeyword("Hnorm")) {
      SetAtmosHnorm(algorithm["Hnorm"]);
    }

    if(algorithm.hasKeyword("Iord")) {
      SetAtmosIord(algorithm["Iord"][0] == "YES");
    }
    else {
      p_atmosAddOffset = false;
    }

    if(algorithm.hasKeyword("EstTau")) {
      SetAtmosEstTau(algorithm["EstTau"][0] == "YES");
    }
    else {
      p_atmosEstTau = false;
    }
  }

  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) {
    //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(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(IException::Unknown,
                       "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) {
    //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((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(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(IException::Unknown,
                       "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(IException::Unknown,
                       "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,
                                 double *transs) {

    // 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;
    *transs = p_transs;
  }

  void AtmosModel::SetStandardConditions(bool standard) {
    p_standardConditions = standard;
    if(p_standardConditions) {
      p_atmosTausave = p_atmosTau;
      p_atmosTau = p_atmosTauref;
    }
    else {
      p_atmosTau = p_atmosTausave;
    }
  }

  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::GenerateHahgTablesShadow() {
    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
    int nincl = 19;     // used instead of p_atmosNinc
    double deltaInc;    // limits table size and increment used

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

    p_atmosHahgsb = 0.0;
    NumericalAtmosApprox qromb;
    deltaInc = 90.0/(double)(nincl-1.0);

    for(inccnt = 0; inccnt < nincl; inccnt++) {
      p_atmosInc = deltaInc * (double) inccnt;
      p_atmosMunot = cos((PI / 180.0) * p_atmosInc);
      p_atmosSini = sin((PI / 180.0) * p_atmosInc);

      p_atmosAtmSwitch = 2;

      qromb.Reset();
      if (p_atmosInc >= 90.0) {
        fun_temp = 0.0;
      } else {
        fun_temp = qromb.RombergsMethod(this, sub, 0, 180);
      }

      if((inccnt == 0) || (inccnt == nincl - 1)) {
        factor = 0.5;
      }
      else {
        factor = 1.0;
      }
      p_atmosHahgsb = p_atmosHahgsb + p_atmosSini * factor * fun_temp *
        AtmosWha() / 360.0;
    }

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

  void AtmosModel::SetAtmosAtmSwitch(const int atmswitch) {
    if(atmswitch < 0 || atmswitch > 3) {
      string msg = "Invalid value of atmospheric atmswitch [" + IString(atmswitch) + "]";
      throw IException(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(IException::User, msg, _FILEINFO_);
    }

    p_atmosBha = bha;
  }

  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(IException::User, msg, _FILEINFO_);
    }

    p_atmosHga = hga;
  }

  void AtmosModel::SetAtmosInc(const double inc) {
    if(inc < 0.0 || inc > 90.0) {
      string msg = "Invalid value of atmospheric inc [" + IString(inc) + "]";
      throw IException(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(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(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(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(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(IException::User, msg, _FILEINFO_);
    }

    p_atmosWha = wha;
  }

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

  void AtmosModel::SetAtmosHnorm(const double hnorm) {
    if(hnorm < 0.0) {
      QString msg = "Invalid value of Atmospheric hnorm [" + toString(hnorm) + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    p_atmosHnorm = hnorm;
  }

  void AtmosModel::SetAtmosIord(const string offset) {
    IString temp(offset);
    temp = temp.UpCase();

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

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

  void AtmosModel::SetAtmosEstTau(const string esttau) {
    IString temp(esttau);
    temp = temp.UpCase();

    if(temp != "NO" && temp != "YES") {
      string msg = "Invalid value of Atmospheric optical depth estimation[" + temp + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }

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