source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/common/atmosphere/src/LinsleyAtmosphere.cc @ 117

// $Id: LinsleyAtmosphere.cc 2695 2006-05-27 19:09:42Z moreggia $
// S. Moreggia created 27 October 2003

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: LinsleyAtmosphere                                                           *
 *  Package: atmosphere                                                      *
 *  Coordinator: S. Moreggia                                                 *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// LinsleyAtmosphere
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "LinsleyAtmosphere.hh"
#include "LinsleyAtmosphereData.hh"
#include <TH1F.h>
#include <TROOT.h>
#include <math.h>
#include <TSpline.h>

Atmosphere* LinsleyAtmosphere::fMe = NULL;

using namespace sou;

ClassImp(LinsleyAtmosphere)

//___________________________________________________________________________________
LinsleyAtmosphere::LinsleyAtmosphere(string type) : Atmosphere(), fData(0) {
    //
    // ctor
    //
    fType = type;
    Msg(EsafMsg::Info) << "*** LinsleyAtmosphere ***"<< MsgDispatch;
    Build();
}

//___________________________________________________________________________________
LinsleyAtmosphere::~LinsleyAtmosphere() {
    //
    // dtor
    //
    SafeDelete(fData);
}

//___________________________________________________________________________________
void LinsleyAtmosphere::CreateInstance() {
    //
    // Create the single atmosphere instance as being a LinsleyAtmosphere object
    //
    if(!Atmosphere::fChild) Atmosphere::fChild = new LinsleyAtmosphere("linsley");
    else Atmosphere::fChild->Msg(EsafMsg::Warning) << "Trying to create two atmospheres !!"<< MsgDispatch;
}

//___________________________________________________________________________________
void LinsleyAtmosphere::Build() {
    //
    // Build AtmosphereData object
    //
    if(!fData) fData = new LinsleyAtmosphereData();
    if(!fData) Msg(EsafMsg::Panic) << "Pb of memory allocation in LinsleyAtmosphere" << MsgDispatch ; 
    
    if(fData->GetName() == "Uniform") fTOA = fData->Bottom(1); // tow layers only
    
    // interpolators initialization
    /*
    const Double_t* data = 0;
    data =  fData->GetTemperatureTable();
    TArrayD t_arr(100,data);
    fTempInterpol = new TSpline3("fTempInterpol",0,100*km,t_arr.GetArray(),100);    
    
    data =  fData->GetAir_DensityTable();
    TArrayD rho_arr(100,data);
    fDensInterpol = new TSpline3("fDensInterpol",0,100*km,rho_arr.GetArray(),100);    
    
    data =  fData->GetPressureTable();
    TArrayD press_arr(100,data);
    fPressInterpol = new TSpline3("fPressInterpol",0,100*km,press_arr.GetArray(),100);    
    */
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::Pressure(Double_t h) const {
    //
    // Pressure at a given altitude
    //
    if(h > fTOA)  return 0.;
    return 0;//fPressInterpol->Eval(h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::Temperature(Double_t h) const {
    //
    // Temperature at a given altitude
    //
    if(h > fTOA)  return 0.;
    return 0;//fTempInterpol->Eval(h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::AbsoluteHumidity(Double_t h) const {
    //
    // Absolute humidity at a given altitude (water vapor density)
    //
    return 0;
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::Air_Density(Double_t h) const {
    //
    // Air density at a given altitude
    //
    if(h > fTOA)  return 0.;
        
    return fData->DensityParametrization(h);    
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::O_Density(Double_t h) const {
    //
    // Oxigen density at a given altitude
    //
     if(h > fTOA )
        return 0.;

    // linear interpolation between table values
    string val = "Odens";
    string mod = "expon";
    return 0;//Interpolate(val,mod,h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::O2_Density(Double_t h) const {
    //
    // Di-oxigen density at a given altitude
    //
   if(h > fTOA)
        return 0.;

   // linear interpolation between table values
   string val = "O2dens";
   string mod = "expon";
   return 0;//Interpolate(val,mod,h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::O3_Density(Double_t h) const {
    //
    // Ozone density at a given altitude
    //
    if(h > fTOA)
        return 0.;

    // linear interpolation between table values
    string val = "O3dens";
    string mod = "linear";
    return 0;//Interpolate(val,mod,h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::N2_Density(Double_t h) const {
    //
    // Di-nitrogen density at a given altitude
    //
    if(h > fTOA)
        return 0.;

    // linear interpolation between table values
    string val = "N2dens";
    string mod = "expon";
    return 0;//Interpolate(val,mod,h);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::Aerosols_Density(string& type,Double_t h) const {
    //
    // Aerosols density for a given type at a given altitude
    //
    return 0;
}

//___________________________________________________________________________________
long double LinsleyAtmosphere::Index(Double_t h, Double_t wl) const {
    //
    // Index as fonction of altitude ONLY (formula drawn from Corsika)
    //
    return 1. + Index_Minus1(h,wl);
}

//___________________________________________________________________________________
Double_t LinsleyAtmosphere::Index_Minus1(Double_t h, Double_t wl) const {
    //
    // returns (n-1)
    //
    return 0.000283*Air_Density(h)/Air_Density(0.);
}

//___________________________________________________________________________________
double LinsleyAtmosphere::Interpolate(string& val,string& mod,double h) const {
    //
    // Calculate interpolation (from tables of parameter)
    //
    const double* data = 0;
    if(val == "Odens")        data = fData->GetO_DensityTable();
    else if(val == "O2dens")  data = fData->GetO2_DensityTable();
    else if(val == "N2dens")  data = fData->GetN2_DensityTable();
    else  Msg(EsafMsg::Panic) << "Wrong val argument in LinsleyAtmosphere::Interpolate"  << MsgDispatch ;
    
    if(h<0) return data[0];
    int nb=100;
    double step=1*km;
    if(h > (nb-1)*step) return data[nb-1];
 
    double val1, val2;
    double h1;
    int i = 0;
    i = Int_t(floor(h/step));
    h1 = i * step;
    val1 = data[i];
    val2 = data[i+1];
        
    if(mod == "linear" || val1 == 0 || val2 == 0)  return (val2 - val1)*(h - h1)/step + val1;
    else if(mod == "expon")                        return val1 * exp(-(h - h1) * log(val1/val2)/step);
    else   Msg(EsafMsg::Panic) << "Wrong mod argument in LinsleyAtmosphere::Interpolate"  << MsgDispatch ;
    return 0;
}


//_____________________________________________________________________________
Double_t LinsleyAtmosphere::GetLatitude() const {
    //
    // used to identify the model
    //
    
    return fData->GetLatitude();
}

//_____________________________________________________________________________
Double_t LinsleyAtmosphere::GetLongitude() const {
    //
    // dummy
    //
    
    return fData->GetLongitude();
}

//_____________________________________________________________________________
Double_t LinsleyAtmosphere::GetDate() const {
    //
    // dummy
    //
    
    return fData->GetDate();
}