// ESAF : Euso Simulation and Analysis Framework
// $Id: TestClouds.cc,v 1.12 2005/05/13 11:02:39 moreggia Exp $
// Sylvain Moreggia created Apr, 8 2004
#include "TestClouds.hh"
#include "EarthVector.hh"
#include "Config.hh"
#include "EConst.hh"
#include "EsafRandom.hh"
using EConst::EarthRadius;
// ctor
//______________________________________________________________________________
TestClouds::TestClouds() : Clouds() {
//
// ctor
//
Configure();
Build();
}
//______________________________________________________________________________
TestClouds::~TestClouds() {
//
// dtor
//
}
//______________________________________________________________________________
void TestClouds::Configure() {
//
// read config files
//
// clouds features
fType = Conf()->GetStr("TestClouds.fType");
ConfigFileParser* pConfig3 = Config::Get()->GetCF("Atmosphere","AtmosphereFactory");
fFeat = pConfig3->GetStr("Clouds.fFeat");
}
//______________________________________________________________________________
void TestClouds::Build() {
//
// build the object according to config files
//
TRandom* rndm = EsafRandom::Get();
Double_t ran = 0;
// set by hand
if(fFeat == "perso") {
fOD = Conf()->GetNum("TestClouds.fOD");
fTopAltitude = Conf()->GetNum("TestClouds.fTopAltitude")*km;
fThickness = Conf()->GetNum("TestClouds.fThickness")*km;
}
// set according redbook
Bool_t CS = false;
if(fFeat == "redbook" || fFeat == "redbook_noCS") {
// if clearsky allowed, CS set to true
if(fFeat == "redbook") CS = true;
do {
ran = rndm->Rndm();
if(ran <= 0.041) {
CS = true;
fOD = ran/0.041 * 0.1;
ran = rndm->Rndm() * 0.041;
if(ran <= 0.013) fTopAltitude = ran/0.013 * 3*km;
else if(ran <= 0.02) fTopAltitude = (ran - 0.013)/(0.02 - 0.013) * (7 - 3)*km + 3*km;
else if(ran <= 0.025) fTopAltitude = (ran - 0.02)/(0.025 - 0.02) * (10 - 7)*km + 7*km;
else fTopAltitude = (ran - 0.025)/(0.041 - 0.025) * (12 - 10)*km + 10*km; // 12 (Liou p.180)
}
else if (ran <= 0.334) {
CS = true;
fOD = (ran - 0.041)/(0.334 - 0.041) * (1 - 0.1) + 0.1;
ran = rndm->Rndm() * 0.293;
if(ran <= 0.056) fTopAltitude = ran/0.056 * 3*km;
else if(ran <= 0.078) fTopAltitude = (ran - 0.056)/(0.078 - 0.056) * (7 - 3)*km + 3*km;
else if(ran <= 0.114) fTopAltitude = (ran - 0.078)/(0.114 - 0.078) * (10 - 7)*km + 7*km;
else fTopAltitude = (ran - 0.114)/(0.293 - 0.114) * (12 - 10)*km + 10*km; // 12 (Liou p.180)
}
else if(ran <= 0.479) {
CS = true;
fOD = (ran - 0.334)/(0.479 - 0.334) * (2 - 1) + 1;
ran = rndm->Rndm() * 0.145;
if(ran <= 0.051) fTopAltitude = ran/0.051 * 3*km;
else if(ran <= 0.08) fTopAltitude = (ran - 0.051)/(0.08 - 0.051) * (7 - 3)*km + 3*km;
else if(ran <= 0.115) fTopAltitude = (ran - 0.08)/(0.115 - 0.08) * (10 - 7)*km + 7*km;
else fTopAltitude = (ran - 0.115)/(0.145 - 0.115) * (12 - 10)*km + 10*km; // 12 (Liou p.180)
}
else if(ran <= 0.822) {
CS = true;
fOD = (ran - 0.479)/(0.822 - 0.479) * (3 - 2) + 2; // 3 is arbitrary <=> 95% absorption
ran = rndm->Rndm() * 0.343;
if(ran <= 0.161) fTopAltitude = ran/0.161 * 3*km;
else if(ran <= 0.219) fTopAltitude = (ran - 0.161)/(0.219 - 0.161) * (7 - 3)*km + 3*km;
else if(ran <= 0.282) fTopAltitude = (ran - 0.219)/(0.282 - 0.219) * (10 - 7)*km + 7*km;
else fTopAltitude = (ran - 0.282)/(0.343 - 0.282) * (12 - 10)*km + 10*km; // 12 (Liou p.180)
}
else {
fOD = 0; // means no clouds
fTopAltitude = 0;
}
// thickness according to Liou (rough estimation)
ran = rndm->Rndm();
if(fTopAltitude == 0) fThickness = 0;
else if(fTopAltitude <= 3*km) fThickness = ran * 2*km;
else if(fTopAltitude <= 7*km) fThickness = ran * 4*km;
else fThickness = 2*km;
} while(!CS);
}
// checks
ConfigFileParser* pConfig1 = Config::Get()->GetCF("RadiativeTransfer","TestGround");
ConfigFileParser* pConfig2 = Config::Get()->GetCF("RadiativeTransfer","RadiativeFactory");
if(pConfig2->GetStr("RadiativeFactory.Ground") != "test")
Msg(EsafMsg::Warning) << "if TESTGROUND not used, be carefull of ground altitude w.r.t clouds" << MsgDispatch;
Double_t altg = pConfig1->GetNum("TestGround.fAltitude")*km;
if((fTopAltitude - fThickness) < altg) fThickness = fTopAltitude - altg - TOLERANCE;
if(fThickness < 0) {
// means no clouds
Msg(EsafMsg::Info) << "Clouds and Ground are too close -> No clouds" << MsgDispatch;
fOD = 0;
fTopAltitude = altg - TOLERANCE;
fThickness = 0;
}
// red book formula for clouds albedo as a function of OD
fAlbedo = 0.11*fOD / (1 + 0.11*fOD);
Msg(EsafMsg::Info) << "*** TestClouds built *** OD = " << fOD <<" -- TopAlt = " << fTopAltitude <<" -- Thick = "<< fThickness << MsgDispatch;
}
//______________________________________________________________________________
Bool_t TestClouds::IsInClouds(const EarthVector& pos) const {
//
// true if position is within clouds, false otherwise
//
Bool_t rtn = false;
if(!fOD) return rtn;
if(pos.IsUnderSeaLevel()) return false;
if(pos.Zv() <= fTopAltitude && pos.Zv() >= (fTopAltitude - fThickness)) rtn = true;
return rtn;
}
//______________________________________________________________________________
Double_t TestClouds::TotalOD(const EarthVector& pos,const EarthVector& dir) const {
//
// Returns the total OD along a track defined by pos and dir
// Case with pos within clouds handled
// NB : vertical OD is homogeneously distributed
//
if(!fOD) return 0;
EarthVector exit, enter, direc;
Double_t rtn(0), mag(0);
direc = dir.Unit();
enter = GetCloudImpact(pos,direc);
// if no clouds along the track
if(enter.Z() == HUGE) return 0;
exit = GetCloudImpact(pos,direc,"outgoing");
if(exit.Z() == HUGE) Msg(EsafMsg::Warning) << "if enter exists, exit should as well" << MsgDispatch;
mag = (exit - enter).Mag();
rtn = fOD * mag/fThickness;
return rtn;
}
//_____________________________________________________________________________________________________
EarthVector TestClouds::GetCloudImpact(const EarthVector& pos, const EarthVector& dir, string opt) const {
//
// Impact on clouds top surface form a given (position,direction) pair in atmosphere
// if no impact, (0,0,HUGE) is returned
// upward bunch is handled, but no scattering to the detector, only absorption
//
// if opt == "default" -> incoming impact returned
// if pos is within clouds, pos is returned
// if opt == "outgoing" -> outgoing impact returned
//
// init
EarthVector rtn(0,0,HUGE);
if(!fOD) return rtn;
Bool_t isinclouds, downward;
Double_t Altimpact;
EarthVector direc = dir.Unit();
// says if dir is going up/downward
downward = GoingDownward(pos,direc);
// says if pos is within the layer (altitude round-off effects < 1mm are handled)
isinclouds = IsInClouds(pos);
// set the impact surface altitude according to bunch present situation (up/downward, above/below, enter/exit clouds)
if(isinclouds) {
if(opt == "default") return pos;
else if(opt == "outgoing") {
if(!downward) Altimpact = fTopAltitude;
else Altimpact = (fTopAltitude - fThickness);
}
else Msg(EsafMsg::Panic) << "<GetCloudImpact()> wrong argument (internal pos)" << MsgDispatch;
}
else {
// if not going toward clouds
if(pos.Zv() < (fTopAltitude - fThickness) && downward) return rtn;
if(pos.Zv() > fTopAltitude && !downward) return rtn;
if(opt == "outgoing" && !downward || opt == "default" && downward) Altimpact = fTopAltitude;
else if(opt == "outgoing" && downward || opt == "default" && !downward) Altimpact = (fTopAltitude - fThickness);
else Msg(EsafMsg::Panic) << "<GetCloudImpact()> wrong argument (external pos)" << MsgDispatch;
}
// now calculate impact
Double_t mag;
Double_t b = pos*direc + direc.Z()*EarthRadius();
Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - Altimpact) - pow(Altimpact,2);
pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
if(p.first == 0) {
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"NO IMPACT ON Clouds" << MsgDispatch;
#endif
return rtn;
}
if(p.first == 1) mag = p.second[0];
else if(p.first ==2) {
if(p.second[0]*p.second[1] < 0) mag = max(p.second[0],p.second[1]);
else mag = min(p.second[0],p.second[1]);
}
if(mag < 0) Msg(EsafMsg::Warning) << "PB in impact calculation GetCloudImpact()" << MsgDispatch;
rtn = pos + mag*direc;
return rtn;
}
//_____________________________________________________________________________
Bool_t TestClouds::GoingDownward(const EarthVector& pos,const EarthVector& direc) const {
//
// to know if locally track is going up/downward
// locally means at pos nadir
// if perfectly horizontal (i.e. local angle of halfpi) returns false (i.e. upward)
//
if(!fOD) return false;
Bool_t rtn;
EarthVector dir = direc.Unit();
// find local zenith angle
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
// returns accordingly
if(angle <= halfpi) rtn = false;
else rtn = true;
return rtn;
}
//_____________________________________________________________________________
void TestClouds::Reset() {
//
// reset clodus features
//
fOD = 0;
fThickness = 0;
fTopAltitude = 0;
fAlbedo = 0;
Build();
}
ClassImp(TestClouds)