ParametricXMLRayleighModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <atm/ParametricXMLRayleighModel.h>
#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ScatteringResult.h>
#include <atm/AttenuationResult.h>
#include <atm/ProfileResult.h>

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

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

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

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


inline
double
RefractiveIndex(const double wl)
{
  // A. Bucholtz, Applied Optics 34, 2765 (1995), eq. 5
  // Refractive index at STP: requires input wavelength in microns
  const double invWl2 = pow(wl / micrometer, -2);

  return 1 + 1e-8 * (5791817. / (238.0185 - invWl2) + 167909. / (57.362 - invWl2));
}


ParametricXMLRayleighModel::ParametricXMLRayleighModel() :
  fIntegratedGDistance(false),
  fIntegrationStepWidth(0)
{ }


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

  topB.GetChild("RhoN").GetData(fRhoN);
  fIntegratedGDistance = bool(topB.GetChild("SlantDepthIntegrationStepWidth"));
  if (fIntegratedGDistance) {
    topB.GetChild("SlantDepthIntegrationStepWidth").GetData(fIntegrationStepWidth);

    ostringstream info;
    info << "do grammage integration in curved geometry using deltaDist="
         << fIntegrationStepWidth/m
         << " m";
    INFO(info);
  }
}


ScatteringResult
ParametricXMLRayleighModel::EvaluateRayleighScattering(const utl::Point& xA,
                                                       const utl::Point& xB,
                                                       const double angle,
                                                       const double distance,
                                                       const std::vector<double>& wLength)
  const
{
  AttenuationResult attenuation = EvaluateRayleighAttenuation(xA, xB, wLength);
  return EvaluateRayleighScattering(xA, xB, angle, distance, attenuation);
}


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

  const unsigned int nWl = attenuation.GetNPoints();
  for (unsigned int iWl = 0; iWl < nWl; ++iWl) {
    const double wl = attenuation.GetX(iWl);
    const double sct = EvaluateRayleighScattering(xA, xB, angle, distance, wl, attenuation.GetY(iWl));
    // Errors to be implemented
    frac.PushBack(wl, wError, sct, fractionError);
  }

  return ScatteringResult(frac, angle);
}


double
ParametricXMLRayleighModel::EvaluateRayleighScattering(const utl::Point& xA,
                                                       const utl::Point& /*xB*/,
                                                       const double angle,
                                                       const double distance,
                                                       double wLength,
                                                       const double raylAtt)
  const
{
  return (1 - raylAtt)
         * EvaluateScatteringAngle(xA, angle, wLength)
         / (distance * distance);
}


double
ParametricXMLRayleighModel::EvaluateRayleighScattering(const utl::Point&  xA,
                                                       const utl::Point&  xB,
                                                       const double angle,
                                                       const double distance,
                                                       double wLength)
  const
{
  const double raylAtt = EvaluateRayleighAttenuation(xA, xB, wLength);
  return EvaluateRayleighScattering(xA, xB, angle, distance, wLength, raylAtt);
}


AttenuationResult
ParametricXMLRayleighModel::EvaluateRayleighAttenuation(const utl::Point& xInit,
                                                        const utl::Point& xFinal,
                                                        const vector<double>& wLength)
  const
{
  utl::TabulatedFunctionErrors attenuation;
  const double rAttError = 0;
  const double wError = 0;
  const double slantDistance = GDistance(xInit, xFinal);

  for (std::vector<double>::const_iterator it = wLength.begin();
       it != wLength.end(); ++it) {
    const double attenuationLength = GetAttenuationLength(xInit, *it);
    const double transAtt = exp(-slantDistance/attenuationLength);
    attenuation.PushBack(*it, wError, transAtt, rAttError);
  }

  return AttenuationResult(attenuation);
}


double
ParametricXMLRayleighModel::EvaluateRayleighAttenuation(const utl::Point& xInit,
                                                        const utl::Point& xFinal,
                                                        double wLength)
  const
{
  return exp(-GDistance(xInit, xFinal) / GetAttenuationLength(xInit, wLength));
}


double
ParametricXMLRayleighModel::GDistance(const utl::Point& xA,
                                      const utl::Point& xB)
  const
{
  Detector& theDetector = Detector::GetInstance();
  const Atmosphere& atmosphere = theDetector.GetAtmosphere();

  double xDistance = 0;

  if (fIntegratedGDistance) {
    xDistance = atmosphere.IntegratedGrammage(xA, xB, fIntegrationStepWidth);
  } else {
    ReferenceEllipsoid wgs84(ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84));
    double zA = (wgs84.PointToLatitudeLongitudeHeight(xA)).get<2>();
    double zB = (wgs84.PointToLatitudeLongitudeHeight(xB)).get<2>();

    const ProfileResult& dvh = atmosphere.EvaluateDepthVsHeight();
    double atmHeightMin = dvh.MinX();
    double atmHeightMax = dvh.MaxX();
    double depthA, depthB;

    if (zA < atmHeightMin)
      depthA =  dvh.Y(atmHeightMin);
    else if (zA > atmHeightMax)
      depthA =  dvh.Y(atmHeightMax);
    else
      depthA = dvh.Y(zA);

    if (zB < atmHeightMin)
      depthB =  dvh.Y(atmHeightMin);
    else if (zB > atmHeightMax)
      depthB =  dvh.Y(atmHeightMax);
    else
      depthB = dvh.Y(zB);

    double distAB = (xA - xB).GetMag();

    if (zA != zB)
      xDistance = distAB * abs((depthA - depthB) / (zA - zB));
    else if (zA <= 12.0*km)
      xDistance = distAB * abs(depthA / (0.19 * (44340 * m - zB)));
    else if (zA <= 45.470*km)
      xDistance = distAB * abs(depthA / (6.340*km));
    else
      xDistance = 0.0;
  }

  return xDistance;
}


double
ParametricXMLRayleighModel::GetAttenuationLength(const utl::Point& /*p*/,
                                                 const double wLength)
  const
{
  // N_s = molecules per cubic centimetre at T=15deg, 760mm Hg (Bucholtz)
  const double N_s = 2.54743e19 / cm3;
  const double factorKing = (6 + 3.*fRhoN) / (6 - 7.*fRhoN);

  const double refIndex = RefractiveIndex(wLength);
  const double refIndex2 = refIndex * refIndex;

  const double sigma_molecule =
      24 * pow(kPi, 3)
         * pow(((refIndex2 - 1) / (refIndex2 + 2)), 2)
         /  pow(wLength, 4) / (N_s * N_s) * factorKing;

  return kDryAirMolarMass / (sigma_molecule * kAvogadro);
}


double
ParametricXMLRayleighModel::EvaluateScatteringAngle(const Point& /*p*/,
                                                    const double angle,
                                                    const double /*wLength*/)
  const
{
  const double gamRay = fRhoN / (2. - fRhoN);
  const double factor = 3 / (16 * kPi);
  return factor / (1 + 2*gamRay) *
    ((1 + 3*gamRay) + (1 - gamRay) * pow(cos(angle), 2));
}