source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/radiativetransfer/src/RadiativeProcessesCalculator.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: RadiativeProcessesCalculator.cc 2527 2006-03-02 13:40:43Z moreggia $
// Sylvain Moreggia created Jan,  9 2004

#include "RadiativeProcessesCalculator.hh"
#include "EarthVector.hh"
#include "Atmosphere.hh"
#include <math.h>
#include "BunchOfPhotons.hh"
#include "SinglePhoton.hh"
#include "EsafMsgSource.hh"
#include "EsafRandom.hh"
#include <TF1.h>

ClassImp(RadiativeProcessesCalculator)

using namespace TMath;


//_____________________________________________________________________________
RadiativeProcessesCalculator::RadiativeProcessesCalculator() : EsafMsgSource() {
    //
    // ctor
    //
}

//_____________________________________________________________________________
RadiativeProcessesCalculator::~RadiativeProcessesCalculator() {
    //
    // dtor
    //
}

//_____________________________________________________________________________
void RadiativeProcessesCalculator::RandomDir(string type,EarthVector& dir, Double_t g) const {
    //
    // 'dir' is the incoming direction -> this method set dir to the scattering outgoing direction
    // get a random direction, sampling relevant scattering phase functions
    //
    
    Double_t theta(0);
    Double_t u1, u2;
    TRandom* rndm = EsafRandom::Get();
    
    // get random theta, assuming isotropic distribution in phi
    // theta is the angle between incoming and outgoing directions
    if(type == "rayleigh") {
        // von neumann sampling -> much faster than using TF1
        do {
            u1 = -1. + rndm->Rndm()*2.;
            u2 = rndm->Rndm() * 3./4.;
        } while(u2 >= (3./8. * (1. + u1*u1)));
        theta = acos(u1);
    }
    else if(type == "HG_mie") {
        Msg(EsafMsg::Panic) << "<RandomTheta> Mie scattering not implemented yet"<<MsgDispatch;
        if(g == -10) Msg(EsafMsg::Panic) << "<RandomTheta> Asymmetry parameter should be defined" <<MsgDispatch; 
    }
    else Msg(EsafMsg::Panic) << "<RandomTheta> Wrong type of distribution : "<<type<<MsgDispatch;
    
    // 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);
}

//_____________________________________________________________________________
Double_t RadiativeProcessesCalculator::RayleighPhaseFunction(const EarthVector& bunchdir,const EarthVector& dir) const {
    //
    // rayleigh scattering phase function - Normalized when integrated over full solid angle
    //
    Double_t rtn;
    Double_t theta = dir.Angle(bunchdir);
    rtn = 3 * (1. + cos(theta)*cos(theta)) / (16.*Pi());

    return rtn;
}

//_____________________________________________________________________________
Double_t RadiativeProcessesCalculator::Mie_HG_PhaseFunction(const EarthVector& bunchdir,const EarthVector& dir,Double_t g) const {
    //
    // Mie scattering phase function - Henyey-Greenstein parametrization
    // Normalized when integrated over full solid angle
    //
    Double_t rtn;
    Double_t theta = dir.Angle(bunchdir);
    rtn = 1./(4*Pi()) * (1 - g*g) / pow(1 - 2*g*cos(theta) + g*g,1.5);

    return rtn;
}
//_____________________________________________________________________________
Double_t RadiativeProcessesCalculator::GetMaxPhaseFunction(string type) const {
    //
    // returns maximum value of specified phase function (phase function is normalized when integrated over full solid angle)
    //
    Double_t rtn;
    if(type == "rayleigh") rtn = 3./(8.*Pi());
    else Msg(EsafMsg::Panic) << "<GetMaxPhaseFunction> Wrong type of phase function : "<<type<<MsgDispatch;
    
    return rtn;
}