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 | |
---|
16 | ClassImp(AtmosphereData) |
---|
17 | |
---|
18 | using EConst::Avogadro; |
---|
19 | using EConst::Loschmidt; |
---|
20 | using EConst::R_ideal; |
---|
21 | using EConst::AtomicMass; |
---|
22 | using EConst::STP_Pressure; |
---|
23 | using EConst::STP_Temperature; |
---|
24 | using namespace sou; |
---|
25 | |
---|
26 | //extern void InitLOWTRAN(); |
---|
27 | |
---|
28 | //______________________________________________________________________________ |
---|
29 | AtmosphereData::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 | //______________________________________________________________________________ |
---|
48 | AtmosphereData::~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 | //______________________________________________________________________________ |
---|
80 | void 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 | //______________________________________________________________________________ |
---|
269 | void 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 | //_____________________________________________________________________________ |
---|
422 | Double_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 |
---|
453 | void 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 | /* |
---|
548 | Double_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 | */ |
---|