source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/common/atmosphere/src/LinsleyAtmosphereData.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: LinsleyAtmosphereData.cc 2679 2006-05-03 08:28:58Z moreggia $
// Sylvain Moreggia created May, 17 2004

#include "LinsleyAtmosphereData.hh"
#include "NumbersFileParser.hh"
#include <math.h>

using namespace sou;

ClassImp(LinsleyAtmosphereData)

//_____________________________________________________________________________
LinsleyAtmosphereData::LinsleyAtmosphereData() : AtmosphereData() {
    //
    // ctor
    //
    fNb = (Int_t)Conf()->GetNum("LinsleyAtmosphereData.fNb");
    fName = Conf()->GetStr("LinsleyAtmosphereData.fName");

    fBottom = new Double_t[fNb];
    fA = new Double_t[fNb];
    fB = new Double_t[fNb];
    fC = new Double_t[fNb];
    if(!(fBottom && fA && fB && fC)) Msg(EsafMsg::Panic) <<"Pb for memory allocation in LinsleyAtmosphereData ctor"<<MsgDispatch;
    ParamTables();
    SetTables();
}

//_____________________________________________________________________________
LinsleyAtmosphereData::LinsleyAtmosphereData(const Int_t&, const vector<Double_t>&, const Int_t&)
    : AtmosphereData() {
    //
    // ctor
    //
    fBottom = 0;
    fA      = 0;
    fB      = 0;
    fC      = 0;
}

//_____________________________________________________________________________
LinsleyAtmosphereData::~LinsleyAtmosphereData() {
    //
    // dtor
    //
}

//_____________________________________________________________________________
void LinsleyAtmosphereData::ParamTables() {
    //
    // fill data members (tables)
    //
    string fn = fName + ".data";
    NumbersFileParser nf("config/Atmosphere/LinsleyData/" + fn, fNb);
    Msg(EsafMsg::Info)<< "  Linsley Atmosphere - config :" << fName << MsgDispatch;

    vector<Double_t> param;
    for(Int_t i=0; i<fNb; i++) {
        param = nf.GetCol(i);
        fBottom[i] = param[0]*km;
        fA[i] = param[2]*g/cm2;
        fB[i] = param[3]*g/cm2;
        fC[i] = param[4]*cm;
    }
}

//_____________________________________________________________________________
void LinsleyAtmosphereData::SetTables() {
    //
    // fill data members (tables)
    //
    Int_t fNbL=100;
    fAltitudeTable = new Double_t[fNbL];
    fAir_DensityTable = new Double_t[fNbL];
    
    for(Int_t j=0; j<fNbL; j++){
      Double_t h = Double_t(j);
      fAltitudeTable[j] = h*km;
      for(Int_t i=0; i<fNb; i++) {
        if(fBottom[i]/km <= h && h < fBottom[i+1]/km) {
          fAir_DensityTable[j]= fB[i]/fC[i] * exp(-h*km/fC[i]); }

      }
    }
    string atmod = "linsley";
    Msg(EsafMsg::Warning)<< "<SetTables> GetDefault() and WriteUserModelFile() disabled in LinseyAtmosphere mode" << MsgDispatch;
    //GetDefault(atmod,fNbL);
    //WriteUserModelFile(fNbL);   // interface with Lowtran needs to be fixed
#ifdef DEBUG
    //DebugPlots(fNbL);    
#endif
}

//_____________________________________________________________________________
Double_t LinsleyAtmosphereData::DensityParametrization(Double_t h) const {
    //
    // Linsley parametrization
    //
    Int_t i;
    if(h < -kAltitudeTolerance)  {
        Msg(EsafMsg::Warning) << "Negative altitude in LinsleyAtmosphereData, h = "<<h << MsgDispatch;
        Msg(EsafMsg::Warning) << "ZERO density returned"<< MsgDispatch;
        return 0.;
    }
    if(h < 0) h = 0.;
    
    if(fName == "Uniform") {
        if(h < fBottom[1]) return fB[0] / fC[0];
        else return fB[1] / fC[1];
    }
    
    for(i=0; i<fNb-1; i++) {
        if(h < fBottom[i+1]) return fB[i]/fC[i] * exp(-h/fC[i]);
    }
    return fB[i]/fC[i];
}