#ifndef MARS_MAtmosphere
#define MARS_MAtmosphere

#include <TROOT.h>

class TGraph;
class MPhotonData;
class MCorsikaRunHeader;

class MAtmRayleigh
{
private:
    static const Double_t fgMeanFreePath; // [g/cm^2] Mean free path for scattering Rayleigh XR
    static const Double_t fgEarthRadiush; // [cm] Default Earth Radius

    Double_t fR;         // [cm] Earth radius to be used

    Float_t fHeight[5]; // Layer boundaries (atmospheric layer)

    // Parameters of the different atmospheres. We use the same parametrization
    // shape as in Corsika atmospheric models (see Corsika manual, appendix D).
    // The values here can be/are obtained from the Corsika output
    // fAtmA is not required as it is an offset which cancels out
    //Float_t  fAtmA[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
    Float_t  fAtmB[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)
    Float_t  fAtmC[4];   // The index refers to the atmospheric layer (starting from sea level and going upwards)

    Double_t fRho[5];    // Precalculated integrals for rayleigh scatterning

    void PreCalcRho();

protected:
    Float_t fObsLevel; // [cm] observation level a.s.l.

public:
    // Default constructor
    MAtmRayleigh();

    // Init an atmosphere from the data stored in MCorsikaRunHeader
    MAtmRayleigh(const MCorsikaRunHeader &h)
    {
        Init(h);
    }

    // Check if the ovservation level has been correctly initialized
    // Used as a marker for correct initialization
    Bool_t IsValid() const { return fObsLevel>=0; }

    // Get the Earth radius to be used
    Double_t R() const { return fR; }

    void Init(Float_t obs);
    void Init(Float_t obs, const Float_t *atmb, const Float_t *atmc, const Float_t *height);
    void Init(const MCorsikaRunHeader &h);

    Double_t GetVerticalThickness(Double_t height) const;
    Double_t CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const;

    static Double_t Fcn(const Double_t *xx, const Double_t *par)
    {
        const Double_t *A = par;
        const Double_t *B = par+5;
        const Double_t *C = par+10;

        double H[6] = { 0, 0, 0, 0, 0, 1e100 };

        const Double_t x = xx[0]*1e5; // km -> cm

        for (int i=0; i<4; i++)
        {
            H[i+1] = log(C[i+1]/C[i] * B[i]/B[i+1]) / (1./C[i] - 1./C[i+1]);

            if (x>=H[i] && x<H[i+1])
                return log10(A[i] + B[i] * exp(-x/C[i]));
        }

        return -1;
    }
};

// ==========================================================================

class MAtmosphere : public MAtmRayleigh
{
private:
    static const Double_t STEPTHETA; // aprox. 1 degree

    // Aerosol number density for 31 heights a.s.l., from 0 to 30 km,
    // in 1 km steps (units: cm^-3)
    static const Double_t aero_n[31];

    // Ozone concentration for 51 heights a.s.l., from 0 to 50 km,
    // in 1 km steps (units: cm/km)
    static const Double_t oz_conc[51];

    // aerosol_path contains the path integrals for the aerosol number
    // density (relative to the number density at sea level) between the
    // observation level and a height h for different zenith angles. The
    // first index indicate height above sea level in units of 100m, the
    // second is the zenith angle in degrees.
    float aerosol_path[301][90];

    // ozone_path contains the path integrals for the ozone concentration
    // between the observation level and a height h for different zenith
    // angles. The first index indicate height above sea level in units
    // of 100m, the second is the zenith angle in degrees.
    float ozone_path[501][90];

    // Interpolate the graph at wavelength
    Double_t GetBeta(Double_t wavelength, const TGraph &g) const;

    //MSpline3 *fAbsCoeffOzone;
    //MSpline3 *fAbsCoeffAerosols;

    TGraph *fAbsCoeffOzone;
    TGraph *fAbsCoeffAerosols;

public:
    MAtmosphere(const MCorsikaRunHeader &h) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
    {
        Init(h);//, "ozone.txt", "aerosols.txt");
    }
    MAtmosphere(Float_t obs=0, const char *name1=0, const char *name2=0) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
    {
        MAtmRayleigh::Init(obs);

        if (name1)
            InitOzone(name1);
        if (name2)
            InitAerosols(name2);
    }
    MAtmosphere(const char *name1, const char *name2) : fAbsCoeffOzone(0), fAbsCoeffAerosols(0)
    {
        MAtmRayleigh::Init(0);

        if (name1)
            InitOzone(name1);
        if (name2)
            InitAerosols(name2);
    }

    ~MAtmosphere();

    Float_t GetWavelengthMin() const;
    Float_t GetWavelengthMax() const;

    Bool_t HasValidOzone() const;
    Bool_t HasValidAerosol() const;

    Bool_t IsAllValid() const { return IsValid() && HasValidOzone() && HasValidAerosol(); }

    void PreCalcOzone();
    void PreCalcAerosol();
    Bool_t InitOzone(const TString name="");
    Bool_t InitAerosols(const TString name="");

    void Init(const MCorsikaRunHeader &h, const char *name1=0, const char *name2=0)
    {
        MAtmRayleigh::Init(h);

        InitOzone(name1);
        InitAerosols(name2);
    }

    Double_t CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const;
    Double_t CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const;
    Double_t GetTransmission(const MPhotonData &ph) const;
};

#endif