#ifndef _atm_Atmosphere_h_
#define _atm_Atmosphere_h_

#include <atm/VModel.h>
#include <atm/MolecularIds.h>
#include <utl/AugerUnits.h>

#include <vector>
#include <map>

static const char CVSId_atm_Atmosphere[] =
  "$Id$";


namespace det {
  class Detector;
}

namespace fdet {
  class Pixel;
}

namespace fevt {
  class Pixel;
}

namespace utl {
  class Point;
  class Vector;
  class TabulatedFunction;
}

namespace atm {

  class AttenuationResult;
  class ScatteringResult;
  class CloudResult;
  class ProfileResult;
  class VRayleighModel;
  class VProfileModel;
  class VMieModel;
  class VFluorescenceModel;
  class VCherenkovModel;
  class VCloudModel;
  class InclinedAtmosphericProfile;
  class AerosolDB;
  class MolecularDB;
  class OverallQualityDB;
  class LidarDB;
  class GOESDB;

  class Atmosphere {

  public:
    enum EmissionMode {
      /// Fluorescence emission model
      eFluorescence = 0,
      /// Cerenkov emission model
      eCerenkov
    };

    Atmosphere();

    ~Atmosphere() { Clean(); }

    void Update() { Clean(); }

    /** @name Rayleigh models */
    //@{

    /// Compute Rayleigh attenuation between points
    /** Evaluate the attenuation  between two points
        - \param xInit   : initial point
        - \param xFinal  : final point
        - \param wLength : spectrum of inital photons
    */
    atm::AttenuationResult EvaluateRayleighAttenuation(const utl::Point& xInit,
                                                       const utl::Point& xFinal,
                                                       const std::vector<double>& wLength) const;

    double EvaluateRayleighAttenuation(const utl::Point& xInit,
                                       const utl::Point& xFinal,
                                       const double wLength) const;

    atm::ScatteringResult EvaluateRayleighScattering(const utl::Point& xA,
                                                     const utl::Point& xB,
                                                     const double angle, const double distance,
                                                     const std::vector<double>& xLength) const;

    double EvaluateRayleighScattering(const utl::Point& xA,
                                      const utl::Point& xB,
                                      const double angle, const double distance,
                                      const double xLength) const;

    atm::ScatteringResult EvaluateRayleighScattering(const utl::Point& xA,
                                                     const utl::Point& xB,
                                                     const double angle, const double distance,
                                                     const atm::AttenuationResult& rayleighAttenuation) const;

    double EvaluateRayleighScattering(const utl::Point& xA,
                                      const utl::Point& xB,
                                      const double angle, const double distance,
                                      const double xLength,
                                      const double rayleighAttenuation) const;

    double GetRayleighAttenuationLength(const utl::Point& xA,
                                        const double wLength) const;

    //@}

    /** @name Mie models */
    //@{

    atm::AttenuationResult EvaluateMieAttenuation(const utl::Point& xInit,
                                                  const utl::Point& xFinal,
                                                  const std::vector<double>& wLength) const;

    double EvaluateMieAttenuation(const utl::Point& xInit,
                                  const utl::Point& xFinal,
                                  const double wLength) const;

    atm::ScatteringResult EvaluateMieScattering(const utl::Point& xA,
                                                const utl::Point& xB,
                                                const double angle, const double distance,
                                                const std::vector<double>& xLength) const;

    double EvaluateMieScattering(const utl::Point& xA,
                                 const utl::Point& xB,
                                 const double angle, const double distance,
                                 const double xLength) const;

    atm::ScatteringResult EvaluateMieScattering(const utl::Point& xA,
                                                const utl::Point& xB,
                                                const double angle, const double distance,
                                                const atm::AttenuationResult& mieAttenuation) const;

    double EvaluateMieScattering(const utl::Point& xA,
                                 const utl::Point& xB,
                                 const double angle, const double distance,
                                 const double xLength,
                                 const double mieAttenuation) const;

