source: trunk/Mars/msim/MAtmosphere.cc@ 19766

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appears in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MSimAtmosphere
//
//  Task to calculate wavelength or incident angle dependent absorption
//
//  Input Containers:
//   MPhotonEvent
//   MCorsikaRunHeader
//
//  Output Containers:
//   MPhotonEvent
//
//////////////////////////////////////////////////////////////////////////////
#include "MAtmosphere.h"

#include <fstream>

#include <TGraph.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MParList.h"

#include "MCorsikaRunHeader.h"
#include "MPhotonEvent.h"
#include "MPhotonData.h"

ClassImp(MAtmosphere);
ClassImp(MAtmRayleigh);

using namespace std;

// ==========================================================================
//
// January 2002, A. Moralejo: We now precalculate the slant paths for the
// aerosol and Ozone vertical profiles, and then do an interpolation in
// wavelength for every photon to get the optical depths. The parameters
// used, defined below, have been taken from "Atmospheric Optics", by
// L. Elterman and R.B. Toolin, chapter 7 of the "Handbook of geophysics
// and Space environments". (S.L. Valley, editor). McGraw-Hill, NY 1965.
// 
// WARNING: the Mie scattering and the Ozone absorption are implemented
// to work only with photons produced at a height a.s.l larger than the
// observation level. So this is not expected to work well for simulating
// the telescope pointing at theta > 90 deg (for instance for neutrino
// studies. Rayleigh scattering works even for light coming from below.
// 
// Fixed bugs (of small influence) in Mie absorption implementation: there
// were errors in the optical depths table, as well as a confusion:
// height a.s.l. was used as if it was height above the telescope level.
// The latter error was also present in the Ozone aborption implementation.
// 
// On the other hand, now we have the tables AE_ABI and OZ_ABI with optical
// depths between sea level and a height h (before it was between 2km a.s.l
// and a height h). So that we can simulate also in the future a different
// observation level.
// 
// AM: WARNING: IF VERY LARGE zenith angle simulations are to be done (say
// above 85 degrees, for neutrino primaries or any other purpose) this code
// will have to be adapted accordingly and checked, since in principle it has
// been written and tested to simulate the absorption of Cherenkov photons
// arriving at the telescope from above.
// 
// AM: WARNING 2: not to be used for wavelengths outside the range 250-700 nm
// 
// January 2003, Abelardo Moralejo: found error in Ozone absorption treatment.
// At large zenith angles, the air mass used was the one calculated for
// Rayleigh scattering, but since the Ozone distribution is rather different
// from the global distribution of air molecules, this is not a good
// approximation. Now I have left in this code only the Rayleigh scattering,
// and moved to atm.c the Mie scattering and Ozone absorption calculated in
// a better way.
// 
// A. Moralejo, January 2003: added some parameters for Mie scattering
// and Ozone absorption derived from the clear standard atmosphere model
// in "Atmospheric Optics", by L. Elterman and R.B. Toolin, chapter 7 of
// the "Handbook of geophysics and Space environments". S.L. Valley,
// editor. McGraw-Hill, NY 1965.
// 
// AM, Jan 2003: Changed the meaning of the argument height: now it is the
// true height above sea level at which a photon has been emitted, before
// it was the height given by Corsika, measured in the vertical of the
// observer and not in the vertical of the emitting particle.
//
//
// MAGIC-Winter and MAGIC-Summer by M. Haffke,
// parametrizing profiles obtained with MSIS:
// http://uap-www.nrl.navy.mil/models_web/msis/msis_home.htm
//
//
// The MAGIC-Winter and MAGIC-Summer parametrisations reproduce the MSIS
// profiles for the 3 atmospheric layers from 0 up to 40 km height. Beyond
// that height, it was not possible to achieve a good fit, but the amount
// of residual atmosphere is so small that light absorption would be
// negligible anyway. Showers develop well below 40 km.
//
//
// The mass overburden is given by T = AATM + BATM * exp(-h/CATM)
// The last layer of the US standard atmosphere (in which T varies
// linearly with h) is above 100 km and has not been included here
// because it should not matter.
//

