// $Id: LowtranAerosol.cc,v 1.2 2006/05/03 08:28:58 moreggia Exp $
// Author: Sylvain Moreggia 2006/04/24
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: LowtranAerosol *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// LowtranAerosol
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// WARNING : ground altitude is not handled internally. It must be checked in the other parts of the code, from where herebelow methods are called
//
#include "LowtranAerosol.hh"
#include "EarthVector.hh"
#include "Config.hh"
#include "EConst.hh"
#include "EsafRandom.hh"
#include "NumbersFileParser.hh"
#include "ConfigFileParser.hh"
using EConst::EarthRadius;
using namespace sou;
using namespace TMath;
ClassImp(LowtranAerosol)
//_____________________________________________________________________________
LowtranAerosol::LowtranAerosol(string name) : Aerosol(name) {
//
// Constructor
//
Configure();
Build();
}
//_____________________________________________________________________________
LowtranAerosol::~LowtranAerosol() {
//
// Destructor
//
}
//______________________________________________________________________________
void LowtranAerosol::Configure() {
//
// read config files
//
// init
fWlrange[0] = 200*nm;
fWlrange[1] = 300*nm;
fWlrange[2] = 337.1*nm;
fWlrange[3] = 550*nm;
fWlrange[4] = 694.3*nm;
fWlrange[5] = 1060*nm;
fWlrange[6] = 1536*nm;
// aerosol model
fType = Conf()->GetStr("LowtranAerosol.fType");
}
//______________________________________________________________________________
void LowtranAerosol::Build() {
//
// build the object according to config files
// Properties at 70% humidity are tabulated : extinction, absorption, DHG coeff. for phase function
// DHG coeff. come from fits on tabulated lowtran phase functions
// All are wavelength dependent (7 values tabulated)
// For Altitude profile, 3 values tabulated at 0, 1 and 2km
//
// build aerosol models
if(fType == "rural23") {
fSz[0] = 0.158/km;
fSz[1] = 0.0991/km;
fSz[2] = 0.0621/km;
// values for 70% rel. humidity
fKwl_ext[0] = 2.09544; // 200 nm
fKwl_ext[1] = 1.74165; // 300 nm
fKwl_ext[2] = 1.59981; // 337.1 nm
fKwl_ext[3] = 1.; // 550 nm
fKwl_ext[4] = 0.75316; // 694.3 nm
fKwl_ext[5] = 0.42171; // 1060 nm
fKwl_ext[6] = 0.24323; // 1536 nm
//DEBUG : to test w.r.t lowtran, 0% humidity values
//fKwl_ext[0] = 2.09291; // 200 nm
//fKwl_ext[1] = 1.74582; // 300 nm
//fKwl_ext[2] = 1.60500; // 337.1 nm
//fKwl_ext[3] = 1.; // 550 nm
//fKwl_ext[4] = 0.75203; // 694.3 nm
//fKwl_ext[5] = 0.41943; // 1060 nm
//fKwl_ext[6] = 0.24070; // 1536 nm
// values for 70% rel. humidity
fKwl_abs[0] = 0.62968; // 200 nm
fKwl_abs[1] = 0.10816; // 300 nm
fKwl_abs[2] = 0.07671; // 337.1 nm
fKwl_abs[3] = 0.05380; // 550 nm
fKwl_abs[4] = 0.04684; // 694.3 nm
fKwl_abs[5] = 0.05335; // 1060 nm
fKwl_abs[6] = 0.04614; // 1536 nm
// values for 70% rel. humidity
fG1[0] = 0.7632;
fG1[1] = 0.6928;
fG1[2] = 0.6885; // not available in lowtran, linear interpolation used here
fG1[3] = 0.6638;
fG1[4] = 0.6638;
fG1[5] = 0.6638;
fG1[6] = 0.6590;
fG2[0] = -0.6488;
fG2[1] = -0.6202;
fG2[2] = -0.6076; // not available in lowtran, linear interpolation used here
fG2[3] = -0.5351;
fG2[4] = -0.5351;
fG2[5] = -0.5351;
fG2[6] = -0.3744;
fF[0] = 0.7185;
fF[1] = 0.985;
fF[2] = 0.9837; // not available in lowtran, linear interpolation used here
fF[3] = 0.9765;
fF[4] = 0.9765;
fF[5] = 0.9765;
fF[6] = 0.9622;
}
else if(fType == "rural5") {
fSz[0] = 0.77/km;
fSz[1] = 0.77/km;
fSz[2] = 0.0621/km;
// values for 70% rel. humidity
fKwl_ext[0] = 2.09544; // 200 nm
fKwl_ext[1] = 1.74165; // 300 nm
fKwl_ext[2] = 1.59981; // 337.1 nm
fKwl_ext[3] = 1.; // 550 nm
fKwl_ext[4] = 0.75316; // 694.3 nm
fKwl_ext[5] = 0.42171; // 1060 nm
fKwl_ext[6] = 0.24323; // 1536 nm
//DEBUG : to test w.r.t lowtran, 0% humidity values
//fKwl_ext[0] = 2.09291; // 200 nm
//fKwl_ext[1] = 1.74582; // 300 nm
//fKwl_ext[2] = 1.60500; // 337.1 nm
//fKwl_ext[3] = 1.; // 550 nm
//fKwl_ext[4] = 0.75203; // 694.3 nm
//fKwl_ext[5] = 0.41943; // 1060 nm
//fKwl_ext[6] = 0.24070; // 1536 nm
// values for 70% rel. humidity
fKwl_abs[0] = 0.62968; // 200 nm
fKwl_abs[1] = 0.10816; // 300 nm
fKwl_abs[2] = 0.07671; // 337.1 nm
fKwl_abs[3] = 0.05380; // 550 nm
fKwl_abs[4] = 0.04684; // 694.3 nm
fKwl_abs[5] = 0.05335; // 1060 nm
fKwl_abs[6] = 0.04614; // 1536 nm
// values for 70% rel. humidity
fG1[0] = 0.7632;
fG1[1] = 0.6928;
fG1[2] = 0.6885; // not available in lowtran, linear interpolation used here
fG1[3] = 0.6638;
fG1[4] = 0.6638;
fG1[5] = 0.6638;
fG1[6] = 0.6590;
fG2[0] = -0.6488;
fG2[1] = -0.6202;
fG2[2] = -0.6076; // not available in lowtran, linear interpolation used here
fG2[3] = -0.5351;
fG2[4] = -0.5351;
fG2[5] = -0.5351;
fG2[6] = -0.3744;
fF[0] = 0.7185;
fF[1] = 0.985;
fF[2] = 0.9837; // not available in lowtran, linear interpolation used here
fF[3] = 0.9765;
fF[4] = 0.9765;
fF[5] = 0.9765;
fF[6] = 0.9622;
}
else if(fType == "urban5") {
fSz[0] = 0.77/km;
fSz[1] = 0.77/km;
fSz[2] = 0.0621/km;
// values for 70% rel. humidity
fKwl_ext[0] = 1.95582; // 200 nm
fKwl_ext[1] = 1.64994; // 300 nm
fKwl_ext[2] = 1.53070; // 337.1 nm
fKwl_ext[3] = 1.; // 550 nm
fKwl_ext[4] = 0.77614; // 694.3 nm
fKwl_ext[5] = 0.46639; // 1060 nm
fKwl_ext[6] = 0.29487; // 1536 nm
//DEBUG : to test w.r.t lowtran, 0% humidity values
//fKwl_ext[0] = 1.88816; // 200 nm
//fKwl_ext[1] = 1.63316; // 300 nm
//fKwl_ext[2] = 1.51867; // 337.1 nm
//fKwl_ext[3] = 1.; // 550 nm
//fKwl_ext[4] = 0.77785; // 694.3 nm
//fKwl_ext[5] = 0.47095; // 1060 nm
//fKwl_ext[6] = 0.30006; // 1536 nm
// values for 70% rel. humidity
fKwl_abs[0] = 0.69032; // 200 nm
fKwl_abs[1] = 0.49367; // 300 nm
fKwl_abs[2] = 0.45165; // 337.1 nm
fKwl_abs[3] = 0.29741; // 550 nm
fKwl_abs[4] = 0.24070; // 694.3 nm
fKwl_abs[5] = 0.17399; // 1060 nm
fKwl_abs[6] = 0.13146; // 1536 nm
// values for 70% rel. humidity
fG1[0] = 0.7906;
fG1[1] = 0.7632;
fG1[2] = 0.7478; // not available in lowtran, linear interpolation used here
fG1[3] = 0.6592;
fG1[4] = 0.6592;
fG1[5] = 0.6536;
fG1[6] = 0.6590;
fG2[0] = -0.0;
fG2[1] = -0.6488;
fG2[2] = -0.6797; // not available in lowtran, linear interpolation used here
fG2[3] = -0.8568;
fG2[4] = -0.8568;
fG2[5] = 0.5032;
fG2[6] = -0.3744;
fF[0] = 1.;
fF[1] = 0.7185;
fF[2] = 0.7602; // not available in lowtran, linear interpolation used here
fF[3] = 0.9998;
fF[4] = 0.9998;
fF[5] = 0.7134;
fF[6] = 0.9622;
}
else if(fType == "maritime23") {
fSz[0] = 0.158/km;
fSz[1] = 0.0991/km;
fSz[2] = 0.0621/km;
// values for 70% rel. humidity
fKwl_ext[0] = 1.36924; // 200 nm
fKwl_ext[1] = 1.25443; // 300 nm
fKwl_ext[2] = 1.20835; // 337.1 nm
fKwl_ext[3] = 1.; // 550 nm
fKwl_ext[4] = 0.91367; // 694.3 nm
fKwl_ext[5] = 0.77089; // 1060 nm
fKwl_ext[6] = 0.64987; // 1536 nm
//DEBUG : to test w.r.t lowtran, 0% humidity values
//fKwl_ext[0] = 1.47576; // 200 nm
//fKwl_ext[1] = 1.32614; // 300 nm
//fKwl_ext[2] = 1.26171; // 337.1 nm
//fKwl_ext[3] = 1.; // 550 nm
//fKwl_ext[4] = 0.88133; // 694.3 nm
//fKwl_ext[5] = 0.70297; // 1060 nm
//fKwl_ext[6] = 0.56487; // 1536 nm
// values for 70% rel. humidity
fKwl_abs[0] = 0.23367; // 200 nm
fKwl_abs[1] = 0.03127; // 300 nm
fKwl_abs[2] = 0.02070; // 337.1 nm
fKwl_abs[3] = 0.01297; // 550 nm
fKwl_abs[4] = 0.01063; // 694.3 nm
fKwl_abs[5] = 0.01285; // 1060 nm
fKwl_abs[6] = 0.01190; // 1536 nm
// values for 70% rel. humidity
fG1[0] = 0.8;
fG1[1] = 0.75;
fG1[2] = 0.7458; // not available in lowtran, linear interpolation used here
fG1[3] = 0.7214;
fG1[4] = 0.75;
fG1[5] = 0.75;
fG1[6] = 0.7445;
fG2[0] = -0.5396;
fG2[1] = -0.75;
fG2[2] = -0.7302; // not available in lowtran, linear interpolation used here
fG2[3] = -0.6163;
fG2[4] = -0.5840;
fG2[5] = -0.5840;
fG2[6] = -0.6509;
fF[0] = 0.9492;
fF[1] = 0.9662;
fF[2] = 0.9644; // not available in lowtran, linear interpolation used here
fF[3] = 0.9539;
fF[4] = 0.9643;
fF[5] = 0.9643;
fF[6] = 0.9926;
}
else Msg(EsafMsg::Panic) << "<Build> Wrong type of LowtranAerosol model = "<<fType << MsgDispatch;
Msg(EsafMsg::Info) << "*** LowtranAerosol built ***"<< MsgDispatch;
}
//______________________________________________________________________________
Double_t LowtranAerosol::Kwl(Double_t wl, Bool_t abs) const {
//
// interpolate fKwl coeff.
// if abs=kTRUE, return value for absorption, else return value for extinction
// Linear interpolation as Lowtran does
//
Double_t rtn(0.), k1(0.), k2(0.), wl1(0.), wl2(0.);
CheckWlRange(wl);
// extinction coeff.
if(!abs) {
for(Int_t i=0; i < 6; i++) {
if(wl >= fWlrange[i] && wl <= fWlrange[i+1]) {
k1 = fKwl_ext[i];
k2 = fKwl_ext[i+1];
wl1 = fWlrange[i];
wl2 = fWlrange[i+1];
break;
}
}
}
// absorption coeff.
else {
for(Int_t i=0; i < 6; i++) {
if(wl >= fWlrange[i] && wl <= fWlrange[i+1]) {
k1 = fKwl_abs[i];
k2 = fKwl_abs[i+1];
wl1 = fWlrange[i];
wl2 = fWlrange[i+1];
break;
}
}
}
rtn = (k2 - k1) / (wl2 - wl1) * (wl - wl2) + k2;
return rtn;
}
//______________________________________________________________________________
Bool_t LowtranAerosol::IsInAerosol(const EarthVector& pos) const {
//
// true if position is within aerosol, false otherwise
//
return IsInAerosolZone(pos,ALL);
}
//______________________________________________________________________________
Bool_t LowtranAerosol::IsInAerosolZone(const EarthVector& pos, AerZone zone) const {
//
// true if position is within aerosol, false otherwise
// deals with aerosol zones
//
// init
Bool_t rtn = false;
Double_t alt_low = 0*km;
Double_t alt_high = 2*km;
if(pos.IsUnderSeaLevel()) return false;
// depends on aerosol zone
if(zone == BOUND1) alt_high = 1*km;
else if(zone == BOUND2) alt_low = 1*km;
alt_low -= kAltitudeTolerance;
alt_high += kAltitudeTolerance;
if( (pos.Zv() < alt_high) && (pos.Zv() > alt_low ) ) rtn = true;
return rtn;
}
//______________________________________________________________________________
Bool_t LowtranAerosol::CheckWlRange(Double_t wl) const {
//
// true if wl within [200,1536]nm
//
// init
Bool_t rtn = kTRUE;
if(wl < 200*nm) {
wl = 200*nm;
rtn = kFALSE;
}
else if(wl > 1536*nm) {
wl = 1536*nm;
rtn = kFALSE;
}
if(!rtn) {
Msg(EsafMsg::Warning) << "<CheckWlRange> Required wavelength is out of range [200,1536nm] : wl = "<<wl << MsgDispatch;
Msg(EsafMsg::Warning) << "---> thus set to = "<<wl << MsgDispatch;
}
return rtn;
}
//______________________________________________________________________________
Double_t LowtranAerosol::Trans(const EarthVector& startpos,const EarthVector& finalposi,Double_t wl, Double_t& scat) const {
//
// Returns the total transmission along the track defined by startpos and finalpos
// Aerosol zones handled
// Altitude profile is exponential, excepted for low visibility where profile is uniform between [0-1km]
// Local earth-sphericity not taken into account for transmission calculation, to fasten simulation (but used for geometry (impact..))
// Negligeable anyway as boundary layers are pretty thin (< 2km)
//
// init
EarthVector exit(0.,0.,0.), exit1(0.,0.,0.), enter(0.,0.,0.), enter1(0.,0.,0.), direc(0.,0.,0.);
EarthVector finalpos(finalposi);
Double_t exponext(0.), exponscat(0.), mag(0.), magmiddle(0.);
Double_t alt1(0.), alt2(0.), altmiddle(0.), h0(0.);
direc = (finalpos - startpos).Unit();
if(direc.Mag() == 0) return 1.;
Int_t split = -1;
// 1. check if startpos is a valid position
if(startpos.IsUnderSeaLevel()) {
Msg(EsafMsg::Warning) << "<Trans> Starting position is underground, zero returned" << MsgDispatch;
return 0.;
}
// 2. get entering position
enter = GetImpact(startpos,direc);
// if no aerosol along the track
if(enter.Z() == HUGE) return 0;
alt1 = enter.Zv();
// 2.1. get outgoing position for AerZone=ALL
exit = GetZoneImpact(startpos,direc,ALL,"outgoing");
if(exit.Z() == HUGE) Msg(EsafMsg::Warning) << "if enter exists, exit should as well" << MsgDispatch;
alt2 = exit.Zv();
// 3. check if finalpos is a valid position
if(finalpos.IsUnderSeaLevel()) {
Msg(EsafMsg::Warning) << "<Trans> Final position is underground, set to exit position" << MsgDispatch;
finalpos = exit;
}
// 4. look if needs to split aerosol layer in two parts (profile is not exponential in [0-1km] for low visibility)
if(fType == "urban5" || fType == "rural5") {
if(startpos.Zv() >= 1*km && finalpos.Zv() <= 1*km) split = 1; // finalpos is inside BOUND1
else if(finalpos.Zv() >= 1*km && startpos.Zv() <= 1*km) split = 2; // startpos is inside BOUND1
else if(finalpos.Zv() <= 1*km && startpos.Zv() <= 1*km) split = 3; // both pos are inside BOUND1
else split = 0; // both positions are outside BOUND1
}
// 5. get wl dependent factor
Double_t kwlext = Kwl(wl);
Double_t kwlabs = Kwl(wl,kTRUE);
// 6. define path length and related positions
mag = (exit - enter).Mag();
if(mag > (finalpos - enter).Mag()) {
exit = finalpos;
mag = (exit - enter).Mag();
alt2 = exit.Zv();
}
// 7. if needed, get intermediate enter/exit position on top of AerZone=BOUND1
if(split > 0) {
if(split == 1) {
enter1 = GetZoneImpact(startpos,direc,BOUND1);
if(enter1.Z() == HUGE) Msg(EsafMsg::Warning) << "<Trans> In this case (split==1), enter position should exist" << MsgDispatch;
altmiddle = enter1.Zv();
magmiddle = (enter1 - enter).Mag();
}
else if(split == 2) {
exit1 = GetZoneImpact(startpos,direc,BOUND1,"outgoing");
if(exit1.Z() == HUGE) Msg(EsafMsg::Warning) << "<Trans> In this case (split==2), exit position should exist" << MsgDispatch;
altmiddle = exit1.Zv();
magmiddle = (exit1 - enter).Mag();
}
}
// 8. find other parameters for integration : scale height, local zenith angle
h0 = 2.14378*km; // for high visibility (>= 23km)
if(split >= 0) h0 = 0.3972*km; // for visibility = 5km
Double_t angle;
EarthVector temp(startpos);
temp(2) += EarthRadius(); // pos expressed in earth-centered frame -> gives local vertical direction (expressed in MES frame)
angle = direc.Angle(temp); // angle between dir and local vertical
// 9. integrate BETAext*dl along the track
// 9.1. if horizontal track
if(fabs(angle - PiOver2()) < kTolerance) {
if(split == 3) {
exponscat = fSz[0] * (kwlext - kwlabs) * mag;
exponext = fSz[0] * kwlext * mag;
}
else { // should occur rarely
exponscat = fSz[0] * (kwlext - kwlabs) * mag * Exp(- 0.5*(alt1 + alt2)/h0);
exponext = fSz[0] * kwlext * mag * Exp(- 0.5*(alt1 + alt2)/h0);
}
}
// 9.2. else integrate exponential
else {
if(angle > PiOver2()) {
EarthVector dirinvert = -direc;
angle = dirinvert.Angle(temp);
}
if(split == 0) {
exponscat = fSz[1] * (kwlext - kwlabs) * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-alt1/h0) - Exp(-alt2/h0) );
exponext = fSz[1] * kwlext * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-alt1/h0) - Exp(-alt2/h0) );
}
else if(split == 1) {
exponscat = fSz[1] * (kwlext - kwlabs) * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-alt1/h0) - Exp(-altmiddle/h0) );
exponext = fSz[1] * kwlext * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-alt1/h0) - Exp(-altmiddle/h0) );
exponscat += fSz[0] * (kwlext - kwlabs) * (mag - magmiddle);
exponext += fSz[0] * kwlext * (mag - magmiddle);
}
else if(split == 2) {
exponscat = fSz[0] * (kwlext - kwlabs) * magmiddle;
exponext = fSz[0] * kwlext * magmiddle;
exponscat += fSz[1] * (kwlext - kwlabs) * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-altmiddle/h0) - Exp(-alt2/h0) );
exponext += fSz[1] * kwlext * h0 / Cos(angle) * Exp(1*km/h0) * fabs( Exp(-altmiddle/h0) - Exp(-alt2/h0) );
}
// for low visibility and both positions inside BOUND1 : aerosol density is constant
else if(split == 3) {
exponscat = fSz[0] * (kwlext - kwlabs) * mag;
exponext = fSz[0] * kwlext * mag;
}
else {
exponscat = fSz[0] * (kwlext - kwlabs) * h0 / Cos(angle) * fabs( Exp(-alt1/h0) - Exp(-alt2/h0) );
exponext = fSz[0] * kwlext * h0 / Cos(angle) * fabs( Exp(-alt1/h0) - Exp(-alt2/h0) );
}
}
scat = Exp(-exponscat);
return Exp(-exponext);
}
//______________________________________________________________________________
Bool_t LowtranAerosol::RandomScatPos(const EarthVector& posi,const EarthVector& dir,Double_t wl,EarthVector& scatpos) const {
//
// return kFALSE if absorpion, kTRUE if scattering
// get a scattering position along a given track
// handle aerosols zones
// Altitude profile is exponential, excepted for low visibility where profile is uniform between [0-1km]
// Local earth-sphericity not taken into account for transmission calculation, to fasten simulation (but used for geometry (impact..))
// Negligeable anyway as boundary layers are pretty thin (< 2km)
//
// - if scatpos is outside aerosol layer : (0,0,HUGE) returned
// - if scatpos is underground : (0,0,HUGE) returned
//
// init
Bool_t rtn(kFALSE);
TRandom* rndm = EsafRandom::Get();
EarthVector exit1(0.,0.,0.), enter(0.,0.,0.), enter1(0.,0.,0.), outgopos(0.,0.,0.);
EarthVector pos(posi);
Double_t mag(0.), magmiddle(0.);
Double_t alt1(0.), alt2(0.), altmiddle(0.), h0(0.);
EarthVector direc = dir.Unit();
Int_t split = -1; // to deal with zones
// check if pos is a valid position
if(pos.IsUnderSeaLevel()) {
Msg(EsafMsg::Warning) << "<RandomScatPos> Starting position is underground, HUGE vector returned" << MsgDispatch;
scatpos.SetXYZ(0,0,HUGE);
return kTRUE;
}
// check if pos is in aerosol layer (it must be)
// if not, pos is set to enter position
if(!IsInAerosol(pos)) {
Msg(EsafMsg::Warning) << "<RandomScatPos()> Starting position SHOULD be in aerosol layers, pos is set to entering position in layer" << MsgDispatch;
// get entering position
enter = GetImpact(pos,direc);
// if no aerosol along the track
if(enter.Z() == HUGE) {
Msg(EsafMsg::Warning) << "<RandomScatPos()> STRANGE CASE : randomscatpos asked, but track does not seem to intersect aerosol layer" << MsgDispatch;
scatpos.SetXYZ(0,0,HUGE);
return kTRUE;
}
pos = enter;
}
// 0. look if needs to split aerosol layer in two parts (profile is not exponential in [0-1km] for low visibility)
outgopos = GetImpact(pos,dir,"outgoing");
if(fType == "urban5" || fType == "rural5") {
if(pos.Zv() >= 1*km && outgopos.Zv() <= 1*km) split = 1; // outgopos is inside BOUND1
else if(outgopos.Zv() >= 1*km && pos.Zv() <= 1*km) split = 2; // startpos is inside BOUND1
else if(outgopos.Zv() <= 1*km && pos.Zv() <= 1*km) split = 3; // both pos are inside BOUND1
else split = 0; // both positions are outside BOUND1
}
// 1. determine if scattering process or absorption process
Double_t kwlext = Kwl(wl);
Double_t kwlabs = Kwl(wl,kTRUE);
Double_t abs = kwlabs / kwlext;
Double_t scat = 1 - abs;
Double_t K = kwlabs;
if(scat > rndm->Rndm()) {
K = kwlext - kwlabs;
rtn = kTRUE;
}
// 2. if needed, get intermediate enter/exit position on top of AerZone=BOUND1
if(split > 0) {
if(split == 1) {
enter1 = GetZoneImpact(pos,direc,BOUND1);
if(enter1.Z() == HUGE) Msg(EsafMsg::Warning) << "<RandomScatPos> In this case (split==1), enter position should exist" << MsgDispatch;
altmiddle = enter1.Zv();
magmiddle = (enter1 - pos).Mag();
}
else if(split == 2) {
exit1 = GetZoneImpact(pos,direc,BOUND1,"outgoing");
if(exit1.Z() == HUGE) {
Msg(EsafMsg::Warning) << "<RandomScatPos> In this case (split==2), exit position should exist" << MsgDispatch;
Msg(EsafMsg::Warning) << "pos = "<<pos << MsgDispatch;
Msg(EsafMsg::Warning) << "direc = "<<direc << MsgDispatch;
// find local zenith angle
Double_t anglee;
EarthVector tempp(pos);
tempp(2) += EarthRadius(); // pos expressed in earth-centered frame -> gives local vertical direction (expressed in MES frame)
anglee = direc.Angle(tempp); // angle between dir and local vertical
Msg(EsafMsg::Warning) << "local zenith angle/halfpi = "<<anglee/PiOver2()<<" rad, angle = "<<anglee*RadToDeg()<<" deg" << MsgDispatch;
Msg(EsafMsg::Warning) << "<RandomScatPos> Fake aerosol scat position returned" << MsgDispatch;
scatpos.SetXYZ(0,0,HUGE);
return rtn;
}
altmiddle = exit1.Zv();
magmiddle = (exit1 - pos).Mag();
}
}
// 3. find other parameters for integration
// scale height, local zenith angle, pos altitude
h0 = 2.14378*km; // for high visibility (>= 23km)
if(split >= 0) h0 = 0.3972*km; // for visibility = 5km
Double_t angle;
EarthVector temp(pos);
temp(2) += EarthRadius(); // pos expressed in earth-centered frame -> gives local vertical direction (expressed in MES frame)
angle = direc.Angle(temp); // angle between dir and local vertical
alt1 = pos.Zv();
// 4. generate random number to sample distribution function : Integral(Beta(s).ds) = -ln(1 - RNDM)
Double_t RNDM = -Log(1 - rndm->Rndm());
// 5. find path length, thus determining position of interaction
Bool_t upward = kTRUE;
// 5.1. if horizontal track
if(fabs(angle - PiOver2()) < kTolerance) {
if(split == 3) mag = RNDM / (fSz[0] * K);
else mag = RNDM / (fSz[0] * K * Exp(-alt1/h0)); // should occur rarely
}
// 5.2. else deals with zones
else {
if(angle > PiOver2()) {
EarthVector dirinvert = -direc;
angle = dirinvert.Angle(temp);
upward = kFALSE;
}
if(split == 0) {
if(upward) {
if(RNDM > (fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * Exp(-alt1/h0))) {
// means no valid position found
scatpos.SetXYZ(0,0,HUGE);
return rtn;
}
mag = -h0 / Cos(angle) * Log(1 - RNDM / (fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * Exp(-alt1/h0)));
}
else mag = h0 / Cos(angle) * Log(1 + RNDM / (fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * Exp(-alt1/h0)));
}
else if(split == 1) { // upward is false
// get a position assuming everything occurs in BOUND2
mag = h0 / Cos(angle) * Log(1 + RNDM / (fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * Exp(-alt1/h0)));
// if found position altitude is below altmiddle, needs to consider BOUND1 too
alt2 = alt1 - mag * Cos(angle);
if(alt2 < altmiddle) {
// substract first part of integration until altmiddle (in BOUND2)
RNDM -= fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * (Exp(-altmiddle/h0) - Exp(-alt1/h0));
mag = magmiddle;
mag += RNDM / (fSz[0] * K);
}
}
else if(split == 2) { // upward is true
// get a position assuming everything occurs in BOUND1
mag = RNDM / (fSz[0] * K);
// if found position altitude is above altmiddle, needs to consider BOUND2 too
alt2 = mag * Cos(angle) + alt1;
if(alt2 > altmiddle) {
// substract first part of integration (in BOUND1)
RNDM -= fSz[0] * K * magmiddle;
mag = magmiddle;
mag += -h0 / Cos(angle) * Log(1 - RNDM / (fSz[1] * K * h0 / Cos(angle) * Exp(1*km/h0) * Exp(-altmiddle/h0)));
}
}
// for low visibility and both positions inside BOUND1 : aerosol density is constant
else if(split == 3) mag = RNDM / (fSz[0] * K);
else {
if(upward) {
if(RNDM > (fSz[0] * K * h0 / Cos(angle) * Exp(-alt1/h0))) {
// means no valid position found
scatpos.SetXYZ(0,0,HUGE);
return rtn;
}
mag = -h0 / Cos(angle) * Log(1 - RNDM / (fSz[0] * K * h0 / Cos(angle) * Exp(-alt1/h0)));
}
else mag = h0 / Cos(angle) * Log(1 + RNDM / (fSz[0] * K * h0 / Cos(angle) * Exp(-alt1/h0)));
}
}
// 6. set position of interaction and check if it a valid position
scatpos = pos + mag*direc;
if(!IsInAerosol(scatpos)) scatpos.SetXYZ(0,0,HUGE); // handle "under sea level" as well
return rtn;
}
//______________________________________________________________________________
Double_t LowtranAerosol::PhaseFunction(Double_t theta,Double_t wl) const {
//
// Mie scattering phase function
// Normalized when integrated over full solid angle
//
// check if wl is valid
CheckWlRange(wl);
// find bounds for interpolation
Double_t G1(0.), G2(0.), F(0.);
Int_t i=0;
for(i=0; i < 6; i++) {
if(wl >= fWlrange[i] && wl <= fWlrange[i+1]) break;
}
G1 = (fG1[i+1] - fG1[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fG1[i];
G2 = (fG2[i+1] - fG2[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fG2[i];
F = (fF[i+1] - fF[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fF[i];
// calculate value
Double_t rtn(0.);
Double_t HG1 = (1 - G1*G1) / pow(1 - 2*G1*Cos(theta) + G1*G1,1.5);
Double_t HG2 = (1 - G2*G2) / pow(1 - 2*G2*Cos(theta) + G2*G2,1.5);
rtn = 1./(4*Pi()) * (F*HG1 + (1 - F)*HG2);
return rtn;
}
//______________________________________________________________________________
Double_t LowtranAerosol::PhaseFunction(const EarthVector& incom_dir,const EarthVector& outgo_dir,Double_t wl) const {
//
// Mie scattering phase function
// Normalized when integrated over full solid angle
//
Double_t theta = outgo_dir.Angle(incom_dir);
return PhaseFunction(theta,wl);
}
//______________________________________________________________________________
Double_t LowtranAerosol::GetMaxPhaseFunction(Double_t wlmin,Double_t wlmax) const {
//
// returns maximum value of currently used scattering phase function
// (phase function is normalized when integrated over full solid angle)
//
Double_t rtn(0.), temp(0.);
for(Int_t i=0; i < 7; i++) {
temp = PhaseFunction(0.,fWlrange[i]);
if(temp > rtn && fWlrange[i] > wlmin && fWlrange[i] < wlmax) rtn = temp;
}
temp = PhaseFunction(0.,wlmin);
if(temp > rtn) rtn = temp;
temp = PhaseFunction(0.,wlmax);
if(temp > rtn) rtn = temp;
return rtn;
}
//______________________________________________________________________________
void LowtranAerosol::RandomDir(EarthVector& dir,Double_t wl) const {
//
// 'dir' is the incoming direction -> this method set dir to the scattering outgoing direction
// get a random direction, sampling clouds scattering phase function
//
Double_t theta(0);
Double_t u1(0.), u2(0.);
TRandom* rndm = EsafRandom::Get();
// check if wl is valid
CheckWlRange(wl);
// find bounds for interpolation
Double_t G1(0.), G2(0.), F(0.);
Int_t i=0;
for(i=0; i < 6; i++) {
if(wl >= fWlrange[i] && wl <= fWlrange[i+1]) break;
}
G1 = (fG1[i+1] - fG1[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fG1[i];
G2 = (fG2[i+1] - fG2[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fG2[i];
F = (fF[i+1] - fF[i]) / (fWlrange[i+1] - fWlrange[i]) * (wl - fWlrange[i]) + fF[i];
// get random theta, assuming isotropic distribution in phi
// theta is the angle between incoming and outgoing directions
// von neumann sampling -> much faster than using TF1
Double_t HG1 = 0.;
Double_t HG2 = 0.;
Double_t pfvalue = 0.;
do {
u1 = -1 + rndm->Rndm()*2;
HG1 = (1 + G1) / (1 - G1) / (1 - G1);
HG2 = (1 + G2) / (1 - G2) / (1 - G2);
u2 = rndm->Rndm() * 0.5 * (F*HG1 + (1 - F)*HG2);
HG1 = (1 - G1*G1) / pow(1 - 2*G1*u1 + G1*G1,1.5);
HG2 = (1 - G2*G2) / pow(1 - 2*G2*u1 + G2*G2,1.5);
pfvalue = 0.5 * (F*HG1 + (1 - F)*HG2);
} while(u2 >= pfvalue);
theta = acos(u1);
// then get a random direction
EarthVector v = dir.Orthogonal();
v.Rotate(rndm->Rndm()*2*Pi(),dir);
EarthVector axis = v.Cross(dir);
dir.Rotate(theta,axis);
}
//_____________________________________________________________________________________________________
EarthVector LowtranAerosol::GetZoneImpact(const EarthVector& pos, const EarthVector& dir, AerZone zone, string opt) const {
//
// Impact on clouds top surface form a given (position,direction) pair in atmosphere
// if no impact, (0,0,HUGE) is returned
// upward direction is handled
//
// Aerosol considered geometry is defined using "zone" argument
//
// if opt == "default" -> incoming impact returned
// if pos is within clouds, pos is returned
// if opt == "outgoing" -> outgoing impact returned
//
// init
EarthVector rtn(0,0,HUGE);
Bool_t isinaerosol, downward;
Double_t Altimpact(0);
EarthVector direc = dir.Unit();
Double_t topaltitude = 2*km;
Double_t bottomaltitude = 0.;
if(zone == BOUND1) topaltitude = 1*km;
if(zone == BOUND2) bottomaltitude = 1*km;
// says if dir is going up/downward
downward = GoingDownward(pos,direc);
// says if pos is within boundary layers (altitude round-off effects < 1mm are handled)
isinaerosol = IsInAerosolZone(pos,zone);
// set the impact surface altitude according to photon present situation (up/downward, above/below, enter/exit aerosol)
if(isinaerosol) {
if(opt == "default") return pos;
else if(opt == "outgoing") {
if(!downward) Altimpact = topaltitude;
else Altimpact = bottomaltitude;
}
else Msg(EsafMsg::Panic) << "<GetZoneImpact()> wrong argument (internal pos)" << MsgDispatch;
}
else {
// if not going toward aerosol layer
if(pos.Zv() < bottomaltitude && downward) return rtn;
if(pos.Zv() > topaltitude && !downward) return rtn;
if(opt == "outgoing" && !downward || opt == "default" && downward) Altimpact = topaltitude;
else if(opt == "outgoing" && downward || opt == "default" && !downward) Altimpact = bottomaltitude;
else Msg(EsafMsg::Panic) << "<GetZoneImpact()> wrong argument (external pos)" << MsgDispatch;
}
// now calculate impact
Double_t mag(0);
/* flat earth
if(fabs(direc.Z()) < kTolerance) rtn.SetXYZ(0,0,HUGE);
else {
mag = (Altimpact - pos.Z()) / direc.Z();
if(mag < -kAltitudeTolerance) rtn.SetXYZ(0,0,HUGE);
else rtn = pos + mag*direc;
}
if(mag < -kAltitudeTolerance) Msg(EsafMsg::Warning) << "PB in impact calculation GetCloudImpact(), mag = ["<<mag<< MsgDispatch;
*/
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) {
// if direction is nearly horizontal, try the other layer bound
Double_t anglee;
EarthVector tempp(pos);
tempp(2) += EarthRadius();
anglee = dir.Angle(tempp);
if(fabs(anglee - PiOver2()) < 0.1) {
if(Altimpact == topaltitude) Altimpact = bottomaltitude;
else Altimpact = topaltitude;
mag = 0;
b = pos*direc + direc.Z()*EarthRadius();
c = pos.Mag2() + 2*EarthRadius()*(pos.Z() - Altimpact) - pow(Altimpact,2);
p = findRoots(1.,2*b,c);
}
else return rtn;
}
if(p.first == 0) return rtn;
else if(p.first == 1) mag = p.second[0];
else if(p.first ==2) {
if((p.second[0]+kAltitudeTolerance) * (p.second[1]+kAltitudeTolerance) < 0) mag = max(p.second[0],p.second[1]);
else mag = min(p.second[0],p.second[1]);
}
if(mag < -kAltitudeTolerance) Msg(EsafMsg::Warning) << "PB in impact calculation GetZoneImpact(), [min,max] = ["<<mag<<", "<< max(p.second[0],p.second[1])<<"]" << MsgDispatch;
rtn = pos + mag*direc;
return rtn;
}
//_____________________________________________________________________________________________________
EarthVector LowtranAerosol::GetImpact(const EarthVector& pos, const EarthVector& dir, string opt) const {
//
// Impact on aerosol top surface form a given (position,direction) pair in atmosphere
// if no impact, (0,0,HUGE) is returned
// upward direction is handled
//
// if opt == "default" -> incoming impact returned
// if pos is within clouds, pos is returned
// if opt == "outgoing" -> outgoing impact returned
//
return GetZoneImpact(pos,dir,ALL,opt);
}
//_____________________________________________________________________________
Bool_t LowtranAerosol::GoingDownward(const EarthVector& pos,const EarthVector& direc) 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;
/* flat earth
if(direc.Z() > 0) return false;
else return true;
*/
EarthVector dir = direc.Unit();
// 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 <= PiOver2()) rtn = false;
else rtn = true;
return rtn;
}
//_____________________________________________________________________________
void LowtranAerosol::Reset() {
//
// reset clodus features
//
// reset current model tables
for(Int_t i=0; i < 7; i++) {
if(i < 3) fSz[i] = 0;
fKwl_ext[i] = 0;
fKwl_abs[i] = 0;
fG1[i] = 0;
fG2[i] = 0;
fF[i] = 0;
}
// build new model tables
Build();
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.