LunarLambertEmpirical.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#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::SetPhotoPhaseList(PvlKeyword phaseList) {
 
    // If the format is Keyword="1,2,3,4,5" rather than Keyword = (1,2,3,4,5) 
    if (phaseList.size() == 1) {
      SetPhotoPhaseList(QString(phaseList));
      return;
    }
 
    double phaseAngle;
    p_photoPhaseList.clear();
 
    for (int i=0; i< phaseList.size(); i++) {
      phaseAngle = phaseList[i].toDouble();
 
      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::SetPhotoLList(PvlKeyword lstrList) {
 
    // If the format is Keyword="1,2,3,4,5" rather than Keyword = (1,2,3,4,5) 
    if (lstrList.size() == 1) {
      SetPhotoLList(QString(lstrList));
      return;
    }
 
    p_photoLList.clear();
    for (int i=0; i<lstrList.size(); i++) {
      p_photoLList.push_back(lstrList[i].toDouble());
    }
  }
 
 
  void LunarLambertEmpirical::SetPhotoPhaseCurveList(QString phasecurvestrlist) {
    double phasecurve;
    IString strlist(phasecurvestrlist);
    p_photoPhaseCurveList.clear();
 
    while (strlist.length()) {
      phasecurve = strlist.Token(",");
      p_photoPhaseCurveList.push_back(phasecurve);
    }
  }
 
 
  void LunarLambertEmpirical::SetPhotoPhaseCurveList(PvlKeyword photocurvestrList) {
 
    // If the format is Keyword="1,2,3,4,5" rather than Keyword = (1,2,3,4,5) 
    if (photocurvestrList.size() == 1) {
      SetPhotoPhaseCurveList(QString(photocurvestrList));
      return;
    }
 
    p_photoPhaseCurveList.clear();
    for (int i=0; i<photocurvestrList.size(); i++) {
      p_photoPhaseCurveList.push_back(photocurvestrList[i].toDouble());
    }
  }
 
 
  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);
}