Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

LinsleyAtmosphere - source file

// $Id: LinsleyAtmosphere.cc,v 1.24 2005/10/02 14:17:07 thea Exp $
// 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>

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 ; 
}

//___________________________________________________________________________________
 Double_t LinsleyAtmosphere::Pressure(Double_t h) const {
    //
    // Pressure at a given altitude
    //
    if(h > 100*km)
        return 0.;

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

//___________________________________________________________________________________
 Double_t LinsleyAtmosphere::Temperature(Double_t h) const {
    //
    // Temperature at a given altitude
    //
    if(h > 100*km)
        return 0.;

    // linear interpolation between table values
    string val = "temp";
    string mod = "linear";
    return Interpolate(val,mod,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 > 100*km)
        return 0.;
	
    return fData->DensityParametrization(h);    
}

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

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

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

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

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

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

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

    // linear interpolation between table values
    string val = "N2dens";
    string mod = "linear";
    return 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 == "press")
         data = fData->GetPressureTable();
    else if(val == "temp")
        data = fData->GetTemperatureTable();
    else if(val == "Airdens")
        data = fData->GetAir_DensityTable();
    else 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=33;
    double step=1*km;
    if(h > nb*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")
        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;
}




About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005