const Double_t MAtmosphere::STEPTHETA = 1.74533e-2; // aprox. 1 degree

const Double_t MAtmRayleigh::fgMeanFreePath = 2970;

const Double_t MAtmosphere::aero_n[31] = {200, 87, 38, 16, 7.2, 3.1, 1.1, 0.4, 0.14, 5.0e-2, 2.6e-2, 2.3e-2, 2.1e-2, 2.3e-2, 2.5e-2, 4.1e-2, 6.7e-2, 7.3e-2, 8.0e-2, 9.0e-2, 8.6e-2, 8.2e-2, 8.0e-2, 7.6e-2, 5.2e-2, 3.6e-2, 2.5e-2, 2.4e-2, 2.2e-2, 2.0e-2, 1.9e-2};

const Double_t MAtmosphere::oz_conc[51]={0.3556603E-02, 0.3264150E-02, 0.2933961E-02, 0.2499999E-02, 0.2264150E-02, 0.2207546E-02, 0.2160377E-02, 0.2226414E-02, 0.2283018E-02, 0.2811320E-02, 0.3499999E-02, 0.4603772E-02, 0.6207545E-02, 0.8452828E-02, 0.9528299E-02, 0.9905657E-02, 0.1028302E-01, 0.1113207E-01, 0.1216981E-01, 0.1424528E-01, 0.1641509E-01, 0.1839622E-01, 0.1971697E-01, 0.1981131E-01, 0.1933962E-01, 0.1801886E-01, 0.1632075E-01, 0.1405660E-01, 0.1226415E-01, 0.1066037E-01, 0.9028300E-02, 0.7933960E-02, 0.6830187E-02, 0.5820753E-02, 0.4830188E-02, 0.4311319E-02, 0.3613206E-02, 0.3018867E-02, 0.2528301E-02, 0.2169811E-02, 0.1858490E-02, 0.1518867E-02, 0.1188679E-02, 0.9301884E-03, 0.7443394E-03, 0.5764149E-03, 0.4462263E-03, 0.3528301E-03, 0.2792452E-03, 0.2226415E-03, 0.1858490E-03};

MAtmRayleigh::MAtmRayleigh() : fR(MCorsikaRunHeader::fgEarthRadius),
    //fHeight{0, 300000, 1e+06, 5e+06, 1e+07},
    //fAtmA{-144.838, -124.071, 0.360027, -0.000824761, 0.00221589},
    fAtmB{1192.34, 1173.98, 1412.08, 810.682, 1},
    fAtmC{994186, 967530, 636143, 735640, 4.92961e9},
    fObsLevel(-1)
{
}

// --------------------------------------------------------------------------
//
// Precalcalculate the integrals from the observer level to the next
// atmpsheric layer below including the lower boundary. Thus a
// correct calculation is reduced to the calculation of the upper
// boundary.
//
// fRho[0] = B0;
// fRho[1] = B0-A0 + B1;
// fRho[2] = B0-A0 + B1-A1 + B2;
// fRho[3] = B0-A0 + B1-A1 + B2+A2 + B3;
// fRho[4] = B0-A0 + B1-A1 + B2+A2 + B3 - A3;
//
void MAtmRayleigh::PreCalcRho()
{
    // Limits (height in cm) of the four layers in which
    // atmosphere is parametrized.
    // This is a stupid trick giving 0 for the integrals below
    // the observer

    // FIXME: Could be replaced by 0, AtmLay[0]-fAtmLay[3]

    const Double_t h[5] =
    {
        fObsLevel,                     // fObsLevel,                   //   0km
        TMath::Max(fObsLevel, 7.75e5), // TMath::Max(fObsLevel, 4e5),  //   4km
        16.5e5,                        //  10e5,                       //  10km
        50.0e5,                        //  40e5,                       //  40km
        105.0e5,                       // 100e5                        // 100km
    };

    memcpy(fHeight, h, sizeof(Double_t)*5);

    fRho[0] = 0;
    for (int i=0; i<4; i++)
    {
        const Double_t b = fAtmB[i];
        const Double_t c = fAtmC[i];

        const Double_t h1 = h[i+1];
        const Double_t h0 = h[i];

        const Double_t B = b*TMath::Exp(-h0/c);
        const Double_t A = b*TMath::Exp(-h1/c);

        // Calculate rho for the i-th layer from the lower
        // to the higher layer boundary.
        // If height is within the layer only calculate up to height.
        fRho[i]  += B;
        fRho[i+1] = fRho[i] - A;
    }
}

