Photometry.cpp

#include "IException.h"
#include "IString.h"
#include "Pvl.h"
#include "Photometry.h"
#include "PhotoModelFactory.h"
#include "PhotoModel.h"
#include "AtmosModelFactory.h"
#include "AtmosModel.h"
#include "NormModelFactory.h"
#include "NormModel.h"
#include "Plugin.h"
#include "FileName.h"

namespace Isis {
  Photometry::Photometry(Pvl &pvl) {
    p_phtAmodel = NULL;
    p_phtPmodel = NULL;
    p_phtNmodel = NULL;
    if(pvl.hasObject("PhotometricModel")) {
      p_phtPmodel = PhotoModelFactory::Create(pvl);
    } else {
      std::string msg = "A Photometric model must be specified to do any type of photometry";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    if(pvl.hasObject("AtmosphericModel")) {
      p_phtAmodel = AtmosModelFactory::Create(pvl, *p_phtPmodel);
    }
    if (pvl.hasObject("NormalizationModel")) {
      if (p_phtAmodel != NULL) {
        p_phtNmodel = NormModelFactory::Create(pvl, *p_phtPmodel, *p_phtAmodel);
      } else {
        p_phtNmodel = NormModelFactory::Create(pvl, *p_phtPmodel);
      }
    }
  }

  Photometry::~Photometry() {
    if(p_phtAmodel != NULL) {
      delete p_phtAmodel;
      p_phtAmodel = NULL;
    }

    if(p_phtPmodel != NULL) {
      delete p_phtPmodel;
      p_phtPmodel = NULL;
    }

    if(p_phtNmodel != NULL) {
      delete p_phtNmodel;
      p_phtNmodel = NULL;
    }
  }

  void Photometry::SetPhotomWl(double wl) {
    p_phtNmodel->SetNormWavelength(wl);
  }

  void Photometry::Compute(double pha, double inc, double ema,
                           double dn, double &albedo, double &mult,
                           double &base) {

    // Calculate the surface brightness
    p_phtNmodel->CalcNrmAlbedo(pha, inc, ema, dn, albedo, mult, base);
    return;
  }

  void Photometry::Compute(double pha, double inc, double ema,
                           double deminc, double demema, double dn,
                           double &albedo, double &mult, double &base) {

    // Calculate the surface brightness
    p_phtNmodel->CalcNrmAlbedo(pha, inc, ema, deminc, demema, dn, albedo, mult, base);
    return;
  }

  int Photometry::brentsolver(double x_lo, double x_hi, gsl_function *Func, double tolerance, double &root){
    int status;
    int iter=0, max_iter=100;
    const gsl_root_fsolver_type *T;
    gsl_root_fsolver *s;

    T = gsl_root_fsolver_brent;
    s = gsl_root_fsolver_alloc (T);
    gsl_root_fsolver_set(s, Func, x_lo, x_hi);

    do {
      iter++;
      status = gsl_root_fsolver_iterate(s);
      root   = gsl_root_fsolver_x_lower(s);
      x_lo   = gsl_root_fsolver_x_lower(s);
      x_hi   = gsl_root_fsolver_x_upper(s);
      status = gsl_root_test_interval(x_lo, x_hi, 0, tolerance);

    } while (status != GSL_SUCCESS && iter < max_iter);

    gsl_root_fsolver_free(s);
    return status;
  }

  int Photometry::brentminimizer(double x_lower, double x_upper, gsl_function *Func,
      double & x_minimum, double tolerance) {
    int status;
    int iter=0, max_iter=100;

    const gsl_min_fminimizer_type *T;
    gsl_min_fminimizer *s;
    //double m_expected = M_PI;

    T = gsl_min_fminimizer_brent;
    s = gsl_min_fminimizer_alloc(T);

    // This function sets, or resets, an existing minimizer s to use the function Func and
    // the initial search interval [x_lower, x_upper], with a guess for the location of
    // the minimum x_minimum. If the interval given does not contain a minimum, then
    // the function returns an error code of GSL_EINVAL.
    gsl_min_fminimizer_set(s, Func, x_minimum, x_lower, x_upper);

    do {
      iter++;
      status    = gsl_min_fminimizer_iterate(s);
      x_minimum = gsl_min_fminimizer_x_minimum(s);
      x_lower   = gsl_min_fminimizer_x_lower(s);
      x_upper   = gsl_min_fminimizer_x_upper(s);

      status = gsl_min_test_interval(x_lower, x_upper, tolerance, 0.0);
    } while(status == GSL_CONTINUE && iter < max_iter);

    // This function frees all the memory associated with the minimizer s.
    gsl_min_fminimizer_free(s);

    return status;
  }

  void Photometry::minbracket(double &xa, double &xb, double &xc, double &fa, double &fb,
      double &fc, double Func(double par, void *params), void *params) {
    double eps = 1.0e-21;
    double Gold = 1.618034;
    double GrowLimit = 110;
    int maxiter = 1000;

    fa = Func(xa, params);
    fb = Func(xb, params);
    if (fa < fb) {
      double tmp = xa;
      xa = xb;
      xb = tmp;
      tmp = fa;
      fa = fb;
      fb = tmp;
    }
    xc = xb + Gold * (xb - xa);
    fc = Func(xc, params);
    int iter = 0;

    while (fc < fb) {
      double tmp1 = (xb - xa) * (fb - fc);
      double tmp2 = (xb - xc) * (fb - fa);
      double val = tmp2 - tmp1;
      double denom;
      if (fabs(val) < eps) {
        denom = 2.0 * eps;
      } else {
        denom = 2.0 * val;
      }
      double w = xb - ((xb -xc) * tmp2 - (xb - xa) * tmp1) / denom;
      double wlim = xb + GrowLimit * (xc - xb);
      if (iter > maxiter) {
        IString msg = "Maximum iterations exceeded in minimum bracketing ";
        msg += "algorithm (minbracket) - root cannot be bracketed";
        throw IException(IException::User, msg, _FILEINFO_);
      }
      iter = iter + 1;
      double fw;
      if (((w-xc)*(xb-w)) > 0.0) {
        fw = Func(w, params);
        if (fw < fc) {
          xa = xb;
          xb = w;
          fa = fb;
          fb = fw;
          return;
        } else if (fw > fb) {
          xc = w;
          fc = fw;
          return;
        }
        w = xc + Gold * (xc - xb);
        fw = Func(w, params);
      } else if (((w-wlim)*(wlim-xc)) >= 0.0) {
        w = wlim;
        fw = Func(w, params);
      } else if (((w-wlim)*(xc-w)) > 0.0) {
        fw = Func(w, params);
        if (fw < fc) {
          xb = xc;
          xc = w;
          w = xc + Gold * (xc - xb);
          fb = fc;
          fc = fw;
          fw = Func(w, params);
        }
      } else {
        w = xc + Gold * (xc - xb);
        fw = Func(w, params);
      }
      xa = xb;
      xb = xc;
      xc = w;
      fa = fb;
      fb = fc;
      fc = fw;
    }
    return;
  }
}