Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

UniformLayersClouds - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: UniformLayersClouds.cc,v 1.7 2005/10/02 14:17:07 thea Exp $
// Author: Sylvain Moreggia   Aug, 18 2004

#include "UniformLayersClouds.hh"
#include "EConst.hh"

ClassImp(UniformLayersClouds)

using EConst::EarthRadius;
using namespace sou;

//_____________________________________________________________________________
 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 <= TMath::PiOver2()) rtn = false;
    else  rtn = true;
    
    return rtn;
}




About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005