LunarLambertEmpirical.cpp

#include <cmath>
#include "LunarLambertEmpirical.h"
#include "IException.h"

namespace Isis {
  LunarLambertEmpirical::LunarLambertEmpirical(Pvl &pvl) : PhotoModel(pvl) {
    PvlGroup &algo = pvl.findObject("PhotometricModel")
                     .findGroup("Algorithm", Pvl::Traverse);
    // There are no default values for the Lunar Lambert Empirical function; if user
    // does not provide information, then an exception is thrown
    if (algo.hasKeyword("PhaseList")) {
      SetPhotoPhaseList(algo["PhaseList"]);
    } else {
      QString msg = "The empirical Lunar Lambert phase list was not provided by user";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    if (algo.hasKeyword("LList")) {
      SetPhotoLList(algo["LList"]);
    } else {
      QString msg = "The empirical Lunar Lambert l exponent list was not provided by user";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    if (algo.hasKeyword("PhaseCurveList")) {
      SetPhotoPhaseCurveList(algo["PhaseCurveList"]);
    } else {
      QString msg = "The empirical Lunar Lambert phase brightness list was not provided by user";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    // Make sure all the vectors are the same size
    p_photoPhaseAngleCount = (int)p_photoPhaseList.size();

    if (p_photoPhaseAngleCount != (int)p_photoLList.size()) {
      QString msg = "Number of empirical Lunar Lambert l list values must be equal";
      msg += "to number of phase angles provided";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    if (p_photoPhaseAngleCount != (int)p_photoPhaseCurveList.size()) {
      QString msg = "Number of empirical Lunar Lambert phase curve list values must be equal";
      msg += "to number of phase angles provided";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    // Create Cubic piecewise linear splines
    p_photoLSpline.Reset();
    p_photoLSpline.SetInterpType(NumericalApproximation::CubicClamped);
    p_photoLSpline.AddData(p_photoPhaseList,p_photoLList);
    p_photoLSpline.SetCubicClampedEndptDeriv(1.0e30,1.0e30);
    p_photoBSpline.Reset();
    p_photoBSpline.SetInterpType(NumericalApproximation::CubicClamped);
    p_photoBSpline.AddData(p_photoPhaseList,p_photoPhaseCurveList);
    p_photoBSpline.SetCubicClampedEndptDeriv(1.0e30,1.0e30);
  }

  LunarLambertEmpirical::~LunarLambertEmpirical() {
    p_photoLSpline.Reset();
    p_photoBSpline.Reset();
    p_photoPhaseList.clear();
    p_photoLList.clear();
    p_photoPhaseCurveList.clear();
  }

  void LunarLambertEmpirical::SetPhotoPhaseList(QString phasestrlist) {
    double phaseangle;
    IString strlist(phasestrlist);
    p_photoPhaseList.clear();

    while (strlist.length()) {
      phaseangle = strlist.Token(",");
      if (phaseangle < 0.0 || phaseangle > 180.0) {
        QString msg = "Invalid value of empirical Lunar Lambert phase angle list value [" +
                          toString(phaseangle) + "]";
        throw IException(IException::User, msg, _FILEINFO_);
      }
      p_photoPhaseList.push_back(phaseangle);
    }
  }

  void LunarLambertEmpirical::SetPhotoLList(QString lstrlist) {
    double lvalue;
    IString strlist(lstrlist);
    p_photoLList.clear();

    while (strlist.length()) {
      lvalue = strlist.Token(",");
      p_photoLList.push_back(lvalue);
    }
  }

  void LunarLambertEmpirical::SetPhotoPhaseCurveList(QString phasecurvestrlist) {
    double phasecurve;
    IString strlist(phasecurvestrlist);
    p_photoPhaseCurveList.clear();

    while (strlist.length()) {
      phasecurve = strlist.Token(",");
      p_photoPhaseCurveList.push_back(phasecurve);
    }
  }

  double LunarLambertEmpirical::PhotoModelAlgorithm(double phase, double incidence,
                                       double emission) {
    static double pht_lunarlambert_empirical;
    double incrad;
    double emarad;
    double munot;
    double mu;
    double lInterpolated = 0;
    double bInterpolated = 0;

    static double old_phase = -9999;
    static double old_incidence = -9999;
    static double old_emission= -9999;

    // Don't need to do anything if the photometric angles are the same as before
    if (old_phase == phase && old_incidence == incidence && old_emission == emission) {
      return pht_lunarlambert_empirical;
    }

    old_incidence = incidence;
    old_emission = emission;

    incrad = incidence * Isis::PI / 180.0;
    emarad = emission * Isis::PI / 180.0;
    munot = cos(incrad);
    mu = cos(emarad);

    if (phase != old_phase) {
      lInterpolated = p_photoLSpline.Evaluate(phase, NumericalApproximation::Extrapolate);
      bInterpolated = p_photoBSpline.Evaluate(phase, NumericalApproximation::Extrapolate);
      old_phase = phase;
    }

    if(munot <= 0.0 || mu <= 0.0) {
      pht_lunarlambert_empirical = 0.0;
    }
    else if(lInterpolated == 0.0) {
      pht_lunarlambert_empirical = munot * bInterpolated;
    }
    else if(lInterpolated == 1.0) {
      pht_lunarlambert_empirical = bInterpolated * 2.0 * munot / (munot + mu);
    }
    else {
      pht_lunarlambert_empirical = bInterpolated * munot * ((1.0 - lInterpolated) + 2.0 * lInterpolated / (munot + mu));
    }

    return pht_lunarlambert_empirical;
  }
}

extern "C" Isis::PhotoModel *LunarLambertEmpiricalPlugin(Isis::Pvl &pvl) {
  return new Isis::LunarLambertEmpirical(pvl);
}