// ESAF : Euso Simulation and Analysis Framework
// $Id: AtmosphereData.cc,v 1.22 2005/05/13 08:35:13 moreggia Exp $
// Sylvain Moreggia created Dec, 1 2003
#include "AtmosphereData.hh"
#include "Atmosphere.hh"
#include "LinsleyAtmosphere.hh"
#include "NumbersFileParser.hh"
#include "ConfigFileParser.hh"
#include "Config.hh"
#include <TH1F.h>
#include <TROOT.h>
#include "EConst.hh"
ClassImp(AtmosphereData)
using EConst::Avogadro;
using EConst::Loschmidt;
using EConst::R_ideal;
using EConst::AtomicMass;
extern void InitLOWTRAN();
//______________________________________________________________________________
AtmosphereData::AtmosphereData() : EsafConfigurable(), EsafMsgSource(){
//
// ctor
//
fAltitudeTable = 0;
fPressureTable = 0;
fTemperatureTable = 0;
fAir_DensityTable = 0;
fO_DensityTable = 0;
fO2_DensityTable = 0;
fO3_DensityTable = 0;
fN2_DensityTable = 0;
fCO2_DensityTable = 0;
fH2O_DensityTable = 0;
fIndexTable = 0;
}
//______________________________________________________________________________
AtmosphereData::~AtmosphereData() {
//
// dtor
//
if(fAltitudeTable) delete[] fAltitudeTable;
fAltitudeTable = 0;
if(fPressureTable) delete[] fPressureTable;
fPressureTable = 0;
if(fTemperatureTable) delete[] fTemperatureTable;
fTemperatureTable = 0;
if(fAir_DensityTable) delete[] fAir_DensityTable;
fAir_DensityTable = 0;
if(fO_DensityTable) delete[] fO_DensityTable;
fO_DensityTable = 0;
if(fO2_DensityTable) delete[] fO2_DensityTable;
fO2_DensityTable = 0;
if(fO3_DensityTable) delete[] fO3_DensityTable;
fO3_DensityTable = 0;
if(fN2_DensityTable) delete[] fN2_DensityTable;
fN2_DensityTable = 0;
if(fCO2_DensityTable) delete[] fCO2_DensityTable;
fCO2_DensityTable = 0;
if(fH2O_DensityTable) delete[] fH2O_DensityTable;
fH2O_DensityTable = 0;
if(fIndexTable) delete[] fIndexTable;
fIndexTable = 0;
}
//______________________________________________________________________________
void AtmosphereData::GetDefault(string& atmod, Int_t fNbL) {
//
// fill missing data with default values
//
if (!fAltitudeTable){
fAltitudeTable = new Double_t[fNbL];}
if (!fPressureTable){
fPressureTable = new Double_t[fNbL];}
if (!fTemperatureTable){
fTemperatureTable = new Double_t[fNbL];}
if (!fAir_DensityTable){
fAir_DensityTable = new Double_t[fNbL];}
if (!fO_DensityTable){
fO_DensityTable = new Double_t[fNbL];}
if (!fO2_DensityTable){
fO2_DensityTable = new Double_t[fNbL];}
if (!fO3_DensityTable){
fO3_DensityTable = new Double_t[fNbL];}
if (!fN2_DensityTable){
fN2_DensityTable = new Double_t[fNbL];}
if (!fCO2_DensityTable){
fCO2_DensityTable = new Double_t[fNbL];}
if (!fH2O_DensityTable){
fH2O_DensityTable = new Double_t[fNbL];}
if (!fIndexTable){
fIndexTable = new Double_t[fNbL];}
string gaslaw = "lowtran"; // can be ideal for debugging
/* If Linsley : complete data with USStandard and MSISE */
/* profiles */
/* If Msise : complete data with lowtran models, */
/* according to latitude, longitude */
/* If Lowtran : complete data with Msise, according to model*/
const Double_t R=R_ideal()/joule;
const Double_t mAir=AtomicMass("Air");
const Double_t mO=AtomicMass("Oxygen");
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");
if(atmod =="linsley"){
Msg(EsafMsg::Info) << "The present model is Linsley" << MsgDispatch;
Msg(EsafMsg::Info) << "Missing profiles are provided by default files" <<MsgDispatch;
Msg(EsafMsg::Info) << "US Standard and MSISE lat=45, long=270, day=1"<<MsgDispatch;
NumbersFileParser nf("config/Atmosphere/DefaultData/Default_6",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> co2 = nf.GetCol(6);
NumbersFileParser ng("config/Atmosphere/DefaultData/lat_45_long_270_day_1",4);
// vector<Double_t> alt = nf.GetCol(0);
vector<Double_t> n2 = ng.GetCol(1);
vector<Double_t> o2 = ng.GetCol(2);
vector<Double_t> o = ng.GetCol(3);
Double_t pss,tss,convert;
for(Int_t i=0; i<fNbL; i++) {
fTemperatureTable[i] = temp[i];
fPressureTable[i] = fAir_DensityTable[i] / (mAir*g/mole) *
(R * joule/kelvin/mole) * fTemperatureTable[i];
pss = fPressureTable[i]/STP_Pressure;
tss = STP_Temperature/fTemperatureTable[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] = n2[i]*mN2*convert;
fCO2_DensityTable[i] = co2[i]*mCO2*convert;
fH2O_DensityTable[i] = h2o[i]*mH2O*convert;
fO_DensityTable[i] = o[i]*convert*mO;
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);}
}
}
if(atmod =="msise"){
ConfigFileParser* pConf =
Config::Get()->GetCF("Atmosphere","MSISE_00Atmosphere");
Double_t day = (Int_t)pConf->GetNum("MSISE_00Atmosphere.doy");
Double_t g_lat = pConf->GetNum("MSISE_00Atmosphere.g_lat");
Double_t g_long = pConf->GetNum("MSISE_00Atmosphere.g_long");
string fmodel;
if (-25.<= g_lat && g_lat <= 25.){
fmodel ="Default_1"; //Tropical
}
else{
if (80 <=day && day<= 266) {
fmodel="Default_2";} //MidLatSummer
else {
fmodel="Default_3";}} //MidLatWInt_ter
NumbersFileParser nf("config/Atmosphere/DefaultData/" + fmodel,7);
vector<Double_t> alt = nf.GetCol(0);
vector<Double_t> h2o = nf.GetCol(3);
vector<Double_t> ozone= nf.GetCol(4);
vector<Double_t> co2 = nf.GetCol(6);
Double_t pss,tss,convert;
Msg(EsafMsg::Info) << "The present model is Msise : lat = "<<g_lat
<< " long = "<<g_long << " day = " <<day <<MsgDispatch;
Msg(EsafMsg::Info) << "Missing profiles are provided by files" <<MsgDispatch;
Msg(EsafMsg::Info) << "from the corresponding Lowtran Model : "<< fmodel << MsgDispatch;
for(Int_t i=0; i<fNbL; i++) {
pss=fPressureTable[i]/STP_Pressure;
tss=STP_Temperature/fTemperatureTable[i];
convert=1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
fO3_DensityTable[i] = ozone[i]*mO3*convert;
fCO2_DensityTable[i] = co2[i]*mCO2*convert;
fH2O_DensityTable[i] = h2o[i]*mH2O*convert;
if (gaslaw=="lowtran"){
fO3_DensityTable[i] = fO3_DensityTable[i]*pow(pss,0.420)*pow(tss,1.3903);
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);}
}
}
if(atmod =="lowtran"){
Msg(EsafMsg::Info) << "The present model is Lowtran" << MsgDispatch ;
ConfigFileParser* pConf =
Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
Int_t Model = (Int_t)pConf->GetNum("LowtranAtmosphere.model");
string fmodel;
if (Model==1){
// Tropical model
fmodel ="lat_0_long_0_day_1";
}
else{
if (Model==2) {
//"MidLatSummer"
fmodel="lat_45_long_270_day_200";}
else {
//MidLatWInt_ter or US Standard
fmodel="lat_45_long_270_day_1";}
}
Msg(EsafMsg::Info) << "Missing profiles are provided by files" <<MsgDispatch;
Msg(EsafMsg::Info) << "from the corresponding MSISE Model : "<< fmodel << MsgDispatch;
NumbersFileParser nf("config/Atmosphere/DefaultData/" + fmodel,4);
vector<Double_t> alt = nf.GetCol(0);
vector<Double_t> n2 = nf.GetCol(1);
vector<Double_t> o2 = nf.GetCol(2);
vector<Double_t> o = nf.GetCol(3);
Double_t pss,tss,convert;
for(Int_t i=0; i<fNbL; i++) {
pss=fPressureTable[i]/STP_Pressure;
tss=STP_Temperature/fTemperatureTable[i];
convert= 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
for(Int_t j=0; j<101; j++) {
if (alt[j]==fAltitudeTable[i]/km ||
(Int_t)alt[j]==(Int_t)(fAltitudeTable[i]/km)){
fO2_DensityTable[i] = o2[j]*convert*mO2;
fO_DensityTable[i] = o[j]*convert*mO;}}
}
}
#ifdef DEBUG
Msg(EsafMsg::Debug)
<< "densities computation using " << gaslaw << " gas law "<< MsgDispatch;
#endif
return;
}
//______________________________________________________________________________
void AtmosphereData::WriteUserModelFile(Int_t fNbL) {
/* prepare and write a file necessary to use lowtran radiative
process calculator. It correspond to tape5 */
FILE* tape5=fopen("config/RadiativeTransfer/Lwtrn_Tape5.dat","w");
ConfigFileParser* pConf =
Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
/* prepare CARD1 */
/* the model has been defined by the user */
const Int_t Model=7,M1=0,M2=0,M3=0,M4=6,M5=6,M6=6;
const Int_t Mdef =1; //
const Int_t Itype =2; // vertical or slant path between two altitudes
const Int_t Iemsct=0; // program execution in transmittance mode
const Int_t Imult =0; // no multiple scattering in transmittance mode
const Double_t Tbound=0; // FOR NORMAL OPERATION OF PROGRAM.
Int_t Noprt = (Int_t)pConf->GetNum("LowtranAtmosphere.Noprt");
Int_t Im = (Int_t)pConf->GetNum("LowtranAtmosphere.Im");
Double_t Salb = pConf->GetNum("LowtranAtmosphere.Salb");
Int_t Ihaze = (Int_t)pConf->GetNum("LowtranAtmosphere.Ihaze");
Int_t Iseasn = (Int_t)pConf->GetNum("LowtranAtmosphere.Iseasn");
Int_t Ivulcn = (Int_t)pConf->GetNum("LowtranAtmosphere.Ivulcn");
Int_t Icstl = (Int_t)pConf->GetNum("LowtranAtmosphere.Icstl");
Int_t Icld = (Int_t)pConf->GetNum("LowtranAtmosphere.Icld");
Int_t Ivsa = (Int_t)pConf->GetNum("LowtranAtmosphere.Ivsa");
Double_t Vis = pConf->GetNum("LowtranAtmosphere.Vis");
Double_t Wss = pConf->GetNum("LowtranAtmosphere.Wss");
Double_t Whh = pConf->GetNum("LowtranAtmosphere.Whh");
Double_t Rainrt = pConf->GetNum("LowtranAtmosphere.Rainrt");
Double_t Gndalt = pConf->GetNum("LowtranAtmosphere.Gndalt");
/* */
fprintf(tape5,"%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%8.3f%7.2fn",
Model,Itype,Iemsct,Imult,M1,M2,M3,M4,M5,M6,Mdef,Im,Noprt,Tbound,Salb);
fprintf(tape5,"%5d%5d%5d%5d%5d%5d%10.3f%10.3f%10.3f%10.3f%10.3fn",
Ihaze,Iseasn,Ivulcn,Icstl,Icld,Ivsa,Vis,Wss,Whh,Rainrt,Gndalt);
Int_t Nmax=34,count=0;
fprintf(tape5,"%5d 0 0 user mode ln",Nmax);
Double_t pss,tss,convert;
Double_t Xco2,Xh2o,Xo3;
string gaslaw = "lowtran"; //can be ideal for debugging
const Double_t mO3=AtomicMass("TriOxygen");
const Double_t mH2O=AtomicMass("H2O");
const Double_t mCO2=AtomicMass("CO2");
for(Int_t i=0; i<fNbL; i++) {
pss=fPressureTable[i]/STP_Pressure;
tss=STP_Temperature/fTemperatureTable[i];
convert= 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
Double_t alt=fAltitudeTable[i]/km;
Xh2o= fH2O_DensityTable[i]/(convert*mH2O);
Xco2= fCO2_DensityTable[i]/(convert*mCO2);
Xo3 = fO3_DensityTable[i]/(convert*mO3);
if (gaslaw=="lowtran"){
Xo3=Xo3/(pow(pss,0.420)*pow(tss,1.3903));
Xco2=Xco2/(pow(pss,0.6705)*pow(tss,-2.2560));
Xh2o=Xh2o/(pow(pss,0.981)*pow(tss,0.3324));}
if ((alt<26 || (Int_t)alt%5==0) && count<Nmax ){
fprintf(tape5,"%10.3f %9.3f %9.3f %9.3e %9.3e %9.3en",
fAltitudeTable[i]/km,
fPressureTable[i]/(bar*0.001),
fTemperatureTable[i],
Xh2o,Xco2,Xo3);
count++;
}
}
/* Card 3 */
// Double_t H1=0.,H2=400.,Angle=0.,Range=0.,Beta=0.,R0=0.;
Double_t H1=0.,H2=400.,Angle=45.,Range=0.,Beta=0.,R0=0.;
Int_t Len=0;
fprintf(tape5,"%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%5dn",
H1,H2,Angle,Range,Beta,R0,Len);
/* Card 4 */
Double_t V1=22200.,V2=40000,DV=180;
// Double_t V1=25000.,V2=34500,DV=95;
fprintf(tape5,"%10.3f%10.3f%10.3fn",V1,V2,DV);
/* Card 5 */
Int_t Irpt=0;
fprintf(tape5,"%5dn",Irpt);
fclose(tape5);
#ifdef DEBUG
FILE* Debug=fopen("output/atmosphere.lis","w");
for(Int_t i=0; i<Nmax; i++) {
fprintf(Debug,
"%10.3f %9.3f %9.3e %9.3e %9.3e %9.3e %9.3e %9.3e\n",
fAltitudeTable[i]/km ,fTemperatureTable[i],
fPressureTable[i] ,fAir_DensityTable[i],
fO3_DensityTable[i],fH2O_DensityTable[i],
fO_DensityTable[i] ,fCO2_DensityTable[i]);}
fclose(Debug);
#endif
Msg(EsafMsg::Info) << "==> lowtran RT input file (tape5) closed " << MsgDispatch;
// lowtran initialization - done only once - function defined in LowtranRadiativeProcessesCalculator.cc
InitLOWTRAN();
Msg(EsafMsg::Info) << "Lowtran routines initialized" << MsgDispatch;
return;
}
// Plot the implemented tables
void AtmosphereData::DebugPlots(Int_t fNb) const {
Char_t nb [50];
Char_t name[50];
Msg(EsafMsg::Debug) << "AtmosphereData::DebugPlots : fNb= "<<fNb <<MsgDispatch;
static Int_t counter = 0;
sprintf(nb,"%d",counter++);
strcpy(name,"press");
strcat(name,nb);
TH1F* press = (TH1F*)gROOT->FindObject(name);
if(!press) press = new TH1F(name,"Pressure profile",50,0,50);
strcpy(name,"temp");
strcat(name,nb);
TH1F* temp = (TH1F*)gROOT->FindObject(name);
if(!temp) temp = new TH1F(name,"Temperature profile",100,0,100);
strcpy(name,"airdens");
strcat(name,nb);
TH1F* airdens = (TH1F*)gROOT->FindObject(name);
if(!airdens) airdens = new TH1F(name,"Air Density profile",50,0,50);
strcpy(name,"Odens");
strcat(name,nb);
TH1F* Odens = (TH1F*)gROOT->FindObject(name);
if(!Odens) Odens = new TH1F(name,"O Density profile",50,0,50);
strcpy(name,"O2dens");
strcat(name,nb);
TH1F* O2dens = (TH1F*)gROOT->FindObject(name);
if(!O2dens) O2dens = new TH1F(name,"O2 Density profile",50,0,50);
strcpy(name,"O3dens");
strcat(name,nb);
TH1F* O3dens = (TH1F*)gROOT->FindObject(name);
if(!O3dens) O3dens = new TH1F(name,"O3 Density profile",70,0,70);
strcpy(name,"N2dens");
strcat(name,nb);
TH1F* N2dens = (TH1F*)gROOT->FindObject(name);
if(!N2dens) N2dens = new TH1F(name,"N2 Density profile",50,0,50);
strcpy(name,"H2Odens");
strcat(name,nb);
TH1F* H2Odens = (TH1F*)gROOT->FindObject(name);
if(!H2Odens) H2Odens = new TH1F(name,"H2O Density profile",20,0,20);
strcpy(name,"CO2dens");
strcat(name,nb);
TH1F* CO2dens = (TH1F*)gROOT->FindObject(name);
if(!CO2dens) CO2dens = new TH1F(name,"CO2 Density profile",50,0,50);
/* for(Int_t i=0; i<fNb; i++) {
Double_t alt=fAltitudeTable[i]/km;
press->Fill(alt,fPressureTable[i]/bar);
temp->Fill(alt,fTemperatureTable[i]/kelvin);
airdens->Fill(alt,fAir_DensityTable[i]*m3/kg);
Odens->Fill(alt,fO_DensityTable[i]*m3/kg);
O2dens->Fill(alt,fO2_DensityTable[i]*m3/kg);
O3dens->Fill(alt,fO3_DensityTable[i]*m3/kg);
N2dens->Fill(alt,fN2_DensityTable[i]*m3/kg); */
string val;
for(Int_t i=0; i<fNb; i++) {
Double_t alt =fAltitudeTable[i]/km;
Double_t alt1=fAltitudeTable[i+1]/km;
press->Fill(alt,fPressureTable[i]/bar);
temp->Fill(alt,fTemperatureTable[i]/kelvin);
airdens->Fill(alt,fAir_DensityTable[i]*m3/kg);
Odens->Fill(alt,fO_DensityTable[i]*m3/kg);
O2dens->Fill(alt,fO2_DensityTable[i]*m3/kg);
O3dens->Fill(alt,fO3_DensityTable[i]*m3/kg);
N2dens->Fill(alt,fN2_DensityTable[i]*m3/kg);
CO2dens->Fill(alt,fCO2_DensityTable[i]*m3/kg);
H2Odens->Fill(alt,fH2O_DensityTable[i]*m3/kg);
if (alt1-alt>1){
for(Double_t h=alt+1; h<alt1; h++){
if(h!=(Int_t)alt1) {
val="temp";
temp->Fill(h,Smooth(val,i,h)/kelvin);
val="press";
press->Fill(h,Smooth(val,i,h)/bar);
val="Airdens";
airdens->Fill(h,Smooth(val,i,h)*m3/kg);
val="Odens";
Odens->Fill(h,Smooth(val,i,h)*m3/kg);
val="O2dens";
O2dens->Fill(h,Smooth(val,i,h)*m3/kg);
val="O3dens";
O3dens->Fill(h,Smooth(val,i,h)*m3/kg);
val="N2dens";
N2dens->Fill(h,Smooth(val,i,h)*m3/kg);
val="CO2dens";
CO2dens->Fill(h,Smooth(val,i,h)*m3/kg);
val="H2Odens";
H2Odens->Fill(h,Smooth(val,i,h)*m3/kg);
}}
}
}
}
Double_t AtmosphereData::Smooth(string& val, Int_t i, Double_t h) const {
Double_t val1, val2;
Double_t h1;
Double_t step;
string mod;
h1 = fAltitudeTable[i]/km;
step =(fAltitudeTable[i+1]-fAltitudeTable[i])/km;
if(val == "press"){
val1 = fPressureTable[i];
val2 = fPressureTable[i+1];}
else if(val == "temp"){
val1 = fTemperatureTable[i];
val2 = fTemperatureTable[i+1];}
else if(val == "Airdens"){
val1 = fAir_DensityTable[i];
val2 = fAir_DensityTable[i+1];}
else if(val == "Odens"){
val1 = fO_DensityTable[i];
val2 = fO_DensityTable[i+1];}
else if(val == "O2dens"){
val1 = fO2_DensityTable[i];
val2 = fO2_DensityTable[i+1];}
else if(val == "O3dens"){
val1 = fO3_DensityTable[i];
val2 = fO3_DensityTable[i+1];}
else if(val == "N2dens"){
val1 = fN2_DensityTable[i];
val2 = fN2_DensityTable[i+1];}
else if(val == "H2Odens"){
val1 = fH2O_DensityTable[i];
val2 = fH2O_DensityTable[i+1];}
else if(val == "CO2dens"){
val1 = fCO2_DensityTable[i];
val2 = fCO2_DensityTable[i+1];}
else
Msg(EsafMsg::Panic)<<"Wrong val argument in AtmosphereData::Smooth"<< MsgDispatch;
mod="expon";
if (val=="temp") mod="linear";
if(mod == "linear")
return (val2 - val1)*(h - h1)/step + val1;
else if(mod == "expon")
if (val2==0)
return 0;
else
return val1 * exp(-(h - h1) / (step/log(val1/val2)) );
else
Msg(EsafMsg::Panic)<< "Wrong mod argument in AtmosphereData::Smooth" << MsgDispatch;
return 0;
}