Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

LowtranRadiativeProcessesCalculator - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: LowtranRadiativeProcessesCalculator.cc,v 1.47 2005/05/11 13:42: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;

//____________________________________________________________________________________
    //
    // 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 = 0;
    fFlag = false;
    fInitPos.SetXYZ(0,0,HUGE);
    fNbwl = 0;
    fName = "lowtran";
}

//________________________________________________________________________________________
LowtranRadiativeProcessesCalculator::~LowtranRadiativeProcessesCalculator() {
    //
    // dtor
    //
    Reset();
}

//________________________________________________________________________________________
 Double_t LowtranRadiativeProcessesCalculator::Trans(const SinglePhoton& s, const EarthVector& finalpos) const {
    //
    // Total transmission between photon position and finalpos
    // 
    // 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 rtn, trans1, trans2, coeff1, coeff2;
    UInt_t Wli, Stepi;
    Float_t wl = s.Wl();
    Double_t OD = Atmosphere::Get()->GetClouds()->TotalOD(s.Pos(),s.Dir());
    Double_t travel, entrystep, endstep, ontrack;
    travel = (finalpos - s.Pos()).Mag();
    const ListPhotonsInAtmosphere& track = *fTrack;
    UInt_t tr_size = track.GetNbTrackSteps();
    
    
    // if fToDetectorCoeff card doesn't exist
    if(fToDetectorCoeff.size() == 0) {
        LowtranTransferOutput output;
	CallLowtran100(s.Pos(),finalpos,output);
        Double_t rtn(0);
        // 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 = InterpolWl("linear",wl,output.WaveLength[Wli],output.WaveLength[Wli+1],output.Transmission[Wli],output.Transmission[Wli+1]);
        return rtn * exp(-OD);
    }

    // if transmission card to detector exists :
    // 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)) {
        if(s.Status() != Reflected) Msg(EsafMsg::Debug) <<"scattered photons position further than track end" <<MsgDispatch;
	if(fImpact.Z() == HUGE) Msg(EsafMsg::Warning) <<"SinglePhoton reflected but No track impact" <<MsgDispatch;
	coeff1 = fToDetectorCoeff[fToDetectorCoeff.size()-1][Wli];
        if(Wli == 99) coeff2 = coeff1;
        else coeff2 = fToDetectorCoeff[fToDetectorCoeff.size()-1][Wli+1];
        rtn = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2);   // reference transmission
        endstep = (finalpos - fImpact).Mag();
	return pow(rtn,travel/endstep) * 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();
    
    // interpolation in lambda for the two step bounds
    coeff1 = fToDetectorCoeff[Stepi][Wli];
    if(Wli == 99) coeff2 = coeff1;
    else coeff2 = fToDetectorCoeff[Stepi][Wli+1];
    trans1 = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2);
    //
    coeff1 = fToDetectorCoeff[Stepi+1][Wli];
    if(Wli == 99) coeff2 = coeff1;
    else coeff2 = fToDetectorCoeff[Stepi+1][Wli+1];
    trans2 = InterpolWl("linear",wl,fWlArray[Wli],fWlArray[Wli+1],coeff1,coeff2);

    rtn = (trans2 - trans1)/(endstep - entrystep) * (mag_projec - entrystep) + trans1;  // reference transmission
    ontrack = (finalpos - projec).Mag();
    
    return pow(rtn,travel/ontrack) * 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 = pConfig->GetNum("BunchRadiativeTransfer.fStep")*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);
    if(Nb == 0) Msg(EsafMsg::Warning) << "BunchPreProcess() : Nb should not be equal to zero" << MsgDispatch;
    if(Nb > 1000) Nb = 1000;
    fStep = 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);
	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::PreProcessSingles(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()
    //
#ifdef DEBUG
    Msg(EsafMsg::Debug) <<"Preprocessing Singles propagation" << MsgDispatch;
#endif
    
    // 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);
    Double_t* Transmis = 0;
    
    // loop makes calculations then fills the map
    // an array for each step
    for(UInt_t k=0; k < tr_size; k++) {
	Transmis = new Double_t[100];
        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][i] = output.Transmission[i];
	}
    }
    // last step with init point at the track impact
    Transmis = new Double_t[100];
    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][i] = output.Transmission[i];
    }
    else {
        for(UInt_t i=0; i<100; i++) fToDetectorCoeff[tr_size][i] = 0;        
	Msg(EsafMsg::Info) << "PreProcessSingles: NO IMPACT on ground -> Trans=0" << MsgDispatch;
    }
}

//________________________________________________________________________________________
 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 == Double_t(fAlongTrackCoeff.size())) return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
    if(mag > fAlongTrackCoeff.size()*fStep) {
        if(fabs(mag - fAlongTrackCoeff.size()*fStep) > 1*mm) {
#ifdef DEBUG
	    //Msg(EsafMsg::Debug) << "<Interpolate> mag is too large -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep) << MsgDispatch;
            if(fabs(mag - fAlongTrackCoeff.size()*fStep) > 5e4*mm)
	        Msg(EsafMsg::Warning) << "<Interpolate> mag is TOO large (>50m) -> " <<fabs(mag - fAlongTrackCoeff.size()*fStep) << MsgDispatch;
#endif
	    mag = fAlongTrackCoeff.size()*fStep;
	}
        else return (*fAlongTrackCoeff[fAlongTrackCoeff.size() - 1])[i][j];
    }
    Double_t iter = floor(mag/fStep);
    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)/fStep * log(Tf/Ti));
    // linear interpolation
    //rtn = Ti + (Tf - Ti) * (mag - n*fStep)/fStep;

    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;
    h1 = floor(h1*1000 + 0.5)/1000;  // to avoid lowtran round-off bugs
    h2 = floor(h2*1000 + 0.5)/1000;
    Float_t theta = (Float_t)(pos2 - pos1).Theta()/deg;
    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];
        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(type == "rayleigh") {
        Double_t expon = (w*w*w*w - w1*w1*w1*w1) / 
                         (w2*w2*w2*w2 - w1*w1*w1*w1)* w2*w2*w2*w2;
        rtn = T1*exp(expon * log(T2/T1));
    }
    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++) {
		if(fToDetectorCoeff[i]) delete[] fToDetectorCoeff[i];
		fToDetectorCoeff[i] = 0;
            }
            fToDetectorCoeff.clear();
	}
    }
    else Msg(EsafMsg::Warning) << "<CleanBuffer> Wrong argument -> POTENTIAL MEMORY LEAK" << MsgDispatch;
}

//________________________________________________________________________________________
 void LowtranRadiativeProcessesCalculator::Reset() {
    //
    // get ready for next event
    //
    CleanBuffer("alongtrack");
    CleanBuffer("todetector");
    fFlag = false;
    fInitPos.SetXYZ(0,0,HUGE);
    fStep = 0;
    fNbwl = 0;
}

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31