ParametricXMLMieModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

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

#include <atm/ParametricXMLMieModel.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/TabulatedFunction.h>
#include <utl/RandomSamplerFromPDF.h>

#include <algorithm>
#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#include <limits>

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


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

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

  profileB.GetChild("attLength").GetData(fLMix);
  profileB.GetChild("attLengthError").GetData(fLMixError);

  profileB.GetChild("mixHeight").GetData(fHMix);
  profileB.GetChild("mixHeightError").GetData(fHMixError);

  profileB.GetChild("scaleHeight").GetData(fHScl);
  profileB.GetChild("scaleHeightError").GetData(fHSclError);

  // Get phase function.
  Branch pfB = topB.GetChild("phaseFunction");
  pfB.GetChild("f").GetData(fPhaseFunction.f);
  pfB.GetChild("fError").GetData(fPhaseFunction.fError);
  pfB.GetChild("g").GetData(fPhaseFunction.g);
  pfB.GetChild("gError").GetData(fPhaseFunction.gError);

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

  // Set prevailing ground height (use PampaAmarilla as reference coordinate system).
  // obtaining the reference coordinate system from the detector does not work,
  // because in this Init() funciton there should not be any GetInstance() of the detector.
  const utl::CoordinateSystemPtr referenceCS = fwk::CoordinateSystemRegistry::Get("PampaAmarilla");

  const utl::Point origin(0,0,0, referenceCS);
  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  fHGnd = wgs84.PointToLatitudeLongitudeHeight(origin).get<2>();
}


ScatteringResult
ParametricXMLMieModel::EvaluateMieScattering(const utl::Point& xA,
                                             const utl::Point& xB,
                                             const double angle,
                                             const double distance,
                                             const vector<double>& wLength)
  const
{
  const AttenuationResult& mieAttenuation = EvaluateMieAttenuation(xA, xB, wLength);
  return EvaluateMieScattering(xA, xB, angle, distance, mieAttenuation);
}


ScatteringResult
ParametricXMLMieModel::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
ParametricXMLMieModel::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
ParametricXMLMieModel::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
ParametricXMLMieModel::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
ParametricXMLMieModel::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>() - fHGnd;
  double hFinal = wgs84.PointToLatitudeLongitudeHeight(xFinal).get<2>() - fHGnd;

  const double alphaM = 1. / fLMix;
  const double distance = (xInit-xFinal).GetMag();
  const double deltaHeight = abs(hFinal - hInit);
  const double cosTheta = deltaHeight / distance;

  double vaodSlant;

  if (deltaHeight > 1*m) {
    // Integration limits expect hInit < hFinal.
    if (hInit > hFinal)
      swap(hInit, hFinal);

    // Integrate extinction coefficient; be careful of integration limits.
    if (hInit < fHMix && hFinal < fHMix)
      vaodSlant = alphaM * (hFinal - hInit);
    else if (hInit < fHMix && hFinal > fHMix)
      vaodSlant = alphaM * (fHMix - hInit +
                  fHScl * (1. - exp(-(hFinal - fHMix) / fHScl)));
    else
      vaodSlant = alphaM * fHScl *
                  (exp(-(hInit - fHMix)/fHScl) - exp(-(hFinal - fHMix)/fHScl));

    // Convert vertical depth to slant depth.
    vaodSlant /= cosTheta;
  } else {
    // Assume attenuation at a fixed altitude.
    if (hInit < fHMix)
      vaodSlant = alphaM * distance;
    else
      vaodSlant = alphaM * exp(-(hInit - fHMix) / fHScl) * distance;
  }

  // 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
ParametricXMLMieModel::GetVerticalAerosolOpticalDepth(const unsigned int eyeId,
                                                      const double altitude)
  const
{
  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  const utl::Point& eyePos =
    det::Detector::GetInstance().GetFDetector().GetEye(eyeId).GetPosition();
  const double hEye =
    wgs84.PointToLatitudeLongitudeHeight(eyePos).get<2>() - fHGnd;
  const double hFinal = altitude - fHGnd;

  const double alphaM = 1. / fLMix;

  // Ignore ground level differences across the array and fix the mixing
  // layer height with respect to the prevailing ground height.
  if ((hFinal - hEye) > 1*m) {
    // Integrate extinction coefficient; be careful of integration limits.
    if (hEye < fHMix && hFinal < fHMix)
      return alphaM * (hFinal - hEye);
    else if (hEye < fHMix && hFinal > fHMix)
      return alphaM * (fHMix - hEye +
                       fHScl * (1. - exp(-(hFinal - fHMix) / fHScl)));
    else
      return alphaM * fHScl *
                  (exp(-(hEye - fHMix)/fHScl) - exp(-(hFinal - fHMix)/fHScl));
  } else
    return 0;
}


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

  double alpha;

  if (height < fHMix)
    alpha = 1. / fLMix;
  else
    alpha = 1. / (fLMix * exp((height - fHMix) / fHScl));

  return 1. / (alpha * pow(wLength / 355. / nanometer, -fGamma));
}


double
ParametricXMLMieModel::EvaluateScatteringAngle(const Point& /*p*/,
                                               const double angle,
                                               const double /*wLength*/)
  const
{
  const double g2 = fPhaseFunction.g * fPhaseFunction.g;
  const double mu = cos(angle);
  const double a = 0.25 * (1 - g2) / kPi;
  const double b = pow(1 + g2 - 2 * fPhaseFunction.g * mu, -1.5);
  const double c = fPhaseFunction.f * 0.5 * (3*mu*mu - 1) * pow(1 + g2, -1.5);

  return a * (b + c);
}