    double GetMieAttenuationLength(const utl::Point& xA,
                                   const double wLength) const;


    /// Retrieve vertical optical depth at some location and altitude.
    double GetVerticalAerosolOpticalDepth(const unsigned int eyeId,
                                          const double altitude = 4.5*utl::km) const;

    //@}

    /** @name Cloud obscuration models */
    //@{

    atm::CloudResult EvaluateCloudCoverage(const fdet::Pixel& pix,
                                           const utl::Point& x) const;

    atm::CloudResult EvaluateCloudCoverage(const fevt::Pixel& pix,
                                           const utl::Point& x) const;

    atm::CloudResult EvaluateCloudCoverage(const unsigned int eyeId,
                                           const unsigned int telId,
                                           const unsigned int pixId,
                                           const utl::Point& x) const;

    //@}

    /** @name Molecular profile models */
    //@{

    /// Tabulated function giving depth as a function of height
    const atm::ProfileResult& EvaluateDepthVsHeight() const;

    /// Tabulated function giving height as a function of depth
    const atm::ProfileResult& EvaluateHeightVsDepth() const;

    /// Tabulated function giving air pressure as a function of height
    const atm::ProfileResult& EvaluatePressureVsHeight() const;

    /// Tabulated function giving H20 vapor pressure as a function of height
    const atm::ProfileResult& EvaluateVaporPressureVsHeight() const;

    /// Tabulated function giving temperature as a function of height
    const atm::ProfileResult& EvaluateTemperatureVsHeight() const;

    /// Tabulated function giving density as a function of height
    const atm::ProfileResult& EvaluateDensityVsHeight() const;

    /// Tabulated function giving refraction index as a function of height
    const atm::ProfileResult& EvaluateRefractionIndexVsHeight() const;

    /// Tabulated function giving refraction index as a function of height and wavelength
    const atm::ProfileResult& EvaluateRefractionIndexVsHeight(const double wavelength) const;

    //@}

    /** @name Slant profile model (non-vertical geometry)
        In order to use these methods, first call the InitSlantProfileModel method.
        Otherwise, an exception will be thrown.
        Note also that there is no method to get the distance given a height.
        This is on purpose.
     */
    //@{
    /**
       Initialize the slant depth model for simulations
       @param[in] core The point from which distances are measured.
       @param[in] dir The direction along the path, starting at the core. This can be the negativ of the shower direction, for example.
       @param[in] deltaX The step size to do the integration in depth.
     */
    void InitSlantProfileModel (const utl::Point& core,
                                const utl::Vector& dir,
                                double deltaX) const;

    /**
       Integral of the density along a line connecting two points.
       @param[in] delta The step size to do the integration in depth.
     */
    double IntegratedGrammage(const utl::Point& pStart,
                              const utl::Point& pStop,
                              double delta) const;

    /**
       Table of slant depth as a function of distance
       \attention{The distance is positive in the \b opposite direction given by the direction vector provided when initializing the profile model.}
     */
    const atm::ProfileResult& EvaluateSlantDepthVsDistance() const;

    /// Table of distance as a function of slant depth. \attention{The distance is positive in the \b opposite direction given by the direction vector provided when initializing the profile model.}
    const atm::ProfileResult& EvaluateDistanceVsSlantDepth() const;

    /// Table of height as a function of slant depth
    const atm::ProfileResult& EvaluateHeightVsSlantDepth() const;

    /// Table of height as a function of distance
    const atm::ProfileResult& EvaluateHeightVsDistance() const;
    //@}

    /** @name Fluorescence models */
    //@{

    /// Evaluated Fluorescence Yield for a specific model
    /** The fluorescence yield is usually given as
      * y= dEdX/dEdX0 * f(T,P)
      * This function will return f(T,P) */
    const utl::TabulatedFunction&
    EvaluateFluorescenceYield(const double heightAboveSeaLevel) const;

