// ESAF : Euso Simulation and Analysis Framework
// $Id: TestGround.cc,v 1.35 2005/10/06 08:09:55 moreggia Exp $
// S. Moreggia created 27 October 2003
#include "euso.hh"
#include "TestGround.hh"
#include "utils.hh"
#include <stdlib.h>
#include "EConst.hh"
#include "EsafRandom.hh"
#include "SinglePhoton.hh"
#include "ListPhotonsInAtmosphere.hh"
#include <TF1.h>
ClassImp(TestGround)
using EConst::EarthRadius;
using namespace TMath;
using namespace sou;
//_____________________________________________________________________________________________________
TestGround::TestGround() {
//
// ctor
//
Configure();
Msg(EsafMsg::Info) << "*** TestGround built ***" << MsgDispatch;
}
//_____________________________________________________________________________________________________
TestGround::~TestGround(){
//
// dtor
//
}
//_____________________________________________________________________________________________________
void TestGround::Configure() {
//
// Configure the object
//
fType = Conf()->GetStr("TestGround.fType");
fAlbedo = Conf()->GetNum("TestGround.fAlbedo");
fAltitude = Conf()->GetNum("TestGround.fAltitude")*km;
}
//_____________________________________________________________________________________________________
void TestGround::RandomDir(const EarthVector& pos, EarthVector& incom_dir) const {
//
// set incom_dir to the outgoing direction of scattering, sampling Outgoing phase function
//
// get random LOCAL theta (not in the MES)
TRandom* rndm = EsafRandom::Get();
Double_t theta = acos(sqrt(rndm->Rndm())); // 2pi * phase_func(theta) * sin(theta) --- cos(theta) = x -----> dProba/dx = 2x
// then get a random direction in the MES
EarthVector local_Uz(pos);
local_Uz(2) += EarthRadius(); // now local_Uz is the local vertical direction
local_Uz.Unit();
EarthVector v = local_Uz.Orthogonal();
v.Rotate(rndm->Rndm()*2*Pi(),local_Uz);
EarthVector axis = v.Cross(local_Uz);
local_Uz.Rotate(theta,axis);
incom_dir = local_Uz; // outgoing direction
#ifdef DEBUG
if(IsUnderGround(pos + (100*m)*incom_dir)) Msg(EsafMsg::Warning) << "<RandomDir> Outgoing direction is not UPWARD, it should"<<MsgDispatch;
#endif
}
//_____________________________________________________________________________________________________
Double_t TestGround::Outgoing_phase_function(const EarthVector& pos1, const EarthVector& pos2) const {
//
// Phase function of the outgoing light, emitted from pos1 toward pos2
// earth sphericity taken into account
// Phase function is normalized when integrated over an hemisphere, i.e. theta within [0,Pi()/2]
//
Double_t rtn(0);
Double_t angle;
EarthVector temp(pos1);
temp(2) += EarthRadius(); // pos1 expressed in earth-centered frame -> gives local vertical in MES frame
angle = (pos2 - pos1).Angle(temp); // angle between outgoing direction and local vertical
if(fType == "isotropic")
rtn = 1./TwoPi();
else if(fType == "lambertian")
rtn = cos(angle)/Pi();
else invalid_argument("unvalid ground type for brdf");
return rtn;
}
//_____________________________________________________________________________________________________
EarthVector TestGround::GetImpact(const EarthVector& pos, const EarthVector& dir) const {
//
// Impact on ground form a given (position,direction) pair in atmosphere
//
// if pos is underground
EarthVector rtn;
if((pos.Zv() + 1*mm) < fAltitude) {
Msg(EsafMsg::Debug) << "Starting position is underground, so no ground impact calculated" << MsgDispatch;
rtn.SetXYZ(0,0,-HUGE);
return rtn;
}
#ifdef DEBUG
// Msg(EsafMsg::Debug) <<"Impact on ground = " << MsgDispatch;
#endif
// to know if dir is locally going upward
Double_t angle;
EarthVector temp(pos);
temp(2) += EarthRadius(); // pos expressed in earth-centered frame -> gives local vertical direction expressed in MES frame
angle = dir.Angle(temp); // angle between dir and local vertical
if(angle < (PiOver2() - kTolerance)) {
//Msg(EsafMsg::Debug) << MsgDispatch;
//Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND, track is locally going upward" << MsgDispatch;
rtn.SetXYZ(0,0,HUGE);
return rtn;
}
// now impact calculation
Double_t mag(0);
EarthVector direc = dir.Unit();
// spherical earth
Double_t b = pos*direc + direc.Z()*EarthRadius();
Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - fAltitude) - fAltitude*fAltitude;
pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
if(p.first == 0) {
#ifdef DEBUG
//Msg(EsafMsg::Debug) <<"NO IMPACT ON GROUND" << MsgDispatch;
#endif
rtn.SetXYZ(0,0,HUGE);
return rtn;
}
if(p.first == 1) mag = p.second[0];
else if(p.first ==2) mag = min(p.second[0],p.second[1]);
rtn = pos + mag*direc;
#ifdef DEBUG
//Msg(EsafMsg::Debug) << rtn << MsgDispatch;
#endif
return rtn;
}