NonParametricXMLMieModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <fdet/FDetector.h>
#include <fdet/Eye.h>

#include <atm/NonParametricXMLMieModel.h>
#include <atm/ScatteringResult.h>
#include <atm/AttenuationResult.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/Point.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/Vector.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/RandomSamplerFromPDF.h>

#include <limits>

using namespace utl;
using namespace atm;
using namespace fwk;
using namespace std;


void
NonParametricXMLMieModel::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("NonParametricXMLMieModel");

  // Get model parameters for vertical distribution of aerosols.
  Branch profileB = topB.GetChild("aerosolProfile");

  vector<double> height;
  vector<double> vaod;

  profileB.GetChild("height").GetData(height);
  profileB.GetChild("opticalDepth").GetData(vaod);

  // Get phase function and wavelength dependence.
  topB.GetChild("phaseFunction").GetData(fPF);

  // Get wavelength dependence (assume Angstrom parameterization)
  Branch waveB = topB.GetChild("wavelengthDependence");
  waveB.GetChild("gamma").GetData(fGamma);
  waveB.GetChild("gammaError").GetData(fGammaError);

  // Fill the VAOD and extinction coefficient tables
  vector<double> alpha;
  const double dh = height[1] - height[0];
  const int N = height.size();

  // Define extinction @ top of altitude bin (CLF convention)
  if (vaod[0] > 0) {
    alpha.push_back(vaod[0] / dh);
    for (int i = 1; i < N; ++i)
      alpha.push_back((vaod[i] - vaod[i-1]) / dh);
  }
  // Ignore lowest bin where VAOD=0 (force the CLF convention)
  else {
    for (int i = 1; i < N; ++i)
      alpha.push_back((vaod[i] - vaod[i-1]) / dh);
    vaod.erase(vaod.begin());
    height.erase(height.begin());
  }

  fVAODVsHeight = new TabulatedFunction(height, vaod);
  fAlphaVsHeight = new TabulatedFunction(height, alpha);
}


void
NonParametricXMLMieModel::SetOpticalDepth(const utl::TabulatedFunction& vaod,
                                          const utl::TabulatedFunction& alpha)
{
  // Old pointers are automatically deleted when we do this
  fVAODVsHeight = new TabulatedFunction(vaod);
  fAlphaVsHeight = new TabulatedFunction(alpha);
}


void
NonParametricXMLMieModel::SetPhaseFunction(const PhaseFunction& phaseFunc)
{
  fPF = phaseFunc;
}


ScatteringResult
NonParametricXMLMieModel::EvaluateMieScattering(const utl::Point& xA,
                                             const utl::Point& xB,
                                             const double angle,
                                             const double distance,
                                             const vector<double>& wLength)
  const
{
  utl::TabulatedFunctionErrors frac;
  const double fractionError = 0;
  const double wError = 0;

  for (std::vector<double>::const_iterator it = wLength.begin();
       it != wLength.end(); ++it) {
    const double value = EvaluateMieScattering(xA, xB, angle, distance, *it);
    frac.PushBack(*it, wError, value, fractionError);
  }

  return ScatteringResult(frac, angle);
}


ScatteringResult
NonParametricXMLMieModel::EvaluateMieScattering(const utl::Point& xA,
                                             const utl::Point& xB,
                                             const double angle,
                                             const double distance,
                                             const AttenuationResult& mieAttenuation)
  const
{
  const utl::TabulatedFunctionErrors& attenuation = mieAttenuation.GetTransmissionFactor();
  utl::TabulatedFunctionErrors frac;
  const double fractionError = 0;
  const double wError = 0;

  const unsigned int nWl = attenuation.GetNPoints();
  for (unsigned int iWl = 0; iWl < nWl; ++iWl) {
    const double wl = attenuation.GetX(iWl);
    const double value = EvaluateMieScattering(xA, xB, angle, distance, wl, attenuation.GetY(iWl));
    frac.PushBack(wl, wError, value, fractionError);
  }

  return ScatteringResult(frac, angle);
}


double
NonParametricXMLMieModel::EvaluateMieScattering(const utl::Point& xA,
                                             const utl::Point& xB,
                                             const double angle,
                                             const double distance,
                                             const double wLength)
  const
{
  // Attenuation is needed in the scattering calculations.
  const double mAtt = EvaluateMieAttenuation(xA, xB, wLength);
  return EvaluateMieScattering(xA, xB, angle, distance, wLength, mAtt);
}


double
NonParametricXMLMieModel::EvaluateMieScattering(const utl::Point& xA,
                                                const utl::Point& /*xB*/,
                                                const double angle,
                                                const double distance,
                                                const double wLength,
                                                const double mieAttenuation)
  const
{
  const double xphs = EvaluateScatteringAngle(xA, angle, wLength);
  return (1 - mieAttenuation) * xphs / pow(distance, 2);
}


