Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

Atmosphere - source file

// $Id: Atmosphere.cc,v 1.42 2005/05/13 08:12:10 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 <TF1.h>

ClassImp(Atmosphere)

using EConst::EarthRadius;

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 = ((EConst::EarthRadius() + Altitude)/(EConst::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) 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;
    
    if(opt == "pos") {
	EarthVector Vmin;
	EarthVector Vmax;
	if(V1.Zv() <= V2.Zv()) {
            Vmin = V1;
	    Vmax = V2;
	}
	else {
            Vmin = V2;
	    Vmax = V1;
	}

	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();
	    rho2 = Air_Density(h2);
	    dX = (rho1 + rho2)/2 * dl;
	    if(dX == 0) break;
	    X += dX;
	    rho1 = rho2;
        }
    }
	
	    
    else if(opt == "dir") {
        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;
	    if(dX == 0) break;
	    X += dX;
	    rho1 = rho2;
        }
   }
   
    else  Msg(EsafMsg::Warning) << "Wrong option for grammage calculation" << MsgDispatch;

    return X;
}
//______________________________________________________________________________
 EarthVector Atmosphere::InvertGrammage(const EarthVector& pos1, const EarthVector& direc, Double_t depth) const {
    //
    // returns final position for given pos1 and air depth
    // precision achieved (in air depth) <1e-3 ( <=> <10m for 10km vertical track)
    //
    
    // init
    EarthVector max(pos1);
    EarthVector min(pos1);
    EarthVector middle(pos1);
    EarthVector dir = direc.Unit();
    Double_t rho0 = Air_Density(0);
    Double_t deltaL = depth / rho0;
    max = ImpactASL(pos1,dir);
    
    // if sea level reached before foreseen air depth been travelled
    if((max-pos1).Mag() < deltaL) {
	return max;
    }
    min = pos1 + deltaL*dir;    
    
    // if no impact
    if(max.Z() == HUGE) max = pos1 + 500e6*dir;
    while((Grammage(pos1,max) - Grammage(pos1,min)) > 0.001*depth) {
       middle = 0.5*(max+min);
       if ( fabs(Grammage(pos1,middle) - depth) < 0.001*depth) return middle;
       if ( (Grammage(pos1,middle) - depth)  < 0) min = middle;
       else max = middle;
    }
    
    return min;
}

//______________________________________________________________________________
 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 < halfpi) {
        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;    
}

//______________________________________________________________________________
 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};
	return depthIntegral->Integral(eps,1-eps,pars,fDepthCalculationPrecision);
    }
    else { 
	Double_t Hstar = ( EConst::EarthRadius() + H)*TMath::Sin(TMath::Pi() - Theta ) - 
	    EConst::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 14-May-2005 21:31