USStdADBProfileModel.cc

Go to the documentation of this file.

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/TabulatedFunction.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/CentralConfig.h>

#include <atm/USStdADBProfileModel.h>
#include <atm/MolecularDB.h>
#include <atm/ProfileResult.h>
#include <atm/MolecularLayer.h>
#include <atm/MolecularZone.h>
#include <atm/ProfileResult.h>
#include <det/Detector.h>

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

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


// 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


bool
USStdADBProfileModel::HasData()
  const
{
  if (fTablesFilled)
    return true;
  FillTables();
  if (fTablesFilled)
    return true;
  return false;
}


const ProfileResult&
USStdADBProfileModel::EvaluatePressureVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabLogPressureVsHeight;
}


const ProfileResult&
USStdADBProfileModel::EvaluateTemperatureVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabTemperatureVsHeight;
}


const ProfileResult&
USStdADBProfileModel::EvaluateVaporPressureVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabLogVaporPressureVsHeight;
}


const ProfileResult&
USStdADBProfileModel::EvaluateDensityVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabLogDensityVsHeight;
}


const ProfileResult&
USStdADBProfileModel::EvaluateDepthVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabLogDepthVsHeight;
}


const ProfileResult&
USStdADBProfileModel::EvaluateHeightVsDepth()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabHeightVsLogDepth;
}


const ProfileResult&
USStdADBProfileModel::EvaluateRefractionIndexVsHeight()
  const
{
  if (!fTablesFilled)
    FillTables();
  return *fTabRIVsHeight;
}


void
USStdADBProfileModel::FillTables()
  const
{
  const MolecularDB& mDB =
    det::Detector::GetInstance().GetAtmosphere().GetMolecularDB(MolecularIds::eUSStdA);

  TabulatedFunctionErrors tabLogDensityVsHeight;
  TabulatedFunctionErrors tabRIVsHeight;
  TabulatedFunctionErrors tabLogPressureVsHeight;
  TabulatedFunctionErrors tabTemperatureVsHeight;
  TabulatedFunctionErrors tabLogVaporPressureVsHeight;
  TabulatedFunctionErrors tabLogDepthVsHeight;
  TabulatedFunctionErrors tabHeightVsLogDepth;

  // Fill tabulated functions.
  // Perform conversions into Auger base units.
  int nZonesCheck = 0;
  const double minLog = kLn10 * numeric_limits<double>::min_exponent10;

  for (MolecularDB::ZoneIterator zIt = mDB.ZonesBegin();
       zIt != mDB.ZonesEnd(); ++zIt) {

    for (MolecularZone::LayerIterator lIt = zIt->LayersBegin();
         lIt != zIt->LayersEnd(); ++lIt) {

      // TEMPORARY.. in case some entries in database are not ordered
      // by ordinate, skip the unorder entries for now.  This should be
      // fixed automatically when TabulatedFunction is redesigned + reimplemented,
      // at which time code up to 'END TEMPORARY' should be removed.

      if (lIt != zIt->LayersBegin()) {
        MolecularZone::LayerIterator newIt = lIt;
         if (newIt->GetHeight() <= (--newIt)->GetHeight())
           break;
      }
      // END TEMPORARY patch to compensate shitty TabulatedFunction class

      tabLogDepthVsHeight.PushBack(
        lIt->GetHeight(), lIt->GetHeightError(),
        log(lIt->GetDepth()), lIt->GetDepthError() / lIt->GetDepth()
      );

      const double refracIndex =
        utl::RefractionIndex::LorentzLorentz(lIt->GetPressure() * (kOverburdenSeaLevel/atmosphere));
      tabRIVsHeight.PushBack(
        lIt->GetHeight(), 0,  // ERRORS need to be computed
        refracIndex, 0
      );

      tabLogPressureVsHeight.PushBack(
        lIt->GetHeight(), lIt->GetHeightError(),
        log(lIt->GetPressure()), lIt->GetPressureError() / lIt->GetPressure()
      );

      // Calculate vapor pressure from relative humidity
      const double pSat = utl::SaturationVaporPressure(lIt->GetTemperature());
      const double pVapor = lIt->GetHumidity() * pSat;
      const double dpVapor = lIt->GetHumidityError() * pSat;  // Ignore dT...

      tabLogVaporPressureVsHeight.PushBack(
        lIt->GetHeight(), lIt->GetHeightError(),
        (pVapor > 0) ? log(pVapor) : minLog, (pVapor > 0 && dpVapor > 0) ? pVapor / dpVapor : minLog
      );

      tabTemperatureVsHeight.PushBack(
        lIt->GetHeight(), lIt->GetHeightError(),
        lIt->GetTemperature(), lIt->GetTemperatureError()
      );

      tabLogDensityVsHeight.PushBack(
        lIt->GetHeight(), lIt->GetHeightError(),
        log(lIt->GetAirDensity()),
        lIt->GetAirDensityError() / lIt->GetAirDensity()
      );

    }

    MolecularZone::LayerIterator lIt = zIt->LayersEnd();
    while (lIt != zIt->LayersBegin()) {
      --lIt;

      // TEMPORARY.. in case some entries in database are not ordered
      // by ordinate, skip the unorder entries for now.  This should be
      // fixed automatically when TabulatedFunction is redesigned + reimplemented,
      // at which time code up to 'END TEMPORARY' should be removed.
      if (lIt != --(zIt->LayersEnd())) {
        MolecularZone::LayerIterator newIt = lIt;
        if (newIt->GetHeight() >= (++newIt)->GetHeight())
          break;
      }
      // END TEMPORARY patch to compensate shitty TabulatedFunction class

      tabHeightVsLogDepth.PushBack(
        log(lIt->GetDepth()), lIt->GetDepthError() / lIt->GetDepth(),
        lIt->GetHeight(), lIt->GetHeightError()
      );
    }
    ++nZonesCheck;
  }
  if (nZonesCheck != 1) {
    ostringstream msg;
    msg << "Expected to find one zone in DB for US Std Atmosphere, but found"
        << nZonesCheck;
    DEBUGLOG(msg);
    return;
  }

  // create profile results
  *fTabLogPressureVsHeight =
    ProfileResult(tabLogPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogVaporPressureVsHeight =
    ProfileResult(tabLogVaporPressureVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabTemperatureVsHeight =
    ProfileResult(tabTemperatureVsHeight, ProfileResult::eLinear, ProfileResult::eLinear);

  *fTabLogDensityVsHeight =
    ProfileResult(tabLogDensityVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabLogDepthVsHeight =
    ProfileResult(tabLogDepthVsHeight, ProfileResult::eLinear, ProfileResult::eLog);

  *fTabHeightVsLogDepth =
    ProfileResult(tabHeightVsLogDepth, ProfileResult::eLog, ProfileResult::eLinear);

  *fTabRIVsHeight =
    ProfileResult(tabRIVsHeight, ProfileResult::eLinear, ProfileResult::eLinear);

  // Extend all tables to a height of 100 km
  ExtendProfilesTo100km();

  fTablesFilled = true;
}