AttenuationResult
NonParametricXMLMieModel::EvaluateMieAttenuation(const utl::Point& xInit,
                                              const utl::Point& xFinal,
                                              const std::vector<double>& wLength)
  const
{
  utl::TabulatedFunctionErrors attenuation;
  const double wError = 0;
  const double mieAttError = 0;

  for (std::vector<double>::const_iterator wlIter = wLength.begin();
       wlIter != wLength.end(); ++wlIter) {
    const double mieAtt = EvaluateMieAttenuation(xInit, xFinal, *wlIter);
    attenuation.PushBack(*wlIter, wError, mieAtt, mieAttError);
  }

  return AttenuationResult(attenuation);
}


double
NonParametricXMLMieModel::EvaluateMieAttenuation(const utl::Point& xInit,
                                              const utl::Point& xFinal,
                                              const double wLength)
  const
{
  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  double hInit = wgs84.PointToLatitudeLongitudeHeight(xInit).get<2>();
  double hFinal = wgs84.PointToLatitudeLongitudeHeight(xFinal).get<2>();

  const double hMin = fVAODVsHeight->XFront();
  const double hMax = fVAODVsHeight->XBack();

  const double distance = (xFinal - xInit).GetMag();
  const double deltaHeight = abs(hFinal - hInit);
  const double cosTheta = deltaHeight / distance;
  double averageHeight = min(hInit, hFinal) + 0.5 * deltaHeight;

  // Use closest measured VAOD point
  if (hInit > hFinal)
    swap(hInit, hFinal);

  if (hInit < hMin)
    hInit = hMin;
  if (hInit > hMax)
    hInit = hMax;
  if (hFinal < hMin)
    hFinal = hMin;
  if (hFinal > hMax)
    hFinal = hMax;
  if (averageHeight < hMin)
    averageHeight = hMin;
  if (averageHeight > hMax)
    averageHeight = hMax;

  double vaodSlant;

  // Estimate optical depth of slant path
  if (deltaHeight > 1 * m) {
    const double vaodInit = fVAODVsHeight->InterpolateY(hInit,3);
    const double vaodFinal = fVAODVsHeight->InterpolateY(hFinal,3);
    vaodSlant = abs(vaodFinal - vaodInit) / cosTheta;
  } else {
    const double alpha = fAlphaVsHeight->Y(hInit);
    vaodSlant = distance*alpha;
  }

  // Add the wavelength dependence.
  const double arg = -vaodSlant * pow(wLength/355./nanometer, -fGamma);

  if (arg < std::numeric_limits<double>::max_exponent10 * kLn10)
    return exp(arg);

  return std::numeric_limits<double>::max();
}


double
NonParametricXMLMieModel::GetVerticalAerosolOpticalDepth(const unsigned int eyeId,
                                                      const double altitude)
  const
{
  const double hMax = fVAODVsHeight->XBack();
  double VAOD;

  if (altitude >= hMax) {
    ostringstream warn;
    warn << "Requested VAOD at " << altitude / km << " km at eye " << eyeId
         << ", but VAOD profile only goes to " << hMax / km << "km."
         << " Returning heighest measured VAOD.";
    WARNING(warn);

    VAOD = fVAODVsHeight->YBack();
  }
  else
    VAOD =  fVAODVsHeight->InterpolateY(altitude, 3);

  return VAOD;
}


double
NonParametricXMLMieModel::GetAttenuationLength(const utl::Point& p,
                                               const double wLength)
  const
{
  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  double height = wgs84.PointToLatitudeLongitudeHeight(p).get<2>();

  const double hMin = fAlphaVsHeight->XFront();
  const double hMax = fAlphaVsHeight->XBack();

  // Use closest measured VAOD point
  if (height < hMin)
    height = hMin;
  if (height > hMax)
    height = hMax;

  const double alpha = fAlphaVsHeight->Y(height);
  return 1. / (alpha * pow(wLength / 355. / nanometer, -fGamma));
}


double
NonParametricXMLMieModel::EvaluateScatteringAngle(const Point& /*p*/,
                                               const double angle,
                                               const double /*wLength*/)
  const
{
  const int nbin = fPF.size();
  const double binSize = kPi / (nbin - 1);
  int ibin = int(angle / binSize);

  if (ibin < 0)
    ibin = 0;
  else if (ibin > nbin)
    ibin = nbin - 1;

  return exp(log(fPF[ibin]) +
    (angle - binSize * ibin) * log(fPF[ibin+1] / fPF[ibin]) / binSize);
}