    const std::vector<double>& GetWavelengths(const EmissionMode mode = eFluorescence) const;

    /// get reference energy deposit for fluorescence yield model
    double GetdEdX0() const;

    /// get fluorescence de-excitation time
    double GetDeExcitationTime(const double height) const;

    /// Evaluate light signal time-of-flight between two altitudes
    double GetVerticalTimeOfFlight(const double height1,
                                   const double height2) const;

    //@}

    /** @name Cherenkov models */
    //@{
    void SetCherenkovEnergyCutoff(const double eCut) const;

    /// cumulative of angular Cherenkov distribution from 0 to theta
    double AngularCherenkovCDF(const double theta, const double verticalDepth,
                               const double showerAge) const;

    /// angular Cherenkov light distribution
    double AngularCherenkovPDF(const double theta, const double verticalDepth,
                               const double showerAge) const;

    const utl::TabulatedFunction&
    EvaluateCherenkovPhotons(const utl::Point& xA,
                             const utl::Point& xB,
                             const double showerAge) const;

    const utl::TabulatedFunction&
    EvaluateCherenkovDirect(const utl::Point& xA,
                            const utl::Point& xB,
                            const utl::Point& xEye,
                            const double showerAge) const;

    double
    EvaluateDirectCherenkovProbability(const utl::Point& xA,
                                       const utl::Point& xB,
                                       const utl::Point& xEye,
                                       const double showerAge) const;
    //@}

    enum ModelWithUncertainty {
      eMie
    };

    /// alter Model "model" by "nSigma" standard deviations
    void SetUncertaintyBound(const ModelWithUncertainty model,
                             const double nSigma) const;

    /** @name Low-level database interfaces */
    //@{

    /// Low-level interface to portion of the database with aerosol information
    /** For an example of using this interface, see the ReadAerosolLLNS::ReadAerosolLL example. */
    const AerosolDB& GetAerosolDB() const;

    /// Low-level interface to portion of the database with molecular information
    /** For an example of using this interface, see the ReadMolecularLLNS::ReadMolecularLL example. */
    const MolecularDB& GetMolecularDB(const MolecularIds::ProfileId& id) const;

    /// low-level interface to portion of the database summarizing overall quality
    /** For an example of using this interface, see the ReadAtmQualityLLNS::ReadAtmQualityLL example. */
    const OverallQualityDB& GetOverallQualityDB() const;

    /// low-level interface to portion of the database with cloud data from lidar analysis
    /** For an example of using this interface, see the ReadLidarLLNS::ReadLidarLL example. */
    const LidarDB& GetLidarDB() const;

    /// low-level interface to the cloud information from the GOES database
    /** For an example of using this interface, see the ReadGOESLL example. */
    const GOESDB& GetGOESDB() const;

    //@}

  private:
    Atmosphere(const Atmosphere&);
    Atmosphere& operator=(const Atmosphere&);

    void Init();

    void Clean();

    VRayleighModel* fRayleighModel;
    VProfileModel* fProfileModel;
    VMieModel* fMieModel;
    VFluorescenceModel* fFluorescenceModel;
    VCherenkovModel* fCherenkovModel;
    VCloudModel* fCloudModel;
    mutable InclinedAtmosphericProfile* fInclinedProfileModel;

    bool HasModel(VModel* const model) const;

    mutable AerosolDB* fAerosolDB;
    mutable OverallQualityDB* fOverallQualityDB;
    mutable LidarDB* fLidarDB;
    mutable GOESDB* fGOESDB;

    typedef std::map<MolecularIds::ProfileId, MolecularDB*> MolecularDBMap;
    mutable MolecularDBMap fMolecularDBMap;

    friend class det::Detector;

  };

} // atm


#endif // _atm_Atmosphere_h_

// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. Atmosphere/Atmosphere.o  -k"
// End: