// ESAF : Euso Simulation and Analysis Framework
// $Id: LowtranAtmosphereData.cc,v 1.11 2005/05/02 16:51:54 moreggia Exp $
// Corinne Berat created Mar, 30 2004
//
// Handle LOWTRAN atmospheric parameter tables
#include "LowtranAtmosphereData.hh"
#include "Config.hh"
#include "NumbersFileParser.hh"
#include "ConfigFileParser.hh"
#include <TProfile.h>
#include <TROOT.h>
#include "EConst.hh"
ClassImp(LowtranAtmosphereData)
using EConst::Avogadro;
using EConst::R_ideal;
using EConst::Loschmidt;
using EConst::AtomicMass;
//_____________________________________________________________________________
LowtranAtmosphereData::LowtranAtmosphereData() : AtmosphereData() {
//
// ctor
//
fNb = 50;
SetTables();
}
//_____________________________________________________________________________
LowtranAtmosphereData::~LowtranAtmosphereData() {
//
// dtor
//
}
//_____________________________________________________________________________
void LowtranAtmosphereData::Init() {
//
// Initialize lowtran parameters to work out tables claculations
//
ConfigFileParser* pConf = Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
fInput.model = (Int_t)pConf->GetNum("LowtranAtmosphere.model");
Msg(EsafMsg::Info)<< "LOWTRAN model number " << fInput.model << " named as " << fModelName << MsgDispatch;
switch(fInput.model) {
case 1: fModelName = "Tropical"; break;
case 2: fModelName = "MidLatSummer"; break;
case 3: fModelName = "MidLatWInt_ter"; break;
case 6: fModelName = "USStandard"; break;
default: Msg(EsafMsg::Warning)<< "incorrect LOWTRAN model name " << MsgDispatch;
}
}
//_____________________________________________________________________________
void LowtranAtmosphereData::SetTables() {
//
// Set tables using data file
//
Init();
fAltitudeTable = new Double_t[fNb];
fAir_DensityTable = new Double_t[fNb];
fO2_DensityTable = new Double_t[fNb];
fCO2_DensityTable = new Double_t[fNb];
fO3_DensityTable = new Double_t[fNb];
fN2_DensityTable = new Double_t[fNb];
fH2O_DensityTable = new Double_t[fNb];
fTemperatureTable = new Double_t[fNb];
fPressureTable = new Double_t[fNb];
NumbersFileParser nf("config/Atmosphere/LowtranData/" + fModelName,7);
vector<Double_t> alt = nf.GetCol(0);
vector<Double_t> temp = nf.GetCol(1);
vector<Double_t> pres = nf.GetCol(2);
vector<Double_t> h2o = nf.GetCol(3);
vector<Double_t> ozone= nf.GetCol(4);
vector<Double_t> o2 = nf.GetCol(5);
vector<Double_t> co2 = nf.GetCol(6);
const Double_t mAir=AtomicMass("Air");
const Double_t mO2=AtomicMass("DiOxygen");
const Double_t mO3=AtomicMass("TriOxygen");
const Double_t mN2=AtomicMass("DiNitrogen");
const Double_t mH2O=AtomicMass("H2O");
const Double_t mCO2=AtomicMass("CO2");
Double_t mAir_over_R=mAir/(R_ideal()/joule);
const Double_t p0=STP_Pressure/(bar*0.001),t0=STP_Temperature;
Double_t pss,tss,convert;
#ifdef DEBUG
//ofstream Debug("output/lowtran.out");
//Debug << mAir << " " << R_ideal()/joule << " " << mAir_over_R << " " << endl;
//Debug << Loschmidt() << " " << Loschmidt()*cm3 << " " << endl;
#endif
string gaslaw = "lowtran"; // can be "ideal" for debugging
for(Int_t i=0; i<fNb; i++) {
fAltitudeTable[i] = alt[i]*km;
fTemperatureTable[i] = temp[i];
fPressureTable[i] = pres[i]*bar*0.001;
fAir_DensityTable[i] = mAir_over_R*1.e-4*(pres[i]/temp[i])*g/cm3;
pss=pres[i]/p0;
tss=t0/temp[i];
convert = 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
fO2_DensityTable[i] = o2[i]*mO2*convert;
fO3_DensityTable[i] = ozone[i]*mO3*convert;
fN2_DensityTable[i] = 0.781*1.e6*mN2*convert;
fCO2_DensityTable[i] = co2[i]*mCO2*convert;
fH2O_DensityTable[i] = h2o[i]*mH2O*convert;
if (gaslaw=="lowtran"){
fO2_DensityTable[i] = fO2_DensityTable[i]*pow(pss,1.1879)*pow(tss,2.9738);
fO3_DensityTable[i] = fO3_DensityTable[i]*pow(pss,0.420)*pow(tss,1.3903);
fN2_DensityTable[i] = 0.781*1.e6*mN2*convert *pow(pss,0.5)*pow(tss,0.5);
fCO2_DensityTable[i] = fCO2_DensityTable[i]*pow(pss,0.6705)*pow(tss,-2.2560);
fH2O_DensityTable[i] = fH2O_DensityTable[i]*pow(pss,0.981)*pow(tss,0.3324);
}
#ifdef DEBUG
/*
Debug << " *** " <<endl;
Debug << " lowtran " << alt[i] <<" "
<< temp[i] <<" "
<< pres[i] <<endl;
Debug << " " << " N2 O2 ozone co2 h2o\n";
Debug << " " <<0.781*(pss*pss)*tss*sqrt(tss)<<" "
<<o2[i]*0.1*pow(pss,2.1879)*pow(tss,3.9738)<< " "
<<ozone[i]*0.1*pow(pss,1.420)*pow(tss,2.3908)<< " "
<<co2[i]*0.1*pow(pss,1.6705)*pow(tss,-1.256)<< " "
<<h2o[i]*0.1*pow(pss,1.981)*pow(tss,1.3324)
*mH2O*Loschmidt()*cm3/Avogadro() <<endl;
Debug << " converted in ESAF units " <<endl;
Debug << " " <<0.781*(pss*pss)*tss*sqrt(tss)
* (mN2*g) /(cm2*km) << " "
<<o2[i]*1.e-6*pow(pss,2.1879)*pow(tss,3.9738)
* (mO2*g) /(cm2*km) << " "
<<ozone[i]*1.e-6*pow(pss,1.420)*pow(tss,2.3908)
* (mO3*g) /(cm2*km) <<" "
<<co2[i]*1.e-6*pow(pss,1.6705)*pow(tss,-1.256)
* (mCO2*g) /(cm2*km) << " "
<<h2o[i]*1.e-6*pow(pss,1.981)*pow(tss,1.3324)
* (mH2O*g) /(cm2*km) <<endl;
Debug << " computed and stored in tables " <<endl;
Debug << " "
<< fTemperatureTable[i] << " "
<< fAir_DensityTable[i] << " "
<< fPressureTable[i] << endl;
Debug << " "
<< fN2_DensityTable[i] << " "
<< fO2_DensityTable[i] << " "
<< fO3_DensityTable[i] << " "
<< fCO2_DensityTable[i] << " "
<< fH2O_DensityTable[i] << endl;
*/
#endif
}
string atmod = "lowtran";
GetDefault(atmod,fNb);
WriteUserModelFile(fNb);
#ifdef DEBUG
//DebugPlots(fNb);
#endif
}