// ESAF : Euso Simulation and Analysis Framework
// $Id: MSISE_00AtmosphereData.cc,v 1.13 2005/05/02 15:06:59 moreggia Exp $
// Sylvain Moreggia created Dec, 12 2003
#include "MSISE_00AtmosphereData.hh"
#include "Config.hh"
#include <TH1F.h>
#include <TROOT.h>
#include "EConst.hh"
ClassImp(MSISE_00AtmosphereData)
using EConst::Avogadro;
using EConst::AtomicMass;
using EConst::R_ideal;
//________________________________________________________________________________________
MSISE_00AtmosphereData::MSISE_00AtmosphereData() : AtmosphereData(), fTool(0) {
//
// ctor
//
fNb = 100;
fStep = 100*km / fNb;
SetTables();
}
//________________________________________________________________________________________
MSISE_00AtmosphereData::~MSISE_00AtmosphereData() {
//
// dtor
//
if(fTool) delete fTool;
fTool = 0;
}
//________________________________________________________________________________________
void MSISE_00AtmosphereData::Init() {
//
// Initialize msise parameters to work out tables claculations
//
fTool = new MSISEtool();
if(!fTool) Msg(EsafMsg::Panic) << "No new MSISEtool, memory pb" <<
MsgDispatch ;
ConfigFileParser* pConf =
Config::Get()->GetCF("Atmosphere","MSISE_00Atmosphere");
fInput.year = (Int_t)pConf->GetNum("MSISE_00Atmosphere.year");
fInput.doy = (Int_t)pConf->GetNum("MSISE_00Atmosphere.doy");
fInput.sec = pConf->GetNum("MSISE_00Atmosphere.sec");
fInput.g_lat = pConf->GetNum("MSISE_00Atmosphere.g_lat");
fInput.g_long = pConf->GetNum("MSISE_00Atmosphere.g_long");
fInput.lst = fInput.sec/3600 + fInput.g_long/15;
if(fInput.lst > 24.) fInput.lst -= 24.;
if(fInput.lst < 0.) fInput.lst += 24.;
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"LAT. = " <<fInput.g_lat <<MsgDispatch;
Msg(EsafMsg::Debug) <<"LONG. = " <<fInput.g_long <<MsgDispatch;
Msg(EsafMsg::Debug) <<"Local Solar Time = " <<fInput.lst <<" h"<<MsgDispatch;
#endif
fInput.f107A = pConf->GetNum("MSISE_00Atmosphere.f107A");
fInput.f107 = pConf->GetNum("MSISE_00Atmosphere.f107");
fInput.ap = pConf->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..
}
//________________________________________________________________________________________
void MSISE_00AtmosphereData::SetTables() {
//
// Set tables using MSISEtool
//
Init();
fAltitudeTable = new Double_t[fNb];
fAir_DensityTable = new Double_t[fNb];
fO_DensityTable = new Double_t[fNb];
fO2_DensityTable = new Double_t[fNb];
fN2_DensityTable = new Double_t[fNb];
fTemperatureTable = new Double_t[fNb];
fPressureTable = new Double_t[fNb];
struct nrlmsise_output _out;
for(Int_t i=0; i<fNb; i++) {
fInput.alt = i * fStep/km;
fTool->gtd7(&fInput,&fTuning,&_out);
fO_DensityTable[i] = _out.d[1] * (AtomicMass("Oxygen")*g/(mole*Avogadro())) / cm3;
fN2_DensityTable[i] = _out.d[2] * (AtomicMass("DiNitrogen")*g/(mole*Avogadro())) / cm3;
fO2_DensityTable[i] = _out.d[3] * (AtomicMass("DiOxygen")*g/(mole*Avogadro())) / cm3;
fAir_DensityTable[i] = _out.d[5] * g/cm3;
fTemperatureTable[i] = _out.t[1] * kelvin;
fPressureTable[i] = fAir_DensityTable[i] / (AtomicMass("Air")*g/mole) *
(R_ideal()/kelvin/mole) * fTemperatureTable[i];
fAltitudeTable[i] = i * fStep;
}
string atmod = "msise";
GetDefault(atmod,fNb);
// int fNbL=33;
WriteUserModelFile(fNb);
#ifdef DEBUG
DebugPlots(fNb);
#endif
}