// ESAF : Euso Simulation and Analysis Framework
// $Id: UniformLayersClouds.cc,v 1.6 2005/04/14 15:48:49 moreggia Exp $
// Author: Sylvain Moreggia Aug, 18 2004
#include "UniformLayersClouds.hh"
#include "EConst.hh"
ClassImp(UniformLayersClouds)
using EConst::EarthRadius;
//_____________________________________________________________________________
UniformLayersClouds::UniformLayersClouds() : Clouds(), fBox(0) {
//
// Constructor
//
Build();
}
//_____________________________________________________________________________
UniformLayersClouds::~UniformLayersClouds() {
//
// Destructor
//
}
//_____________________________________________________________________________
EarthVector UniformLayersClouds::GetCloudImpact(const EarthVector& pos, const EarthVector& dir,string opt) const {
//
// Return impact for given position and direction
// if opt is "outgoing", returns the track exit point
// Convention : - if already within clouds, returns pos
// - if no impact, returns (0,0,HUGE)
//
if(IsInClouds(pos)) return pos;
Int_t layerID = 0;
Double_t top = fBox->GetTopAlt();
Double_t bottom = fBox->GetBottomAlt();
// impact on the bottom
if(pos.Zv() < bottom && opt == "default" || pos.Zv() > top && opt == "outgoing") layerID = fBox->GetNumberOfLayers() - 1;
// impact on the top
else if(pos.Zv() > top && opt == "default" || pos.Zv() < bottom && opt == "outgoing") layerID = 0;
return GetLayerImpact(layerID,pos,dir,opt);
}
//_____________________________________________________________________________
Double_t UniformLayersClouds::TotalOD(const EarthVector& pos,const EarthVector& dir) const {
//
// returns total Optical Depth along the given track (defined by pos and dir)
//
Double_t rtn = 0;
Int_t i;
Bool_t downward = GoingDownward(pos,dir);
EarthVector incom, outgo;
Int_t layerID = LayerNumber(pos);
// downward track
if(downward) {
// if pos outside of clouds
if(layerID == -1) {
if(!PointToClouds(pos,dir)) return 0;
layerID = 0;
rtn = 0;
}
// if pos inside clouds
else {
outgo = GetLayerImpact(layerID,pos,dir);
rtn = InclinedODloc(layerID,pos,outgo);
layerID++;
}
// loop over layers
for(i=layerID; i<fBox->GetNumberOfLayers() - 1; i++) {
incom = GetLayerImpact(i,pos,dir);
outgo = GetLayerImpact(i,pos,dir,"outgoing");
rtn += InclinedODloc(i,incom,outgo);
}
}
// upward track
else {
// if pos outside of clouds
if(layerID == -1) {
if(!PointToClouds(pos,dir)) return 0;
layerID = fBox->GetNumberOfLayers() - 1;
rtn = 0;
}
// if pos inside clouds
else {
outgo = GetLayerImpact(layerID,pos,dir);
rtn = InclinedODloc(layerID,pos,outgo);
layerID--;
}
for(i=layerID; i>-1; i--) {
incom = GetLayerImpact(i,pos,dir);
outgo = GetLayerImpact(i,pos,dir,"outgoing");
rtn += InclinedODloc(i,incom,outgo);
}
}
return rtn;
}
//_____________________________________________________________________________
Double_t UniformLayersClouds::InclinedODloc(Int_t layerID,const EarthVector& incom,const EarthVector& outgo) const {
//
// gives layer OD between incoming and outgoing points
//
Double_t mag = (outgo - incom).Mag();
return fBox->GetODlocal(layerID) * mag/fBox->GetLayerWidth();
}
//_____________________________________________________________________________
EarthVector UniformLayersClouds::GetLayerImpact(Int_t layerID,const EarthVector& pos,const EarthVector& dir,string opt) const {
//
// Returns the impact of a track on the specified CloudLayer
// If pos is external, returned impact is incoming by default or outgoing if specified
// if no impact, returns (0,0,HUGE)
//
EarthVector rtn(0,0,HUGE);
Bool_t isinlayer, downward;
Double_t Altimpact;
Double_t top = fBox->GetLayerTop(layerID);
Double_t bottom = top - fBox->GetLayerWidth();
// to know if locally dir is going up/downward
downward = GoingDownward(pos,dir);
// to know if pos within the layer (altitude round-off effects < 1mm are handled)
isinlayer = IsInLayer(layerID,pos,dir);
if(!isinlayer) {
if(PointToLayer(layerID,pos,dir)) throw runtime_error("Wrong call to UniformLayersClouds::GetLayerImpact()");
if(opt == "outgoing" && !downward || opt == "default" && downward) Altimpact = top;
else if(opt == "outgoing" && downward || opt == "default" && !downward) Altimpact = bottom;
else throw runtime_error("UniformLayersClouds::GetLayerImpact -> this option is not handled");
}
else {
if(downward) Altimpact = bottom;
else Altimpact = top;
}
// now calculate impact
Double_t mag;
EarthVector direc = dir.Unit();
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
cerr <<"NO IMPACT ON CloudLayer number = " << layerID << endl;
#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) throw runtime_error("PB in impact calculation UniformLayersClouds::GetLayerImpact()");
rtn = pos + mag*direc;
return rtn;
}
//_____________________________________________________________________________
void UniformLayersClouds::Build() {
//
// build the CloudsBox object
// Standard clouds types drawn from Liou, Radiation and Cloud Processes in the Atmosphere, 1992
// It gives : - typical values for BetaE (= sigmaE * N) for wl = 0.5 microns
// - altitudes (influences of geographical coo. not considered), not given precisely
// clouds splitted into 100 layers
//
string cloudstype = Conf()->GetStr("UniformLayersClouds.fCloudsType");
Double_t BetaE, Hbottom, Htop;
// water clouds
if(cloudstype == "FairWeatherCumulus") {
BetaE = 165.9*pow(micrometer,2) * 293/cm3;
Hbottom = 1.5*km;
Htop = 2.2*km;
}
else if(cloudstype == "CongestusCumulus") {
BetaE = 468.3*pow(micrometer,2) * 207/cm3;
Hbottom = 1.5*km;
Htop = 2.2*km;
}
else if(cloudstype == "CumuloNimbus") {
BetaE = 481.6*pow(micrometer,2) * 72/cm3;
Hbottom = 1.5*km;
Htop = 4.5*km;
}
else if(cloudstype == "StratoCumulus") {
BetaE = 120.4*pow(micrometer,2) * 350/cm3;
Hbottom = 1*km;
Htop = 1.2*km;
}
else if(cloudstype == "NimboStratus") {
BetaE = 410.1*pow(micrometer,2) * 330/cm3;
Hbottom = 1*km;
Htop = 3*km;
}
else if(cloudstype == "Stratus") {
BetaE = 473.1*pow(micrometer,2) * 464/cm3; // over ocean (for land 464 -> 260)
Hbottom = 0.5*km;
Htop = 0.7*km;
}
// ice clouds
else if(cloudstype == "CirroStratus") {
BetaE = 0.3865/km;
Hbottom = 8*km;
Htop = 10*km;
}
else if(cloudstype == "UncinusCirrus") {
BetaE = 2.6058/km;
Hbottom = 8*km;
Htop = 10*km;
}
else if(cloudstype == "WarmCirrus") {
BetaE = 0.6525/km;
Hbottom = 6*km;
Htop = 7.8*km;
}
else if(cloudstype == "ColdCirrus") {
BetaE = 0.1662/km;
Hbottom = 10*km;
Htop = 11.8*km;
}
fBox = new CloudsBox(exp(-BetaE*(Htop-Hbottom)),Hbottom,Htop,100);
if(!fBox) throw runtime_error("UniformLayersClouds::Build() -> new command returned NIL pointer for fBox construction");
}
//_____________________________________________________________________________
Bool_t UniformLayersClouds::IsInLayer(Int_t layerID,const EarthVector& pos, EarthVector dir) const {
//
// says if pos is in the specified layer
// if dir is specified, altitude round-off effects < 1mm are handled (cf. cases with pos on a layer edge)
//
Bool_t rtn;
Double_t top = fBox->GetLayerTop(layerID);
Double_t bottom = top - fBox->GetLayerWidth();
if(bottom < pos.Zv() && pos.Zv() < top) rtn = true;
else if(dir.Mag()) {
Bool_t downward = GoingDownward(pos,dir);
if(fabs(bottom - pos.Zv()) < 1 && !downward) rtn = true;
else if(fabs(top - pos.Zv()) < 1 && downward) rtn = true;
else rtn = false;
}
return rtn;
}
//_____________________________________________________________________________
Bool_t UniformLayersClouds::IsInClouds(const EarthVector& pos) const {
//
// says if pos is in clouds or not
//
Bool_t rtn;
if(fBox->GetBottomAlt() <= pos.Zv() && pos.Zv() <= fBox->GetTopAlt()) rtn = true;
else if(fabs(pos.Zv() - fBox->GetBottomAlt()) < 1 || fabs(pos.Zv() - fBox->GetTopAlt()) < 1) rtn = true;
else rtn = false;
return rtn;
}
//_____________________________________________________________________________
Int_t UniformLayersClouds::PointToLayer(Int_t layerID,const EarthVector& pos,const EarthVector& dir) const {
//
// says if track points toward clouds
// returns -1 if is inside the layer
//
if(IsInLayer(layerID,pos)) return -1;
Int_t rtn;
Bool_t downward = GoingDownward(pos,dir);
Double_t top = fBox->GetLayerTop(layerID);
if(downward && pos.Zv() < top - fBox->GetLayerWidth() || !downward && pos.Zv() > top) rtn = 0;
else rtn = 1;
return rtn;
}
//_____________________________________________________________________________
Int_t UniformLayersClouds::PointToClouds(const EarthVector& pos,const EarthVector& dir) const {
//
// says if track points toward clouds
// returns -1 if is inside clouds
//
Int_t rtn;
if(IsInClouds(pos)) rtn = -1;
else if(PointToLayer(fBox->GetNumberOfLayers() - 1,pos,dir) || PointToLayer(0,pos,dir)) rtn = 1;
else rtn = 0;
return rtn;
}
//_____________________________________________________________________________
Int_t UniformLayersClouds::LayerNumber(const EarthVector& pos) const {
//
// gives the layer number where pos is
// returns -1 if is not in the cloud
//
Bool_t in = IsInClouds(pos);
if(!in) return -1;
Double_t diff = fBox->GetTopAlt() - pos.Zv();
return (Int_t) floor(diff/fBox->GetLayerWidth());
}
//_____________________________________________________________________________
Bool_t UniformLayersClouds::GoingDownward(const EarthVector& pos,const EarthVector& dir) 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)
//
Bool_t rtn;
// 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;
}