Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

Atmosphere - source file

// $Id: Atmosphere.cc,v 1.51 2005/11/14 11:21:16 moreggia Exp $
// S. Moreggia created 27 October 2003

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: Atmosphere                                                           *
 *  Package: atmosphere                                                      *
 *  Coordinator: S. Moreggia                                                 *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// Atmosphere
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "Atmosphere.hh"
#include "AtmosphereFactory.hh"
#include "LinsleyAtmosphere.hh"
#include "LowtranAtmosphere.hh"
#include "MSISE_00Atmosphere.hh"
#include "Config.hh"
#include "EConst.hh"
#include "utils.hh"

#include <TF1.h>

ClassImp(Atmosphere)

using namespace sou;
using namespace TMath;
using namespace EConst;

Atmosphere* Atmosphere::fChild = NULL;

//______________________________________________________________________________
Double_t densityIntegral(Double_t *x, Double_t *par)
{
    // 
    Double_t Altitude = par[0],  Angle = par[1];
    Double_t H0=8.43*km; // An auxiulary constant, km
    Double_t h = Altitude - H0*TMath::Log(x[0]);
    Double_t f = ((EarthRadius() + Altitude)/(EarthRadius() + h))*TMath::Sin(Angle);
    return H0*Atmosphere::Get()->Air_Density(h)/TMath::Sqrt(1.-TMath::Power(f,2))/x[0];
}

//______________________________________________________________________________
Atmosphere::Atmosphere() : EsafConfigurable(), EsafMsgSource() {
    //
    // ctor
    //
    fClouds = AtmosphereFactory::Get()->GetClouds();
    if(!fClouds) Msg(EsafMsg::Panic) <<"No clouds built : memory problems"<< MsgDispatch;
    fDepthCalculationPrecision = 1.e-5;
}

//______________________________________________________________________________
 Atmosphere::~Atmosphere() {
    //
    // dtor
    //
    if(fClouds) delete fClouds;
    fClouds = 0;
}

//______________________________________________________________________________
 const Atmosphere* Atmosphere::Get() {
    //
    // Static method to get the right atmosphere singleton from anywhere in the code
    //
    if(fChild == 0) {
        ConfigFileParser *pConfig = Config::Get()->GetCF("Atmosphere","AtmosphereFactory");
        string type = pConfig->GetStr("Atmosphere.type");
        
	if(type == "linsley")
            LinsleyAtmosphere::CreateInstance();
        
	else if(type == "msise00")
	    MSISE_00Atmosphere::CreateInstance();

        else if(type == "lowtran")
            LowtranAtmosphere::CreateInstance(); 

        else cout<<"Wrong type of atmospheren";
        fChild->Msg(EsafMsg::Info) << "Atmosphere built" << MsgDispatch;
    }
    
     return fChild;
}
    
//______________________________________________________________________________
 void Atmosphere::Reset() {
    //
    // Static method to delete atmosphere
    //
    if(fChild) {
	delete fChild;
	fChild = 0;
#ifdef DEBUG
        cout << "\n[Debug]    Atmosphere has been deleted\n";
#endif
    }
    else cout << "\n[Warning]   There is no Atmosphere to reset\n";
}

//______________________________________________________________________________
long double Atmosphere::Index(Double_t h, Double_t wl) const {
    //
    // Returns air index value at a given altitude
    // By default, it is taken CONSTANT with wavelength (lambda=350nm)
    // (formula drawn from Reinhard Beer, Hanbook of Optics, Chap.2)
    //
    long double rtn;
    rtn = (long double)Index_Minus1(h,wl);
    return 1 + rtn;
}

//______________________________________________________________________________
 Double_t Atmosphere::Index_Minus1(Double_t h, Double_t wl) const {
    //
    // Returns (n-1) at a given altitude
    // By default, it is taken CONSTANT with wavelength (lambda=350nm)
    // (formula drawn from Reinhard Beer, Hanbook of Optics, Chap.2)
    //
    Double_t rtn;
    rtn = 88.43 + 185.08/(1 - 1/pow(wl/cm*1.14e5,2)) + 4.11/(1 - 1/pow(wl/cm*6.24e4,2));
    rtn *= (Pressure(h) - WaterVaporPartialPressure(h))/Temperature(h);
    rtn *= 296.15*kelvin/atmosphere;
    rtn += (43.49 - 1/pow(wl/cm*1.7e4,2))*WaterVaporPartialPressure(h)/atmosphere;
    return rtn*1e-6; 
}

