Framework/Atmosphere/Atmosphere.cc

Go to the documentation of this file.

#include <atm/Atmosphere.h>

#include <atm/VRayleighModel.h>
#include <atm/VProfileModel.h>
#include <atm/VMieModel.h>
#include <atm/VFluorescenceModel.h>
#include <atm/VCherenkovModel.h>
#include <atm/VCloudModel.h>

#include <atm/ScatteringResult.h>
#include <atm/AttenuationResult.h>
#include <atm/ProfileResult.h>
#include <atm/CloudResult.h>
#include <atm/InclinedAtmosphericProfile.h>

#include <atm/ModelRegister.h>

#include <det/Detector.h>
#include <fdet/Pixel.h>
#include <fdet/Telescope.h>

#include <fevt/Pixel.h>

#include <fwk/CentralConfig.h>

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>

#include <atm/AerosolDB.h>
#include <atm/MolecularDB.h>
#include <atm/OverallQualityDB.h>
#include <atm/LidarDB.h>
#include <atm/GOESDB.h>

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


namespace atm {

  template<class Factory>
  typename Factory::ObjectPtrType
  CreateAndInit(const Branch& topB, const string& modelType)
  {
    typename Factory::ObjectPtrType model = nullptr;
    const Branch branch = topB.GetChild(modelType);
    if (!branch) {
      ostringstream msg;
      msg << "Could not find branch " << modelType << ". "
             "This model will not be available";
      WARNING(msg);
    } else {
      const auto modelName = branch.Get<string>();
      model = Factory::Create(modelName);
      if (model) {
        model->Init();
      } else {
        ostringstream msg;
        msg << "Failed loading model " << modelName << ". "
               "Available models are: ";
        for (typename Factory::Iterator it = Factory::Begin();
             it != Factory::End(); ++it) {
          msg << it->first << ", ";
        }
        WARNING(msg);
      }
    }
    return model;
  }


  AttenuationResult
  Atmosphere::EvaluateRayleighAttenuation(const Point& xInit,
                                          const Point& xFinal,
                                          const vector<double>& wLength)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