void MAtmRayleigh::Init(Double_t obs, const Float_t *atmb, const Float_t *atmc)
{
    // Observation level above earth radius
    fObsLevel = obs;

    //memcpy(fAtmA, (Float_t*)h.GetAtmosphericCoeffA(), sizeof(Float_t)*4);
    if (atmb)
        memcpy(fAtmB, atmb, sizeof(Float_t)*5);
    if (atmc)
        memcpy(fAtmC, atmc, sizeof(Float_t)*5);

    PreCalcRho();
}

// Init an atmosphere from the data stored in MCorsikaRunHeader
// This initialized fObsLevel, fR, fAtmB and fAtmC and
// PreCalcRho
void MAtmRayleigh::Init(const MCorsikaRunHeader &h)
{
    // Use earth radius as defined in Corsika
    fR = h.EarthRadius();

    Init(h.GetObsLevel(), h.GetAtmosphericCoeffB(), h.GetAtmosphericCoeffC());
}

// Return the vertical thickness between the observer and height.
// Therefor the integral of the layers below (including the lower
// boudary) Have been precalculated by PreCalcRho
Double_t MAtmRayleigh::GetVerticalThickness(Double_t height) const
{
    // FIXME: We could store the start point above obs-level
    //        (Does this really gain anything?)
    Int_t i=0;
    while (i<4 && height>fHeight[i+1])
        i++;

    const Double_t b = fAtmB[i];
    const Double_t c = fAtmC[i];

    return fRho[i] - b*TMath::Exp(-height/c);
}

/*
// The "orginal" code for the vertical thickness
Double_t GetVerticalThickness(Double_t obslev, Double_t height) const
{
    // This is a C++-version of the original code from attenu.c

    // Limits (height in cm) of the four layers in which atmosphere is parametrized:
    const double lahg[5] =
    {
        obslev,
        TMath::Max(obslev, 4e5),
        1.0e6,
        4.0e6,
        1.0e7
    };

    Double_t Rho_Tot = 0.0;
    for (int i=0; i<4; i++)
    {
        const Double_t b = fAtmB[i];
        const Double_t c = fAtmC[i];

        const Double_t h0 = TMath::Min(height, lahg[i+1]);
        const Double_t h1 =                    lahg[i];

        // Calculate rho for the i-th layer from the lower
        // to the higher layer boundary.
        // If height is within the layer only calculate up to height.
        Rho_Tot += b*(exp(-h1/c) - exp(-h0/c));

        if (lahg[i+1] > height)
            break;
    }

    return Rho_Tot;
}
*/

