// $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;
}
}