Isis 3.0 Object Programmers' Reference

Home

---

NumericalAtmosApprox.cpp

Go to the documentation of this file.

#include <cmath>
#include "AtmosModel.h"
#include "NumericalAtmosApprox.h"
#include "NumericalApproximation.h"
#include "iException.h"
#include "iString.h"

using namespace std;
namespace Isis {
  double NumericalAtmosApprox::RombergsMethod (AtmosModel *am, IntegFunc sub, double a, double b) throw (iException &){
    // This method was derived from an algorithm in the text 
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 4.3 by Flannery, Press, Teukolsky, and Vetterling
    int maxits = 20;       // maximium number of iterations allowed to converge
    double dss = 0;        // error estimate for 
    double h[maxits+1];    // relative stepsizes for trap
    double trap[maxits+1]; // successive trapeziodal approximations
    double epsilon;        // desired fractional accuracy
    double epsilon2;       // desired fractional accuracy
    double ss;             // result

    epsilon = 1.0e-4;
    epsilon2 = 1.0e-6;

    h[0] = 1.0;
    try{
      NumericalApproximation interp(NumericalApproximation::PolynomialNeville);
      for (int i=0; i<maxits; i++) { 
        // i will determine number of trapezoidal partitions of area  
        // under curve for "integration" using refined trapezoidal rule
        trap[i] = RefineExtendedTrap(am,sub,a,b,trap[i],i+1); // validates data here
        if (i >= 4) {
          interp.AddData(5, &h[i-4], &trap[i-4]);
          ss = interp.Evaluate(0.0,NumericalApproximation::Extrapolate);
          dss = interp.PolynomialNevilleErrorEstimate()[0];
          interp.Reset();
          // we work only until our necessary accuracy is achieved
          if (fabs(dss) <= epsilon*fabs(ss)) return ss;
          if (fabs(dss) <= epsilon2) return ss;
        }
        trap[i+1] = trap[i];
        h[i+1] = 0.25 * h[i];
        // This is a key step:  the factor is 0.25d0 even though    
        // the stepsize is decreased by 0.5d0.  This makes the      
        // extraplolation a polynomial in h-squared as allowed      
        // by the equation from Numerical Recipes 4.2.1 pg.132, 
        // not just a polynomial in h.
      }
    }
    catch (iException e){ // catch error from RefineExtendedTrap, Constructor, Evaluate, PolynomialNevilleErrorEstimate
      throw e.Message(e.Type(),
                      "NumericalAtmosApprox::RombergsMethod() - Caught the following error: ",
                      _FILEINFO_);
    }
    throw iException::Message(iException::Programmer, 
                              "NumericalAtmosApprox::RombergsMethod() - Failed to converge in " 
                              + iString(maxits) + " iterations.", 
                              _FILEINFO_);
  }

  double NumericalAtmosApprox::RefineExtendedTrap(AtmosModel *am, IntegFunc sub, double a, double b, double s, unsigned int n) throw (iException &){
    // This method was derived from an algorithm in the text 
    // Numerical Recipes in C: The Art of Scientific Computing
    // Section 4.2 by Flannery, Press, Teukolsky, and Vetterling
    try{
      if (n == 1) {
        double begin;
        double end;
        if (sub == InnerFunction) {
          begin = InrFunc2Bint(am,a);
          end = InrFunc2Bint(am,b);
        }
        else {
          begin = OutrFunc2Bint(am,a);
          end = OutrFunc2Bint(am,b);
        }
        return (0.5 * (b - a) * (begin + end));
      }
      else {
        int it;
        double delta,tnm,x,sum;
        it = (int)(pow(2.0,(double)(n-2)));
        tnm = it;
        delta = (b - a) / tnm; // spacing of the points to be added
        x = a + 0.5 * delta;
        sum = 0.0;
        for (int i=0; i<it; i++) {
          if (sub == InnerFunction) {
            sum = sum + InrFunc2Bint(am,x);
          }
          else {
            sum = sum + OutrFunc2Bint(am,x);
          }
          x = x + delta;
        }
        return (0.5 * (s + (b - a) * sum / tnm));// replace s with refined value
      }
    }
    catch (iException e){ // catch exception from Evaluate()
      throw e.Message(e.Type(),
                      "NumericalAtmosApprox::RefineExtendedTrap() - Caught the following error: ",
                      _FILEINFO_);
    }
  }