     if (!fRayleighModel->HasData())
       throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighAttenuation(xInit, xFinal, wLength);
  }


  double
  Atmosphere::EvaluateRayleighAttenuation(const Point& xInit,
                                          const Point& xFinal,
                                          const double wLength)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

     if (!fRayleighModel->HasData())
       throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighAttenuation(xInit, xFinal, wLength);
  }


  AttenuationResult
  Atmosphere::EvaluateMieAttenuation(const Point& xInit,
                                     const Point& xFinal,
                                     const vector<double>& wLength)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieAttenuation(xInit, xFinal, wLength);
  }


  double
  Atmosphere::EvaluateMieAttenuation(const Point& xInit,
                                     const Point& xFinal,
                                     const double wLength)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieAttenuation(xInit, xFinal, wLength);
  }


  double
  Atmosphere::GetMieAttenuationLength(const Point& xA,
                                      const double wLength)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->GetAttenuationLength(xA, wLength);
  }


  double
  Atmosphere::GetVerticalAerosolOpticalDepth(const unsigned int eyeId,
                                             const double altitude)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->GetVerticalAerosolOpticalDepth(eyeId, altitude);
  }


  CloudResult
  Atmosphere::EvaluateCloudCoverage(const fdet::Pixel& pixel,
                                    const Point& x)
    const
  {
    if (!HasModel(fCloudModel))
      throw NonExistentComponentException();

    if (!fCloudModel->HasData())
      throw NoDataForModelException("Cloud model could not locate data.");

    return fCloudModel->EvaluateCloudCoverage(pixel, x);
  }


  CloudResult
  Atmosphere::EvaluateCloudCoverage(const fevt::Pixel& pixel,
                                    const Point& x)
    const
  {
    if (!HasModel(fCloudModel))
      throw NonExistentComponentException();

    if (!fCloudModel->HasData())
      throw NoDataForModelException("Cloud model could not locate data.");

    return fCloudModel->EvaluateCloudCoverage(pixel, x);
  }


  CloudResult
  Atmosphere::EvaluateCloudCoverage(const unsigned int eyeId,
                                    const unsigned int telId,
                                    const unsigned int pixId,
                                    const Point& x)
    const
  {
    if (!HasModel(fCloudModel))
      throw NonExistentComponentException();

    if (!fCloudModel->HasData())
      throw NoDataForModelException("Cloud model could not locate data.");

    return fCloudModel->EvaluateCloudCoverage(eyeId, telId, pixId, x);
  }


  const ProfileResult&
  Atmosphere::EvaluateDepthVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateDepthVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateHeightVsDepth()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateHeightVsDepth();
  }


  const ProfileResult&
  Atmosphere::EvaluatePressureVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluatePressureVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateTemperatureVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateTemperatureVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateVaporPressureVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateVaporPressureVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateDensityVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateDensityVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateRefractionIndexVsHeight()
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateRefractionIndexVsHeight();
  }


  const ProfileResult&
  Atmosphere::EvaluateRefractionIndexVsHeight(const double wavelength)
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->EvaluateRefractionIndexVsHeight(wavelength);
  }


  void
  Atmosphere::InitSlantProfileModel(const Point& core,
                                    const Vector& dir,
                                    const double deltaX)
    const
  {
    delete fInclinedProfileModel;
    try {
      fInclinedProfileModel = new InclinedAtmosphericProfile(core, dir, deltaX);
    } catch (...) {
      fInclinedProfileModel = nullptr;
      throw;
    }
  }


  double
  Atmosphere::IntegratedGrammage(const Point& pStart,
                                 const Point& pStop,
                                 const double delta)
    const
  {
    return InclinedAtmosphericProfile::IntegratedGrammage(pStart, pStop, delta);
  }


  const ProfileResult&
  Atmosphere::EvaluateSlantDepthVsDistance()
    const
  {
    if (!fInclinedProfileModel)
      throw InclinedAtmosphericProfile::InclinedAtmosphereModelException("First initialize the inclined atmosphere model!");

    return fInclinedProfileModel->EvaluateSlantDepthVsDistance();
  }


  const ProfileResult&
  Atmosphere::EvaluateDistanceVsSlantDepth()
    const
  {
    if (!fInclinedProfileModel)
      throw InclinedAtmosphericProfile::InclinedAtmosphereModelException("First initialize the inclined atmosphere model!");

    return fInclinedProfileModel->EvaluateDistanceVsSlantDepth();
  }


  const ProfileResult&
  Atmosphere::EvaluateHeightVsSlantDepth()
    const
  {
    if (!fInclinedProfileModel)
      throw InclinedAtmosphericProfile::InclinedAtmosphereModelException("First initialize the inclined atmosphere model!");

    return fInclinedProfileModel->EvaluateHeightVsSlantDepth();
  }


  const ProfileResult&
  Atmosphere::EvaluateHeightVsDistance()
    const
  {
    if (!fInclinedProfileModel)
      throw InclinedAtmosphericProfile::InclinedAtmosphereModelException("First initialize the inclined atmosphere model!");

    return fInclinedProfileModel->EvaluateHeightVsDistance();
  }


  ScatteringResult
  Atmosphere::EvaluateRayleighScattering(const Point& x1,
                                         const Point& x2,
                                         const double angle,
                                         const double distance,
                                         const vector<double>& wLength)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

    if (!fRayleighModel->HasData())
      throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighScattering(x1, x2, angle, distance, wLength);
  }


  ScatteringResult
  Atmosphere::EvaluateRayleighScattering(const Point& x1,
                                         const Point& x2,
                                         const double angle,
                                         const double distance,
                                         const AttenuationResult& rayleighAttenuation)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

    if (!fRayleighModel->HasData())
      throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighScattering(x1, x2, angle, distance, rayleighAttenuation);
  }


  double
  Atmosphere::EvaluateRayleighScattering(const Point& x1,
                                         const Point& x2,
                                         const double angle,
                                         const double distance,
                                         const double wLength)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

    if (!fRayleighModel->HasData())
      throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighScattering(x1, x2, angle, distance, wLength);
  }


  double
  Atmosphere::EvaluateRayleighScattering(const Point& x1,
                                         const Point& x2,
                                         const double angle,
                                         const double distance,
                                         const double wLength,
                                         const double rayleighAttenuation)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

    if (!fRayleighModel->HasData())
      throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->EvaluateRayleighScattering(x1, x2, angle, distance, wLength, rayleighAttenuation);
  }


  double
  Atmosphere::GetRayleighAttenuationLength(const Point& xA,
                                           const double wLength)
    const
  {
    if (!HasModel(fRayleighModel))
      throw NonExistentComponentException();

    if (!fRayleighModel->HasData())
      throw NoDataForModelException("Rayleigh model could not locate data");

    return fRayleighModel->GetAttenuationLength(xA, wLength);
  }


  ScatteringResult
  Atmosphere::EvaluateMieScattering(const Point& x1,
                                    const Point& x2,
                                    const double angle,
                                    const double distance,
                                    const vector<double>& wLength)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieScattering(x1, x2, angle, distance, wLength);
  }


  double
  Atmosphere::EvaluateMieScattering(const Point& x1,
                                    const Point& x2,
                                    const double angle,
                                    const double distance,
                                    const double wLength)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieScattering(x1, x2, angle, distance, wLength);
  }


  ScatteringResult
  Atmosphere::EvaluateMieScattering(const Point& x1,
                                    const Point& x2,
                                    const double angle,
                                    const double distance,
                                    const AttenuationResult& mieAttenuation)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieScattering(x1, x2, angle, distance, mieAttenuation);
  }


  double
  Atmosphere::EvaluateMieScattering(const Point& x1,
                                    const Point& x2,
                                    const double angle,
                                    const double distance,
                                    const double wLength,
                                    const double mieAttenuation)
    const
  {
    if (!HasModel(fMieModel))
      throw NonExistentComponentException();

    if (!fMieModel->HasData())
      throw NoDataForModelException("Mie model could not locate data.");

    return fMieModel->EvaluateMieScattering(x1, x2, angle, distance, wLength, mieAttenuation);
  }


  const TabulatedFunction&
  Atmosphere::EvaluateFluorescenceYield(const double heightAboveSeaLevel)
    const
  {
    if (!HasModel(fFluorescenceModel))
      throw NonExistentComponentException();

    return fFluorescenceModel->EvaluateFluorescenceYield(heightAboveSeaLevel);
  }


  const vector<double>&
  Atmosphere::GetWavelengths(const EmissionMode mode)
    const
  {
    switch (mode) {
    case eFluorescence:
      if (!HasModel(fFluorescenceModel))
        throw NonExistentComponentException();
      return fFluorescenceModel->GetWavelengths();
    case eCerenkov:
      if (!HasModel(fCherenkovModel))
        throw NonExistentComponentException();
      return fCherenkovModel->GetWavelengths();
    default:
      throw NonExistentComponentException("Asked for not-implemeted EmissionMode");
    }
  }


  double
  Atmosphere::GetdEdX0()
    const
  {
    if (!HasModel(fFluorescenceModel))
      throw NonExistentComponentException();

    return fFluorescenceModel->GetdEdX0();
  }


  double
  Atmosphere::GetDeExcitationTime(const double height)
    const
  {
    if (!HasModel(fFluorescenceModel))
      throw NonExistentComponentException();

    return fFluorescenceModel->GetDeExcitationTime(height);
  }


  double
  Atmosphere::GetVerticalTimeOfFlight(const double height1,
                                      const double height2)
    const
  {
    if (!HasModel(fProfileModel))
      throw NonExistentComponentException();

    if (!fProfileModel->HasData())
      throw NoDataForModelException("Profile model could not locate data");

    return fProfileModel->GetVerticalTimeOfFlight(height1, height2);
  }


  void
  Atmosphere::SetCherenkovEnergyCutoff(const double eCut)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->SetEnergyCutoff(eCut);
  }


  const TabulatedFunction&
  Atmosphere::EvaluateCherenkovPhotons(const Point& xA,
                                       const Point& xB,
                                       const double meanShowerAge)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->EvaluateCherenkovPhotons(xA, xB,
                                                     meanShowerAge);
  }


  const TabulatedFunction&
  Atmosphere::EvaluateCherenkovDirect(const Point& xA,
                                      const Point& xB,
                                      const Point& xEye,
                                      const double meanShowerAge)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->EvaluateCherenkovDirect(xA, xB, xEye,
                                                    meanShowerAge);
  }


  double
  Atmosphere::EvaluateDirectCherenkovProbability(const Point& xA,
                                                 const Point& xB,
                                                 const Point& xEye,
                                                 const double meanShowerAge)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->EvaluateDirectCherenkovProbability(xA, xB, xEye,
                                                               meanShowerAge);
  }

  double
  Atmosphere::EvaluateDirectCherenkovProbability(const Point& xA,
                                                 const Point& xB,
                                                 const Point& xEye,
                                                 const double meanShowerAge,
                                                 const double wavelength)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->EvaluateDirectCherenkovProbability(xA, xB, xEye,
                                                               meanShowerAge, wavelength);
  }

  const AerosolDB&
  Atmosphere::GetAerosolDB()
    const
  {
    if (!fAerosolDB)
      fAerosolDB = new AerosolDB;
    return *fAerosolDB;
  }


  const OverallQualityDB&
  Atmosphere::GetOverallQualityDB()
    const
  {
    if (!fOverallQualityDB)
      fOverallQualityDB = new OverallQualityDB;
    return *fOverallQualityDB;
  }


  const LidarDB&
  Atmosphere::GetLidarDB()
    const
  {
    if (!fLidarDB)
      fLidarDB = new LidarDB;
    return *fLidarDB;
  }

  const GOESDB&
  Atmosphere::GetGOESDB()
    const
  {
    if (!fGOESDB)
      fGOESDB = new GOESDB;
    return *fGOESDB;
  }


  const MolecularDB&
  Atmosphere::GetMolecularDB(const MolecularIds::ProfileId& id)
    const
  {
    // Check if requested MolecularDB is in the map.  If not, create it
    // and put it in the map.  This creation could be done in
    // the VAtmosphere constructor instead.

    MolecularDBMap::iterator it = fMolecularDBMap.find(id);
    if (it != fMolecularDBMap.end())
      return *it->second;
    else {
      MolecularDB* const m = new MolecularDB(id);
      fMolecularDBMap.insert(make_pair(id, m));
      return *m;
    }
  }


  bool
  Atmosphere::HasModel(VModel* const model)
  {
    if (!model) {
      ERROR("Call to a not registered or not configured model");
      return false;
    }
    return true;
  }


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

    if (!topB) {
      DEBUGLOG("Could not find branch Atmosphere; "
               "Atmosphere will not be available");
      return;
    }

    delete fProfileModel;
    fProfileModel = CreateAndInit<VProfileModelFactory>(topB, "ProfileModel");

    delete fRayleighModel;
    fRayleighModel = CreateAndInit<VRayleighModelFactory>(topB, "RayleighModel");

    delete fMieModel;
    fMieModel = CreateAndInit<VMieModelFactory>(topB, "MieModel");

    delete fFluorescenceModel;
    fFluorescenceModel = CreateAndInit<VFluorescenceModelFactory>(topB, "FluorescenceModel");

    delete fCherenkovModel;
    fCherenkovModel = CreateAndInit<VCherenkovModelFactory>(topB, "CherenkovModel");

    delete fCloudModel;
    fCloudModel = CreateAndInit<VCloudModelFactory>(topB, "CloudModel");
  }


  void
  Atmosphere::Clean()
  {
    delete fOverallQualityDB;
    fOverallQualityDB = nullptr;

    delete fLidarDB;
    fLidarDB = nullptr;

    delete fAerosolDB;
    fAerosolDB = nullptr;

    for (MolecularDBMap::iterator it = fMolecularDBMap.begin();
         it != fMolecularDBMap.end(); ++it)
      delete it->second;

    fMolecularDBMap.clear();

    delete fInclinedProfileModel;
    fInclinedProfileModel = nullptr;

    delete fGOESDB;
    fGOESDB = nullptr;
  }


  void
  Atmosphere::SetUncertaintyBound(const ModelWithUncertainty model,
                                  const double nSigma)
    const
  {
    switch (model) {
    case eMie:
      fMieModel->SetUncertaintyBound(nSigma);
      break;
    default:
      ERROR(" SetUncertaintyBound not implemented ...");
    }
  }


  double
  Atmosphere::AngularCherenkovCDF(const double theta,
                                  const double verticalDepth,
                                  const double showerAge)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->AngularCDF(theta, verticalDepth, showerAge);
  }


  double
  Atmosphere::AngularCherenkovPDF(const double theta,
                                  const double verticalDepth,
                                  const double showerAge)
    const
  {
    if (!HasModel(fCherenkovModel))
      throw NonExistentComponentException();

    return fCherenkovModel->AngularPDF(theta, verticalDepth, showerAge);
  }

}