MonthlyAvgDBProfileModel.cc

Go to the documentation of this file.

#include <iostream>
#include <string>
#include <sstream>
#include <vector>
#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/UTCDateTime.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>

#include <fwk/CentralConfig.h>

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


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


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


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


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


const ProfileResult&
MonthlyAvgDBProfileModel::EvaluateDepthVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No monthly average data found for depth.");
  return *fTabLogDepthVsHeight;
}


const ProfileResult&
MonthlyAvgDBProfileModel::EvaluateHeightVsDepth()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No monthly average data found for height vs depth.");
  return *fTabHeightVsLogDepth;
}


const ProfileResult&
MonthlyAvgDBProfileModel::EvaluateRefractionIndexVsHeight()
  const
{
  if (!CheckForUpdates())
    throw NoDataForModelException("No monthly average data found for refraction index.");
  return *fTabRIVsHeight;
}


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


bool
MonthlyAvgDBProfileModel::CheckForUpdates()
  const
{
  // Check if the detector has been updated to a new month,
  // and fill verious lookup tables if so
  if (fCurrentTime &&
      UTCDateTime(fCurrentTime).GetMonth() ==
      UTCDateTime(det::Detector::GetInstance().GetTime()).GetMonth())
    return true;

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

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

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

  // determine (seasonal average) profileType based on month
  const unsigned int iMonth = UTCDateTime(fCurrentTime).GetMonth();

  MolecularIds::ProfileId profileId = MolecularIds::eNone;
  switch (iMonth) {
  case 1:
    profileId = MolecularIds::eMalargueJanuaryModel;
    break;
  case 2:
    profileId = MolecularIds::eMalargueFebruaryModel;
    break;
  case 3:
    profileId = MolecularIds::eMalargueMarchModel;
    break;
  case 4:
    profileId = MolecularIds::eMalargueAprilModel;
    break;
  case 5:
    profileId = MolecularIds::eMalargueMayModel;
    break;
  case 6:
    profileId = MolecularIds::eMalargueJuneModel;
    break;
  case 7:
    profileId = MolecularIds::eMalargueJulyModel;
    break;
  case 8:
    profileId = MolecularIds::eMalargueAugustModel;
    break;
  case 9:
    profileId = MolecularIds::eMalargueSeptemberModel;
    break;
  case 10:
    profileId = MolecularIds::eMalargueOctoberModel;
    break;
  case 11:
    profileId = MolecularIds::eMalargueNovemberModel;
    break;
  case 12:
    profileId = MolecularIds::eMalargueDecemberModel;
    break;
  default:
    ostringstream msg;
    msg << "Invalid month number : " << iMonth;
    ERROR(msg);
  }

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

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

      // Move calculation of refraction index to base class once
      // TabulatedFunction is fixed
      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()
      );

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

    // Fill the log depth vs height profile, starting from the top layer
    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) {
    ostringstream msg;
    msg << "Monthly average data requested from database, "
           "but no data found for month "
        << UTCDateTime(fCurrentTime).GetMonth();
    DEBUGLOG(msg);
    return false;
  }

  if (nZonesCheck != 1) {
    ERROR("some kind of error here");
  }

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