VProfileModel.cc

Go to the documentation of this file.

#include <atm/VProfileModel.h>
#include <atm/ProfileResult.h>

#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/TabulatedFunctionErrors.h>

#include <iostream>
#include <limits>

using namespace utl;
using namespace std;


// For some reason, CLHEP SystemOfUnits.h #defines pascal to hep_pascal,
// putting it in the global namespace and screwing up our own definition of
// pascal from AugerUnits.h.  It appears the CLHEP units are being pulled in
// via the geometry package.  The following line undoes the CLHEP sabotage.
#undef pascal


namespace atm {

  VProfileModel::VProfileModel() :
    fTabLogDepthVsHeight(new ProfileResult),
    fTabHeightVsLogDepth(new ProfileResult),
    fTabLogPressureVsHeight(new ProfileResult),
    fTabLogVaporPressureVsHeight(new ProfileResult),
    fTabTemperatureVsHeight(new ProfileResult),
    fTabLogDensityVsHeight(new ProfileResult),
    fTabRIVsHeight(new ProfileResult)
  { }


  VProfileModel::~VProfileModel()
  {
    delete fTabLogDepthVsHeight;
    delete fTabHeightVsLogDepth;
    delete fTabLogPressureVsHeight;
    delete fTabLogVaporPressureVsHeight;
    delete fTabTemperatureVsHeight;
    delete fTabLogDensityVsHeight;
    delete fTabRIVsHeight;

    CleanRIVsWavelength();
  }


  const ProfileResult&
  VProfileModel::EvaluateRefractionIndexVsHeight()
    const
  {
    // eventually, all calculations of refraction indices could be moved
    // here (since calculation is same for all the derived classes).
    // Unfortunately doing so is best left until after we have
    // TabulatedFunction and TabulatedFunctionErrors classes which
    // do not suck so severely.

    // Where the relation : (n^2-1)/(n^2+2) = const.*density
    //   double d = (((n0*n0) - 1)/((n0*n0) + 2))/overburdenSeaLevel;

    //   ProfileResult rProf = EvaluateDepthVsHeight();

    //   for (TabulatedFunction::Iterator it = fTabLogDepthVsHeight->Begin();
    //        it != fTabLogDepthVsHeight->End() ; ++it){

    //     //    it->Y() = sqrt

    //   }
    return *fTabRIVsHeight;
  }


  const
  ProfileResult&
  VProfileModel::EvaluateRefractionIndexVsHeight(const double waveLength)
    const
  {
    // Find wavelength profile if already stored
    const AltWLFunction::const_iterator profIt = fTabRIVsHeightAndWaveLength.find(waveLength);
    if (profIt != fTabRIVsHeightAndWaveLength.end())
      return *profIt->second;

    // Else, evaluate the profile, assuming other tables have been filled already
    TabulatedFunctionErrors tabRefIndex;
    for (TabulatedFunction::Array::const_iterator heightIt = fTabTemperatureVsHeight->fProfile.XBegin();
         heightIt != fTabTemperatureVsHeight->fProfile.XEnd(); ++heightIt) {
      const double temperature = fTabTemperatureVsHeight->Y(*heightIt);
      const double pressure = fTabLogPressureVsHeight->Y(*heightIt);
      const double vaporPressure = fTabLogVaporPressureVsHeight->Y(*heightIt);

      const double nAir = utl::RefractionIndex::Ciddor95(waveLength, temperature, pressure, vaporPressure);

      tabRefIndex.PushBack(*heightIt, 0, nAir, 0);  // Ignore errors for now
    }

    ProfileResult* const refIndex =
      new ProfileResult(tabRefIndex, ProfileResult::eLinear, ProfileResult::eLinear);
    fTabRIVsHeightAndWaveLength[waveLength] = refIndex;

    return *refIndex;
  }


  double
  VProfileModel::GetVerticalTimeOfFlight(const double h1, const double h2)
    const
  {
    const double height1 = min(h1, h2);
    const double height2 = max(h1, h2);

    // Average Malargue atmosphere parameterization
    // (here we use the April atmosphere from GAP-2005-021)
    const int nLayers = 6;
    const double borders[nLayers + 1] = {
      -numeric_limits<double>::max(),
      0*m,
      9900*m,
      12100*m,
      38000*m,
      100000*m,
      numeric_limits<double>::max()
    };

    const double bParams[nLayers] = {
      0,
      1174.1385*g/cm2,
      1215.7355*g/cm2,
      1425.6919*g/cm2,
      503.93170*g/cm2,
      0
    };

    const double cParams[nLayers] = {
      0,
      967164.96*cm,
      768155.10*cm,
      621852.25*cm,
      785601.62*cm,
      0
    };

    enum Layers {
      eBelowSeaLevel = 0,
      eSeaLevel      = 1,
      eVacuum        = nLayers-1
    };

    const double rho0 = bParams[eSeaLevel] / cParams[eSeaLevel];

    // difference to c_vaccum at sea level
    const double deltaN = kRefractiveIndexSeaLevel - 1;

    // integrate from height1 to height2
    // (see GAP-2007-99, Eq. (6.25), but note the wrong sign!)

    double verticalLength = 0;

    for (int i = 0; i < nLayers; ++i) {

      if (height1 > borders[i] && height1 <= borders[i+1]) {

        double startHeight = height1;

        for (int j = i; j < nLayers; ++j) {

          const double stopHeight =
            (height2 > borders[j] && height2 <= borders[j+1]) ? height2 : borders[j+1];

          const double deltaH = stopHeight - startHeight;

          if (j == eBelowSeaLevel)
            verticalLength += kRefractiveIndexSeaLevel * deltaH;
          else if (j == eVacuum)
            verticalLength += deltaH;
          else
            verticalLength += deltaH +
              deltaN*bParams[j]/rho0 * (exp(-startHeight / cParams[j]) - exp(-stopHeight / cParams[j]));

          if (stopHeight == height2)
            break;

          startHeight = stopHeight;

        }
        break;

      }

    }

    return verticalLength / kSpeedOfLight;
  }


