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