// ESAF : Euso Simulation and Analysis Framework
// $Id: MSISE_00Atmosphere.cc,v 1.25 2005/05/02 15:06:59 moreggia Exp $
// S. Moreggia created 18 November 2003
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: MSISE_00Atmosphere *
* Package: atmosphere *
* Coordinator: S. Moreggia *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// MSISE_00Atmosphere
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include <math.h>
#include "MSISE_00Atmosphere.hh"
#include "MSISEtool.hh"
#include "Config.hh"
#include "MSISE_00AtmosphereData.hh"
#include "EConst.hh"
ClassImp(MSISE_00Atmosphere)
using EConst::Avogadro;
using EConst::R_ideal;
using EConst::AtomicMass;
//____________________________________________________________________________
MSISE_00Atmosphere::MSISE_00Atmosphere(string type) : Atmosphere(), fData(0), fTool(0) {
//
// ctor
//
fType = type;
Msg(EsafMsg::Info) << "*** MSISE_00Atmosphere ***"<< MsgDispatch;
Build();
}
//____________________________________________________________________________
MSISE_00Atmosphere::~MSISE_00Atmosphere() {
//
// dtor
//
SafeDelete(fTool);
SafeDelete(fData);
}
//___________________________________________________________________________________
void MSISE_00Atmosphere::CreateInstance() {
//
// Create the single atmosphere instance as being a MSISE_00Atmosphere object
//
string type = "msise00";
if(!Atmosphere::fChild) Atmosphere::fChild = new MSISE_00Atmosphere(type);
else Atmosphere::fChild->Msg(EsafMsg::Warning) << "Trying to create two atmospheres !!"<< MsgDispatch;
}
//____________________________________________________________________________
void MSISE_00Atmosphere::Build() {
//
// Build atmosphere according to config parameters
// Arrays data always built. No arrays kept for historical reasons
//
string f = "yes";
if(f == "yes")
fFlag = true;
else if(f == "no")
fFlag = false;
else
Msg(EsafMsg::Panic) << "Wrong configuration of MSISE_00Atmosphere.fFlag" << MsgDispatch ;
if(fFlag) {
fData = new MSISE_00AtmosphereData();
if(!fData) Msg(EsafMsg::Panic)<<"No new MSISEData, memory pb"<<MsgDispatch;
}
else {
fTool = new MSISEtool();
if(!fTool) Msg(EsafMsg::Panic)<< "No new MSISEtool, memory pb"<<MsgDispatch;
fInput.year = (Int_t)Conf()->GetNum("MSISE_00Atmosphere.year");
fInput.doy = (Int_t)Conf()->GetNum("MSISE_00Atmosphere.doy");
fInput.sec = Conf()->GetNum("MSISE_00Atmosphere.sec");
fInput.g_lat = Conf()->GetNum("MSISE_00Atmosphere.g_lat");
fInput.g_long = Conf()->GetNum("MSISE_00Atmosphere.g_long");
#ifdef DEBUG
Msg(EsafMsg::Info)<<"LAT. = " <<fInput.g_lat <<MsgDispatch ;
Msg(EsafMsg::Info)<<"LONG. = " <<fInput.g_long << MsgDispatch;
#endif
fInput.lst = fInput.sec/3600 + fInput.g_long/15;
if(fInput.lst > 24.) fInput.lst -= 24.;
if(fInput.lst < 0.) fInput.lst += 24.;
fInput.f107A = Conf()->GetNum("MSISE_00Atmosphere.f107A");
fInput.f107 = Conf()->GetNum("MSISE_00Atmosphere.f107");
fInput.ap = Conf()->GetNum("MSISE_00Atmosphere.ap");
fTuning.switches[0] = 0;
for(Int_t i=1; i<24; i++)
fTuning.switches[i] = 1; //TOFIX : need to allow better tuning..
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::Interpolate(string& val,string& mod,Double_t h) const {
//
// Calculate interpolation
//
const Double_t* 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 MSISE_00Atmosphere::Interpolate" << MsgDispatch;
if(h < 0) return data[0];
Int_t nb = fData->NumberOfElements();
Double_t step = fData->LayerSize();
if(h > step*nb) return data[nb-1];
Double_t val1, val2;
Double_t h1;
Int_t 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 MSISE_00Atmosphere::Interpolate"
<< MsgDispatch ;
return 0;
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::Pressure(Double_t h) const {
//
// Pressure at a given altitude
//
if(h > 100*km)
return 0.;
if(fFlag) {
// linear interpolation between table values
string val = "press";
string mod = "expon";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.d[5]/(AtomicMass("Air")*g/mole) * (R_ideal()/kelvin/mole) * _out.t[1]*kelvin;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::Temperature(Double_t h) const {
//
// Temperature at a given altitude
//
if(h > 100*km)
return 0.;
if(fFlag) {
// linear interpolation between table values
string val = "temp";
string mod = "linear";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.t[1]*kelvin;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::AbsoluteHumidity(Double_t h) const {
//
// Absolute humidity at a given altitude (water vapor density)
//
return 0;
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::Air_Density(Double_t h) const {
//
// Atmospheric air density at a given altitude
//
if(h >= 100*km)
return 0.;
if(fFlag) {
// exponential interpolation between table values
string val = "Airdens";
string mod = "expon";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.d[5] * g/cm3;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::O_Density(Double_t h) const {
//
// Oxygen density at a given altitude
//
if(h > 100*km)
return 0.;
if(fFlag) {
// exponential interpolation between table values
string val = "Odens";
string mod = "linear";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.d[1] * (AtomicMass("Oxygen")*g/(mole*Avogadro())) / cm3;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::O2_Density(Double_t h) const {
//
// Di-oxygen density at a given altitude
//
if(h > 100*km)
return 0.;
if(fFlag) {
// exponential interpolation between table values
string val = "O2dens";
string mod = "linear";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.d[3] * (AtomicMass("DiOxygen")*g/(mole*Avogadro())) / cm3;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::O3_Density(Double_t h) const {
//
// Ozone density at a given altitude
//
Msg(EsafMsg::Panic)<<"O3_Density not available for MSISE_00Atmosphere"
<< MsgDispatch ;
return 0;
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::N2_Density(Double_t h) const {
//
// Di-nitrogen density at a given altitude
//
if(h > 100*km)
return 0.;
if(fFlag) {
// exponential interpolation between table values
string val = "N2dens";
string mod = "linear";
return Interpolate(val,mod,h);
}
else {
// msise calculation
struct nrlmsise_output _out;
fInput.alt = h/km;
fTool->gtd7(&fInput,&fTuning,&_out);
return _out.d[2] * (AtomicMass("DiNitrogen")*g/(mole*Avogadro())) / cm3;
}
}
//____________________________________________________________________________
Double_t MSISE_00Atmosphere::Aerosols_Density(string& type,Double_t h) const {
//
// Aerosols density for a given type of aerosols at a given altitude
//
Msg(EsafMsg::Panic)<<"Aerosol_Density not available for MSISE_00Atmosphere"
<< MsgDispatch ;
return 0;
}