source: JEM-EUSO/esaf_cc_at_lal/packages/common/atmosphere/src/LowtranAerosol.cc @ 114

// $Id: LowtranAerosol.cc 2922 2011-06-12 14:21:23Z mabl $
// 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
        
//TOFIX : there is a pb in the tabulated phase function coeff. for urban5 model
// so far, they have been replaced by rural5 coeff.
        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;

/*
        // values for 70% rel. humidity (urban5 coeff)
        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) {
        scat = 1.;
        return 1.;
    }
    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. this change of variable comes from the old version of the code, which was bugged
    Double_t K = Kwl(wl);  // = extinction coeff. !!
    // 1(bis). tell the kind of interaction : scatt. or abs.
    Double_t kwlabs = Kwl(wl,kTRUE);
    Double_t abs = kwlabs / K;
    Double_t scat = 1 - abs;
    if(scat > rndm->Rndm()) 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();
}