MinnaertEmpirical.cpp

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

namespace Isis {
  MinnaertEmpirical::MinnaertEmpirical(Pvl &pvl) : PhotoModel(pvl) {
    PvlGroup &algo = pvl.findObject("PhotometricModel")
                     .findGroup("Algorithm", Pvl::Traverse);
    // There are no default values for the Minnaert Empirical function; if user
    // does not provide information, then an exception is thrown
    if (algo.hasKeyword("PhaseList")) {
      SetPhotoPhaseList(algo["PhaseList"]);
    } else {
      std::string msg = "The empirical Minnaert phase list was not provided by user";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    if (algo.hasKeyword("KList")) {
      SetPhotoKList(algo["KList"]);
    } else {
      std::string msg = "The empirical Minnaert k exponent list was not provided by user";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    if (algo.hasKeyword("PhaseCurveList")) {
      SetPhotoPhaseCurveList(algo["PhaseCurveList"]);
    } else {
      std::string msg = "The empirical Minnaert 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_photoKList.size()) {
      std::string msg = "Number of empirical Minnaert k 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()) {
      std::string msg = "Number of empirical Minnaert 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_photoKSpline.Reset();
    p_photoKSpline.SetInterpType(NumericalApproximation::CubicClamped);
    p_photoKSpline.AddData(p_photoPhaseList,p_photoKList);
    p_photoKSpline.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);
  }

  MinnaertEmpirical::~MinnaertEmpirical() {
    p_photoKSpline.Reset();
    p_photoBSpline.Reset();
    p_photoPhaseList.clear();
    p_photoKList.clear();
    p_photoPhaseCurveList.clear();
  }

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

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

  void MinnaertEmpirical::SetPhotoKList(QString kstrlist) {
    double kvalue;
    IString strlist(kstrlist);
    p_photoKList.clear();

    while (strlist.length()) {
      kvalue = strlist.Token(",");
      if (kvalue < 0.0) {
        std::string msg = "Invalid value of Minnaert k list value [" +
                          IString(kvalue) + "]";
        throw IException(IException::User, msg, _FILEINFO_);
      }
      p_photoKList.push_back(kvalue);
    }
  }

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

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

  double MinnaertEmpirical::PhotoModelAlgorithm(double phase, double incidence,
                                       double emission) {
    static double pht_minnaert_empirical;
    double incrad;
    double emarad;
    double munot;
    double mu;
    double kInterpolated = 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_minnaert_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) {
      kInterpolated = p_photoKSpline.Evaluate(phase, NumericalApproximation::Extrapolate);
      bInterpolated = p_photoBSpline.Evaluate(phase, NumericalApproximation::Extrapolate);
      old_phase = phase;
    }

    if(munot <= 0.0 || mu <= 0.0 || incidence == 90.0 ||
        emission == 90.0) {
      pht_minnaert_empirical = 0.0;
    }
    else if(kInterpolated == 1.0) {
      pht_minnaert_empirical = munot * bInterpolated;
    }
    else {
      pht_minnaert_empirical = bInterpolated * munot * pow((munot * mu), (kInterpolated - 1.0));
    }

    return pht_minnaert_empirical;
  }
}

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