source: JEM-EUSO/esaf_cc_at_lal/packages/common/atmosphere/src/MSISE_00Atmosphere.cc @ 114

Last change on this file since 114 was 114, checked in by moretto, 11 years ago

actual version of ESAF at CCin2p3
File size: 10.1 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: MSISE_00Atmosphere.cc 2794 2008-02-11 07:56:34Z naumov $
3// S. Moreggia created 18 November 2003
4
5/*****************************************************************************
6 * ESAF: Euso Simulation and Analysis Framework                              *
7 *                                                                           *
8 *  Id: MSISE_00Atmosphere                                                   *
9 *  Package: atmosphere                                                      *
10 *  Coordinator: S. Moreggia                                                 *
11 *                                                                           *
12 *****************************************************************************/
13
14//_____________________________________________________________________________
15//
16// MSISE_00Atmosphere
17//
18// <extensive class description>
19//
20//   Config file parameters
21//   ======================
22//
23//   <parameter name>: <parameter description>
24//   -Valid options: <available options>
25//
26
27#include <math.h>
28#include "MSISE_00Atmosphere.hh"
29#include "MSISEtool.hh"
30#include "Config.hh"
31#include "MSISE_00AtmosphereData.hh"
32#include "EConst.hh"
33#include <TSpline.h>
34#include <TArrayD.h>
35ClassImp(MSISE_00Atmosphere)
36
37using EConst::Avogadro;
38using EConst::R_ideal;
39using EConst::AtomicMass;
40using namespace sou;
41using namespace TMath;
42
43//____________________________________________________________________________
44MSISE_00Atmosphere::MSISE_00Atmosphere(string type) : Atmosphere(), fData(0) {
45    //
46    // ctor
47    //
48    fType = type;
49    Msg(EsafMsg::Info) << "*** MSISE_00Atmosphere ***"<< MsgDispatch;
50    Build();
51}
52
53//____________________________________________________________________________
54MSISE_00Atmosphere::~MSISE_00Atmosphere() {
55    //
56    // dtor
57    //
58    SafeDelete(fData);
59}
60
61//___________________________________________________________________________________
62void MSISE_00Atmosphere::CreateInstance() {
63    //
64    // Create the single atmosphere instance as being a MSISE_00Atmosphere object
65    //
66    string type = "msise00";
67    if(!Atmosphere::fChild) Atmosphere::fChild = new MSISE_00Atmosphere(type);
68    else Atmosphere::fChild->Msg(EsafMsg::Panic) << "Trying to create two atmospheres !!"<< MsgDispatch;
69}
70
71//___________________________________________________________________________________
72void MSISE_00Atmosphere::ResetInstance() {
73    //
74    // in random mode, reset tabulated profiles and randomly build a new atmosphere
75    //
76   
77    SafeDelete(fData);
78    Msg(EsafMsg::Info) << "*** Build a new MSISE_00Atmosphere ***"<< MsgDispatch;
79    Build();
80}
81
82//____________________________________________________________________________
83void MSISE_00Atmosphere::Build() {
84    //
85    // Build atmosphere according to config parameters
86    // Data arrays are built
87    //
88
89    if(!fData) {
90        fData =  new MSISE_00AtmosphereData(IsRandomMode());
91        if(!fData) Msg(EsafMsg::Panic)<<"No new MSISEData, memory pb"<<MsgDispatch;
92    }
93
94    const Double_t* data = 0;
95    // interpolators initialization
96    data =  fData->GetTemperatureTable();
97    TArrayD t_arr(fData->NumberOfElements(),data);
98    fTempInterpol = new TSpline3("fTempInterpol",0,fData->NumberOfElements()*km,t_arr.GetArray(),fData->NumberOfElements());   
99
100    data =  fData->GetAir_DensityTable();
101    TArrayD rho_arr(fData->NumberOfElements(),data);
102    fDensInterpol = new TSpline3("fDensInterpol",0,fData->NumberOfElements()*km,rho_arr.GetArray(),fData->NumberOfElements());   
103
104    data =  fData->GetPressureTable();
105    TArrayD press_arr(fData->NumberOfElements(),data);
106    fPressInterpol = new TSpline3("fPressInterpol",0,fData->NumberOfElements()*km,press_arr.GetArray(),fData->NumberOfElements());   
107}
108
109//____________________________________________________________________________
110Double_t MSISE_00Atmosphere::Interpolate(string& val,string& mod,Double_t h) const {
111    //
112    // Calculate interpolation
113    //   
114    const Double_t* data = 0;
115    if(val == "Odens")            data = fData->GetO_DensityTable();
116    else if(val == "O2dens")      data = fData->GetO2_DensityTable();
117    else if(val == "H2Odens")     data = fData->GetH2O_DensityTable();
118    else if(val == "N2dens")      data = fData->GetN2_DensityTable();
119    else if(val == "O3dens")      data = fData->GetO3_DensityTable();
120    else if(val == "O3densPPMV")  data = fData->GetO3_DensityPPMVTable();
121    else  Msg(EsafMsg::Panic)<< "Wrong val argument in MSISE_00Atmosphere::Interpolate ==> "<<val <<  MsgDispatch;
122   
123    if(h < 0) return data[0];
124    Int_t nb = fData->NumberOfElements();
125    Double_t step = fData->LayerSize();
126    if(h > step*(nb-1)) return data[nb-1];
127    Double_t val1, val2;
128    Double_t h1;
129    Int_t i = 0;
130    i = Int_t(floor(h/step));
131    h1 = i * step;
132    val1 = data[i];
133    val2 = data[i+1];
134   
135    if(mod == "linear" || val1 == 0 || val2 == 0)  return (val2 - val1)*(h - h1)/step + val1;
136    else if(mod == "expon")                        return val1 * exp(-(h - h1) * log(val1/val2)/step);
137    else Msg(EsafMsg::Panic)<<"Wrong mod argument in MSISE_00Atmosphere::Interpolate" << MsgDispatch ;
138   
139    return 0;
140}
141
142//____________________________________________________________________________
143Double_t MSISE_00Atmosphere::Pressure(Double_t h) const {
144    //
145    // Pressure at a given altitude
146    //
147    if(h > fTOA)  return 0.;
148
149    // may occur with interfaces with external codes
150    if(IsNaN(h)) Msg(EsafMsg::Panic)<<"<Pressure> Given altitude is NaN" << MsgDispatch ;
151    return fPressInterpol->Eval(h);
152}
153
154//____________________________________________________________________________
155Double_t MSISE_00Atmosphere::Temperature(Double_t h) const {
156    //
157    // Temperature at a given altitude
158    //
159    if(h > fTOA)  return 0.;
160
161    // may occur with interfaces with external codes
162    if(IsNaN(h)) Msg(EsafMsg::Panic)<<"<Pressure> Given altitude is NaN" << MsgDispatch ;
163    return fTempInterpol->Eval(h);
164}
165
166//____________________________________________________________________________
167Double_t MSISE_00Atmosphere::AbsoluteHumidity(Double_t h) const {
168    //
169    // Absolute humidity at a given altitude (water vapor density)
170    //
171   
172    if(h > fTOA)  return 0.;
173    string val = "H2Odens";
174    string mod = "linear";
175   
176    return Interpolate(val,mod,h);
177}
178
179//____________________________________________________________________________
180Double_t MSISE_00Atmosphere::Air_Density(Double_t h) const {   
181    //
182    // Atmospheric air density at a given altitude
183    //
184    if(h > fTOA)  return 0.;
185
186    // may occur with interfaces with external codes
187    if(IsNaN(h)) Msg(EsafMsg::Panic)<<"<Pressure> Given altitude is NaN" << MsgDispatch ;
188    return fDensInterpol->Eval(h);
189}
190
191//____________________________________________________________________________
192Double_t MSISE_00Atmosphere::O_Density(Double_t h) const {
193    //
194    // Oxygen density at a given altitude
195    //
196    if(h > fTOA)
197        return 0.;
198
199    // exponential interpolation between table values
200    string val = "Odens";
201    string mod = "expon";
202    return Interpolate(val,mod,h);
203}
204
205//____________________________________________________________________________
206Double_t MSISE_00Atmosphere::O2_Density(Double_t h) const {
207    //
208    // Di-oxygen density at a given altitude
209    //
210    if(h > fTOA)
211        return 0.;
212
213    // exponential interpolation between table values
214    string val = "O2dens";
215    string mod = "expon";
216    return Interpolate(val,mod,h);
217}
218
219//____________________________________________________________________________
220Double_t MSISE_00Atmosphere::O3_Density(Double_t h) const {
221    //
222    // Ozone density at a given altitude
223    //
224   
225    // exponential interpolation between table values
226    string val = "O3dens";
227    string mod = "linear";
228    return Interpolate(val,mod,h);
229}
230
231//____________________________________________________________________________
232Double_t MSISE_00Atmosphere::O3_DensityPPMV(Double_t h) const {
233    //
234    // Ozone PPMV density at a given altitude
235    //
236   
237    // exponential interpolation between table values
238    string val = "O3densPPMV";
239    string mod = "linear";
240    return Interpolate(val,mod,h);
241}
242
243//____________________________________________________________________________
244Double_t MSISE_00Atmosphere::N2_Density(Double_t h) const {
245    //
246    // Di-nitrogen density at a given altitude
247    //
248    if(h > fTOA)
249        return 0.;
250
251    // exponential interpolation between table values
252    string val = "N2dens";
253    string mod = "expon";
254    return Interpolate(val,mod,h);
255}
256
257//____________________________________________________________________________
258Double_t MSISE_00Atmosphere::Aerosols_Density(string& type,Double_t h) const {
259    //
260    // Aerosols density for a given type of aerosols at a given altitude
261    //
262    Msg(EsafMsg::Panic)<<"Aerosol_Density not available for MSISE_00Atmosphere"
263                         << MsgDispatch ;
264    return 0;
265}
266
267
268/*
269// For MSISEscan Macro --- WARNING : Build() call in ctor MUST BE DISABLED
270//____________________________________________________________________________
271void MSISE_00Atmosphere::Build(Int_t doy,Double_t lst, Double_t lat, Double_t longi) {
272    //
273    // Build atmosphere according to config parameters
274    // Data arrays are built
275    //
276
277    if(!fData) {
278        fData =  new MSISE_00AtmosphereData(IsRandomMode());
279        if(!fData) Msg(EsafMsg::Panic)<<"No new MSISEData, memory pb"<<MsgDispatch;
280    }
281   
282    fData->InitForMacro(doy,lst,lat,longi);
283    fData->SetTables();
284
285    const Double_t* data = 0;
286    // interpolators initialization
287    data =  fData->GetTemperatureTable();
288    TArrayD t_arr(fData->NumberOfElements(),data);
289    fTempInterpol = new TSpline3("fTempInterpol",0,fData->NumberOfElements()*km,t_arr.GetArray(),fData->NumberOfElements());   
290
291    data =  fData->GetAir_DensityTable();
292    TArrayD rho_arr(fData->NumberOfElements(),data);
293    fDensInterpol = new TSpline3("fDensInterpol",0,fData->NumberOfElements()*km,rho_arr.GetArray(),fData->NumberOfElements());   
294
295    data =  fData->GetPressureTable();
296    TArrayD press_arr(fData->NumberOfElements(),data);
297    fPressInterpol = new TSpline3("fPressInterpol",0,fData->NumberOfElements()*km,press_arr.GetArray(),fData->NumberOfElements());   
298}
299*/
300
Note: See TracBrowser for help on using the repository browser.