Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

LowtranAtmosphereData - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: LowtranAtmosphereData.cc,v 1.12 2005/10/02 14:17:07 thea 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;
using EConst::STP_Pressure;
using EConst::STP_Temperature;
using namespace sou;

//_____________________________________________________________________________
 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
}

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005