//______________________________________________________________________________
 Double_t Atmosphere::Grammage(const EarthVector& V1, const EarthVector& V2, string opt, Double_t maxalt) const {
    //
    // Calculate grammage
    // - opt = pos : between two positions V1, V2 in the atmosphere
    // - opt = dir : from a point V1 until atmosphere top, along a track of given angles (V2 thus used for direction)
    //
    Double_t dl = 500*m;
    Double_t dX, rho2, h2;
    EarthVector U, inter;
    Double_t X = 0;
    Double_t rho1;
    Bool_t IsTooHigh = false;
    
    if(opt == "pos") {
	EarthVector Vmin;
	EarthVector Vmax;
	if(V1.IsUnderSeaLevel() || V2.IsUnderSeaLevel()) {
	    Msg(EsafMsg::Warning) << "<Grammage> V1 or/and V2 is under sea level : it SHOULD NOT, O returned"<<MsgDispatch;
	    return 0.;
	}
	if(V1.Zv() <= V2.Zv()) {
            Vmin = V1;
	    Vmax = V2;
	}
	else {
            Vmin = V2;
	    Vmax = V1;
	}

	if(Vmax.Zv() > maxalt) IsTooHigh = true;
	U = (Vmax - Vmin).Unit();
	inter = Vmin;
	rho1 = Air_Density(Vmin.Zv());
	
        while((Vmax - Vmin).Mag() > (inter - Vmin).Mag()) {
	    if( (inter + U*dl - Vmin).Mag() >= (Vmax - Vmin).Mag() ) {
	        rho2 = Atmosphere::Get()->Air_Density(Vmax.Zv());
	        dX = (rho1 + rho2)/2 * (Vmax - inter).Mag();
	        X += dX;
	        break;
	    }
	    inter += U*dl;
	    h2 = inter.Zv();
	    if(h2 > maxalt && IsTooHigh) break;
	    rho2 = Air_Density(h2);
	    dX = (rho1 + rho2)/2 * dl;
	    X += dX;
	    rho1 = rho2;
        }
    }
	
	    
    else if(opt == "dir") {
	if(V1.IsUnderSeaLevel()) Msg(EsafMsg::Warning) << "<Grammage> V1 is under sea level : it SHOULD NOTd"<<MsgDispatch;
        EarthVector TOA = ImpactAtTOA(V1,V2);
	if(TOA.Z() == -HUGE) {
	    Msg(EsafMsg::Warning) << "<Grammage("dir")> : starting position is UNDERSEALEVEL, O returned" << MsgDispatch;
	    return 0.;
	}
	if(TOA.Z() == HUGE) {
	    Msg(EsafMsg::Warning) << "<Grammage("dir")> : track DOES NOT GO OUT atmosphere, O returned" << MsgDispatch;
	    return 0.;
	}
	U = V2.Unit();
	inter = V1;
	rho1 =  Air_Density(V1.Zv());
	while(true) {
	    inter += U*dl;
	    h2 = inter.Zv();
	    rho2 = Air_Density(h2);
	    dX = (rho1 + rho2)/2 * dl;
	    X += dX;
	    rho1 = rho2;
	    if(inter.Mag() > (TOA - V1).Mag()) break;
        }
   }
   
    else  Msg(EsafMsg::Warning) << "Wrong option for grammage calculation" << MsgDispatch;

    return X;
}

