// ESAF : Euso Simulation and Analysis Framework
// $Id: LowtranRadiativeProcessesCalculator.cc,v 1.58 2005/11/03 12:01:33 moreggia Exp $
// Sylvain Moreggia created Jan, 29 2004
#include "LowtranRadiativeProcessesCalculator.hh"
#include <vector>
#include "Config.hh"
#include "SinglePhoton.hh"
#include "BunchOfPhotons.hh"
#include "EarthVector.hh"
#include "EConst.hh"
#include "Atmosphere.hh"
#include "Ground.hh"
#include "DetectorGeometry.hh"
#include "ListPhotonsInAtmosphere.hh"
ClassImp(LowtranRadiativeProcessesCalculator)
using EConst::EarthRadius;
using namespace sou;
using namespace TMath;
//____________________________________________________________________________________
//
// for lowtran interface
//
// lowtran output container
// places of declarations are of importance, don't change them
struct LowtranTransferOutput {
Float_t R,Altitude,Angle;
Float_t WaveLength[100];
Float_t Transmission[100];
Float_t Ozone[100];
Float_t Rayl[100];
Float_t Aerosol[100];
};
// Fortran interfaces
extern "C"{
#define lwtrans lwtrans_
#define lwtrans_v2 lwtrans_v2_
// calculate transmission between a point and EUSO
void lwtrans(Int_t* init, Float_t* altitude, Float_t* angle, Float_t* trans);
// calculate transmission between two points, expressed in lowtran frame
void lwtrans_v2(Int_t* init, Float_t* h1, Float_t* h2, Float_t* theta);
}
extern "C"{
#define lowtrantransfer tuple_
struct LowtranTransferOutput tuple_;
}
void InitLOWTRAN() {
// lowtran initialisation
Float_t trans[100];
Int_t flag(0);
Float_t h(0),theta(0);
lwtrans(&flag,&h,&theta,trans);
}
//________________________________________________________________________________________
LowtranRadiativeProcessesCalculator::LowtranRadiativeProcessesCalculator() : RadiativeProcessesCalculator() {
//
// ctor
//
// Configure
fStep_along_track = 0;
fStep_detector = 0;
fFlag = false;
fInitPos.SetXYZ(0,0,HUGE);
fNbwl = 0;
fName = "lowtran";
fNbSinglesFarFromTrack = 0.;
fDetectorAlt = 0.;
fVertical = false;
}
//________________________________________________________________________________________
LowtranRadiativeProcessesCalculator::~LowtranRadiativeProcessesCalculator() {
//
// dtor
//
Reset();
CleanBuffer("todetector");
fStep_detector = 0.;
fDetectorAlt = 0.;
fVertical = false;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::MakeVerticalDatacard(const Ground* ground, Double_t detector_alt) {
//
// build a VERTICAL transmission map (fToDetectorCoeff) from positions uniformly distributed in altitude until detector
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Preprocessing Singles propagation" << MsgDispatch;
#endif
// datacard built only once for a given run
if(fToDetectorCoeff.size()) {
Msg(EsafMsg::Info) <<"Datacard already built" << MsgDispatch;
return;
}
SetVertical();
// initializations
LowtranTransferOutput output;
CleanBuffer("todetector");
ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","LowtranRadiativeProcessesCalculator");
fStep_detector = pConfig->GetNum("LowtranRadiativeProcessesCalculator.fStep_detector")*km;
Double_t TOA = 50*km;
fToDetectorCoeff.reserve(UInt_t(TOA/fStep_detector) + 1);
EarthVector start(1);
fDetectorAlt = detector_alt;
EarthVector detector(1);
detector.SetXYZ(0,0,fDetectorAlt);
vector<vector<Double_t> >* Transmis;
// loop makes calculations then fills the map
// an array for each altitude
for(UInt_t k=0; k < UInt_t(TOA/fStep_detector) + 1; k++) {
Transmis = new vector<vector<Double_t> >(4);
fToDetectorCoeff.push_back(Transmis);
start.SetXYZ(0,0,k*fStep_detector);
CallLowtran100(start,detector,output);
for(UInt_t i=0; i<100; i++) {
if(!k) fWlArray[i] = output.WaveLength[i];
(*fToDetectorCoeff[k])[0].push_back(output.Transmission[i]);
(*fToDetectorCoeff[k])[1].push_back(output.Rayl[i]);
(*fToDetectorCoeff[k])[2].push_back(output.Ozone[i]);
(*fToDetectorCoeff[k])[3].push_back(output.Aerosol[i]);
}
}
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::MakeTrackToDetectorDatacard(const Ground* ground, const EarthVector& detector, const ListPhotonsInAtmosphere& track) {
//
// Copy the light track and
// build a transmission map (fToDetectorCoeff), from positions along the track until detector
// used in BunchRadiativeTransfer::PropagationOfSingles() in RT 'decoupled' mode
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Preprocessing Singles propagation" << MsgDispatch;
#endif
SetVertical(false);
// initializations
LowtranTransferOutput output;
CleanBuffer("todetector");
fTrack = &track;
UInt_t tr_size = track.GetNbTrackSteps();
fToDetectorCoeff.reserve(tr_size + 1); // +1 for photons on ground
EarthVector start(1);
vector<vector<Double_t> >* Transmis;
// loop makes calculations then fills the map
// an array for each step
for(UInt_t k=0; k < tr_size; k++) {
Transmis = new vector<vector<Double_t> >(4);
fToDetectorCoeff.push_back(Transmis);
start = track.GetTrackStep(k);
CallLowtran100(start,detector,output);
for(UInt_t i=0; i<100; i++) {
if(!k) fWlArray[i] = output.WaveLength[i];
(*fToDetectorCoeff[k])[0].push_back(output.Transmission[i]);
(*fToDetectorCoeff[k])[1].push_back(output.Rayl[i]);
(*fToDetectorCoeff[k])[2].push_back(output.Ozone[i]);
(*fToDetectorCoeff[k])[3].push_back(output.Aerosol[i]);
}
}
// last step with init point at the track impact
Transmis = new vector<vector<Double_t> >(4);
fToDetectorCoeff.push_back(Transmis);
fImpact = ground->GetImpact(track.GetTrackStep(0),track.GetTrackStep(tr_size-1) - track.GetTrackStep(0));
start = fImpact;
if(start.Z() != HUGE) {
CallLowtran100(start,detector,output);
for(UInt_t i=0; i<100; i++) {
(*fToDetectorCoeff[tr_size])[0].push_back(output.Transmission[i]);
(*fToDetectorCoeff[tr_size])[1].push_back(output.Rayl[i]);
(*fToDetectorCoeff[tr_size])[2].push_back(output.Ozone[i]);
(*fToDetectorCoeff[tr_size])[3].push_back(output.Aerosol[i]);
}
}
else {
for(UInt_t i=0; i<100; i++) {
(*fToDetectorCoeff[tr_size])[0].push_back(0);
(*fToDetectorCoeff[tr_size])[1].push_back(0);
(*fToDetectorCoeff[tr_size])[2].push_back(0);
(*fToDetectorCoeff[tr_size])[3].push_back(0);
}
Msg(EsafMsg::Info) << "PreProcessSingles: NO IMPACT on ground -> Trans=0" << MsgDispatch;
}
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::Trans(const SinglePhoton& s, const EarthVector& finalpos, Double_t* coeff) const {
//
// Transmission between photon position and finalpos
Double_t wl = s.Wl();
UInt_t Wli;
Double_t rtn[4];
Double_t OD = Atmosphere::Get()->GetClouds()->TotalOD(s.Pos(),s.Dir());
// if fToDetectorCoeff card doesn't exist
if(fToDetectorCoeff.size() == 0) {
LowtranTransferOutput output;
CallLowtran100(s.Pos(),finalpos,output);
// lambda research (array is in decreasing order)
UInt_t nabove, nbelow, middle;
nabove = 101;
nbelow = 0;
while(nabove-nbelow > 1) {
middle = (nabove+nbelow)/2;
if (wl == output.WaveLength[middle-1]) Wli = middle-1;
if (wl > output.WaveLength[middle-1]) nabove = middle;
else nbelow = middle;
}
Wli = nbelow-1;
rtn[0] = InterpolWl("linear",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Transmission[Wli],output.Transmission[Wli+1]);
rtn[1] = InterpolWl("rayleigh",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Rayl[Wli],output.Rayl[Wli+1]);
rtn[2] = InterpolWl("linear",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Ozone[Wli],output.Ozone[Wli+1]);
rtn[3] = InterpolWl("linear",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Aerosol[Wli],output.Aerosol[Wli+1]);
if(coeff) for(Int_t i=0; i<4; i++) coeff[i] = rtn[i];
return rtn[0] * exp(-OD);
}
else if(fVertical) {
//
// Principle of interpolation from datacard :
// - photon altitude tells the datacard element to be used --> gives an averaged interaction length
// - then apply 1/cos(theta), theta being the angle of (EUSO - photon_pos) with vertical axis
// - VALID ONLY WHEN TRACK TO DETECTOR IS NOT STRONGLY INCLINED
//
// init
Double_t trans1[4], trans2[4], coeff1[4], coeff2[4];
for(Int_t i=0; i<4; i++) {
rtn[i] = 0.;
trans1[i] = 0.;
trans2[i] = 0.;
coeff1[i] = 0.;
coeff2[i] = 0.;
}
UInt_t Stepi, Stepiplus1;
// local direction to go to detector
Double_t theta;
EarthVector temp(s.Pos());
temp(2) += EarthRadius(); // pos expressed in earth-centered frame -> gives local vertical direction expressed in MES frame
theta = (finalpos - s.Pos()).Angle(temp); // angle with local vertical
// warning if theta > 60deg
if(theta > (Pi() / 3.))
Msg(EsafMsg::Warning) <<"<Trans(SinglePhoton)> Track to detector strongly inclined :"<<theta*RadToDeg()<<" deg, Interpolation not valid anymore" << MsgDispatch;
// lambda research (array is in decreasing order)
UInt_t nabove, nbelow, middle;
nabove = 101;
nbelow = 0;
while(nabove-nbelow > 1) {
middle = (nabove+nbelow)/2;
if (wl == fWlArray[middle-1]) Wli = middle-1;
if (wl > fWlArray[middle-1]) nabove = middle;
else nbelow = middle;
}
Wli = nbelow-1;
// step research
Stepi = UInt_t(s.Pos().Zv() / fStep_detector);
if(Stepi >= (fToDetectorCoeff.size() - 1)) {
Stepi = fToDetectorCoeff.size() - 1;
Stepiplus1 = Stepi;
}
else Stepiplus1 = Stepi + 1;
for(UInt_t j=0; j<4; j++) {
// interpolation in lambda for the two step bounds
coeff1[j] = (*fToDetectorCoeff[Stepi])[j][Wli];
if(Wli == 99) coeff2[j] = coeff1[j];
else coeff2[j] = (*fToDetectorCoeff[Stepi])[j][Wli+1];
if(j == 1) trans1[j] = InterpolWl("rayleigh",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
else trans1[j] = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
//
coeff1[j] = (*fToDetectorCoeff[Stepiplus1])[j][Wli];
if(Wli == 99) coeff2[j] = coeff1[j];
else coeff2[j] = (*fToDetectorCoeff[Stepiplus1])[j][Wli+1];
if(j == 1) trans2[j] = InterpolWl("rayleigh",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
else trans2[j] = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
if(fDetectorAlt == 0) rtn[j] = trans1[j] * pow(trans2[j] / trans1[j], (s.Pos().Zv() - Stepi*fStep_detector) / fStep_detector); //TOFIX : to be checked when relevant
else rtn[j] = trans2[j] * pow(trans1[j] / trans2[j], (Stepiplus1*fStep_detector - s.Pos().Zv()) / fStep_detector );
rtn[j] = pow(rtn[j],1./fabs(cos(theta))); // fabs() for ground detector case
}
for(Int_t i=0; i<4; i++) coeff[i] = rtn[i];
return rtn[0] * exp(-OD);
}
else {
//
// Principle of interpolation from datacard :
// - find a reference transmission (calcuated from ground for reflected -- interpolated within a step for scattered)
// - then get the interaction length and apply it to the real travel
//
// init
Double_t trans1[4], trans2[4], coeff1[4], coeff2[4];
for(Int_t i=0; i<4; i++) {
rtn[i] = 0.;
trans1[i] = 0.;
trans2[i] = 0.;
coeff1[i] = 0.;
coeff2[i] = 0.;
}
UInt_t Stepi;
Double_t travel, entrystep, endstep, ontrack;
travel = (finalpos - s.Pos()).Mag();
const ListPhotonsInAtmosphere& track = *fTrack;
UInt_t tr_size = track.GetNbTrackSteps();
// lambda research (array is in decreasing order)
UInt_t nabove, nbelow, middle;
nabove = 101;
nbelow = 0;
while(nabove-nbelow > 1) {
middle = (nabove+nbelow)/2;
if (wl == fWlArray[middle-1]) Wli = middle-1;
if (wl > fWlArray[middle-1]) nabove = middle;
else nbelow = middle;
}
Wli = nbelow-1;
// step research
EarthVector projec;
Double_t mag_projec = 0;
Stepi = track.TrackStepIndex(s.Pos(),projec);
mag_projec = (projec - track.GetTrackStep(0)).Mag();
// reflected photons are treated differently
// scattered photons with position further than the end of last fTrack step are treated here as reflected ones
if(s.Status() == Reflected || Stepi == (tr_size-1)) {
#ifdef DEBUG
if(s.Status() != Reflected) {
Msg(EsafMsg::Debug) <<"scattered photons position further than track end" <<MsgDispatch;
fNbSinglesFarFromTrack++;
Msg(EsafMsg::Debug) <<"Nb of Singles in that case = "<<fNbSinglesFarFromTrack << MsgDispatch;
}
#endif
if(fImpact.Z() == HUGE) Msg(EsafMsg::Warning) <<"SinglePhoton reflected but No track impact" <<MsgDispatch;
for(UInt_t j=0; j<4; j++) {
coeff1[j] = (*fToDetectorCoeff[fToDetectorCoeff.size()-1])[j][Wli];
if(Wli == 99) coeff2[j] = coeff1[j];
else coeff2[j] = (*fToDetectorCoeff[fToDetectorCoeff.size()-1])[j][Wli+1];
if(j == 1) rtn[j] = InterpolWl("rayleigh",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]); // reference transmission
else rtn[j] = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]); // reference transmission
endstep = (finalpos - fImpact).Mag();
rtn[j] = pow(rtn[j],travel/endstep);
}
if(coeff) for(Int_t i=0; i<4; i++) coeff[i] = rtn[i];
return rtn[0] * exp(-OD);
}
// for backscattered photons : get reference transmission
endstep = (track.GetTrackStep(0) - track.GetTrackStep(Stepi+1)).Mag();
entrystep = (track.GetTrackStep(0) - track.GetTrackStep(Stepi)).Mag();
ontrack = (finalpos - projec).Mag();
for(UInt_t j=0; j<4; j++) {
// interpolation in lambda for the two step bounds
coeff1[j] = (*fToDetectorCoeff[Stepi])[j][Wli];
if(Wli == 99) coeff2[j] = coeff1[j];
else coeff2[j] = (*fToDetectorCoeff[Stepi])[j][Wli+1];
if(j == 1) trans1[j] = InterpolWl("rayleigh",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
else trans1[j] = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
//
coeff1[j] = (*fToDetectorCoeff[Stepi+1])[j][Wli];
if(Wli == 99) coeff2[j] = coeff1[j];
else coeff2[j] = (*fToDetectorCoeff[Stepi+1])[j][Wli+1];
if(j == 1) trans2[j] = InterpolWl("rayleigh",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
else trans2[j] = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1[j],coeff2[j]);
rtn[j] = (trans2[j] - trans1[j])/(endstep - entrystep) * (mag_projec - entrystep) + trans1[j]; // reference transmission
rtn[j] = pow(rtn[j],travel/ontrack);
}
if(coeff) for(Int_t i=0; i<4; i++) coeff[i] = rtn[i];
return rtn[0] * exp(-OD);
}
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::PreProcess(const BunchOfPhotons& b) {
//
// Preprocessing for bunch transmission along the track (for ckov backscattering simulation)
// Initialization of tables (fAlongTrackCoeff) for a step by step transfer process of a given bunchofphotons
//
#ifdef DEBUG
Msg(EsafMsg::Debug) <<"Bunch Preprocessing.." << MsgDispatch;
#endif
// initializations
ConfigFileParser* pConfig = Config::Get()->GetCF("RadiativeTransfer","BunchRadiativeTransfer");
fStep_along_track = pConfig->GetNum("BunchRadiativeTransfer.fStep_along_track")*km;
fNbwl = b.GetWlSpectrum().GetAxis().GetNbins();
LowtranTransferOutput output;
Double_t wl;
Int_t j;
fInitPos = b.GetPos();
CleanBuffer("alongtrack");
vector<vector<Double_t> >* Transmis;
// calculation of lowtran number of steps, and step length tunning
EarthVector pos, impact, increm;
impact = b.GetFinalImpact();
increm = impact - fInitPos;
Double_t Nb = ceil(increm.Mag()/fStep_along_track);
if(Nb == 0) Msg(EsafMsg::Warning) << "BunchPreProcess() : Nb should not be equal to zero" << MsgDispatch;
if(Nb > 1000) Nb = 1000;
fStep_along_track = increm.Mag() / Nb;
// loop makes calculations then fills each vector (one vector for one step)
for(Int_t k=0; k<Nb; k++) {
increm.SetMag((k+1)*fStep_along_track);
pos = fInitPos + increm;
Transmis = new vector<vector<Double_t> >(4);
fAlongTrackCoeff.push_back(Transmis);
CallLowtran100(fInitPos,pos,output);
Int_t lim = 99;
// fill the fAlongTrackCoeff[k]
for(Int_t i=0; i<fNbwl; i++) {
wl = b.GetWlSpectrum().GetLambdaEntry(i);
for(j=0; j<lim; j++){
if(output.WaveLength[j+1] <= wl && wl <= output.WaveLength[j]){
(*fAlongTrackCoeff[k])[0].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Transmission[j],output.Transmission[j+1]));
(*fAlongTrackCoeff[k])[1].push_back(InterpolWl("rayleigh",wl,output.WaveLength[j],output.WaveLength[j+1],output.Rayl[j],output.Rayl[j+1]));
(*fAlongTrackCoeff[k])[2].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Ozone[j],output.Ozone[j+1]));
(*fAlongTrackCoeff[k])[3].push_back(InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Aerosol[j],output.Aerosol[j+1]));
break;
}
}
if(j == lim) Msg(EsafMsg::Panic) << "PB in searching lambda in LowtranRadiativeProcessCalculator" << MsgDispatch;
lim = j+1; // to optimize following lambda researches
}
}
fFlag = true;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::Trans(const BunchOfPhotons& b,const EarthVector& nextpos,vector<Double_t*>& coeff) const {
//
// Total transmission between bunch position and a final position
// NB: for a considered bunch, step for running lowtran is independant of step defined for bunch propagation
// (lowtran step is larger.. cf computing time consumption)
//
Int_t j,k,m;
LowtranTransferOutput output;
// if single step process
if(!fFlag) {
// lowtran called for transfer from bunch position to nextpos (one step process)
if((nextpos - b.GetNextImpact()).Mag() < TOLERANCE) {
CallLowtran100(b.GetPos(),nextpos,output);
// then fill coeff
Int_t lim = 99;
Double_t wl;
for(Int_t i=0; i<b.GetWlSpectrum().GetAxis().GetNbins(); i++) {
wl = b.GetWlSpectrum().GetLambdaEntry(i);
for(j=0; j<lim; j++){
if(output.WaveLength[j+1] <= wl && wl <= output.WaveLength[j]){
coeff[0][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Transmission[j],output.Transmission[j+1]);
coeff[1][i] = InterpolWl("rayleigh",wl,output.WaveLength[j],output.WaveLength[j+1],output.Rayl[j],output.Rayl[j+1]);
coeff[2][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Ozone[j],output.Ozone[j+1]);
coeff[3][i] = InterpolWl("linear",wl,output.WaveLength[j],output.WaveLength[j+1],output.Aerosol[j],output.Aerosol[j+1]);
break;
}
}
if(j == lim) Msg(EsafMsg::Panic) << "PB in searching lambda in LowtranRadiativeProcessCalculator" << MsgDispatch;
lim = j+1; // to optimize following lambda researches
}
}
else {
Msg(EsafMsg::Warning) << "<PreProcess()> method has to be called before step by step transmission" << MsgDispatch;
Msg(EsafMsg::Warning) << "May rarely occur if clouds exist" <<MsgDispatch;
}
}
// return transmission between bunch position and nextpos, using interpolation
for(k=0; k<4; k++) {
for(m=0; m<fNbwl; m++)
coeff[k][m] = Interpolate(k,m,nextpos)/Interpolate(k,m,b.GetPos());
}
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::Interpolate(Int_t i,Int_t j,const EarthVector& pos) const {
//
// Interpolate the specified transmission coefficient fAlongTrackCoeff[][i][j]
// to get the value at pos (on the transfer track)
//
Double_t rtn, Ti, Tf;
Double_t mag = (pos - fInitPos).Mag();
if(mag/fStep_along_track == Double_t(fAlongTrackCoeff.size())) return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
if(mag > fAlongTrackCoeff.size()*fStep_along_track) {
if(fabs(mag - fAlongTrackCoeff.size()*fStep_along_track) > 1*mm) {
#ifdef DEBUG
//Msg(EsafMsg::Debug) << "<Interpolate> mag is too large -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep_along_track) << MsgDispatch;
if(fabs(mag - fAlongTrackCoeff.size()*fStep_along_track) > 5e4*mm)
Msg(EsafMsg::Warning) << "<Interpolate> mag is TOO large (>50m) -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep_along_track) << MsgDispatch;
#endif
mag = fAlongTrackCoeff.size()*fStep_along_track;
}
else return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
}
Double_t iter = floor(mag/fStep_along_track);
UInt_t n = UInt_t(iter);
if(n > (fAlongTrackCoeff.size() - 1)) n = fAlongTrackCoeff.size() - 1;
if(n == 0) {
Ti = 1;
Tf = (*fAlongTrackCoeff[n])[i][j];
}
else {
Ti = (*fAlongTrackCoeff[n-1])[i][j];
Tf = (*fAlongTrackCoeff[n])[i][j];
}
// exponential interpolation
rtn = Ti*exp((mag - n*fStep_along_track)/fStep_along_track * log(Tf/Ti));
// linear interpolation
//rtn = Ti + (Tf - Ti) * (mag - n*fStep_along_track)/fStep_along_track;
return rtn;
}
//________________________________________________________________________________________
EarthVector LowtranRadiativeProcessesCalculator::InLowtranFrame(const EarthVector& v1,const EarthVector& v2) const {
//
// Change of frame to match lowtran (v1 is the reference vector)
// the frame moves, then v2 expression in the new frame is returned
//
EarthVector rtn(v2);
EarthVector new1(v1);
EarthVector vrot;
if(v1.Theta() == 0) return rtn;
// translation -> frame is earth-centered
rtn(2) += EarthRadius();
new1(2) += EarthRadius();
// rotation defined to get earth-centered frame z-axis aligned with position new1
EarthVector Uz(0,0,1);
vrot = new1.Cross(Uz);
rtn.Rotate(new1.Theta(),vrot);
rtn(2) -= EarthRadius();
return rtn;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::CallLowtran100(const EarthVector& pos, const EarthVector& nextpos, LowtranTransferOutput& output) const {
//
// Call lowtran (100 wavelengths processed) to calculate transmission between pos and nextpos, then fill output
//
EarthVector pos1(1);
EarthVector pos2(1);
pos1 = InLowtranFrame(pos,pos);
pos2 = InLowtranFrame(pos,nextpos);
if(pos1.Zv() > pos2.Zv()) { // to get h1 < h2 (then good theta)
pos1 = InLowtranFrame(nextpos,nextpos);
pos2 = InLowtranFrame(nextpos,pos);
}
Float_t h1 = (Float_t)pos1.Zv()/km;
Float_t h2 = (Float_t)pos2.Zv()/km;
Float_t theta = (Float_t)(pos2 - pos1).Theta()/deg;
h1 = floor(h1*1000 + 0.5)/1000; // to avoid lowtran round-off bugs
h2 = floor(h2*1000 + 0.5)/1000;
theta = floor(theta*1000 + 0.5)/1000;
Int_t flag(1);
lwtrans_v2(&flag,&h1,&h2,&theta);
// copy lowtran tuple_ in output
output.R = lowtrantransfer.R * km;
output.Altitude = lowtrantransfer.Altitude * km;
output.Angle = lowtrantransfer.Angle;
for(Int_t i=0; i<100; i++) {
output.WaveLength[i] = lowtrantransfer.WaveLength[i] * micrometer;
output.Transmission[i] = lowtrantransfer.Transmission[i];
if(isnan(output.Transmission[i])) Msg(EsafMsg::Panic) <<"LOWTRAN returns a NaN" << MsgDispatch;
output.Rayl[i] = lowtrantransfer.Rayl[i];
output.Ozone[i] = lowtrantransfer.Ozone[i];
output.Aerosol[i] = lowtrantransfer.Aerosol[i];
}
}
//________________________________________________________________________________________
Double_t LowtranRadiativeProcessesCalculator::InterpolWl(string type,Double_t w,Double_t w1,Double_t w2,Double_t T1,Double_t T2) const {
//
// Interpolation of transmission coefficients as function of wavelength
//
Double_t rtn(0);
if(w1 == w2) return w1;
if( (w - w1) * (w - w2) > 0 )
Msg(EsafMsg::Warning) << "<InterpolWl> Wl searched in not within bounds : (wl,wl1,wl2) = ("<<w/nm<<", "<<w1/nm<<", "<<w2/nm<<")" << MsgDispatch;
if(type == "rayleigh") {
Double_t expon = log(T2/T1) / (1./(w2*w2*w2*w2) - 1./(w1*w1*w1*w1));
rtn = T1*exp(expon * (1./(w*w*w*w) - 1./(w1*w1*w1*w1)));
}
else if(type == "linear") {
rtn = T1 + (w - w1)/(w2 - w1) * (T2 - T1);
}
else Msg(EsafMsg::Panic) << "Unvalid type of interpolation (LowtranRadiativeProcessesCalculator)" << MsgDispatch;
return rtn;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::CleanBuffer(string DM) {
//
// Clean the indicated internal trans map (fAlong.. or fToDetect..)
//
if(DM == "alongtrack") {
if(fAlongTrackCoeff.size()) {
for(size_t i=0; i<fAlongTrackCoeff.size(); i++) {
SafeDelete(fAlongTrackCoeff[i]);
}
fAlongTrackCoeff.clear();
}
}
else if(DM == "todetector") {
if(fToDetectorCoeff.size()) {
for(size_t i=0; i<fToDetectorCoeff.size(); i++) {
SafeDelete(fToDetectorCoeff[i]);
}
fToDetectorCoeff.clear();
}
}
else Msg(EsafMsg::Warning) << "<CleanBuffer> Wrong argument -> POTENTIAL MEMORY LEAK" << MsgDispatch;
}
//________________________________________________________________________________________
void LowtranRadiativeProcessesCalculator::Reset() {
//
// get ready for next event
//
CleanBuffer("alongtrack");
fFlag = false;
fInitPos.SetXYZ(0,0,HUGE);
fStep_along_track = 0;
fNbwl = 0;
fNbSinglesFarFromTrack = 0.;
}