source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/common/atmosphere/src/Atmosphere.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: 25.2 KB
Line 
1// $Id: Atmosphere.cc 3002 2011-12-06 08:09:02Z fenu $
2// S. Moreggia created 27 October 2003
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: Atmosphere                                                           *
8 *  Package: atmosphere                                                      *
9 *  Coordinator: S. Moreggia                                                 *
10 *                                                                           *
11 *****************************************************************************/
12
13//_____________________________________________________________________________
14//
15// Atmosphere
16//
17// <extensive class description>
18//
19//   Config file parameters
20//   ======================
21//
22//   <parameter name>: <parameter description>
23//   -Valid options: <available options>
24//
25
26#include "Atmosphere.hh"
27#include "AtmosphereFactory.hh"
28#include "LinsleyAtmosphere.hh"
29#include "LowtranAtmosphere.hh"
30#include "MSISE_00Atmosphere.hh"
31#include "Config.hh"
32#include "EConst.hh"
33#include "utils.hh"
34#include "NumbersFileParser.hh"
35
36#include <TF1.h>
37#include <TSpline.h>
38
39ClassImp(Atmosphere)
40
41using namespace sou;
42using namespace TMath;
43using namespace EConst;
44
45Atmosphere* Atmosphere::fChild = NULL;
46
47//______________________________________________________________________________
48Double_t densityIntegral(Double_t *x, Double_t *par)
49{
50    //
51    Double_t Altitude = par[0],  Angle = par[1];
52    Double_t H0=8.43*km; // An auxiliary constant, km
53    Double_t h = Altitude - H0*TMath::Log(x[0]);
54    Double_t f = ((EarthRadius() + Altitude)/(EarthRadius() + h))*TMath::Sin(Angle);
55    return H0*Atmosphere::Get()->Air_Density(h)/TMath::Sqrt(1.-TMath::Power(f,2))/x[0];
56}
57
58//______________________________________________________________________________
59Atmosphere::Atmosphere() : EsafConfigurable(), EsafMsgSource() {
60    //
61    // ctor
62    //
63   
64    // builds clouds and aerosol if required
65    fClouds = AtmosphereFactory::Get()->GetClouds();
66    if(!fClouds) Msg(EsafMsg::Panic) <<"No clouds built : memory problems"<< MsgDispatch;
67    fAerosol = AtmosphereFactory::Get()->GetAerosol();
68    if(!fAerosol) Msg(EsafMsg::Panic) <<"No aerosol built : memory problems"<< MsgDispatch;
69   
70    // init
71    fDepthCalculationPrecision = 1.e-5;
72    fTempInterpol = 0;
73    fDensInterpol = 0;
74    fPressInterpol = 0;
75    fMSISErandom = kFALSE;
76    fTOA = 100*km;
77   
78    // read ozone coeff.
79    fWvNBmin = 27370.;
80    fWvNBmax = 40800.;
81    NumbersFileParser nf0("config/Atmosphere/O3Crossec/O3CrossecT0",1);
82    vector<Double_t> temp0   = nf0.GetCol(0);
83    NumbersFileParser nf1("config/Atmosphere/O3Crossec/O3CrossecT1",1);
84    vector<Double_t> temp1   = nf1.GetCol(0);
85    NumbersFileParser nf2("config/Atmosphere/O3Crossec/O3CrossecT2",1);
86    vector<Double_t> temp2   = nf2.GetCol(0);
87    for(Int_t i=0; i < 2690; i++) {
88        if(i < 2687) fC0[i] = temp0[i];
89        fC1[i] = temp1[i];
90        fC2[i] = temp2[i];
91    }
92}
93
94//______________________________________________________________________________
95Atmosphere::~Atmosphere() {
96    //
97    // dtor
98    //
99    SafeDelete(fClouds);
100    SafeDelete(fAerosol);
101    SafeDelete(fTempInterpol);
102    SafeDelete(fDensInterpol);
103    SafeDelete(fPressInterpol);
104}
105
106//______________________________________________________________________________
107void Atmosphere::Delete() {
108    //
109    // staic method to call dtor
110    //
111   
112    SafeDelete(fChild);
113}
114
115//______________________________________________________________________________
116const Atmosphere* Atmosphere::Get() {
117    //
118    // Static method to get the right atmosphere singleton from anywhere in the code
119    //
120    if(fChild == 0) {
121        ConfigFileParser *pConfig = Config::Get()->GetCF("Atmosphere","AtmosphereFactory");
122        string type = pConfig->GetStr("AtmosphereFactory.Atmosphere.type");
123       
124        if(type == "linsley")
125            LinsleyAtmosphere::CreateInstance();
126       
127        else if(type == "msise00")
128            MSISE_00Atmosphere::CreateInstance();
129
130        else if(type == "RandomMsise00") {
131            MSISE_00Atmosphere::CreateInstance();
132            fChild->SetRandomMode();
133        }
134
135        else if(type == "lowtran")
136            LowtranAtmosphere::CreateInstance(); 
137
138        else cout<<"Wrong type of atmosphere\n";
139        fChild->Msg(EsafMsg::Info) << "Atmosphere built" << MsgDispatch;
140    }
141   
142     return fChild;
143}
144   
145//______________________________________________________________________________
146void Atmosphere::Reset() {
147    //
148    // Static method to reset atmosphere : useful only in MSISErandom mode
149    //
150   
151    if(fChild) {
152        if(fChild->IsRandomMode()) {
153            MSISE_00Atmosphere* atmo = (MSISE_00Atmosphere*)fChild;
154            atmo->ResetInstance();
155        }
156    }
157}
158
159//______________________________________________________________________________
160Double_t Atmosphere::GetOzoneCoeff(Double_t wl, UInt_t i) const {
161    //
162    // get Ozone C coeff. using fC tables (i=0 C0, i=1 C1, i=2 C2)
163    // Lowtran tables have been generated using a WvNB stepwidth of 5
164    //
165   
166    Double_t rtn = 0.;
167    Double_t WvNB = TranslateInWvNB(wl);
168    if(WvNB < fWvNBmin || WvNB > fWvNBmax) return 0.;
169    Int_t index = Int_t((WvNB - fWvNBmin) / 5.);
170    if(i == 0)      rtn = fC0[index];
171    else if(i == 1) rtn = fC1[index];
172    else if(i == 2) rtn = fC2[index];
173    else Msg(EsafMsg::Panic) << "<GetOzoneCoeff> Wrong argument asked" << MsgDispatch;
174   
175    return rtn;
176}
177
178//______________________________________________________________________________
179Double_t Atmosphere::Index(Double_t h, Double_t wl) const {
180    //
181    // Returns air index value at a given altitude
182    // By default, it is taken CONSTANT with wavelength (lambda=350nm)
183    // (formula drawn from Reinhard Beer, Hanbook of Optics, Chap.2)
184    //
185    Double_t rtn;
186    rtn = Index_Minus1(h,wl);
187    return 1 + rtn;
188}
189
190//______________________________________________________________________________
191Double_t Atmosphere::Index_Minus1(Double_t h, Double_t wl) const {
192    //
193    // Returns (n-1) at a given altitude
194    // By default, it is taken CONSTANT with wavelength (lambda=350nm)
195    // (formula drawn from Reinhard Beer, Hanbook of Optics, Chap.2)
196    //
197    Double_t rtn;
198    rtn = 88.43 + 185.08/(1 - 1/pow(wl/cm*1.14e5,2)) + 4.11/(1 - 1/pow(wl/cm*6.24e4,2));
199    rtn *= (Pressure(h) - WaterVaporPartialPressure(h))/Temperature(h);
200    rtn *= 296.15*kelvin/atmosphere;
201    rtn += (43.49 - 1/pow(wl/cm*1.7e4,2))*WaterVaporPartialPressure(h)/atmosphere;
202    return rtn*1e-6; 
203}
204
205//______________________________________________________________________________
206Double_t Atmosphere::Grammage(const EarthVector& V_1, const EarthVector& V_2, string opt, Double_t maxalt) const {
207    //
208    // Calculate grammage
209    // - opt = pos : between two positions V1, V2 in the atmosphere
210    // - opt = dir : from a point V1 until atmosphere top, along a track of given angles (V2 thus used for direction)
211    //
212    Double_t dl = 500*m;
213    Double_t dX, rho2, h2;
214    EarthVector U, inter;
215    EarthVector V1(V_1), V2(V_2);
216    Double_t X = 0;
217    Double_t rho1;
218    Bool_t IsTooHigh = false;
219    if(maxalt == -HUGE) maxalt = fTOA;
220   
221    if(opt == "pos") {
222        EarthVector Vmin;
223        EarthVector Vmax;
224        if(V1.IsUnderSeaLevel() || V2.IsUnderSeaLevel()) {
225            if(V1.IsUnderSeaLevel() && V2.IsUnderSeaLevel()) {
226                Msg(EsafMsg::Warning) << "<Grammage> Both V1 and V2 are under sea level --> SHOULD NOT OCCUR, zero returned"<<MsgDispatch;
227                return 0.;
228            }
229            Msg(EsafMsg::Warning) << "<Grammage> V1 or V2 is under sea level --> search for the true track impact"<<MsgDispatch;
230            EarthVector impact(1);
231            impact = ImpactASL(V1,(V2 - V1).Unit());
232            if(impact.Z() == -HUGE) {
233                impact = ImpactASL(V2,(V1 - V2).Unit());
234                if(impact.Z() == HUGE || impact.Z() == -HUGE) {
235                    Msg(EsafMsg::Warning) << "<Grammage> Failed to find an impact, zero returned"<<MsgDispatch;
236                    return 0.;
237                }
238                else V1 = impact;
239            }
240            else V2 = impact;
241        }
242       
243        if(V1.Zv() <= V2.Zv()) {
244            Vmin = V1;
245            Vmax = V2;
246        }
247        else {
248            Vmin = V2;
249            Vmax = V1;
250        }
251
252        if(Vmax.Zv() > maxalt) IsTooHigh = true;
253        U = (Vmax - Vmin).Unit();
254        inter = Vmin;
255        rho1 = Air_Density(Vmin.Zv());
256       
257        while((Vmax - Vmin).Mag() > (inter - Vmin).Mag()) {
258            if( (inter + U*dl - Vmin).Mag() >= (Vmax - Vmin).Mag() ) {
259                rho2 = Air_Density(Vmax.Zv());
260                dX = 0.5*(rho1 + rho2) * (Vmax - inter).Mag();
261                X += dX;
262                break;
263            }
264            inter += U*dl;
265            h2 = inter.Zv();
266            if(h2 > maxalt && IsTooHigh) break;
267            rho2 = Air_Density(h2);
268            dX = 0.5*(rho1 + rho2) * dl;
269            X += dX;
270            rho1 = rho2;
271        }
272    }
273       
274           
275    else if(opt == "dir") {
276        if(V1.IsUnderSeaLevel()) Msg(EsafMsg::Warning) << "<Grammage> V1 is under sea level : it SHOULD NOT"<<MsgDispatch;
277        EarthVector TOA = ImpactAtTOA(V1,V2,fTOA);
278        if(TOA.Z() == -HUGE) {
279            Msg(EsafMsg::Warning) << "<Grammage(\"dir\")> : starting position is UNDERSEALEVEL, O returned" << MsgDispatch;
280            return 0.;
281        }
282        if(TOA.Z() == HUGE) {
283            Msg(EsafMsg::Warning) << "<Grammage(\"dir\")> : track DOES NOT GO OUT atmosphere, O returned" << MsgDispatch;
284            return 0.;
285        }
286        U = V2.Unit();
287        inter = V1;
288        rho1 =  Air_Density(V1.Zv());
289        while(true) {
290            inter += U*dl;
291            h2 = inter.Zv();
292            rho2 = Air_Density(h2);
293            dX = 0.5*(rho1 + rho2) * dl;
294            X += dX;
295            rho1 = rho2;
296            if((inter - V1).Mag() > (TOA - V1).Mag()) break;
297        }
298   }
299   
300    else  Msg(EsafMsg::Warning) << "Wrong option for grammage calculation" << MsgDispatch;
301
302    return X;
303}
304
305//______________________________________________________________________________
306Double_t Atmosphere::OzoneAmountAlongPath(const EarthVector& V_1, const EarthVector& V_2, Double_t wl) const {
307    //
308    // Get ozone amount along the path, at given wavelength
309    // Integrate    X_path(s,wl) = Crossec(s,wl) * OzDens(s) * ds     along the defined path
310    //
311   
312    // init
313    Double_t dl = 500*m;
314    Double_t maxalt = 60*km;  // Ozone amount negligible above this altitude
315    Double_t dX, rho2, h2;
316    EarthVector U, inter;
317    EarthVector V1(V_1), V2(V_2);
318    Double_t X = 0;
319    Double_t rho1;
320    Bool_t IsTooHigh = false;
321   
322    // find wavelength dependent Ozone coeff.
323    Double_t C0 = GetOzoneCoeff(wl,0);
324    Double_t C1 = GetOzoneCoeff(wl,1);
325    Double_t C2 = GetOzoneCoeff(wl,2);
326   
327    EarthVector Vmin;
328    EarthVector Vmax;
329    if(V1.IsUnderSeaLevel() || V2.IsUnderSeaLevel()) {
330        if(V1.IsUnderSeaLevel() && V2.IsUnderSeaLevel()) {
331            Msg(EsafMsg::Warning) << "<OzoneAmountAlongPath> Both V1 and V2 are under sea level --> SHOULD NOT OCCUR, zero returned"<<MsgDispatch;
332            return 0.;
333        }
334        Msg(EsafMsg::Warning) << "<OzoneAmountAlongPath> V1 or V2 is under sea level --> search for the true track impact"<<MsgDispatch;
335        EarthVector impact(1);
336        impact = ImpactASL(V1,(V2 - V1).Unit());
337        if(impact.Z() == -HUGE) {
338            impact = ImpactASL(V2,(V1 - V2).Unit());
339            if(impact.Z() == HUGE || impact.Z() == -HUGE) {
340                Msg(EsafMsg::Warning) << "<OzoneAmountAlongPath> Failed to find an impact, zero returned"<<MsgDispatch;
341                return 0.;
342            }
343            else V1 = impact;
344        }
345        else V2 = impact;
346    }
347
348    if(V1.Zv() <= V2.Zv()) {
349        Vmin = V1;
350        Vmax = V2;
351    }
352    else {
353        Vmin = V2;
354        Vmax = V1;
355    }
356
357    if(Vmax.Zv() > maxalt) IsTooHigh = true;
358    U = (Vmax - Vmin).Unit();
359    inter = Vmin;
360    Double_t DeltaT = Temperature(Vmin.Zv()) - STP_Temperature();
361    rho1 = O3_DensityPPMV(Vmin.Zv()) * C0 *(1 + C1*DeltaT + C2*DeltaT*DeltaT) * Pressure(Vmin.Zv()) / Temperature(Vmin.Zv());
362
363
364    while((Vmax - Vmin).Mag() > (inter - Vmin).Mag()) {
365        if( (inter + U*dl - Vmin).Mag() >= (Vmax - Vmin).Mag() ) {
366            DeltaT = Temperature(Vmax.Zv()) - STP_Temperature();
367            rho2 = O3_DensityPPMV(Vmax.Zv()) * C0 *(1 + C1*DeltaT + C2*DeltaT*DeltaT) * Pressure(Vmax.Zv()) / Temperature(Vmax.Zv());
368            dX = 0.5*(rho1 + rho2) * (Vmax - inter).Mag()/km;
369            X += dX;
370            break;
371        }
372        inter += U*dl;
373        h2 = inter.Zv();
374        if(h2 > maxalt && IsTooHigh) break;
375        DeltaT = Temperature(h2) - STP_Temperature();
376        rho2 = O3_DensityPPMV(h2) * C0 *(1 + C1*DeltaT + C2*DeltaT*DeltaT) * Pressure(h2) / Temperature(h2);
377        dX = 0.5*(rho1 + rho2) * dl/km;
378        X += dX;
379        rho1 = rho2;
380    }
381       
382    return X*0.0269*STP_Temperature()/STP_Pressure();
383}
384
385//______________________________________________________________________________
386Int_t Atmosphere::InvertGrammage(const EarthVector& pos1, const EarthVector& direc, Double_t depth, EarthVector& rtn, Double_t TOA_alt, Double_t maxtof) const {
387    //
388    // calculate final position for given pos1,direc and air depth
389    //
390    // if 0 < TOA_alt < fTOA  --> TOA_alt is used to make present calculation (to optimize CPU time)
391    //
392    // it returns : -1 if pos1 is under sea level or if none impact (latter should not occur)
393    //               0 if position found
394    //               1 if TOF cut or infinite loop
395    //               2 if TOA reached before
396    //               3 if sea level reached before
397    //
398   
399    // init
400    Int_t status(-1);
401    if(pos1.IsUnderSeaLevel()) return status;
402    if(rtn.Mag()) rtn.SetMag(0);
403    EarthVector dir = direc.Unit();
404    EarthVector impact(1);
405    Double_t TOA(0.);
406   
407    // TOA settings
408    if( (0 < TOA_alt) && (TOA_alt < fTOA) ) TOA = TOA_alt;
409    else TOA = fTOA;
410   
411    // in case depth is null
412    if(!depth) {
413        rtn = pos1;
414        return 0;
415    }
416   
417    // check if sea level reached before foreseen air depth been travelled
418    // check if TOA reached before the foreseen air depth has been travelled
419    impact = ImpactASL(pos1,dir);
420    // if no impact at sea level
421    if(impact.Z() == HUGE) {
422        impact = ImpactAtTOA(pos1,dir,TOA);
423        if(impact.Z() == HUGE) {
424            Msg(EsafMsg::Warning) << "<InvertGrammage> track reaches NEITHER sea level NOR TOA -> SHOUD NOT HAPPEN"<<MsgDispatch;
425            return -1;
426        }
427        if(Grammage(pos1,impact) < depth) {
428            rtn = impact;
429            return 2;
430        }
431    }
432    else {
433        if(Grammage(pos1,impact) < depth) {
434            rtn = impact;
435            return 3;
436        }
437    }
438   
439   
440    // step by step process until it reaches depth value
441    Double_t dl = 500*m;
442    Double_t dX, rho2, h2;
443    EarthVector U, inter;
444    Double_t X = 0.;
445    Double_t rho1;
446    U = dir;
447    inter = pos1;
448    rho1 = Air_Density(pos1.Zv());
449    if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();
450
451    Int_t counter(0), tofcut(0);
452    if(maxtof < 0) tofcut = Int_t(2*(EarthRadius() + fTOA) / dl) + 1;
453    else tofcut = Int_t(maxtof*Clight() / dl);
454    while(true) {
455        if(counter++ > tofcut) {
456#ifdef DEBUG
457            //Msg(EsafMsg::Debug) <<"<atmoInvertGrammage> TOF cut reached : photon dumped here "<<MsgDispatch;
458#endif
459            if(rtn.Mag() == 0) rtn = pos1;
460            return 1;
461        }
462        if(counter > Int_t(2*(EarthRadius() + fTOA) / dl)) {
463            Msg(EsafMsg::Warning) <<"<atmoInvertGrammage> \"infinite\" loop (>2*EarthRadius travelled) broken by hand "<<MsgDispatch;
464            if(rtn.Mag() == 0) rtn = pos1;
465            return 1;
466        }
467        inter += U*dl;
468        h2 = inter.Zv();
469        rho2 = Air_Density(h2);
470        dX = 0.5*(rho1 + rho2) * dl;
471        if((X + dX) >= depth) break;
472        X += dX;
473        rho1 = rho2;
474        if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();
475    }
476    // come back to last iteration for fine tuning
477    rtn = inter - U*dl;
478    Double_t slope = (rho2 - rho1) / dl;  // density profile considered linear locally
479    if(slope) dl = sign(slope) * sqrt( (rho1/slope) * (rho1/slope) + 2.*(depth - X) / slope) - rho1/slope;
480    else dl = (depth - X) / rho1;  // if constant density
481    rtn += U*dl;
482
483    if(rtn.Mag() == 0) rtn = pos1;
484   
485    return 0;
486}
487
488//______________________________________________________________________________
489Int_t Atmosphere::InvertOzoneAmount(const EarthVector& pos1, const EarthVector& direc, Double_t OzAmount, Double_t wl, EarthVector& rtn, Double_t TOA_alt, Double_t maxtof) const {
490    //
491    // calculate position corresponding to ozone amount, using given pos1,direc
492    //
493    // if 0 < TOA_alt < fTOA  --> TOA_alt is used to make present calculation (to optimize CPU time)
494    //
495    // it returns : -1 if pos1 is under sea level or if none impact (latter should not occur)
496    //               0 if position found
497    //               1 if TOF cut or infinite loop
498    //               2 if TOA reached before
499    //               3 if sea level reached before
500    //
501   
502    // init
503    Int_t status(-1);
504    if(pos1.IsUnderSeaLevel()) return status;
505    if(rtn.Mag()) rtn.SetMag(0);
506    EarthVector dir = direc.Unit();
507    EarthVector impact(1);
508    Double_t TOA(0.);
509   
510    // TOA settings
511    if( (0 < TOA_alt) && (TOA_alt < fTOA) ) TOA = TOA_alt;
512    else TOA = fTOA;
513   
514    // in case ozone amount is null
515    if(OzAmount == 0.) {
516        rtn = pos1;
517        return 0;
518    }
519
520    // find wavelength dependent Ozone coeff.
521    Double_t C0 = GetOzoneCoeff(wl,0);
522    Double_t C1 = GetOzoneCoeff(wl,1);
523    Double_t C2 = GetOzoneCoeff(wl,2);
524   
525   
526    // check if sea level reached before foreseen ozone amount been traveled
527    // check if TOA reached before the foreseen ozone amount has been traveled
528    impact = ImpactASL(pos1,dir);
529    // if no impact at sea level
530    if(impact.Z() == HUGE) {
531        impact = ImpactAtTOA(pos1,dir,TOA);
532        if(impact.Z() == HUGE) {
533            Msg(EsafMsg::Warning) << "<InvertOzoneAmount> track reaches NEITHER sea level NOR TOA -> SHOUD NOT HAPPEN"<<MsgDispatch;
534            return -1;
535        }
536        if(OzoneAmountAlongPath(pos1,impact,wl) < OzAmount) {
537            rtn = impact;
538            return 2;
539        }
540    }
541    else {
542        if(OzoneAmountAlongPath(pos1,impact,wl) < OzAmount) {
543            rtn = impact;
544            return 3;
545        }
546    }
547   
548   
549    // step by step process until it reaches OzAmount value
550    Double_t dl = 500*m;
551    Double_t dX, rho2, h2;
552    EarthVector U, inter;
553    Double_t X = 0.;
554    Double_t rho1;
555    U = dir;
556    inter = pos1;
557    Double_t DeltaT = Temperature(pos1.Zv()) - STP_Temperature();
558    rho1 = O3_DensityPPMV(pos1.Zv()) * C0 *(1 + C1*DeltaT + C2*DeltaT*DeltaT) * Pressure(pos1.Zv()) / Temperature(pos1.Zv());
559    if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();
560
561    Int_t counter(0), tofcut(0);
562    Int_t stopvalue = Int_t(2*(EarthRadius() + TOA) / dl);
563    if(maxtof < 0) tofcut = stopvalue + 1;
564    else tofcut = Int_t(maxtof*Clight() / dl);
565    while(true) {
566        if(counter++ > tofcut) {
567#ifdef DEBUG
568            //Msg(EsafMsg::Debug) <<"<InvertOzoneAmount> TOF cut reached : photon dumped here "<<MsgDispatch;
569#endif
570            if(rtn.Mag() == 0) rtn = pos1;
571            return 1;
572        }
573        if(counter > stopvalue) {
574            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> \"infinite\" loop (>2*EarthRadius travelled) broken by hand"<<MsgDispatch;
575            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> counter = "<<counter <<MsgDispatch;
576            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> Stop value =  "<<stopvalue <<MsgDispatch;
577            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> maxtof =  "<<maxtof <<MsgDispatch;
578            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> MCRT TOA = " << TOA <<MsgDispatch;
579            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> initpos = " << pos1 <<MsgDispatch;
580            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> initpos.Zv = " << pos1.Zv() <<MsgDispatch;
581            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> finalpos = " << EarthVector(pos1+inter) <<MsgDispatch;
582            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> finalpos.Zv = " << EarthVector(pos1+inter).Zv() <<MsgDispatch;
583            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> direc = " << direc <<MsgDispatch;
584            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> direc.theta = " << direc.Theta()*RadToDeg() <<MsgDispatch;
585            Msg(EsafMsg::Warning) <<"<InvertOzoneAmount> direc.phi = " << direc.Phi()*RadToDeg() <<MsgDispatch;
586            if(rtn.Mag() == 0) rtn = pos1;
587            return 1;
588        }
589        inter += U*dl;
590        h2 = inter.Zv();
591        DeltaT = Temperature(h2) - STP_Temperature();
592        rho2 = O3_DensityPPMV(h2) * C0 *(1 + C1*DeltaT + C2*DeltaT*DeltaT) * Pressure(h2) / Temperature(h2);
593        dX = 0.5*(rho1 + rho2) * dl/km * 0.0269*STP_Temperature()/STP_Pressure();
594        if((X + dX) >= OzAmount) break;
595        X += dX;
596        rho1 = rho2;
597        if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();
598    }
599
600    // come back to last iteration for fine tuning
601    rtn = inter - U*dl;
602    Double_t slope = (rho2 - rho1) / dl;  // ozone amount profile considered linear locally
603    if(slope) dl = sign(slope) * sqrt( (rho1/slope) * (rho1/slope) + 2.*(OzAmount - X) / (0.0269/km*STP_Temperature()/STP_Pressure()) / slope) - rho1/slope;
604    else dl = (OzAmount - X) / (0.0269/km*STP_Temperature()/STP_Pressure()) / rho1;  // if locally constant ozone amount profile
605    rtn += U*dl;
606
607
608    if(rtn.Mag() == 0) rtn = pos1;
609   
610    return 0;
611}
612
613//______________________________________________________________________________
614EarthVector Atmosphere::ImpactASL(const EarthVector& pos, const EarthVector& dir) const {
615    //
616    // returns impact at sea level of a track defined by starting position and direction
617    //
618
619    // if pos is under sea level
620    EarthVector rtn;
621    if(pos.IsUnderSeaLevel()) {
622        rtn.SetXYZ(0,0,-HUGE);
623        return rtn;
624    }
625
626
627    // to know if dir is locally going upward
628    Double_t angle;
629    EarthVector temp(pos);
630    temp(2) += EarthRadius();  // pos expressed in earth-centered frame -> gives local vertical direction expressed in MES frame
631    angle = dir.Angle(temp);   // angle between dir and local vertical
632    if(angle < (PiOver2() - kTolerance)) {
633        rtn.SetXYZ(0,0,HUGE);
634        return rtn;
635    }
636
637    // now, impact calculation
638    Double_t mag(0);
639    EarthVector direc = dir.Unit();
640
641/*
642//  flat earth
643if(fabs(direc.Z()) < kTolerance) rtn.SetXYZ(0,0,HUGE);
644else {
645    mag = -pos.Z() / direc.Z();
646    if(mag < -kAltitudeTolerance) rtn.SetXYZ(0,0,HUGE);
647    else rtn = pos + mag*direc;
648}
649*/
650    // spherical earth
651    Double_t b = pos*direc + direc.Z()*EarthRadius();
652    Double_t c = pos.Mag2() + 2*EarthRadius()*pos.Z();
653    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
654    if(p.first == 0) {
655        rtn.SetXYZ(0,0,HUGE);
656        return rtn;
657    }
658    if(p.first == 1) mag = p.second[0];
659    else if(p.first ==2) mag = min(p.second[0],p.second[1]);
660   
661    if(mag < -kAltitudeTolerance) Msg(EsafMsg::Debug) <<"SHOULD NOT HAVE IMPACT AT SEA LEVEL, mag = "<<mag << MsgDispatch;
662   
663    rtn = pos + mag*direc;
664    return rtn;   
665}
666
667//______________________________________________________________________________
668EarthVector Atmosphere::ImpactAtTOA(const EarthVector& pos, const EarthVector& dir, Double_t TOA_alt) const {
669    //
670    // returns impact at Top Of Atmosphere of a track defined by starting position and direction
671    // if 0 < TOA_alt < fTOA  --> TOA_alt is used to make present calculation (to optimize CPU time, cf. InvertGrammage())
672    //
673
674    // if pos is under sea level
675    EarthVector rtn;
676    if(pos.IsUnderSeaLevel()) {
677        rtn.SetXYZ(0,0,-HUGE);
678        return rtn;
679    }
680   
681    // TOA settings
682    Double_t altTOA = 0.;
683    if( (0 < TOA_alt) && (TOA_alt < fTOA) ) altTOA = TOA_alt;
684    else altTOA = fTOA;
685
686    // now, impact calculation
687    Double_t mag(0);
688    EarthVector direc = dir.Unit();
689
690/*
691//  flat earth
692if(fabs(direc.Z()) < kTolerance) rtn.SetXYZ(0,0,HUGE);
693else {
694    mag = (altTOA - pos.Z()) / direc.Z();
695    if(mag < -kAltitudeTolerance) rtn.SetXYZ(0,0,HUGE);
696    else rtn = pos + mag*direc;
697}
698*/
699
700    // spherical earth
701    Double_t b = pos*direc + direc.Z()*EarthRadius();
702    Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - altTOA) - pow(altTOA,2);
703    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
704    if(p.first == 0) {
705        rtn.SetXYZ(0,0,HUGE);
706        return rtn;
707    }
708    if(p.first == 1) mag = p.second[0];
709    else if(p.first == 2) {
710        if((p.second[0]+kAltitudeTolerance) * (p.second[1]+kAltitudeTolerance) < 0) mag = max(p.second[0],p.second[1]);
711        else mag = min(p.second[0],p.second[1]);
712    }
713    if(mag < -kAltitudeTolerance) {
714        Msg(EsafMsg::Warning) << "PB in impact calculation ImpactAtTOA(), mag[min,max] = ["<<mag<<", "<< max(p.second[0],p.second[1])<<"]" << MsgDispatch;
715        Msg(EsafMsg::Warning) << "pos = "<<pos << MsgDispatch;
716        Msg(EsafMsg::Warning) << "dir = "<<direc << MsgDispatch;
717    }
718    rtn = pos + mag*direc;
719    return rtn;   
720}
721
722
723
724// used in reco only : obsolete, should use Grammage method
725//______________________________________________________________________________
726Double_t Atmosphere::Depth(const Double_t H, const Double_t Theta) const {
727    // Compute the atmosphere depth between the given point with (h,theta) coordinates and infinity
728    //
729    Double_t eps=1.e-5;
730    TF1 *depthIntegral = new TF1("depthIntegral",densityIntegral,0.,1.,2);
731   
732    if (Theta <= TMath::PiOver2()) {
733        Double_t pars[2] = {H,Theta};
734        Double_t depth = depthIntegral->Integral(eps,1-eps,pars,fDepthCalculationPrecision);
735        delete depthIntegral;
736        return depth;
737    }
738    else { 
739        Double_t Hstar = ( EarthRadius() + H)*TMath::Sin(TMath::Pi() - Theta ) - 
740            EarthRadius();
741        Double_t pars1[2] = {Hstar,TMath::PiOver2()};
742        Double_t pars2[2] = {H,TMath::Pi() - Theta};
743        Double_t depth=2*depthIntegral->Integral(eps,1-eps,pars1) - depthIntegral->Integral(eps,1-eps,pars2,
744            fDepthCalculationPrecision);
745        delete depthIntegral;
746        return depth;
747    }
748}
749
Note: See TracBrowser for help on using the repository browser.