Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

AtmosphereData - source file

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

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(0), val2(0);
    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;
}
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