Double_t MAtmRayleigh::CalcTransmission(Double_t height, Double_t wavelength, Double_t sin2) const
{
    // sin2: sin(theta)^2
    // height is the true height a.s.l.

    // LARGE ZENITH ANGLE FACTOR (AIR MASS FACTOR):
    // Air mass factor "airmass" calculated using a one-exponential
    // density profile for the atmosphere,
    //
    //     rho = rho_0 exp(-height/hscale) with hscale = 7.4 km
    //
    // The air mass factor is defined as I(theta)/I(0), the ratio of
    // the optical paths I (in g/cm2) traversed between two given
    // heights, at theta and at 0 deg z.a.

    const Double_t H = height-fObsLevel;
    const Double_t h = 2*H;

    // Scale-height (cm) for Exponential density profile
    const Double_t hscale = 7.4e5;
    const Double_t f      = 2*hscale;

    // Precalc R*cos(theta)^2 (FIXME: Is ph.GetCosW2 more precise?)
    const Double_t Rcos2 = fR * (1-sin2); // cos2 = 1 - sin2

    const Double_t x1 = TMath::Sqrt((Rcos2      ) / f);
    const Double_t x2 = TMath::Sqrt((Rcos2 + 2*h) / f);
    const Double_t x3 = TMath::Sqrt((fR         ) / f);
    const Double_t x4 = TMath::Sqrt((fR    + 2*h) / f);

    // Return a -1 transmittance in the case the photon comes
    // exactly from the observation level, because in that case the
    // "air mass factor" would become infinity and the calculation
    // is not valid!
    if (fabs(x3-x4) < 1.e-10)
        return  -1.;

    const Double_t e12 = erfc(x1) - erfc(x2);
    const Double_t e34 = erfc(x3) - erfc(x4);

    const Double_t airmass = TMath::Exp(-fR*sin2 / f) * e12/e34;

    // Calculate the traversed "vertical thickness" of air using the
    // US Standard atmosphere:
    const Double_t Rho_tot = GetVerticalThickness(/*fObsLevel,*/ height);

    // We now convert from "vertical thickness" to "slanted thickness"
    // traversed by the photon on its way to the telescope, simply
    // multiplying by the air mass factor m:
    const Double_t Rho_Fi = airmass * Rho_tot;

    // Finally we calculate the transmission coefficient for the Rayleigh
    // scattering:
    // AM Dec 2002, introduced ABS below to account (in the future) for
    // possible photons coming from below the observation level.

    const Double_t wl = 400./wavelength;

    // Mean free path for scattering Rayleigh XR (g/cm^2)
    return TMath::Exp(-TMath::Abs(Rho_Fi/fgMeanFreePath)*wl*wl*wl*wl);
}

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

// Interpolate the graph at wavelength
Double_t MAtmosphere::GetBeta(Double_t wavelength, const TGraph &g) const
{
    // FIXME: This might not be the fastest because range
    // checks are done for each call!
    return g.GetN()==0 ? 0 : g.Eval(wavelength)*1e-5; // from km^-1 to cm^-1
/*
        // Linear interpolation
        // (FIXME: Is it faster to be replaced with a binary search?)
        // (       This might be faster because we have more photons
        //         with smaller wavelengths)
        //int index;
        //for (index = 1; index <g.GetN()-1; index++)
        //    if (wavelength < g.GetX()[index])
        //        break;
        const Int_t index = TMath::BinarySearch(g.GetN(), g.GetX(), wavelength)+1;

        const Double_t t0 = g.GetY()[index-1];
        const Double_t t1 = g.GetY()[index];

        const Double_t w0 = g.GetX()[index-1];
        const Double_t w1 = g.GetX()[index];

        const Double_t beta0 = t0+(t1-t0)*(wavelength-w0)/(w1-w0);

        return beta0 * 1e-5; // from km^-1 to cm^-1
*/
}

MAtmosphere::~MAtmosphere()
{
    if (fAbsCoeffOzone)
        delete fAbsCoeffOzone;
    if (fAbsCoeffAerosols)
        delete fAbsCoeffAerosols;
}

Float_t MAtmosphere::GetWavelengthMin() const
{
    return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Max(fAbsCoeffOzone->GetX()[0], fAbsCoeffAerosols->GetX()[0]) : -1;
}

Float_t MAtmosphere::GetWavelengthMax() const
{
    return fAbsCoeffOzone && fAbsCoeffAerosols ? TMath::Min(fAbsCoeffOzone->GetX()[fAbsCoeffOzone->GetN()-1], fAbsCoeffAerosols->GetX()[fAbsCoeffAerosols->GetN()-1]) : -1;
}

Bool_t MAtmosphere::HasValidOzone() const
{
    return fAbsCoeffOzone && fAbsCoeffOzone->GetN()>0;
}

Bool_t MAtmosphere::HasValidAerosol() const
{
    return fAbsCoeffAerosols && fAbsCoeffAerosols->GetN()>0;
}