//______________________________________________________________________________
 Int_t Atmosphere::InvertGrammage(const EarthVector& pos1, const EarthVector& direc, Double_t depth, EarthVector& rtn, Double_t maxtof) const {
    //
    // calculate final position for given pos1,direc and air depth
    // 
    // it returns : -1 if pos1 is under sea level or if none impact (latter should not occur)
    //               0 if position found
    //               1 if TOF cut
    //               2 if TOA reached before
    //               3 if sea level reached before
    //
    
    // init
    Int_t status(-1);
    if(pos1.IsUnderSeaLevel()) return status;
    rtn.SetMag(0);
    EarthVector dir = direc.Unit();
    if(maxtof < 0) maxtof = 1*second;

    EarthVector impact(1);
    
    // check if sea level reached before foreseen air depth been travelled
    // check if TOA reached before the foreseen air depth has been travelled
    impact = ImpactASL(pos1,dir);
    // if no impact at sea level
    if(impact.Z() == HUGE) {
        impact = ImpactAtTOA(pos1,dir,50*km); // (should be same as in singlepropagator::DoOnePhoton() method) TOA at 50km to avoid pbic cases at high altitudes
	if(impact.Z() == HUGE) {
	    Msg(EsafMsg::Warning) << "<InvertGrammage> track reaches NEITHER sea level NOR TOA -> SHOUD NOT HAPPEN"<<MsgDispatch;
            return -1;
	}
        if(Grammage(pos1,impact) < depth) {
	    rtn = impact;
	    return 2;
        }
    }
    else {
        if(Grammage(pos1,impact) < depth) {
	    rtn = impact;
	    return 3;
        }
    }
    
    
    // step by step process until it reaches depth value
    Double_t dl = 500*m;
    Double_t dX, rho2, h2;
    EarthVector U, inter;
    Double_t X = 0.;
    Double_t rho1;
    U = dir;
    inter = pos1;
    rho1 = Air_Density(pos1.Zv());
    if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();

    Int_t counter = 0;
    Int_t tofcut = Int_t(maxtof*Clight() / dl) + 1;
    while(true) {
        if(counter++ > tofcut) {
#ifdef DEBUG
            //Msg(EsafMsg::Debug) <<"<atmoInvertGrammage> TOF cut reached : photon dumped here "<<MsgDispatch;
#endif
            if(rtn.Mag() == 0) rtn = pos1;
	    return 1;
	}
        if(counter++ > 5000) {
            Msg(EsafMsg::Warning) <<"<atmoInvertGrammage> \"infinite\" loop (>5000 iterations) broken by hand "<<MsgDispatch;
            if(rtn.Mag() == 0) rtn = pos1;
            return 1;
        }
	inter += U*dl;
	h2 = inter.Zv();
	rho2 = Air_Density(h2);
	dX = (rho1 + rho2)/2 * dl;
	if((X + dX) > depth) break;
	X += dX;
	rho1 = rho2;
	if((inter - impact).Mag() < dl) dl = (inter - impact).Mag();
    }
    
    // come back to last iteration for fine tuning
    rtn = inter - U*dl;
    Double_t slope = (rho2 - rho1) / dl;
    dl = sign(slope) * sqrt( (rho1/slope) * (rho1/slope) + 2.*(depth - X) / slope) - rho1/slope;
    rtn += U*dl;

    if(rtn.Mag() == 0) rtn = pos1;
    
    return 0;
}

//______________________________________________________________________________
 EarthVector Atmosphere::ImpactASL(const EarthVector& pos, const EarthVector& dir) const {
    //
    // returns impact at sea level of a track defined by starting position and direction
    //

    // if pos is under sea level
    EarthVector rtn;
    if(pos.IsUnderSeaLevel()) {
        rtn.SetXYZ(0,0,-HUGE);
        return rtn;
    }

    // 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()) {
        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();
    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
    if(p.first == 0) {
        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;

    return rtn;    
}

//______________________________________________________________________________
 EarthVector Atmosphere::ImpactAtTOA(const EarthVector& pos, const EarthVector& dir, Double_t TOA_alt) const {
    //
    // returns impact at Top Of Atmosphere of a track defined by starting position and direction
    // TOA set at 100km altitude by default, but can also be set by dvper (cf. InvertGrammage() method)
    //

    // if pos is under sea level
    EarthVector rtn;
    if(pos.IsUnderSeaLevel()) {
        rtn.SetXYZ(0,0,-HUGE);
        return rtn;
    }

    // now, impact calculation
    Double_t mag(0);
    EarthVector direc = dir.Unit();
    Double_t altTOA = TOA_alt;

    // spherical earth
    Double_t b = pos*direc + direc.Z()*EarthRadius();
    Double_t c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - altTOA) - pow(altTOA,2);
    pair<Int_t,Double_t*>& p = findRoots(1.,2*b,c);
    if(p.first == 0) {
        rtn.SetXYZ(0,0,HUGE);
        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;    
}



// used in reco only : obsolete, use Grammage method
//______________________________________________________________________________
 Double_t Atmosphere::Depth(const Double_t H, const Double_t Theta) const {
    // Compute the atmosphere depth between the given point with (h,theta) coordinates and infinity
    //
    Double_t eps=1.e-5;
    TF1 *depthIntegral = new TF1("depthIntegral",densityIntegral,0.,1.,2);
    
    if (Theta <= TMath::PiOver2()) {
	Double_t pars[2] = {H,Theta};
        Double_t depth = depthIntegral->Integral(eps,1-eps,pars,fDepthCalculationPrecision);
        delete depthIntegral;
	return depth;
    }
    else { 
	Double_t Hstar = ( EarthRadius() + H)*TMath::Sin(TMath::Pi() - Theta ) - 
	    EarthRadius();
	Double_t pars1[2] = {Hstar,TMath::PiOver2()};
	Double_t pars2[2] = {H,TMath::Pi() - Theta};
	Double_t depth=2*depthIntegral->Integral(eps,1-eps,pars1) - depthIntegral->Integral(eps,1-eps,pars2,
            fDepthCalculationPrecision);
	delete depthIntegral;
	return depth;
    }
}

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