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

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

ESAF version compilable on mac OS
File size: 22.5 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: AtmosphereData.cc 2918 2011-06-10 22:22:31Z mabl $
3// Sylvain Moreggia created Dec,  1 2003
4
5#include "AtmosphereData.hh"
6#include "Atmosphere.hh"
7#include "LinsleyAtmosphere.hh"
8#include "NumbersFileParser.hh"
9#include "ConfigFileParser.hh"
10#include "Config.hh"
11#include "LowtranManager.hh"
12#include <TH1F.h>
13#include <TROOT.h>
14#include "EConst.hh"
15
16ClassImp(AtmosphereData)
17
18using EConst::Avogadro;
19using EConst::Loschmidt; 
20using EConst::R_ideal;
21using EConst::AtomicMass;
22using EConst::STP_Pressure;
23using EConst::STP_Temperature;
24using namespace sou;
25
26//extern void InitLOWTRAN();
27
28//______________________________________________________________________________
29AtmosphereData::AtmosphereData() : EsafConfigurable(),  EsafMsgSource(){
30    //
31    // ctor
32    //
33    fAltitudeTable = 0;
34    fPressureTable = 0;
35    fTemperatureTable = 0;
36    fAir_DensityTable = 0;
37    fO_DensityTable = 0;
38    fO2_DensityTable = 0;
39    fO3_DensityTable = 0;
40    fO3_DensityPPMVTable = 0;
41    fN2_DensityTable = 0;
42    fCO2_DensityTable = 0;
43    fH2O_DensityTable = 0;
44    fIndexTable = 0;
45}
46
47//______________________________________________________________________________
48AtmosphereData::~AtmosphereData() {
49    //
50    // dtor
51    //
52    if(fAltitudeTable) delete[] fAltitudeTable;
53    fAltitudeTable = 0;
54    if(fPressureTable) delete[] fPressureTable;
55    fPressureTable = 0;
56    if(fTemperatureTable) delete[] fTemperatureTable;
57    fTemperatureTable = 0;
58    if(fAir_DensityTable) delete[] fAir_DensityTable;
59    fAir_DensityTable = 0;
60    if(fO_DensityTable) delete[] fO_DensityTable;
61    fO_DensityTable = 0;
62    if(fO2_DensityTable) delete[] fO2_DensityTable;
63    fO2_DensityTable = 0;
64    if(fO3_DensityTable) delete[] fO3_DensityTable;
65    fO3_DensityTable = 0; 
66    if(fO3_DensityPPMVTable) delete[] fO3_DensityPPMVTable;
67    fO3_DensityPPMVTable = 0; 
68    if(fN2_DensityTable) delete[] fN2_DensityTable;
69    fN2_DensityTable = 0; 
70    if(fCO2_DensityTable) delete[] fCO2_DensityTable;
71    fCO2_DensityTable = 0; 
72    if(fH2O_DensityTable) delete[] fH2O_DensityTable;
73    fH2O_DensityTable = 0; 
74    if(fIndexTable) delete[] fIndexTable;
75    fIndexTable = 0; 
76
77}
78
79//______________________________________________________________________________
80void AtmosphereData::GetDefault(string& atmod, Int_t fNbL) {
81    //
82    // fill missing data with default values
83    //
84    if (!fAltitudeTable){
85        fAltitudeTable = new Double_t[fNbL];}
86    if (!fPressureTable){
87        fPressureTable = new Double_t[fNbL];}
88    if (!fTemperatureTable){
89        fTemperatureTable = new Double_t[fNbL];}
90    if (!fAir_DensityTable){
91        fAir_DensityTable = new Double_t[fNbL];}
92    if (!fO_DensityTable){
93        fO_DensityTable = new Double_t[fNbL];}
94    if (!fO2_DensityTable){
95        fO2_DensityTable = new Double_t[fNbL];}
96    if (!fO3_DensityTable){
97        fO3_DensityTable = new Double_t[fNbL];}
98    if (!fO3_DensityPPMVTable){
99        fO3_DensityPPMVTable = new Double_t[fNbL];}
100    if (!fN2_DensityTable){
101        fN2_DensityTable = new Double_t[fNbL];}
102    if (!fCO2_DensityTable){
103        fCO2_DensityTable = new Double_t[fNbL];}
104    if (!fH2O_DensityTable){
105        fH2O_DensityTable = new Double_t[fNbL];}
106    if (!fIndexTable){
107        fIndexTable = new Double_t[fNbL];}
108
109    string gaslaw = "lowtran";   // can be ideal for debugging
110
111    /* If Linsley : complete data with USStandard  and MSISE    */
112    /*              profiles                                    */
113    /* If Msise   : complete data with lowtran models,          */ 
114    /*              according to latitude, longitude            */
115    /* If Lowtran : complete data with Msise, according to model*/
116
117    const Double_t R=R_ideal()/joule;
118    const Double_t mAir=AtomicMass("Air");
119    const Double_t mO=AtomicMass("Oxygen");
120    const Double_t mO2=AtomicMass("DiOxygen");
121    const Double_t mO3=AtomicMass("TriOxygen");
122    const Double_t mN2=AtomicMass("DiNitrogen");
123    const Double_t mH2O=AtomicMass("H2O");
124    const Double_t mCO2=AtomicMass("CO2");
125
126
127    if(atmod =="linsley"){
128        Msg(EsafMsg::Info) << "The present model is Linsley" << MsgDispatch; 
129        Msg(EsafMsg::Info) << "Missing profiles are provided by default files" <<MsgDispatch;
130        Msg(EsafMsg::Info) << "US Standard and MSISE lat=45, long=270, day=1"<<MsgDispatch;
131
132        NumbersFileParser nf("config/Atmosphere/DefaultData/Default_6",7);
133        vector<Double_t> alt   = nf.GetCol(0);
134        vector<Double_t> temp  = nf.GetCol(1);
135        vector<Double_t> pres  = nf.GetCol(2);
136        vector<Double_t> h2o   = nf.GetCol(3);
137        vector<Double_t> ozone = nf.GetCol(4); 
138        vector<Double_t> co2   = nf.GetCol(6);
139
140        NumbersFileParser ng("config/Atmosphere/DefaultData/lat_45_long_270_day_1",4);
141        //    vector<Double_t> alt   = nf.GetCol(0);
142        vector<Double_t> n2    = ng.GetCol(1);
143        vector<Double_t> o2    = ng.GetCol(2);
144        vector<Double_t> o     = ng.GetCol(3);
145
146        Double_t pss,tss,convert;
147
148        Double_t interpol_temp(0);
149        for(Int_t i=0; i<fNbL; i++) {
150            interpol_temp            = LinearInterpolateLowtran(i,alt,temp);
151            fTemperatureTable[i]     = interpol_temp;
152            fPressureTable[i]        = fAir_DensityTable[i] / (mAir*g/mole) * (R * joule/kelvin/mole) * fTemperatureTable[i];
153            pss                      = fPressureTable[i]/STP_Pressure();
154            tss                      = STP_Temperature()/fTemperatureTable[i];
155            convert                  = 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
156            fO2_DensityTable[i]      = LinearInterpolateLowtran(i,alt,o2)*mO2*convert;
157            fO3_DensityTable[i]      = LinearInterpolateLowtran(i,alt,ozone)*mO3*convert;
158            fO3_DensityPPMVTable[i]  = LinearInterpolateLowtran(i,alt,ozone);
159            fN2_DensityTable[i]      = LinearInterpolateLowtran(i,alt,n2)*mN2*convert;
160            fCO2_DensityTable[i]     = LinearInterpolateLowtran(i,alt,co2)*mCO2*convert;
161            fH2O_DensityTable[i]     = LinearInterpolateLowtran(i,alt,h2o)*mH2O*convert;
162            fO_DensityTable[i]       = LinearInterpolateLowtran(i,alt,o)*convert*mO;
163
164            if (gaslaw=="lowtran"){
165                fO2_DensityTable[i]  = fO2_DensityTable[i]*pow(pss,1.1879)*pow(tss,2.9738);
166                fO3_DensityTable[i]  = fO3_DensityTable[i]*pow(pss,0.420)*pow(tss,1.3903);
167                fN2_DensityTable[i]  = 0.781*1.e6*mN2*convert *pow(pss,0.5)*pow(tss,0.5);
168                fCO2_DensityTable[i] = fCO2_DensityTable[i]*pow(pss,0.6705)*pow(tss,-2.2560);
169                fH2O_DensityTable[i] = fH2O_DensityTable[i]*pow(pss,0.981)*pow(tss,0.3324);}
170        }
171    }
172
173
174    if(atmod =="msise"){
175        Double_t day    = this->GetDoy();
176        Double_t g_lat  = this->GetLatitude();
177        Double_t g_long = this->GetLongitude();
178
179        string fmodel;
180        // get default according to latitude and season
181        if(fabs(g_lat) <= 25.) fmodel ="Default_1";  //Tropical
182        else{ 
183            if(80 <= day && day <= 266) {
184                if(g_lat > 0) fmodel="Default_2"; //MidLatSummer
185                else          fmodel="Default_3"; //MidLatWinter
186            }
187            else {
188                if(g_lat > 0) fmodel="Default_3"; //MidLatWinter
189                else          fmodel="Default_2"; //MidLatSummer
190            }
191        } 
192        Msg(EsafMsg::Info) << "The present model is Msise : lat = "<<g_lat << " long = "<<g_long << " day = " <<day <<MsgDispatch; 
193        Msg(EsafMsg::Info) << "Missing profiles are provided by files" <<MsgDispatch;   
194        Msg(EsafMsg::Info) << "from the corresponding Lowtran Model : "<< fmodel << MsgDispatch;
195       
196
197        NumbersFileParser nf("config/Atmosphere/DefaultData/" + fmodel,7);
198
199        vector<Double_t> alt = nf.GetCol(0);
200        vector<Double_t> h2o = nf.GetCol(3);
201        vector<Double_t> ozone= nf.GetCol(4);
202        vector<Double_t> co2 = nf.GetCol(6);
203
204        Double_t pss,tss,convert;
205
206        for(Int_t i=0; i<fNbL; i++) {
207            pss                      = fPressureTable[i]/STP_Pressure();
208            tss                      = STP_Temperature()/fTemperatureTable[i];
209            convert                  = 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
210            fO3_DensityTable[i]      = LinearInterpolateLowtran(i,alt,ozone)*mO3*convert;
211            fO3_DensityPPMVTable[i]  = LinearInterpolateLowtran(i,alt,ozone);
212            fCO2_DensityTable[i]     = LinearInterpolateLowtran(i,alt,co2)*mCO2*convert;
213            fH2O_DensityTable[i]     = LinearInterpolateLowtran(i,alt,h2o)*mH2O*convert;
214            if (gaslaw=="lowtran"){
215                fO3_DensityTable[i]  = fO3_DensityTable[i]*pow(pss,0.420)*pow(tss,1.3903);
216                fCO2_DensityTable[i] = fCO2_DensityTable[i]*pow(pss,0.6705)*pow(tss,-2.2560);
217                fH2O_DensityTable[i] = fH2O_DensityTable[i]*pow(pss,0.981)*pow(tss,0.3324);}
218        }
219    }
220   
221   
222    if(atmod =="lowtran"){
223        Msg(EsafMsg::Info) << "The present model is Lowtran" << MsgDispatch ;
224        ConfigFileParser* pConf = Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
225        Int_t Model = (Int_t)pConf->GetNum("LowtranAtmosphere.model");
226        string fmodel;
227        if (Model==1)     fmodel ="lat_0_long_0_day_1";        // Tropical model
228        else{ 
229            if (Model==2) fmodel="lat_45_long_270_day_200";    //"MidLatSummer"
230            else          fmodel="lat_45_long_270_day_1";      //MidLatWinter or US Standard
231        }
232
233        Msg(EsafMsg::Info) << "Missing profiles are provided by files" <<MsgDispatch;   
234        Msg(EsafMsg::Info) << "from the corresponding MSISE Model : "<< fmodel << MsgDispatch;
235
236        NumbersFileParser nf("config/Atmosphere/DefaultData/" + fmodel,4);
237
238
239        vector<Double_t> alt = nf.GetCol(0);
240        vector<Double_t> n2 = nf.GetCol(1);
241        vector<Double_t> o2 = nf.GetCol(2);
242        vector<Double_t> o  = nf.GetCol(3);
243
244
245
246        Double_t pss,tss,convert;
247        for(Int_t i=0; i<fNbL; i++) {
248
249            pss=fPressureTable[i]/STP_Pressure();
250            tss=STP_Temperature()/fTemperatureTable[i];
251            convert= 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
252
253            for(Int_t j=0; j<101; j++) {
254                if (alt[j]==fAltitudeTable[i]/km ||
255                        (Int_t)alt[j]==(Int_t)(fAltitudeTable[i]/km)){
256                    fO2_DensityTable[i] = o2[j]*convert*mO2;
257                    fO_DensityTable[i]  = o[j]*convert*mO;}}
258
259        }
260    }
261#ifdef DEBUG
262    Msg(EsafMsg::Debug) 
263        << "densities computation using " << gaslaw << " gas law "<< MsgDispatch; 
264#endif
265    return;
266}
267
268//______________________________________________________________________________
269void AtmosphereData::WriteUserModelFile(Int_t fNbL) {
270    //
271    // prepare and write a file necessary to use lowtran radiative
272    //   process calculator. It correspond to tape5 */
273    //
274
275    LowtranManager *l = LowtranManager::Get();
276    FILE* tape5=fopen(l->GetTape5(),"w");
277    ConfigFileParser* pConf =
278        Config::Get()->GetCF("Atmosphere","LowtranAtmosphere");
279
280    /* prepare CARD1                          */
281    /* the model has been defined by the user */
282    const Int_t Model=7,M1=0,M2=0,M3=0,M4=6,M5=6,M6=6;
283    const Int_t Mdef  =1;  //
284    const Int_t Itype =2;  // vertical or slant path between two altitudes
285    const Int_t Iemsct=0;  // program execution in transmittance mode
286    const Int_t Imult =0;  // no multiple scattering in transmittance mode
287    const Double_t Tbound=0;  // FOR NORMAL OPERATION OF PROGRAM. 
288
289
290    Int_t Noprt     = (Int_t)pConf->GetNum("LowtranAtmosphere.Noprt"); 
291    Int_t Im        = (Int_t)pConf->GetNum("LowtranAtmosphere.Im");
292    Double_t Salb   = pConf->GetNum("LowtranAtmosphere.Salb");
293
294    Int_t Ihaze     = (Int_t)pConf->GetNum("LowtranAtmosphere.Ihaze");
295    Int_t Iseasn    = (Int_t)pConf->GetNum("LowtranAtmosphere.Iseasn");
296    Int_t Ivulcn    = (Int_t)pConf->GetNum("LowtranAtmosphere.Ivulcn");
297    Int_t Icstl     = (Int_t)pConf->GetNum("LowtranAtmosphere.Icstl");
298    Int_t Icld      = (Int_t)pConf->GetNum("LowtranAtmosphere.Icld");
299    Int_t Ivsa      = (Int_t)pConf->GetNum("LowtranAtmosphere.Ivsa");
300
301    Double_t Vis    = pConf->GetNum("LowtranAtmosphere.Vis");
302    Double_t Wss    = pConf->GetNum("LowtranAtmosphere.Wss");
303    Double_t Whh    = pConf->GetNum("LowtranAtmosphere.Whh");
304    Double_t Rainrt = pConf->GetNum("LowtranAtmosphere.Rainrt");
305    Double_t Gndalt = pConf->GetNum("LowtranAtmosphere.Gndalt");
306
307    Double_t Linf     = pConf->GetNum("LowtranAtmosphere.Linf");
308    Double_t Lsup     = pConf->GetNum("LowtranAtmosphere.Lsup");
309
310    /*  */
311
312    fprintf(tape5,"%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%8.3f%7.2f\n",
313            Model,Itype,Iemsct,Imult,M1,M2,M3,M4,M5,M6,Mdef,Im,Noprt,Tbound,Salb);
314    fprintf(tape5,"%5d%5d%5d%5d%5d%5d%10.3f%10.3f%10.3f%10.3f%10.3f\n",
315            Ihaze,Iseasn,Ivulcn,Icstl,Icld,Ivsa,Vis,Wss,Whh,Rainrt,Gndalt);
316
317    Int_t Nmax=34,count=0;
318    fprintf(tape5,"%5d    0    0    user mode l\n",Nmax); 
319    Double_t pss,tss,convert;
320    Double_t Xco2,Xh2o,Xo3;
321    string gaslaw = "lowtran";   //can be ideal for debugging
322
323
324    const Double_t mO3=AtomicMass("TriOxygen");
325    const Double_t mH2O=AtomicMass("H2O");
326    const Double_t mCO2=AtomicMass("CO2");
327
328    for(Int_t i=0; i<fNbL; i++) {
329
330        pss=fPressureTable[i]/STP_Pressure();
331        tss=STP_Temperature()/fTemperatureTable[i];
332        convert= 1.e-6*(Loschmidt()*cm3/Avogadro())*pss*tss*g/cm3;
333        Double_t alt=fAltitudeTable[i]/km;
334        Xh2o= fH2O_DensityTable[i]/(convert*mH2O);
335        Xco2= fCO2_DensityTable[i]/(convert*mCO2);
336        Xo3 = fO3_DensityTable[i]/(convert*mO3);
337        if (gaslaw=="lowtran"){
338            Xo3=Xo3/(pow(pss,0.420)*pow(tss,1.3903));
339            Xco2=Xco2/(pow(pss,0.6705)*pow(tss,-2.2560));
340            Xh2o=Xh2o/(pow(pss,0.981)*pow(tss,0.3324));}
341
342            // WARNING : this way of writing TAPE5 (used for lowtran RT inputs) entails that
343            // LOWTRAN RT code "see" the atmosphere until 80 km ONLY
344            // NOT A PB, because RT properties in UV not affected above
345            if(alt < 26 && count<Nmax ){
346                fprintf(tape5,"%10.3f %9.3f %9.3f %9.3e %9.3e %9.3e\n",   
347                        fAltitudeTable[i]/km,
348                        fPressureTable[i]/(bar*0.001),
349                        fTemperatureTable[i],
350                        Xh2o,Xco2,Xo3);
351                count++;
352            }
353            else if((alt < 51) && (Int_t(alt)%5==0) && (count<Nmax) ) {
354                fprintf(tape5,"%10.3f %9.3f %9.3f %9.3e %9.3e %9.3e\n",   
355                        fAltitudeTable[i]/km,
356                        fPressureTable[i]/(bar*0.001),
357                        fTemperatureTable[i],
358                        Xh2o,Xco2,Xo3);
359                count++;
360            }
361            else if((alt < 100) && (Int_t(alt-50)%10==0) && (count<Nmax) ) {
362                fprintf(tape5,"%10.3f %9.3f %9.3f %9.3e %9.3e %9.3e\n",   
363                        fAltitudeTable[i]/km,
364                        fPressureTable[i]/(bar*0.001),
365                        fTemperatureTable[i],
366                        Xh2o,Xco2,Xo3);
367                count++;
368            }
369           
370           
371           
372    }
373    /* Card 3 */
374    //    Double_t H1=0.,H2=400.,Angle=0.,Range=0.,Beta=0.,R0=0.;
375    Double_t H1=0.,H2=400.,Angle=45.,Range=0.,Beta=0.,R0=0.;
376    Int_t Len=0;
377    fprintf(tape5,"%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%5d\n",
378            H1,H2,Angle,Range,Beta,R0,Len);
379
380    /* Card 4 */
381    //    Double_t V1=22200.,V2=40000,DV=180;
382    //    Double_t V1=25000.,V2=34500,DV=95;
383    Double_t V1=floor(1e7/Lsup);
384    Double_t V2=floor(1e7/Linf);
385    Double_t DV = floor((V2-V1)/100.);
386    if (fmod(DV,5)!=0) DV=DV+5.;
387    // this last line to avoid two many lambda compared to the table
388    // dimension, since lowtran computes DV to be a multiple of 5
389    fprintf(tape5,"%10.3f%10.3f%10.3f\n",V1,V2,DV);
390
391    /* Card 5 */
392    Int_t Irpt=0;
393    fprintf(tape5,"%5d\n",Irpt);
394
395    fclose(tape5);
396
397#ifdef DEBUG
398    FILE* Debug=fopen("output/atmosphere.lis","w");
399
400    for(Int_t i=0; i<Nmax; i++) {
401        fprintf(Debug,
402                "%10.3f %9.3f %9.3e %9.3e %9.3e %9.3e %9.3e %9.3e\n",   
403                fAltitudeTable[i]/km  ,fTemperatureTable[i],
404                fPressureTable[i]  ,fAir_DensityTable[i],
405                fO3_DensityTable[i],fH2O_DensityTable[i],
406                fO_DensityTable[i] ,fCO2_DensityTable[i]);}
407    fclose(Debug);
408#endif
409
410    Msg(EsafMsg::Info) << "==>  lowtran RT input file (tape5) closed    " << MsgDispatch;
411
412    // lowtran initialization - done only once - function defined in LowtranRadiativeProcessesCalculator.cc
413//    InitLOWTRAN();
414    //l->InitLowtran();
415    Msg(EsafMsg::Info) << "Lowtran routines initialized" << MsgDispatch;
416
417    return;
418}
419
420
421//_____________________________________________________________________________
422Double_t AtmosphereData::LinearInterpolateLowtran(Int_t h, const vector<Double_t>& alt, const vector<Double_t>& data) {
423    //
424    // To translate Lowtran atmosphere profile varying step to esaf constant step
425    // Achieve constant steps when SetTables is called
426    //
427    Double_t val1, val2;
428    Double_t h1;
429    Double_t step;
430    Int_t i = 0;
431   
432    if(h <= 25) return data[h];
433    else if(h <= 50) {
434        step = 2.5;
435        i = 25 + Int_t( floor(Double_t(h - 25)/step) );
436    }
437    else {
438        step = 5;
439        i = 35 + Int_t( floor(Double_t(h - 50)/step) );
440    }
441   
442    h1 = alt[i];
443    val1 = data[i];
444    val2 = data[i+1];
445    return (val2 - val1)*(Double_t(h) - h1)/step + val1;
446}
447
448
449
450//_____________________________________________________________________________
451/*
452// Plot the implemented tables
453void AtmosphereData::DebugPlots(Int_t fNb) const {
454    Char_t nb [50];
455    Char_t name[50];
456    Msg(EsafMsg::Debug) << "AtmosphereData::DebugPlots : fNb=  "<<fNb <<MsgDispatch;
457    static Int_t counter = 0;
458    sprintf(nb,"%d",counter++);
459    strcpy(name,"press");
460    strcat(name,nb);
461    TH1F* press = (TH1F*)gROOT->FindObject(name);
462    if(!press) press = new TH1F(name,"Pressure profile",50,0,50);
463    strcpy(name,"temp");
464    strcat(name,nb);
465    TH1F* temp = (TH1F*)gROOT->FindObject(name);
466    if(!temp) temp = new TH1F(name,"Temperature profile",100,0,100);
467    strcpy(name,"airdens");
468    strcat(name,nb);   
469    TH1F* airdens = (TH1F*)gROOT->FindObject(name);
470    if(!airdens) airdens = new TH1F(name,"Air Density profile",50,0,50);
471    strcpy(name,"Odens");
472    strcat(name,nb); 
473    TH1F* Odens = (TH1F*)gROOT->FindObject(name);
474    if(!Odens) Odens = new TH1F(name,"O Density profile",50,0,50);
475    strcpy(name,"O2dens");
476    strcat(name,nb);   
477    TH1F* O2dens = (TH1F*)gROOT->FindObject(name);
478    if(!O2dens) O2dens = new TH1F(name,"O2 Density profile",50,0,50);
479    strcpy(name,"O3dens");
480    strcat(name,nb);   
481    TH1F* O3dens = (TH1F*)gROOT->FindObject(name);
482    if(!O3dens) O3dens = new TH1F(name,"O3 Density profile",70,0,70);
483    strcpy(name,"N2dens");
484    strcat(name,nb);   
485    TH1F* N2dens = (TH1F*)gROOT->FindObject(name);
486    if(!N2dens) N2dens = new TH1F(name,"N2 Density profile",50,0,50);
487    strcpy(name,"H2Odens");
488    strcat(name,nb);   
489    TH1F* H2Odens = (TH1F*)gROOT->FindObject(name);
490    if(!H2Odens) H2Odens = new TH1F(name,"H2O Density profile",20,0,20);
491    strcpy(name,"CO2dens");
492    strcat(name,nb);   
493    TH1F* CO2dens = (TH1F*)gROOT->FindObject(name);
494    if(!CO2dens) CO2dens = new TH1F(name,"CO2 Density profile",50,0,50);
495*/   
496    /*    for(Int_t i=0; i<fNb; i++) {
497        Double_t alt=fAltitudeTable[i]/km;
498        press->Fill(alt,fPressureTable[i]/bar);
499        temp->Fill(alt,fTemperatureTable[i]/kelvin);
500        airdens->Fill(alt,fAir_DensityTable[i]*m3/kg);
501        Odens->Fill(alt,fO_DensityTable[i]*m3/kg);
502        O2dens->Fill(alt,fO2_DensityTable[i]*m3/kg);
503        O3dens->Fill(alt,fO3_DensityTable[i]*m3/kg);
504        N2dens->Fill(alt,fN2_DensityTable[i]*m3/kg); */
505/*
506        string val;
507
508    for(Int_t i=0; i<fNb; i++) {
509        Double_t alt =fAltitudeTable[i]/km;
510        Double_t alt1=fAltitudeTable[i+1]/km;
511        press->Fill(alt,fPressureTable[i]/bar);
512        temp->Fill(alt,fTemperatureTable[i]/kelvin);
513        airdens->Fill(alt,fAir_DensityTable[i]*m3/kg);
514        Odens->Fill(alt,fO_DensityTable[i]*m3/kg);
515        O2dens->Fill(alt,fO2_DensityTable[i]*m3/kg);
516        O3dens->Fill(alt,fO3_DensityTable[i]*m3/kg);
517        N2dens->Fill(alt,fN2_DensityTable[i]*m3/kg);
518        CO2dens->Fill(alt,fCO2_DensityTable[i]*m3/kg);
519        H2Odens->Fill(alt,fH2O_DensityTable[i]*m3/kg);
520        if (alt1-alt>1){
521          for(Double_t h=alt+1; h<alt1; h++){
522            if(h!=(Int_t)alt1) {
523          val="temp";
524          temp->Fill(h,Smooth(val,i,h)/kelvin);
525          val="press";
526          press->Fill(h,Smooth(val,i,h)/bar);
527          val="Airdens";
528          airdens->Fill(h,Smooth(val,i,h)*m3/kg);
529          val="Odens";
530          Odens->Fill(h,Smooth(val,i,h)*m3/kg);
531          val="O2dens";
532          O2dens->Fill(h,Smooth(val,i,h)*m3/kg);
533          val="O3dens";
534          O3dens->Fill(h,Smooth(val,i,h)*m3/kg);
535          val="N2dens";
536          N2dens->Fill(h,Smooth(val,i,h)*m3/kg);
537          val="CO2dens";
538          CO2dens->Fill(h,Smooth(val,i,h)*m3/kg);
539          val="H2Odens";
540          H2Odens->Fill(h,Smooth(val,i,h)*m3/kg);
541            }}
542        }
543    }
544}
545*/
546
547/*
548Double_t AtmosphereData::Smooth(string& val, Int_t i, Double_t h) const {
549   
550    Double_t val1(0), val2(0);
551    Double_t h1;
552    Double_t step;
553    string mod;
554
555    h1   = fAltitudeTable[i]/km;
556    step =(fAltitudeTable[i+1]-fAltitudeTable[i])/km;
557 
558   
559    if(val == "press"){
560       val1 = fPressureTable[i];
561       val2 = fPressureTable[i+1];}
562    else if(val == "temp"){
563       val1 = fTemperatureTable[i];
564       val2 = fTemperatureTable[i+1];}
565    else if(val == "Airdens"){
566       val1 = fAir_DensityTable[i];
567       val2 = fAir_DensityTable[i+1];}
568    else if(val == "Odens"){
569       val1 = fO_DensityTable[i];
570       val2 = fO_DensityTable[i+1];}
571    else if(val == "O2dens"){
572       val1 = fO2_DensityTable[i];
573       val2 = fO2_DensityTable[i+1];}
574    else if(val == "O3dens"){
575       val1 = fO3_DensityTable[i];
576       val2 = fO3_DensityTable[i+1];}
577    else if(val == "N2dens"){
578       val1 = fN2_DensityTable[i];
579       val2 = fN2_DensityTable[i+1];}
580    else if(val == "H2Odens"){
581       val1 = fH2O_DensityTable[i];
582       val2 = fH2O_DensityTable[i+1];}
583    else if(val == "CO2dens"){
584       val1 = fCO2_DensityTable[i];
585       val2 = fCO2_DensityTable[i+1];}
586    else
587      Msg(EsafMsg::Panic)<<"Wrong val argument in AtmosphereData::Smooth"<< MsgDispatch;
588
589    mod="expon";
590    if (val=="temp") mod="linear";
591
592    if(mod == "linear")
593        return (val2 - val1)*(h - h1)/step + val1;
594    else if(mod == "expon")
595      if (val2==0)
596        return 0;
597      else
598        return val1 * exp(-(h - h1) / (step/log(val1/val2)) );
599    else
600      Msg(EsafMsg::Panic)<< "Wrong mod argument in AtmosphereData::Smooth" << MsgDispatch;   
601    return 0;
602}
603*/
Note: See TracBrowser for help on using the repository browser.