Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestClouds - source file

// 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)
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31