HumidAirRayleighModel.cc

Go to the documentation of this file.

#include <det/Detector.h>

#include <atm/HumidAirRayleighModel.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;


HumidAirRayleighModel::HumidAirRayleighModel() :
  fIntegrationStepWidth(0)
{
}


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

  topB.GetChild("IntegrationStepWidth").GetData(fIntegrationStepWidth);
}


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


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

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

  return ScatteringResult(frac, angle);
}


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


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



AttenuationResult
HumidAirRayleighModel::EvaluateRayleighAttenuation(const Point& xInit,
                                                   const Point& xFinal,
                                                   const vector<double>& wLength)
  const
{
  TabulatedFunctionErrors attenuation;
  const double transAttError = 0;
  const double wError = 0;

  for (std::vector<double>::const_iterator it = wLength.begin();
       it != wLength.end(); ++it) {
    const double transAtt = EvaluateRayleighAttenuation(xInit, xFinal, *it);
    attenuation.PushBack(*it, wError, transAtt, transAttError);
  }

  return AttenuationResult(attenuation);
}


double
HumidAirRayleighModel::EvaluateRayleighAttenuation(const Point& xInit,
                                                   const Point& xFinal,
                                                   double wLength)
  const
{
  const Atmosphere& atmos = Detector::GetInstance().GetAtmosphere();

  const ProfileResult& prVsH = atmos.EvaluatePressureVsHeight();
  const ProfileResult& pvVsH = atmos.EvaluateVaporPressureVsHeight();
  const ProfileResult& tempVsH = atmos.EvaluateTemperatureVsHeight();
  const ProfileResult& riVsH = atmos.EvaluateRefractionIndexVsHeight(wLength);

  const ReferenceEllipsoid& wgs84 = ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  Vector dir = xFinal - xInit;
  const double length = dir.GetMag();
  dir.Normalize();

  double distance = 0;
  double transmittance = 1;

  while (distance < length) {
    const double delta =
      (distance + fIntegrationStepWidth < length) ?
        fIntegrationStepWidth : (length - distance);

    if (delta < 1*m)
      break;

    const Point p1 = xInit + distance * dir;
    distance += delta;
    const Point p2 = xInit + distance * dir;

    const double hP1 = wgs84.PointToLatitudeLongitudeHeight(p1).get<2>();
    const double hP2 = wgs84.PointToLatitudeLongitudeHeight(p2).get<2>();
    const double height = 0.5 * (hP1 + hP2);

    const double pr = prVsH.Y(height);
    const double pv = pvVsH.Y(height);
    const double temp = tempVsH.Y(height);
    const double ri = riVsH.Y(height);

    const double Nair = pr / (kBoltzmann * temp);
    const double sigmaR = RayleighCrossSection(wLength, ri, temp, pr, pv);

    // Combine the optical depth + Beer-Lambert-Bouguer law:
    //   tau_m(0,r; lambda) = \int_0^r N(r') sigmaR(r', lambda) dr',
    //   T(0,r; lambda) = exp{-tau(0,r; lambda)}
    transmittance *= exp(-Nair * sigmaR * delta);
  }

  return transmittance;
}


double
HumidAirRayleighModel::GetAttenuationLength(const Point& p,
                                            const double wLength)
  const
{
  const Atmosphere& atmos = Detector::GetInstance().GetAtmosphere();

  const ReferenceEllipsoid& wgs84 =
      ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);
  const double height = wgs84.PointToLatitudeLongitudeHeight(p).get<2>();

  const double ri = atmos.EvaluateRefractionIndexVsHeight(wLength).Y(height);
  const double pr = atmos.EvaluatePressureVsHeight().Y(height);
  const double pv = atmos.EvaluateVaporPressureVsHeight().Y(height);
  const double temp = atmos.EvaluateTemperatureVsHeight().Y(height);

  const double Nair = pr / (kBoltzmann * temp);
  const double sigmaR = RayleighCrossSection(wLength, ri, temp, pr, pv);

  // Attenuation length = inverse of volume scattering coefficient
  return 1 / (Nair * sigmaR);
}


double
HumidAirRayleighModel::EvaluateScatteringAngle(const Point& p,
                                               const double angle,
                                               const double wLength)
  const
{
  const Atmosphere& atmos = Detector::GetInstance().GetAtmosphere();
  const ReferenceEllipsoid& wgs84 =
    ReferenceEllipsoid::Get(ReferenceEllipsoid::eWGS84);

  const double height = wgs84.PointToLatitudeLongitudeHeight(p).get<2>();
  const double pr = atmos.EvaluatePressureVsHeight().Y(height);
  const double pv = atmos.EvaluateVaporPressureVsHeight().Y(height);

  // Depolarization: use dispersion of the King factor in a humid atmosphere
  const double Fk = KingFactor(wLength, pr, pv);
  const double rhoN = (6. * Fk - 6.) / (7. * Fk + 3.);
  const double costh = cos(angle);

  // Phase function
  return 3 / (8*kPi * (2 + rhoN)) * ((1 + rhoN) + (1 - rhoN) * costh*costh);
}


double
HumidAirRayleighModel::KingFactor(const double wl, const double pressure,
                                  const double vaporPressure)
  const
{
  const double invWl2 = pow(wl / micrometer, -2.);

  // King factors for primary constituents of air.
  // See C. Tomasi et al,. Appl. Opt. 44 (2005) 3320, eqs. 20-22
  const double FkN2  = 1.034 + 3.17e-4 * invWl2;
  const double FkO2  = 1.096 + invWl2 * (1.385e-3 + 1.448e-4 * invWl2);
  const double FkAr  = 1.;
  const double FkCO2 = 1.15;
  const double FkH2O = 1.001;

  const double H2OAirFraction = vaporPressure / pressure;

  // Total King factor is a weighted sum of the constituents of air
  return (FkN2 * kN2AirFraction + FkO2 * kO2AirFraction +
          FkAr * kArAirFraction + FkCO2 * kCO2AirFraction +
          FkH2O * H2OAirFraction) /
         (kN2AirFraction + kO2AirFraction +
          kArAirFraction + kCO2AirFraction +
          H2OAirFraction);
}


double
HumidAirRayleighModel::RayleighCrossSection(const double wl,
                                            const double refIndex,
                                            const double temperature,
                                            const double pressure,
                                            const double vaporPressure)
  const
{
  const double N_air = pressure / (kBoltzmann * temperature);
  const double refIndex2 = refIndex * refIndex;
  const double factorKing = KingFactor(wl, pressure, vaporPressure);

  return 24* pow(kPi, 3)
         * pow(((refIndex2 - 1) / (refIndex2 + 2)), 2)
         * pow(wl, -4) / (N_air * N_air) * factorKing;
}