source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/radiativetransfer/src/TestGround.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: 18.4 KB
Line 
1// ESAF : Euso Simulation and Analysis Framework
2// $Id: TestGround.cc 2733 2006-06-09 12:43:12Z moreggia $
3// S. Moreggia created 27 October 2003
4
5#include "euso.hh"
6#include "TestGround.hh"
7#include "utils.hh"
8#include <stdlib.h>
9#include "EConst.hh"
10#include "EsafRandom.hh"
11#include "SinglePhoton.hh"
12#include "ListPhotonsInAtmosphere.hh"
13#include <TF2.h>
14
15ClassImp(TestGround)
16
17using EConst::EarthRadius;
18using namespace TMath;
19using namespace sou;
20
21
22//______________________________________________________________________________
23Double_t Specular_phase_function(Double_t *x, Double_t *par) {
24    //
25    // for specular reflexion mode
26    // calculate value of 2D phase function using theta-phi
27    //
28    // VECTORS ARE TAKEN LOCALLY, i.e. AT GROUND IMPACT
29    // so that angles are define w.r.t. local vertical
30    //
31    // x[0]    =  theta out
32    // x[1]    =  phi   out
33    // par[0]  =  theta incident  (incident direction has been reversed !!)
34    // par[1]  =  phi   incident  (incident direction has been reversed !!)
35    // par[2]  =  gaussian sigma parameter
36    //
37    //
38   
39    // 1. init
40    EarthVector incomdir, outgodir, Uz;
41    incomdir.SetMagThetaPhi(1.,par[0],par[1]);
42    outgodir.SetMagThetaPhi(1.,x[0],x[1]);
43    Uz.SetXYZ(0,0,1.);
44   
45    // 1. define change of frame : outgodir expressed in the frame where incomdir is the Z-axis (incomdir has been reversed before)
46    EarthVector rotaxis_spec = Uz.Cross(incomdir);
47    Double_t rotangle_spec   = -Uz.Angle(incomdir);
48   
49    // 2. move outgoing direction into new frame and get the out angle
50    outgodir.Rotate(rotangle_spec,rotaxis_spec);
51    Double_t th_out_spec  = outgodir.Theta();
52   
53    // 3. if th_out_spec > 90deg, zero returned
54    if(th_out_spec > PiOver2()) return 0.;
55       
56    return 1. / (Sqrt(TwoPi()) * par[2]) * Exp(-th_out_spec*th_out_spec / (2*par[2]*par[2]));
57}
58
59//_____________________________________________________________________________________________________
60TestGround::TestGround() {
61    //
62    // ctor
63    //
64   
65    Specular_distrib = 0;
66    Configure();
67    Msg(EsafMsg::Info) << "*** TestGround built ***" << MsgDispatch;
68}
69
70//_____________________________________________________________________________________________________
71TestGround::~TestGround(){
72    //
73    // dtor
74    //
75   
76    if(Specular_distrib) SafeDelete(Specular_distrib);
77}
78
79//_____________________________________________________________________________________________________
80void TestGround::Configure() {
81    //
82    // Configure the object
83    //
84   
85    fType = Conf()->GetStr("TestGround.fType");
86    fSpec = Conf()->GetNum("TestGround.fSpec");
87    fSigma = Conf()->GetNum("TestGround.fSigma")*DegToRad();
88    fWindSpeed = Conf()->GetNum("TestGround.fWindSpeed");
89    fAlbedo  = Conf()->GetNum("TestGround.fAlbedo");
90    fAltitude = Conf()->GetNum("TestGround.fAltitude")*km;
91    Msg(EsafMsg::Info) <<"Ground with Albedo = " << fAlbedo << MsgDispatch;
92   
93    // if specular component exists
94    if(fType == "mixed") {
95        if(fSigma < 0) {
96            fSigma = ((fWindSpeed - 2) + 6)*DegToRad();
97            fSpec  = 0.75;
98        }
99        Msg(EsafMsg::Info) <<"Ground with mixed BRDF, specular component = " << fSpec <<", Gaussian width = "<<fSigma*RadToDeg()<< "deg" << MsgDispatch;
100        Specular_distrib = new TF2("Specular_distrib",Specular_phase_function,-PiOver2(),PiOver2(),-Pi(),Pi(),3);
101       
102        // create tables for specular component normalization
103        Msg(EsafMsg::Info) <<"Build specular component normalization tables" << MsgDispatch;
104        NORMmax = 0;
105        for(Int_t k=0; k< 91; k++) {
106            Specular_distrib->SetParameters(k*DegToRad(),0,fSigma);
107            // calculate normalization factor to achieve a self-normalized specular intensity 2D-distribution
108            // use trapezoidal integration (in theta and phi)
109            Double_t norm = 0.;
110            Double_t Integ_phi1(0.), Integ_phi2(0.);
111            Int_t nbsteps = 100;  // allow precision better than 0.001
112            Double_t ttheta(0.), pphi(0.), f1(0.), f2(0.), dOmega1(0.), dOmega2(0.);
113            Double_t dth = PiOver2() / Double_t(nbsteps);
114            Double_t dph = TwoPi() / Double_t(nbsteps);
115
116            // first value for Integ_phi
117            ttheta = 0.;
118            pphi = -Pi();
119            f2 = Specular_distrib->Eval(ttheta,pphi);
120            dOmega2 = 0.;     // zero cause Sin(th=0) term
121            Integ_phi2 = 0.;  // zero cause Sin(th=0) term
122            // loop on theta
123            for(Int_t i=0; i < nbsteps; i++) {
124                ttheta += dth;
125                pphi = -Pi();
126                // first value for f at this theta
127                f2 = Specular_distrib->Eval(ttheta,pphi);
128                Integ_phi1 = Integ_phi2;
129                Integ_phi2 = 0.;
130                // loop on phi
131                for(Int_t j=0; j< nbsteps; j++) {
132                    f1 = f2;
133                    pphi += dph;
134                    f2 = Specular_distrib->Eval(ttheta,pphi);
135                    Integ_phi2 += 0.5 * (f1 + f2);
136                }
137                dOmega1 = dOmega2;
138                dOmega2 = Sin(ttheta) * dth * dph;
139                norm += 0.5 * (Integ_phi1 * dOmega1 + Integ_phi2 * dOmega2);
140            }
141            NORM[k] = norm;
142            if(norm > NORMmax) NORMmax = norm;
143        }
144       
145        Msg(EsafMsg::Info) <<"Tables created" << MsgDispatch;
146    }
147}
148
149//_____________________________________________________________________________________________________
150void TestGround::RandomDir(const EarthVector& pos, EarthVector& incom_dir) const {
151    //
152    // set incom_dir to the outgoing direction of scattering, sampling Outgoing phase function
153    // - pos       = position on ground (MES frame)
154    // - incom_dir = incident direction (MES frame)  :  FINALLY SET TO THE OUTGOING DIRECTION
155    //
156    // Phase function is normalized when integrated over an hemisphere, i.e. theta within [0,Pi()/2]   (SOLID ANGLE integration !!)
157    // Specular pattern is assumed cirular in shape (i.e. independent of phi angle around incomdir_local)
158    // Specular pattern is coded as a gaussian for INTENSITY (NOT for radiance)
159    //
160   
161    // init
162    TRandom* rndm = EsafRandom::Get();
163    Double_t theta_out(0.), th_in(0.), phi_out(0.), phi_in(0.);
164    EarthVector incomdir_local = -incom_dir.Unit();   // WARNING : local incomdir is REVERSED !
165    EarthVector outgodir_local(1,0,0);
166    const EarthVector Uz(0,0,1);
167   
168    // 0. do some checks
169        // 0.1. in weird cases (should not occur)
170    if(incom_dir.Mag() == 0) {
171        Msg(EsafMsg::Warning) << "<RandomDir> NO INCOMING dir, it SHOULD NOT, direction is not changed" << MsgDispatch;
172        return;
173    }
174        // 0.2. check if pos is at ground level
175    if((pos.Zv() - fAltitude) > kAltitudeTolerance) {
176        Msg(EsafMsg::Warning) << "<RandomDir> Ask for BRDF but position is NOT at ground level, delta_h = "<<pos.Zv() - fAltitude<<", direction is not changed" << MsgDispatch;
177        return;
178    }
179
180   
181    // 1. move to local frame
182    //    rotation axis  defined by       Uz x earthlocalaxis(pos)       (expressed in MES)
183    //    rotation angle defined by       Angle(Uz,earthlocalaxis)       (expressed in MES)
184    //
185        // 1.0. check the very specific case where pos is (0,0,0)
186        //      if it is the case, does not apply change of frame
187    EarthVector rotaxis(1,0,0);
188    Double_t rotangle = 0.;
189    if(pos.Perp() > kAltitudeTolerance) {
190        // 1.1. init
191        EarthVector earthlocalaxis(pos);
192        earthlocalaxis(2) +=  EarthRadius();
193        earthlocalaxis.Unit();
194        // 1.2. define rotation
195        rotaxis = Uz.Cross(earthlocalaxis);
196        rotangle   = Uz.Angle(earthlocalaxis);
197        if(rotaxis.Mag() < kTolerance) {      // rotaxis should never be null as pos is not (0,0,0)
198            Msg(EsafMsg::Warning) << "<RandomDir> Rotation axis is NULL (never should!), direction is not changed" << MsgDispatch;
199            return;
200        }
201        // 1.3. move direction into local frame
202        incomdir_local.Rotate(rotangle,rotaxis);
203
204        // 1.4. Define local angles (local incomdir is REVERSED !)
205        th_in   = incomdir_local.Theta();
206        phi_in  = incomdir_local.Phi();
207    }
208    else {
209        th_in   = incomdir_local.Theta();
210        phi_in  = incomdir_local.Phi();
211    }
212
213
214    // 2. checks if incoming local direction is locally going downward (WARNING : local incomdir is REVERSED !)
215    if(th_in > PiOver2() + kTolerance) {
216        Msg(EsafMsg::Warning) << "\n<RandomDir> theta_local/halfpi = "<<th_in/PiOver2()<<", theta_local = "<<th_in*RadToDeg()<<" deg" << MsgDispatch;
217        Msg(EsafMsg::Warning) << "<RandomDir> pos = "<<pos << MsgDispatch;
218        Msg(EsafMsg::Warning) << "<RandomDir> Incoming local direction is going upward (can occur if theta close to 90deg, round-off issue), direction put at 90deg (locally horizontal)" << MsgDispatch;
219        th_in = PiOver2();
220    }
221
222
223    // 3. Get an outgoing direction in sampling Outgoing_phase_function
224    if(fType == "isotropic") {  // theta between [0,pi/2] only
225        theta_out = acos(rndm->Rndm());        // dP/dcos(theta) = 2pi * dP/dphi.dcos(theta) --- dP/dphi.dcos(theta) = 1/2pi, si cos(theta) = x -->  dProba/dx = 1
226        phi_out   = rndm->Rndm()*TwoPi();
227    }
228   
229    else if(fType == "lambertian") {
230        theta_out = acos(sqrt(rndm->Rndm())); // dP/dcos(theta) = 2pi * dP/dphi.dcos(theta)  --- dP/dphi.dcos(theta) = cos(theta)/pi, si cos(theta) = x -->  dProba/dx = 2x
231        phi_out   = rndm->Rndm()*TwoPi();
232    }
233   
234    else if(fType == "mixed") {
235        // photon is reflected diffusely
236        if(fSpec <  rndm->Rndm()) {  // theta between [0,pi/2] only
237            theta_out = acos(sqrt(rndm->Rndm())); // dP/dcos(theta) = 2pi * dP/dphi.dcos(theta)  --- dP/dphi.dcos(theta) = cos(theta)/pi, si cos(theta) = x -->  dProba/dx = 2x
238            phi_out   = rndm->Rndm()*TwoPi();
239        }
240        else {
241            // photon is reflected by specular process
242            // cos(theta) and phi are two independent variables
243            Double_t pfvalue = 0.;
244            Double_t u1(0.), u2(0.);
245            do {
246                //  3.1. specular component value sample uniformly between 0 and max
247                u2 = rndm->Rndm() / (Sqrt(TwoPi()) * fSigma);
248
249                // 3.2. specular component at ASin(u1)
250                    // 3.2.1. init function using incoming angles info
251                Specular_distrib->SetParameters(th_in,phi_in,fSigma);
252                    // 3.2.3. set direction to get specular distribution needed value
253                u1 = rndm->Rndm();  // theta between [0,pi/2] only
254                theta_out = ACos(u1);
255                phi_out = TwoPi()*rndm->Rndm();
256                outgodir_local.SetMagThetaPhi(1.,theta_out,phi_out);
257                pfvalue += Specular_distrib->Eval(theta_out,phi_out);
258            } while(u2 >= pfvalue);
259        }
260    }
261   
262    else Msg(EsafMsg::Panic) << "<RandomDir> Unvalid type of brdf = "<< fType << MsgDispatch;
263   
264    // 4. returns in the MES
265        // check if outgoing direction is upward
266    if(theta_out > PiOver2() + kTolerance) {
267        Msg(EsafMsg::Warning) << "<RandomDir> Outgoing local direction is going downward (never should!),  direction is not changed" << MsgDispatch;
268        return;
269    }
270    outgodir_local.SetMagThetaPhi(1.,theta_out,phi_out);
271    if(pos.Perp() > kAltitudeTolerance) outgodir_local.Rotate(-rotangle,rotaxis);
272    incom_dir = outgodir_local; // outgoing direction
273}
274
275//_____________________________________________________________________________________________________
276Double_t TestGround::Outgoing_phase_function(const EarthVector& pos, const EarthVector& incom, const EarthVector& outgo) const {
277    //
278    // Phase function of the outgoing light (earth sphericity taken into account)
279    // - pos      = position on ground (MES frame)
280    // - incomdir = incident direction (MES frame)
281    // - outgodir = outgoing position  (MES frame)
282    // Phase function is normalized when integrated over an hemisphere, i.e. theta within [0,Pi()/2]   (SOLID ANGLE integration !!)
283    // Specular pattern is assumed cirular in shape (i.e. independent of phi angle around incomdir_local)
284    // Specular pattern is coded as a gaussian for INTENSITY (NOT radiance)
285    //
286   
287    // init
288    Double_t rtn(0);
289    Double_t th_out(0.), th_in(0.), phi_out(0.), phi_in(0.);
290    EarthVector incomdir_local = -incom.Unit();   // WARNING : local incomdir is REVERSED !
291    EarthVector outgodir_local = outgo.Unit();
292    const EarthVector Uz(0,0,1);
293   
294    // 0. do some checks
295        // 0.1. in weird case incomdir comes from a fluo bunchofphotons (dir.mag=0)
296    if(incom.Mag() == 0) {
297        Msg(EsafMsg::Warning) << "<Outgoing_phase_function> NO INCOMING dir, it SHOULD NOT, ZERO RETURNED" << MsgDispatch;
298        return 0.;
299    }
300        // 0.2. check if pos is at ground level
301    if((pos.Zv() - fAltitude) > kAltitudeTolerance) {
302        Msg(EsafMsg::Warning) << "<Outgoing_phase_function> Ask for BRDF but position is NOT at ground level, delta_h = "<<pos.Zv() - fAltitude<<"  ZERO RETURNED" << MsgDispatch;
303        return 0.;
304    }
305
306   
307    // 1. move to local frame
308    //    rotation axis  defined by       Uz x earthlocalaxis(pos)       (expressed in MES)
309    //    rotation angle defined by       Angle(Uz,earthlocalaxis)       (expressed in MES)
310    //
311        // 1.0. check the very specific case where pos is (0,0,0)
312        //      if it is the case, does not apply change of frame
313    if(pos.Perp() > kAltitudeTolerance) {
314        // 1.1. init
315        EarthVector earthlocalaxis(pos);
316        earthlocalaxis(2) +=  EarthRadius();
317        earthlocalaxis.Unit();
318        // 1.2. define rotation
319        EarthVector rotaxis = Uz.Cross(earthlocalaxis);
320        Double_t rotangle   = Uz.Angle(earthlocalaxis);
321        if(rotaxis.Mag() < kTolerance) {      // rotaxis should never be null as pos is not (0,0,0)
322            Msg(EsafMsg::Warning) << "<Outgoing_phase_function> Rotation axis is NULL (never should!), ZERO RETURNED" << MsgDispatch;
323            return 0.;
324        }
325        // 1.3. move directions into local frame
326        incomdir_local.Rotate(rotangle,rotaxis);
327        outgodir_local.Rotate(rotangle,rotaxis);
328
329        // 1.4. Define local angles (local incomdir is REVERSED !)
330        th_out  = outgodir_local.Theta();
331        th_in   = incomdir_local.Theta();
332        phi_out = outgodir_local.Phi();
333        phi_in  = incomdir_local.Phi();
334    }
335    else {
336        th_out  = outgodir_local.Theta();
337        th_in   = incomdir_local.Theta();
338        phi_out = outgodir_local.Phi();
339        phi_in  = incomdir_local.Phi();
340    }
341
342    // 2. checks if :
343    //    - outgoing local direction is locally going upward
344    //    - incoming local direction is locally going downward (WARNING : local incomdir is REVERSED !)
345    if(th_in > PiOver2() + kTolerance) {
346        Msg(EsafMsg::Warning) << "\n<Outgoing_phase_function> theta_local/halfpi = "<<th_in/PiOver2()<<", theta_local = "<<th_in*RadToDeg()<<" deg" << MsgDispatch;
347        Msg(EsafMsg::Warning) << "<Outgoing_phase_function> pos = "<<pos << MsgDispatch;
348        Msg(EsafMsg::Warning) << "<Outgoing_phase_function> Incoming local direction is going upward (can occur if theta close to 90deg, round-off issue), direction put at 90deg (locally horizontal)" << MsgDispatch;
349        th_in = 89.9*DegToRad();
350    }
351    if(th_out > PiOver2()) return 0.;
352
353
354    // 3. Calculate value of reflected intensity (= outgoing phase function)
355    if(fType == "isotropic")
356        rtn = 1./TwoPi();
357   
358    else if(fType == "lambertian")
359        rtn = Cos(th_out)/Pi();
360   
361    else if(fType == "mixed") {
362        // diffuse intensity component (self-normalized)
363        rtn = (1 - fSpec) / Pi() * Cos(th_out);
364        // specular radiance component
365        Specular_distrib->SetParameters(th_in,phi_in,fSigma); // init function using incoming angles info
366        Int_t index = Int_t(th_in*RadToDeg());  // to retrieve factor that achieves a normalized 2D-distribution
367        if(index > 90) {
368            Msg(EsafMsg::Warning) << "<Outgoing_phase_function, norm factor> Int_t(theta) is too high, th = "<<th_in*RadToDeg()<<". Set to 90deg" << MsgDispatch;
369            index = 90;
370        }
371        Double_t norm = (NORM[index+1] - NORM[index]) * (th_in*RadToDeg() - Double_t(index)) + NORM[index];  // linear interpolation
372        rtn += fSpec * Specular_distrib->Eval(th_out,phi_out) / norm;
373    }
374   
375    else Msg(EsafMsg::Panic) << "<Outgoing_phase_function> Unvalid type of brdf = "<< fType << MsgDispatch;
376
377    return rtn;
378}
379
380//______________________________________________________________________________
381Double_t TestGround::GetMaxOutgoing_phase_function() const {
382    //
383    // maximum value of outgoing phase function
384    // NB : Phase function is normalized when integrated over the upper hemisphere, i.e. outoing theta direction within [0,Pi()/2]
385    //
386
387    Double_t rtn = 0.;
388   
389    if(fType == "isotropic")
390        rtn = 1./TwoPi();
391   
392    else if(fType == "lambertian")
393        rtn = 1./Pi();
394   
395    else if(fType == "mixed")
396        rtn = (1 - fSpec)/Pi() + fSpec / NORMmax / (Sqrt(TwoPi()) * fSigma);
397   
398    else Msg(EsafMsg::Panic) << "<Outgoing_phase_function> Unvalid type of brdf = "<< fType << MsgDispatch;
399       
400    return rtn;
401}
402
403//_____________________________________________________________________________________________________
404EarthVector TestGround::GetImpact(const EarthVector& pos, const EarthVector& dir) const {
405    //
406    // Impact on ground form a given (position,direction) pair in atmosphere
407    //
408
409    // if pos is underground
410    EarthVector rtn;
411    if((pos.Zv() + kAltitudeTolerance) < fAltitude) {
412        Msg(EsafMsg::Debug) << "Starting position is underground, so no ground impact calculated" << MsgDispatch;
413        rtn.SetXYZ(0,0,-HUGE);
414        return rtn;
415    }
416
417#ifdef DEBUG
418   // Msg(EsafMsg::Debug) <<"Impact on ground = " << MsgDispatch;
419#endif
420
421    // to know if dir is locally going upward
422    Double_t angle;
423    EarthVector temp(pos);
424    temp(2) += EarthRadius();  // pos expressed in earth-centered frame -> gives local vertical direction expressed in MES frame
425    angle = dir.Angle(temp);   // angle between dir and local vertical
426    if(angle < (PiOver2() - kTolerance)) {
427      //Msg(EsafMsg::Debug) << MsgDispatch;
428        //Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND, track is locally going upward" << MsgDispatch;
429        rtn.SetXYZ(0,0,HUGE);
430        return rtn;
431    }
432
433    // now impact calculation
434    Double_t mag(0);
435    EarthVector direc = dir.Unit();
436
437/*  flat earth
438if(fabs(direc.Z()) < kTolerance) rtn.SetXYZ(0,0,HUGE);
439else {
440    mag = -pos.Z() / direc.Z();
441    if(mag < -kAltitudeTolerance) rtn.SetXYZ(0,0,HUGE);
442    else rtn = pos + mag*direc;
443}
444*/
445
446    // spherical earth
447    Double_t b = pos*direc + direc.Z()*EarthRadius();
448    Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - fAltitude) - fAltitude*fAltitude;
449    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
450    if(p.first == 0) {
451#ifdef DEBUG
452        //Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND" << MsgDispatch;
453#endif
454        rtn.SetXYZ(0,0,HUGE);
455        return rtn;
456    }
457    if(p.first == 1) mag = p.second[0];
458    else if(p.first ==2) mag = min(p.second[0],p.second[1]);
459   
460    if(mag < -kAltitudeTolerance) Msg(EsafMsg::Debug) <<"SHOULD NOT HAVE IMPACT ON GROUND, mag = "<<mag << MsgDispatch;
461   
462    rtn = pos + mag*direc;
463
464#ifdef DEBUG
465    //Msg(EsafMsg::Debug) << rtn << MsgDispatch;
466#endif
467
468
469    return rtn;
470}
471
472
Note: See TracBrowser for help on using the repository browser.