Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

MSISE_00Atmosphere - source file

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

//____________________________________________________________________________
 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 = 0;
    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;
}




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