Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

LowtranRadiativeProcessesCalculator - source file

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

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005