void MAtmosphere::PreCalcOzone()
{
    // It follows a precalculation of the slant path integrals we need
    // for the estimate of the Mie scattering and Ozone absorption:
    Double_t dh = 1.e3;
    //const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree

    // Ozone absorption
    for (Int_t j = 0; j < 90; j++)
    {
        const Double_t theta = j * STEPTHETA;
        const Double_t sin2  = sin(theta)*sin(theta);
        const Double_t H     = R()+fObsLevel;

        Double_t path_slant = 0;
        for (Double_t h=fObsLevel; h<50e5+dh/2; h+=dh)
        {
            // h is the true height vertical above ground
            //if (fmod(h,1e4) == 0)
            {
                ozone_path[TMath::FloorNint(h/1e4)][j] = path_slant;
            }

            const Double_t km  = h/1e5;
            const Int_t    i   = TMath::FloorNint(km);
            const Double_t l   = R()+h;

            const Double_t L   = TMath::Sqrt(l*l - H*H * sin2);
            const Double_t f   = dh * l / L;

            // Linear interpolation at h/1e5
            Double_t interpol = oz_conc[i] + fmod(km, 1) * (oz_conc[i+1]-oz_conc[i]);

            path_slant += f * interpol;
        }
    }
}

void MAtmosphere::PreCalcAerosol()
{
    // It follows a precalculation of the slant path integrals we need
    // for the estimate of the Mie scattering and Ozone absorption:
    Double_t dh = 1.e3;
    //const Double_t STEPTHETA = 1.74533e-2; // aprox. 1 degree

    /* Mie (aerosol): */
    for (Int_t j = 0; j < 90; j++)
    {
        const Double_t theta = j * STEPTHETA;
        const Double_t sin2  = sin(theta)*sin(theta);
        const Double_t H     = R()+fObsLevel;

        Double_t path_slant = 0;
        for (Double_t h=fObsLevel; h<=30e5+dh/2; h+=dh)
        {
            // h is the true height vertical above ground
            //if (fmod(h,1e4) == 0)
            {
                aerosol_path[TMath::FloorNint(h/1e4)][j] = path_slant;
            }

            const Double_t km  = h/1e5;
            const Int_t    i   = TMath::FloorNint(km);
            const Double_t l   = R()+h;

            const Double_t L   = TMath::Sqrt(l*l - H*H * sin2);
            const Double_t f   = dh * l / L;

            // Linear interpolation at h/1e5
            Double_t interpol = aero_n[i] + fmod(km, 1)*(aero_n[i+1]-aero_n[i]);

            path_slant += f * interpol/aero_n[0];    // aero_n [km^-1]
        }
    }
}

Bool_t MAtmosphere::InitOzone(const TString name)
{
    if (!name.IsNull())
    {
        if (fAbsCoeffOzone)
            delete fAbsCoeffOzone;

        fAbsCoeffOzone = new TGraph(name);
        fAbsCoeffOzone->Sort();
    }

    if (!HasValidOzone())
        return kFALSE;

    if (IsValid())
        PreCalcOzone();

    return kTRUE;
}

Bool_t MAtmosphere::InitAerosols(const TString name)
{
    if (!name.IsNull())
    {
        if (fAbsCoeffAerosols)
            delete fAbsCoeffAerosols;

        fAbsCoeffAerosols = new TGraph(name);
        fAbsCoeffAerosols->Sort();
    }

    if (!HasValidAerosol())
        return kFALSE;

    if (IsValid())
        PreCalcAerosol();

    return kTRUE;
}

/*
Double_t GetOz(Double_t height, Double_t theta) const
{
    // Distance between two points D = 1km /cos(theta)
    // Density along y within this km:   f =  (x[i+1]-x[i])/1km * dy
    // Integral of this density  f =  (x[i+1]-x[i])/1km * (y[i+1]-y[i])
    // f(h) = int [ (c1-c0)/1km*(h-h0)*dh + c0 ] dh
    //      = (c-co)*(h-h0)

    Double_t rc = 0;
    int i;
    for (i=0; i<49; i++)
        if (i>=2 && i+1<height/1e5)    // cm -> km
            rc += oz_conc[i] * 1e5/cos(theta);

    rc -= oz_conc[2]*0.2*1e5/cos(theta);
    rc += oz_conc[i+1]*fmod(height/1e5,1)*1e5/cos(theta);

    return rc;
}
*/

