NumericalAtmosApprox.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#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) {
    // 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 IException(e,
                       e.errorType(),
                      "NumericalAtmosApprox::RombergsMethod() - Caught the following error: ",
                      _FILEINFO_);
    }
    throw IException(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) {
    // 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 IException(e,
                       e.errorType(),
                      "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 IException(e,
                       e.errorType(),
                      "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(IException::Programmer, msg, _FILEINFO_);
      }
    }
    return result;
  }
}//end NumericalAtmosApprox.cpp