// $Id: LinsleyAtmosphere.cc,v 1.22 2005/02/10 17:59:22 moreggia 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;
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;
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;
}