Double_t MAtmosphere::CalcOzoneAbsorption(Double_t h, Double_t wavelength, Double_t theta) const
{
    if (!fAbsCoeffOzone)
        return 1;

    //******* Ozone absorption *******

    // Vigroux Ozone absorption coefficient a.s.l. through interpolation:
    //const float oz_vigroux[15]= {1.06e2, 1.01e1, 8.98e-1, 6.40e-2, 1.80e-3, 0, 0, 3.50e-3, 3.45e-2, 9.20e-2, 1.32e-1, 6.20e-2, 2.30e-2, 1.00e-2, 0.00};
    //const Double_t beta0 = getbeta(wavelength, oz_vigroux);
    const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffOzone);

    // Now use the pre-calculated values of the path integral
    // for h and theta
    const UInt_t H = TMath::Min(500, TMath::Nint(h/1e4));
    const UInt_t T = TMath::Min( 89, TMath::Nint(theta/STEPTHETA));

    const Double_t path = ozone_path[H][T];

    return TMath::Exp(-beta0*path);
}

Double_t MAtmosphere::CalcAerosolAbsorption(Double_t h, Double_t wavelength, Double_t theta) const
{
    if (!fAbsCoeffAerosols)
        return 1;

    //******* Mie (aerosol) *******

    //const float aero_betap[15] = {0.27, 0.26, 0.25, 0.24, 0.24, 0.23, 0.20, 0.180, 0.167, 0.158, 0.150, 0.142, 0.135, 0.127, 0.120};
    //const Double_t beta0 = getbeta(wavelength, aero_betap);
    const Double_t beta0 = GetBeta(wavelength, *fAbsCoeffAerosols);

    // Now use the pre-calculated values of the path integral
    // for h and theta
    const UInt_t H = TMath::Min(300, TMath::Nint(h/1e4));
    const UInt_t T = TMath::Min( 89, TMath::Nint(theta/STEPTHETA));

    const Double_t path = aerosol_path[H][T];

    return TMath::Exp(-beta0*path);
}

Double_t MAtmosphere::GetTransmission(const MPhotonData &ph) const
{
    const Double_t wavelength = ph.GetWavelength();
    const Double_t height     = ph.GetProductionHeight();

    // Reduce the necessary number of floating point operations
    // by storing the intermediate results
    const Double_t sin2  = ph.GetSinW2();
    const Double_t cost  = TMath::Sqrt(1-sin2);
    const Double_t theta = TMath::ACos(cost);

    // Path from production height to obslevel
    const Double_t z = height-fObsLevel;

    // Distance of emission point to incident point on ground
    const Double_t d = z/cost;

    // Avoid problems if photon is very close to telescope:
    if (TMath::Abs(d)<1)
        return 1;

    // Earth radius plus observation height (distance of telescope
    // from earth center)
    const Double_t H = R() + fObsLevel;

    // We calculate h, the true height a.s.l.
    // of the photon emission point in cm
    const Double_t h = TMath::Sqrt(H*H + d*d + 2*H*z) - R();

    //**** Rayleigh scattering: *****
    const Double_t T_Ray = CalcTransmission(h, wavelength, sin2);
    if (T_Ray<0)
        return 0;

    //****** Ozone absorption: ******
    const Double_t T_Oz  = CalcOzoneAbsorption(h, wavelength, theta);

    //******** Mie (aerosol) ********
    const Double_t T_Mie = CalcAerosolAbsorption(h, wavelength, theta);

    // FIXME: What if I wanna display these values?

    // Calculate final transmission coefficient
    return T_Ray * T_Oz * T_Mie;
}