// ESAF : Euso Simulation and Analysis Framework
// $Id: LinsleyAtmosphereData.cc,v 1.9 2004/12/05 02:13:14 thea Exp $
// Sylvain Moreggia created May, 17 2004
#include "LinsleyAtmosphereData.hh"
#include "NumbersFileParser.hh"
#include <math.h>
ClassImp(LinsleyAtmosphereData)
//_____________________________________________________________________________
LinsleyAtmosphereData::LinsleyAtmosphereData() : AtmosphereData() {
//
// ctor
//
fNb = (Int_t)Conf()->GetNum("LinsleyAtmosphereData.fNb");
fName = Conf()->GetStr("LinsleyAtmosphereData.fName");
fBottom = new Double_t[fNb];
fA = new Double_t[fNb];
fB = new Double_t[fNb];
fC = new Double_t[fNb];
if(!(fBottom && fA && fB && fC)) Msg(EsafMsg::Panic) <<"Pb for memory allocation in LinsleyAtmosphereData ctor"<<MsgDispatch;
ParamTables();
SetTables();
}
//_____________________________________________________________________________
LinsleyAtmosphereData::LinsleyAtmosphereData(const Int_t&, const vector<Double_t>&, const Int_t&)
: AtmosphereData() {
//
// ctor
//
fBottom = 0;
fA = 0;
fB = 0;
fC = 0;
}
//_____________________________________________________________________________
LinsleyAtmosphereData::~LinsleyAtmosphereData() {
//
// dtor
//
}
//_____________________________________________________________________________
void LinsleyAtmosphereData::ParamTables() {
//
// fill data members (tables)
//
string fn = fName + ".data";
NumbersFileParser nf("config/Atmosphere/LinsleyData/" + fn, fNb);
Msg(EsafMsg::Info)<< " Linsley Atmosphere - config :" << fName << MsgDispatch;
vector<Double_t> param;
for(Int_t i=0; i<fNb; i++) {
param = nf.GetCol(i);
fBottom[i] = param[0]*km;
fA[i] = param[2]*g/cm2;
fB[i] = param[3]*g/cm2;
fC[i] = param[4]*cm;
}
}
//_____________________________________________________________________________
void LinsleyAtmosphereData::SetTables() {
//
// fill data members (tables)
//
Int_t fNbL=100;
fAltitudeTable = new Double_t[fNbL];
fAir_DensityTable = new Double_t[fNbL];
for(Int_t j=0; j<fNbL; j++){
Double_t h=double(j);
fAltitudeTable[j] = h*km;
for(Int_t i=0; i<fNb; i++) {
if(fBottom[i]/km <= h && h < fBottom[i+1]/km) {
fAir_DensityTable[j]= fB[i]/fC[i] * exp(-h*km/fC[i]); }
}
}
string atmod = "linsley";
GetDefault(atmod,fNbL);
WriteUserModelFile(fNbL);
#ifdef DEBUG
DebugPlots(fNbL);
#endif
}
//_____________________________________________________________________________
Double_t LinsleyAtmosphereData::DensityParametrization(Double_t h) const {
//
// Linsley parametrization
//
Int_t i;
if(h < 0) Msg(EsafMsg::Panic) << "Negative altitude not authorized in LinsleyAtmosphereData" << MsgDispatch;
for(i=0; i<fNb-1; i++) {
if(h < fBottom[i+1]) return fB[i]/fC[i] * exp(-h/fC[i]);
}
return fB[i]/fC[i];
}