  double NumericalAtmosApprox::OutrFunc2Bint(AtmosModel *am, double phi) {
    double result;
    NumericalAtmosApprox::IntegFunc sub;
    sub = NumericalAtmosApprox::InnerFunction;

    am->p_atmosPhi = phi;
    am->p_atmosCosphi = cos((PI/180.0)*phi);

    NumericalAtmosApprox qromb;
    try{
      result = qromb.RombergsMethod(am,sub,1.0e-6,1.0);
      return result;
    }
    catch (iException e){ // catch exception from RombergsMethod()
      throw e.Message(e.Type(),
                      "NumericalAtmosApprox::OutrFunc2Bint() - Caught the following error: ",
                      _FILEINFO_);
    }
  }

  double NumericalAtmosApprox::InrFunc2Bint(AtmosModel *am, double mu) {
    double ema;    //pass in mu, get emission angle
    double sine;   //sin(ema)
    double alpha; 
    double phase;  //angle between sun and camera
    double result;
    double phasefn;//Henyey-Greenstein phase fn
    double xx;     //temp var to calculate emunot, emu
    double emunot; //exp(-tau/munot)
    double emu;    //exp(-tau/mu)
    double tfac;   //factor that occurs in the integrals for transmission

    //  calculate the phase angle
    //  also calculate any of the other redundant parameters
    ema = acos(mu) * (180.0/PI);
    sine = sin(ema*(PI/180.0));
    if ((am->p_atmosAtmSwitch == 0) || (am->p_atmosAtmSwitch == 2)) { // Reflection phase <= 90
      alpha = am->p_atmosSini * sine * am->p_atmosCosphi +
          am->p_atmosMunot * mu;
    }
    else { // Transmission phase <= 90
      alpha = am->p_atmosSini * sine * am->p_atmosCosphi -
          am->p_atmosMunot * mu;
    }
    phase = acos(alpha) * (180.0/PI);
    //  now evaluate the integrand; all needed parameters
    //  have been hidden separately and passed to it in COMMON.
    if (am->p_atmosAtmSwitch == 0) {
      // integrand for hemispheric albedo
      result = mu * am->p_atmosPM->CalcSurfAlbedo(phase,am->p_atmosInc,ema);
    }
    else {
      phasefn = (1.0 - am->AtmosHga() * am->AtmosHga()) /pow((1.0+2.0*am->AtmosHga()*alpha+am->AtmosHga()*am->AtmosHga()),1.5);
      xx = -am->AtmosTau() / max(am->p_atmosMunot,1.0e-30);
      if (xx < -69.0) {
        emunot = 0.0;
      }
      else if (xx > 69.0) {
        emunot = 1.0e30;
      }
      else {
        emunot = exp(xx);
      }
      xx = -am->AtmosTau() / max(mu,1.0e-30);

      if (xx < -69.0) {
        emu = 0.0;
      }
      else if (xx > 69.0) {
        emu = 1.0e30;
      }
      else {
        emu = exp(xx);
      }
      if (mu == am->p_atmosMunot) {
        tfac = am->AtmosTau() * emunot / (am->p_atmosMunot * am->p_atmosMunot);
      }
      else {
        tfac = (emunot - emu) / (am->p_atmosMunot - mu);
      }
      if (am->p_atmosAtmSwitch == 1) {
        result = mu * (phasefn - 1.0) * tfac;
      }
      else if (am->p_atmosAtmSwitch == 2) {
        result = am->p_atmosMunot * mu * (phasefn - 1.0) * (1.0 - emunot * emu) / (am->p_atmosMunot + mu);
      }
      else if (am->p_atmosAtmSwitch == 3) {
        result = -sine * am->p_atmosCosphi * (phasefn - 1.0) * tfac;
      }
      else {
        string msg = "NumericalAtmosApprox::InrFunc2Bint() - Invalid value of atmospheric ";
        msg += "switch used as argument to this function";
        throw iException::Message(iException::Programmer,msg,_FILEINFO_);
      }
    }
    return result;
  }
}//end NumericalAtmosApprox.cpp