GDASProfileModel.cc

Go to the documentation of this file.

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

#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/TimeStamp.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/CentralConfig.h>

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

#include <utl/ErrorLogger.h>

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

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

GDASProfileModel::GDASProfileModel() :
  VProfileModel(),
  fCurrentTime(0),
  fDbIsEmptyNow(false)
{
}


void
GDASProfileModel::Init()
{
}


bool
GDASProfileModel::HasData()
  const
{
  return CheckForUpdates();
}


bool
GDASProfileModel::CheckForUpdates()
  const
{
  // check if detector has been updated and fill
  // various lookup tables if so
  //
  if (fCurrentTime == det::Detector::GetInstance().GetTime() &&
      fCurrentTime)
    return !fDbIsEmptyNow;

  fCurrentTime = det::Detector::GetInstance().GetTime();
  fDbIsEmptyNow = false;

  // Atmosphere state is changing: clean up calculated refraction indices
  CleanRIVsWavelength();

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

  const MolecularDB& mDB =
    det::Detector::GetInstance().GetAtmosphere().GetMolecularDB(MolecularIds::eMalargueGDAS);

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

      // Skip everything above 900hPa (below ~1000m a.s.l.)
      // since InclinedAtmosphericProfile bitches otherwise
      //if(lIt->GetPressure() > 900*hPa)
      //  continue;
      // Should be fine with new DB filler

      // TEMPORARY.. skip of unordered entries
      bool skip = false;
      if (lIt != zIt->LayersBegin()) {
        MolecularZone::LayerIterator newIt = lIt;
          if (newIt->GetHeight() <= (--newIt)->GetHeight())
            skip = true;
      }
      if (skip)
        break;
      // END TEMPORARY

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

      const double refracIndex = utl::RefractionIndex::LorentzLorentz(lIt->GetDepth());
      tabRIVsHeight.PushBack(lIt->GetHeight(), lIt->GetHeightError(),
                             refracIndex, 0.);

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

      // Convert relative humidity to vapor pressure
      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());

    }
    ++nZonesCheck;
  }

  if (nZonesCheck == 0) {
    fDbIsEmptyNow = true;

    ostringstream msg;
    msg << "GDAS model requested data from database, "
        << "but no data found for time "
        << det::Detector::GetInstance().GetTime();
    DEBUGLOG(msg);

    return false;
  }

  if (nZonesCheck != 1) {
    // some kind of error here.
  }


  // Extrapolate to ground level
  const double groundHeight = 1. * km;
  const double minHeight = tabTemperatureVsHeight.GetX(0);
  const int nLayers = distance(tabTemperatureVsHeight.Begin(),
                               tabTemperatureVsHeight.End());

  if (nLayers > 2 && minHeight > groundHeight) {
    const double h1   = tabTemperatureVsHeight.GetX(0);
    const double h2   = tabTemperatureVsHeight.GetX(1);
    const double relH = (groundHeight-h1)/(h2-h1);

    // Linear extrapolation
    const double groundLogDepth = tabLogDepthVsHeight.GetY(0) + 
      (tabLogDepthVsHeight.GetY(1) - tabLogDepthVsHeight.GetY(0)) * relH;

    const double groundLogDensity = tabLogDensityVsHeight.GetY(0) + 
      (tabLogDensityVsHeight.GetY(1) - tabLogDensityVsHeight.GetY(0)) * relH;

    const double groundLogPressure = tabLogPressureVsHeight.GetY(0) + 
      (tabLogPressureVsHeight.GetY(1) - tabLogPressureVsHeight.GetY(0)) * relH;

    const double groundLogVaporPressure = tabLogVaporPressureVsHeight.GetY(0) + 
      (tabLogVaporPressureVsHeight.GetY(1) - tabLogVaporPressureVsHeight.GetY(0)) * relH;

    const double groundTemperature = tabTemperatureVsHeight.GetY(0) + 
      (tabTemperatureVsHeight.GetY(1) - tabTemperatureVsHeight.GetY(0)) * relH;

    const double groundRefIndex = tabRIVsHeight.GetY(0) + 
      (tabRIVsHeight.GetY(1) - tabRIVsHeight.GetY(0)) * relH;

    // Insert the ground level values
    tabLogDepthVsHeight.Insert(tabLogDepthVsHeight.Begin(),
          groundHeight, 0., groundLogDepth, minLog);

    tabLogDensityVsHeight.Insert(tabLogDensityVsHeight.Begin(),
          groundHeight, 0., groundLogDensity, minLog);

    tabLogPressureVsHeight.Insert(tabLogPressureVsHeight.Begin(),
          groundHeight, 0., groundLogPressure, minLog);

    tabLogVaporPressureVsHeight.Insert(tabLogVaporPressureVsHeight.Begin(),
          groundHeight, 0., groundLogVaporPressure, minLog);

    tabTemperatureVsHeight.Insert(tabTemperatureVsHeight.Begin(),
          groundHeight, 0., groundTemperature, 0.);

    tabRIVsHeight.Insert(tabRIVsHeight.Begin(),
          groundHeight, 0., groundRefIndex, 0.);
  }

  // If the database contains values below 1000 meters the linear extrapolation
  // is skipped and nLayers is not increased, causing the next loop to crash
  const int nLayersDepth = distance(tabLogDepthVsHeight.Begin(),
                                    tabLogDepthVsHeight.End());

  // Fill the log depth vs height profile, starting from the top layer
  for (int i = nLayersDepth-1; i >= 0; --i)
    tabHeightVsLogDepth.PushBack(
          tabLogDepthVsHeight.GetY(i), tabLogDepthVsHeight.GetYErr(i),
          tabLogDepthVsHeight.GetX(i), tabLogDepthVsHeight.GetXErr(i));

  // 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();

  return true;
}


const ProfileResult&
GDASProfileModel::EvaluatePressureVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data found for pressure");
  return *fTabLogPressureVsHeight;
}


const ProfileResult&
GDASProfileModel::EvaluateTemperatureVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data found for temperature");
  return *fTabTemperatureVsHeight;
}


const ProfileResult&
GDASProfileModel::EvaluateVaporPressureVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data found for humidity");
  return *fTabLogVaporPressureVsHeight;
}


const ProfileResult&
GDASProfileModel::EvaluateDensityVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data found for density");
  return *fTabLogDensityVsHeight;
}


const ProfileResult&
GDASProfileModel::EvaluateDepthVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data for depth");
  return *fTabLogDepthVsHeight;
}


const ProfileResult&
GDASProfileModel::EvaluateHeightVsDepth()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data for height vs depth");
  return *fTabHeightVsLogDepth;
}


const ProfileResult&
GDASProfileModel::EvaluateRefractionIndexVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No GDAS data for refraction index");
  return *fTabRIVsHeight;
}