  void
  VProfileModel::CleanRIVsWavelength() const
  {
    // Remove stored wavelength-dependent refraction index height profiles
    for (const auto& ref : fTabRIVsHeightAndWaveLength)
      delete ref.second;

    fTabRIVsHeightAndWaveLength.clear();
  }


  void
  VProfileModel::ExtendProfilesTo100km()
    const
  {
    // Extend all tables to 100 km
    // DEPRECATED: when Atm_Molecular_0_A is REPLACED, this will be removed...
    TabulatedFunctionErrors& logDepthVsHeight = fTabLogDepthVsHeight->fProfile;
    TabulatedFunctionErrors& tempVsHeight = fTabTemperatureVsHeight->fProfile;
    TabulatedFunctionErrors& logPVsHeight = fTabLogPressureVsHeight->fProfile;
    TabulatedFunctionErrors& logDensityVsHeight = fTabLogDensityVsHeight->fProfile;

    const double maxHeight = logDepthVsHeight.XBack();
    const double minLog = kLn10 * numeric_limits<double>::min_exponent10;

    // Balloon profiles may not extend above the troposphere
    if (maxHeight < 11*km) {
      logDepthVsHeight.PushBack(11*km, 0, log(231.81*g/cm2), minLog);
      tempVsHeight.PushBack(11*km, 0, 216.65*kelvin, 0.);
      logPVsHeight.PushBack(11*km, 0, log(22633*pascal), minLog);
      logDensityVsHeight.PushBack(11*km, 0, log(0.3638*kg/m3), minLog);
    }

    if (maxHeight < 20*km) {
      logDepthVsHeight.PushBack(20*km, 0, log(56.24*g/cm2), minLog);
      tempVsHeight.PushBack(20*km, 0, 216.65*kelvin, 0);
      logPVsHeight.PushBack(20*km, 0, log(5475*pascal), minLog);
      logDensityVsHeight.PushBack(20*km, 0, log(0.088*kg/m3), minLog);
    }

    // Monthly DBs usually extend to 30 km only
    logDepthVsHeight.PushBack( 32*km, 0, log(8.95*g/cm2), minLog);
    logDepthVsHeight.PushBack( 47*km, 0, log(1.15*g/cm2), minLog);
    logDepthVsHeight.PushBack( 51*km, 0, log(0.68*g/cm2), minLog);
    logDepthVsHeight.PushBack( 71*km, 0, log(0.04*g/cm2), minLog);
    logDepthVsHeight.PushBack(100*km, 0, log(1e-3*g/cm2), minLog);

    tempVsHeight.PushBack( 32*km, 0, 228.00*kelvin, 0);
    tempVsHeight.PushBack( 47*km, 0, 270.65*kelvin, 0);
    tempVsHeight.PushBack( 51*km, 0, 270.65*kelvin, 0);
    tempVsHeight.PushBack( 71*km, 0, 214.65*kelvin, 0);
    tempVsHeight.PushBack(100*km, 0, 214.65*kelvin, 0);

    logPVsHeight.PushBack( 32*km, 0, log(868*pascal), minLog);
    logPVsHeight.PushBack( 47*km, 0, log(111*pascal), minLog);
    logPVsHeight.PushBack( 51*km, 0, log( 66*pascal), minLog);
    logPVsHeight.PushBack( 71*km, 0, log(  4*pascal), minLog);
    logPVsHeight.PushBack(100*km, 0, log(0.1*pascal), minLog);

    logDensityVsHeight.PushBack( 32*km, 0, log(0.0132*kg/m3), minLog);
    logDensityVsHeight.PushBack( 47*km, 0, log(0.0014*kg/m3), minLog);
    logDensityVsHeight.PushBack( 51*km, 0, log(0.0009*kg/m3), minLog);
    logDensityVsHeight.PushBack( 71*km, 0, log(0.0001*kg/m3), minLog);
    logDensityVsHeight.PushBack(100*km, 0, log(0.00001*kg/m3), minLog);

    fTabLogVaporPressureVsHeight->fProfile.PushBack(100*km, 0, minLog, minLog);
    fTabRIVsHeight->fProfile.PushBack(100*km, 0, 1 + std::numeric_limits<double>::epsilon(), 0);

    // Fill log(depth) table to height of 100 km (NOTE: use reverse order)
    TabulatedFunctionErrors& heightVsLogDepth = fTabHeightVsLogDepth->fProfile;
    TabulatedFunctionErrIterator hxIt = heightVsLogDepth.Begin();

    if (maxHeight < 11*km)
      hxIt = heightVsLogDepth.Insert(hxIt, log(231.81*g/cm2), minLog, 11*km, 0);
    if (maxHeight < 20*km)
      hxIt = heightVsLogDepth.Insert(hxIt, log(56.24*g/cm2), minLog, 20*km, 0);

    hxIt = heightVsLogDepth.Insert(hxIt, log(8.95*g/cm2), minLog,  32*km, 0);
    hxIt = heightVsLogDepth.Insert(hxIt, log(1.15*g/cm2), minLog,  47*km, 0);
    hxIt = heightVsLogDepth.Insert(hxIt, log(0.68*g/cm2), minLog,  51*km, 0);
    hxIt = heightVsLogDepth.Insert(hxIt, log(0.04*g/cm2), minLog,  71*km, 0);
    hxIt = heightVsLogDepth.Insert(hxIt, log(1e-3*g/cm2), minLog, 100*km, 0);
  }

}