!----------------------------------------------------------------------------------------!
! Author: Benjamin Marchant (2008)
! 
! Outputs: phase function (PF) values as fct of scattering angle for  
! randomly oriented pristine hexagonal crystals. 
! Methods: Ray-Tracing + Maggi-Rubinowicz (Diffraction)

program gom_phm

implicit none

integer, parameter :: pr=selected_real_kind(15,3)

integer, parameter :: nb_teta = 200
integer, parameter :: nb_teta_mr = 200!60
integer, parameter :: nb_phi_mr = 360

real(pr) :: wl = 0.66d0  ! Wavelenght (micron)

! Real and imaginary part of the refractive index
real(pr) :: inr = 1.308d0 , ini = 1.67d-8  

real(pr) :: wk  ! Wave Number (micron).
real(pr) :: wk2

real(pr) :: reff = 300.0d0  ! Effective Radius (micron) 

real(pr) :: qf = 2.5d0 ! Shape Factor 

integer :: nod = 1000 
integer :: nor = 1000000

integer :: ray_coeff = 2 

integer :: nray_min,nray  ! ray-tracing number   

integer :: max_rm

real(pr) :: eps = 1.0d-6 ! Stop condition multiple scattering 

integer :: i,j,k,l,m,n,p

real(pr) :: rcsi,rcs_tot

real(pr) :: abcs

real(pr) :: pi = acos( -1.0d0 )
real(pr) :: pi2 

real(pr) :: cdr 

real(pr) :: const1 

real(pr) :: ray,length,length2

real(pr), dimension(12,3) :: cs   
real(pr), dimension(8,3)  :: cn  

real(pr), dimension(8,2) :: acs,acn

real(pr) :: projected_area,ps1,ps2,ps3,ps4

integer :: face

real(pr) :: cteta,steta,cphi,sphi  ! cosinus directeur.
real(pr) :: t11,t12,t13,t21,t22,t23,t31,t32,t33

real(pr) :: gg,bg,sm  
real(pr) :: al,al1,al2

integer :: bin

real(pr) :: norme

real(pr), dimension(6) :: TAB_temp2,TAB_temp3

real(pr) :: temp1,temp2,temp3,temp4,temp5,temp6,temp7,temp8

real(pr) :: cond,cond2

real(pr) :: xrt,yrt,zrt
real(pr) :: distance,distmin,distmax,dist,dist_tot

real(pr) :: xmax,ymax

real(pr) :: rtcond

real(pr) :: m1,m2,m3,m4
real(pr) :: s21,s23,s24,s31,s34,s41
real(pr) :: d21,d23,d24,d31,d34,d41

real(pr) :: tp,ts,rp,rs
real(pr) :: rpr,rpi,rsr,rsi

real(pr) :: diff_solid_angle

real(pr) :: xsr,ysr

real(pr) :: teta,phi

real(pr), dimension(nb_teta) :: tab_teta
real(pr), dimension(nb_phi_mr) :: tab_phi

real(pr), dimension(nb_teta) :: F11_MR,F12_MR,F13_MR,F14_MR,F21_MR,F22_MR,F23_MR,F24_MR&
&                              ,F31_MR,F32_MR,F33_MR,F34_MR,F41_MR,F42_MR,F43_MR,F44_MR
real(pr), dimension(nb_teta) :: P11_MR,P12_MR,P13_MR,P14_MR,P21_MR,P22_MR,P23_MR,P24_MR&
&                              ,P31_MR,P32_MR,P33_MR,P34_MR,P41_MR,P42_MR,P43_MR,P44_MR

real(pr), dimension(nb_teta) :: F11_RT,F12_RT,F13_RT,F14_RT,F21_RT,F22_RT,F23_RT,F24_RT&
&                              ,F31_RT,F32_RT,F33_RT,F34_RT,F41_RT,F42_RT,F43_RT,F44_RT
real(pr), dimension(nb_teta) :: P11_RT,P12_RT,P13_RT,P14_RT,P21_RT,P22_RT,P23_RT,P24_RT&
&                              ,P31_RT,P32_RT,P33_RT,P34_RT,P41_RT,P42_RT,P43_RT,P44_RT

real(pr), dimension(nb_teta) :: P11,P12,P13,P14,P21,P22,P23,P24&
&                              ,P31,P32,P33,P34,P41,P42,P43,P44

real(pr), dimension(2,2) :: ea1,ea2

real(pr) :: ecs,scs

real(pr) :: albedo

real(pr) :: f_delta

real(pr) :: lc  
real(pr) :: prub,piub
real(pr) :: part1,part2,part3,part4,part5
real(pr) :: fx,fy
real(pr) :: csx1,csx2,csy1,csy2
real(pr) :: coeffdir,oo

real(pr) :: dx,dy

real(pr), dimension(3) :: n1,np

real(pr), dimension(2,2,2) :: tab,tab1,tab2,tab3,tab4,tab5,tab6

real(pr), dimension(2,2) :: gamma1,gammap

real(pr) :: ps

real(pr), dimension(2,2,2) :: fresnel_r
real(pr), dimension(2,2) :: fresnel_t

real(pr) :: absorp

real(pr) :: side,base

real(pr) :: rn,rn1,rn2

real(pr) :: intensity,is

real(pr) :: tetai,tetaf

real(pr), dimension(3) :: ei

real(pr), dimension(3) :: e_i1,e_t1,e_s1
real(pr), dimension(3) :: e_ip,e_rp,e_sp

real(pr), dimension(3) :: beta_i,alpha_i
real(pr), dimension(3) :: beta_s,alpha_s

real(pr), dimension(3) :: beta_i1,alpha_i1
real(pr), dimension(3) :: beta_r1,alpha_r1
real(pr), dimension(3) :: beta_t1,alpha_t1

real(pr), dimension(3) :: beta,alpha

real(pr), dimension(3) :: beta_ip,alpha_ip
real(pr), dimension(3) :: beta_rp,alpha_rp
real(pr), dimension(3) :: beta_tp,alpha_tp

real(pr), dimension(2,2) :: gamma_1,gamma_p
real(pr), dimension(2,2) :: gammas,gamma_i

real(pr), dimension(2,2) :: r_1,t_1,t_p
real(pr), dimension(2,2,2) :: r_p

real(pr), dimension(2,4) :: ma

real(pr), dimension(2,2) :: taba,tabb,tabc,tabd

real(pr) :: min

real(pr) :: norme_rt

real(pr) :: num_dirac_function

cdr = pi / 180.0d0

pi2 = 2.0d0 * pi

const1 = dcos( pi / 6.0d0 )

wk  = pi2 / wl  ! Wavenumber.
wk2 = 2.0d0 * wk

!----------------------------------------------------------------------------------------!

temp1 = pi2 / ( 9.0d0 * const1 * qf )

ray = reff * temp1 ** 0.33d0  ! Rayon du cercle circonscrit à la base de l'hexagone.

length = qf * 2.0d0 * ray  ! Longueur de l'hexagone.

!ray = 80.0d0 
!length = 400.0d0

length2 = length / 2.0d0

side = ray * length
base = 3.0d0 * ray * ray * const1

if(base .lt. side)then
rcsi = base / 10.0
else
rcsi = side / 10.0
end if

rcsi = ( wl * wl ) / ( 4.0d0 * pi )  !##### Cross section of the incident ray #####!

!----------------------------------------------------------------------------------------!
! Coordonnées des sommets du modèle PHM.

cs(1,1)  =   ray
cs(1,2)  =   0.0d0
cs(1,3)  = - length2
cs(2,1)  =   ray * 0.5d0
cs(2,2)  =   ray * const1
cs(2,3)  = - length2
cs(3,1)  = - ray * 0.5d0
cs(3,2)  =   ray * const1
cs(3,3)  = - length2
cs(4,1)  = - ray
cs(4,2)  =   0.0d0
cs(4,3)  = - length2
cs(5,1)  = - ray * 0.5d0
cs(5,2)  = - ray * const1
cs(5,3)  = - length2
cs(6,1)  =   ray * 0.5d0
cs(6,2)  = - ray * const1
cs(6,3)  = - length2

cs(7,1)  =   ray
cs(7,2)  =   0.0d0
cs(7,3)  =   length2
cs(8,1)  =   ray * 0.5d0
cs(8,2)  =   ray * const1
cs(8,3)  =   length2
cs(9,1)  = - ray * 0.5d0
cs(9,2)  =   ray * const1
cs(9,3)  =   length2
cs(10,1) = - ray
cs(10,2) =   0.0d0
cs(10,3) =   length2
cs(11,1) = - ray * 0.5d0
cs(11,2) = - ray * const1
cs(11,3) =   length2
cs(12,1) =   ray * 0.5d0
cs(12,2) = - ray * const1
cs(12,3) =   length2

distmax = sqrt( 4.0d0 * ray * ray + length * length ) 

!----------------------------------------------------------------------------------------!
! Coordonnées des normales (Orientées vers l'interieure).

cn(1,1) = - const1
cn(1,2) = - 0.5d0
cn(1,3) =   0.0d0
cn(2,1) =   0.0d0
cn(2,2) = - 1.0d0
cn(2,3) =   0.0d0
cn(3,1) =   const1
cn(3,2) = - 0.5d0
cn(3,3) =   0.0d0
cn(4,1) =   const1 
cn(4,2) =   0.5d0
cn(4,3) =   0.0d0
cn(5,1) =   0.0d0
cn(5,2) =   1.0d0
cn(5,3) =   0.0d0
cn(6,1) = - const1
cn(6,2) =   0.5d0
cn(6,3) =   0.0d0
cn(7,1) =   0.0d0
cn(7,2) =   0.0d0
cn(7,3) =   1.0d0
cn(8,1) =   0.0d0
cn(8,2) =   0.0d0
cn(8,3) = - 1.0d0

!----------------------------------------------------------------------------------------!
! Choix des angles de diffusion.

do i = 1 , 11
tab_teta(i) = 0.1d0 * float(i-1)
end do

do i = 12, 16
tab_teta(i) = 1.0d0 + 0.2d0 * float(i-11)
end do

do i = 17, 20
tab_teta(i) = 2.0d0 + 0.25d0 * float(i-16)
end do

tab_teta(21) = 3.3
tab_teta(22) = 3.6
tab_teta(23) = 4.0
tab_teta(24) = 4.5
tab_teta(25) = 5.0

do i = 26, nb_teta
tab_teta(i) = 5.0 + float(i-25)
end do

!----------------------------------------------------------------------------------------!
! Initialisation de la matrice de phase 
! pour la diffraction et le lancer de rayon.

do i = 1 , nb_teta

F11_RT(i) = 0.0d0
F12_RT(i) = 0.0d0
F13_RT(i) = 0.0d0
F14_RT(i) = 0.0d0
F21_RT(i) = 0.0d0
F22_RT(i) = 0.0d0
F23_RT(i) = 0.0d0
F24_RT(i) = 0.0d0
F31_RT(i) = 0.0d0
F32_RT(i) = 0.0d0
F33_RT(i) = 0.0d0
F34_RT(i) = 0.0d0
F41_RT(i) = 0.0d0
F42_RT(i) = 0.0d0
F43_RT(i) = 0.0d0
F44_RT(i) = 0.0d0

P11_RT(i) = 0.0d0
P12_RT(i) = 0.0d0
P13_RT(i) = 0.0d0
P14_RT(i) = 0.0d0
P21_RT(i) = 0.0d0
P22_RT(i) = 0.0d0
P23_RT(i) = 0.0d0
P24_RT(i) = 0.0d0
P31_RT(i) = 0.0d0
P32_RT(i) = 0.0d0
P33_RT(i) = 0.0d0
P34_RT(i) = 0.0d0
P41_RT(i) = 0.0d0
P42_RT(i) = 0.0d0
P43_RT(i) = 0.0d0
P44_RT(i) = 0.0d0

F11_MR(i) = 0.0d0
F12_MR(i) = 0.0d0
F13_MR(i) = 0.0d0
F14_MR(i) = 0.0d0
F21_MR(i) = 0.0d0
F22_MR(i) = 0.0d0
F23_MR(i) = 0.0d0
F24_MR(i) = 0.0d0
F31_MR(i) = 0.0d0
F32_MR(i) = 0.0d0
F33_MR(i) = 0.0d0
F34_MR(i) = 0.0d0
F41_MR(i) = 0.0d0
F42_MR(i) = 0.0d0
F43_MR(i) = 0.0d0
F44_MR(i) = 0.0d0

P11_MR(i) = 0.0d0
P12_MR(i) = 0.0d0
P13_MR(i) = 0.0d0
P14_MR(i) = 0.0d0
P21_MR(i) = 0.0d0
P22_MR(i) = 0.0d0
P23_MR(i) = 0.0d0
P24_MR(i) = 0.0d0
P31_MR(i) = 0.0d0
P32_MR(i) = 0.0d0
P33_MR(i) = 0.0d0
P34_MR(i) = 0.0d0
P41_MR(i) = 0.0d0
P42_MR(i) = 0.0d0
P43_MR(i) = 0.0d0
P44_MR(i) = 0.0d0

end do

!----------------------------------------------------------------------------------------!
! Calculation of the normalized phase functions for the part of diffraction.

gg = 1  ! graine du générateur aléatoire.
bg = 2147483647  ! entier le plus grand sur 4 octets.

do n = 1 , nod

!----------------------------------------------------------------------------------------!
! Orientation aleatoire.

sm = mod( gg * 16807 , bg )  ! standard Minimal d'après "Park & Miller".
gg = sm

al    = real( sm ) / real( bg )  ! nombre aléatoire compris entre 0 et 1.             
cteta = al 
steta = sqrt( 1.0d0 - cteta * cteta ) 

sm = mod( gg * 16807 , bg )          
gg = sm

al   = real( sm ) / real( bg )                         
phi  = 30.0d0 * al                       
cphi = cos( phi * cdr )
sphi = sin( phi * cdr )

!----------------------------------------------------------------------------------------!
! Matrice dce changement de base.

t11 =   cphi * cteta
t12 = - sphi
t13 =   steta * cphi
t21 =   sphi * cteta
t22 =   cphi
t23 =   steta * sphi
t31 = - steta
t32 =   0.0d0
t33 =   cteta

ei(1) = t13
ei(2) = t23
ei(3) = t33

!----------------------------------------------------------------------------------------!
! Diffraction: transformation de Maggi-Rubinowicz. 

!----------------------------------------------------------------------------------------!
! Aperture.

acs(1,1) = t11 * cs(9,1) + t21 * cs(9,2) + t31 * cs(9,3)
acs(1,2) = t12 * cs(9,1) + t22 * cs(9,2) + t32 * cs(9,3)
acs(2,1) = t11 * cs(3,1) + t21 * cs(3,2) + t31 * cs(3,3)
acs(2,2) = t12 * cs(3,1) + t22 * cs(3,2) + t32 * cs(3,3)
acs(3,1) = t11 * cs(2,1) + t21 * cs(2,2) + t31 * cs(2,3)
acs(3,2) = t12 * cs(2,1) + t22 * cs(2,2) + t32 * cs(2,3)
acs(4,1) = t11 * cs(1,1) + t21 * cs(1,2) + t31 * cs(1,3)
acs(4,2) = t12 * cs(1,1) + t22 * cs(1,2) + t32 * cs(1,3)
acs(5,1) = t11 * cs(6,1) + t21 * cs(6,2) + t31 * cs(6,3)
acs(5,2) = t12 * cs(6,1) + t22 * cs(6,2) + t32 * cs(6,3)
acs(6,1) = t11 * cs(12,1) + t21 * cs(12,2) + t31 * cs(12,3)
acs(6,2) = t12 * cs(12,1) + t22 * cs(12,2) + t32 * cs(12,3)
acs(7,1) = t11 * cs(11,1) + t21 * cs(11,2) + t31 * cs(11,3)
acs(7,2) = t12 * cs(11,1) + t22 * cs(11,2) + t32 * cs(11,3)
acs(8,1) = t11 * cs(10,1) + t21 * cs(10,2) + t31 * cs(10,3)
acs(8,2) = t12 * cs(10,1) + t22 * cs(10,2) + t32 * cs(10,3)

!----------------------------------------------------------------------------------------!
! Calcul des normales.

do i = 1 , 8  

if( i == 8 ) then  
dy = acs(1,2) - acs(i,2)
dx = acs(1,1) - acs(i,1)
else
dy = acs(i+1,2) - acs(i,2)
dx = acs(i+1,1) - acs(i,1)
end if
temp1 = dx / dy
if( abs( dy ) > 0.001d0 ) then  
acn(i,1) = 1.0d0
acn(i,2) = - temp1
norme = sqrt(1 + temp1 * temp1 )
acn(i,1) = acn(i,1) / norme
acn(i,2) = acn(i,2) / norme
else
acn(i,1) = 0.0d0
acn(i,2) = 1.0d0
end if
temp1 = acn(i,1) * acs(i,1) + acn(i,2) * acs(i,2)
if( temp1 > 0d0 ) then  ! Orientation.
acn(i,1) = - acn(i,1)
acn(i,2) = - acn(i,2)
end if

end do

!----------------------------------------------------------------------------------------!

part1 = 1.0d0 / ( 4.0d0 * pi )

do j = 2 , nb_teta_mr 

teta = tab_teta(j)

do k = 1 , nb_phi_mr

phi = tab_phi(k)

xsr = dsin( teta * cdr ) * dcos( phi * cdr )
ysr = dsin( teta * cdr ) * dsin( phi * cdr )

temp1 = tan( teta * cdr )
temp2 = cos( phi * cdr )
temp3 = sin( phi * cdr )

fx = ( temp1 * temp2 ) / wl
fy = ( temp1 * temp3 ) / wl

prub = 0.0d0
piub = 0.0d0

!----------------------------------------------------------------------------------------!
! Debut de la sommation sur les cotés du polynome.

do i = 1 , 8  

if( i == 8 ) then

csx1 = acs(i,1)
csy1 = acs(i,2)

csx2 = acs(1,1)
csy2 = acs(1,2)

else

csx1 = acs(i,1)
csy1 = acs(i,2)

csx2 = acs(i+1,1)
csy2 = acs(i+1,2)

end if

temp1 = acn(i,1) * dcos( phi * cdr ) + acn(i,2) * dsin( phi * cdr )
temp2 = dtan( ( teta * cdr ) / 2.0d0 )
part2 = temp1 / temp2

dy = csy2 - csy1
dx = csx2 - csx1

if( abs( dx ) > 0.01d0 ) then  ! Cas ou le segment est quelconque

coeffdir = dy / dx  ! Coefficient directeur de la droite. 
oo = csy2 - coeffdir * csx2  ! Ordonnée a l'origine de la droite.

temp1 = abs( dx )
temp2 = ( wk / 2.0d0 ) * ( xsr + coeffdir * ysr ) * temp1
part3 = temp1 * fsc( temp2 )

part4 = sqrt( 1.0d0 + coeffdir * coeffdir )

temp1 = ( wk / 2.0d0 ) * ( xsr + coeffdir * ysr ) * ( csx1 + csx2 )
temp2 = wk * ysr * oo
part5 = temp1 + temp2

prub = prub + part2 * part3 * part4 * cos( part5 )
piub = piub + part2 * part3 * part4 * sin( part5 )

else  ! Cas ou le segment est parallele à l'axe y.

temp1 = abs( dy )
temp2 = ( wk / 2.0d0 ) * ysr * temp1
part3 = temp1 * fsc( temp2 )

temp1 = ( wk / 2.0d0 ) * ysr * ( csy1 + csy2 )
temp2 = wk * xsr * csx1
part4 = temp1 + temp2

prub = prub + part2 * part3 * cos( part4 )
piub = piub + part2 * part3 * sin( part4 )

end if

end do  

! Fin de la sommation sur les cotés du polynome.
!----------------------------------------------------------------------------------------!
! Calcul de l'intensité.

prub = part1 * prub
piub = part1 * piub
temp1 = ( prub * prub + piub * piub ) / ( 16.0d0 * pi * pi )

F11_MR(j) = F11_MR(j) +  wk * wk * temp1

end do

end do

!----------------------------------------------------------------------------------------!
! Calcul de l'intensité en 0.

ps1 = ( ei(1) * cn(3,1) + ei(2) * cn(3,2) ) * side
ps2 = ( ei(1) * cn(4,1) + ei(2) * cn(4,2) ) * side
ps3 = ( ei(1) * cn(5,1) + ei(2) * cn(5,2) ) * side
ps4 =   ei(3) * base

projected_area = ps1 + ps2 + ps3 + ps4

temp1 = ( wk * wk * projected_area ) / pi2
F11_MR(1) = F11_MR(1) + temp1 * temp1 

!----------------------------------------------------------------------------------------!

end do  ! Fin de la boucle sur les orientations.

!----------------------------------------------------------------------------------------!
! Calcul de la matrice de phase correspondant à la diffraction

!----------------------------------------------------------------------------------------!
! Calcul de la norme

norme = 0.0d0

do i = 1 , nb_teta_mr

if(i==1)then

tetai = ( tab_teta(i) + 0.0d0 ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

else

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

end if

norme = norme + 0.5d0 * F11_MR(i) * ( cos( tetai * cdr ) - cos( tetaf * cdr ) )

end do

!----------------------------------------------------------------------------------------!
! Normalisation

do i = 1 , nb_teta

P11_MR(i) = F11_MR(i) / norme
P22_MR(i) = F11_MR(i) / norme
P33_MR(i) = F11_MR(i) / norme
P44_MR(i) = F11_MR(i) / norme

end do

!----------------------------------------------------------------------------------------!
! Ecriture des resultats

open(1,file = 'hexagomal_pristine_crystal_MR.txt')

do i = 1 , nb_teta_mr

write(1,*) tab_teta(i),P11_MR(i)

end do

close(1)

!----------------------------------------------------------------------------------------!
! Calculation of the normalized phase functions for the part of reflection and refraction.

rcs_tot = 0.0d0
abcs = 0.0d0

gg = 1  ! graine du générateur aléatoire.
bg = 2147483647  ! entier le plus grand sur 4 octets.

do n = 1 , nor

!----------------------------------------------------------------------------------------!
! Orientation aleatoire.

sm = mod( gg * 16807 , bg )  ! standard Minimal d'après "Park & Miller".
gg = sm

al    = real( sm ) / real( bg )  ! nombre aléatoire compris entre 0 et 1.             
cteta = al 
steta = sqrt( 1.0d0 - cteta * cteta ) 

sm = mod( gg * 16807 , bg )          
gg = sm

al   = real( sm ) / real( bg )                         
phi  = 30.0d0 * al                       
cphi = cos( phi * cdr )
sphi = sin( phi * cdr )

!----------------------------------------------------------------------------------------!
! changement de base.

t11 =   cphi * cteta
t12 = - sphi
t13 =   steta * cphi
t21 =   sphi * cteta
t22 =   cphi
t23 =   steta * sphi
t31 = - steta
t32 =   0.0
t33 =   cteta

!----------------------------------------------------------------------------------------!
! Ray_tracing.

ei(1) = t13
ei(2) = t23
ei(3) = t33

ps1 = ( ei(1) * cn(3,1) + ei(2) * cn(3,2) ) * side
ps2 = ( ei(1) * cn(4,1) + ei(2) * cn(4,2) ) * side
ps3 = ( ei(1) * cn(5,1) + ei(2) * cn(5,2) ) * side
ps4 =   ei(3) * base

projected_area = ps1 + ps2 + ps3 + ps4

nray_min = int( projected_area / rcsi ) 

nray = 1 ! nray_min * ray_coeff 

!write(6,*) 'nray : ' , nray

do j = 1 , nray

rcs_tot = rcs_tot + rcsi

ea1(1,1) = 1.0d0
ea1(1,2) = 0.0d0
ea1(2,1) = 0.0d0
ea1(2,2) = 0.0d0

ea2(1,1) = 0.0d0
ea2(1,2) = 1.0d0
ea2(2,1) = 0.0d0
ea2(2,2) = 0.0d0

e_i1(1) = ei(1)
e_i1(2) = ei(2)
e_i1(3) = ei(3)

rtcond = 1.0d0

!----------------------------------------------------------------------------------------!
! Calcul du point d'impact.

temp1 = ps1 + ps2
temp2 = temp1 + ps3

sm = mod( gg * 16807 , bg ) 
gg = sm

rn = real( sm ) / real( bg )  
rn = rn * projected_area

!----------------------------------------------------------------------------------------!
if( rn < ps1 ) then  

n1(1) = - cn(3,1)
n1(2) = - cn(3,2)
n1(3) = - cn(3,3)

sm = mod( gg * 16807 , bg ) 
gg = sm

rn1 = real( sm ) / real( bg )  
xrt = ( 2.0d0 - rn1 ) * ( - ray / 2.0d0 ) 

sm = mod( gg * 16807 , bg ) 
gg = sm

rn2 = real( sm ) / real( bg )  
zrt = ( 1.0d0 - 2.0d0 * rn2 ) * length2

yrt = ( -1.0d0 / n1(2) ) * ( xrt * n1(1) + zrt * n1(3) - const1 * ray )

end if

!----------------------------------------------------------------------------------------!
if( ps1 < rn .and. rn < temp1 ) then

n1(1) = - cn(4,1)
n1(2) = - cn(4,2)
n1(3) = - cn(4,3)

sm = mod( gg * 16807 , bg ) 
gg = sm

rn1 = real( sm ) / real( bg )  
xrt = ( 2.0d0 - rn1 ) * ( - ray / 2.0d0 ) 

sm = mod( gg * 16807 , bg ) 
gg = sm

rn2 = real( sm ) / real( bg )  
zrt = ( 1.0d0 - 2.0d0 * rn2 ) * length2

yrt = ( -1.0d0 / n1(2) ) * ( xrt * n1(1) + zrt * n1(3) - const1 * ray )

end if

!----------------------------------------------------------------------------------------!
if( temp1 < rn .and. rn < temp2 ) then

n1(1) = - cn(5,1)
n1(2) = - cn(5,2)
n1(3) = - cn(5,3)

sm = mod( gg * 16807 , bg ) 
gg = sm

rn1 = real( sm ) / real( bg )  
xrt = ( 1.0d0 - 2.0d0 * rn1 ) * ( ray / 2.0d0 )

sm = mod( gg * 16807 , bg ) 
gg = sm

rn2 = real( sm ) / real( bg )  
zrt = ( 1.0d0 - 2.0d0 * rn2 ) * length2

yrt = ( -1.0d0 / n1(2) ) * ( xrt * n1(1) + zrt * n1(3) - const1 * ray )

end if

!----------------------------------------------------------------------------------------!
if( temp2 < rn .and. rn < projected_area ) then

n1(1) = - cn(7,1)
n1(2) = - cn(7,2)
n1(3) = - cn(7,3)

sm = mod( gg * 16807 , bg ) 
gg = sm

rn1 = real( sm ) / real( bg )  
yrt = ( 1.0d0 - 2.0d0 * rn1) * const1 * ray

if( yrt > 0.0 ) then
xmax = ( -1.0d0/ cn(1,1) ) * ( yrt * cn(1,2) + const1 * ray) 
xmax = abs( xmax )
else
xmax = ( -1.0d0 / cn(6,1) ) * ( yrt * cn(6,2) + const1 * ray) 
xmax = abs( xmax )
end if

sm = mod( gg * 16807 , bg ) 
gg = sm

rn2 = real( sm ) / real( bg )  
xrt = ( 1.0d0 - 2.0d0 * rn2) * xmax

zrt = - length2

end if

!----------------------------------------------------------------------------------------!
! Calcul des cosinus directeurs du rayon refracte et reflechi ( air -> glace ).

ps = e_i1(1) * n1(1) + e_i1(2) * n1(2) + e_i1(3) * n1(3)
           
temp1 = dsqrt( inr * inr - 1.0d0 + ps * ps )

e_t1(1) = ( 1.0d0 / inr ) * ( e_i1(1) - ps * n1(1) - temp1 * n1(1) )
e_t1(2) = ( 1.0d0 / inr ) * ( e_i1(2) - ps * n1(2) - temp1 * n1(2) )
e_t1(3) = ( 1.0d0 / inr ) * ( e_i1(3) - ps * n1(3) - temp1 * n1(3) )

e_s1(1) = e_i1(1) - 2.0d0 * ps * n1(1)
e_s1(2) = e_i1(2) - 2.0d0 * ps * n1(2)
e_s1(3) = e_i1(3) - 2.0d0 * ps * n1(3)

!----------------------------------------------------------------------------------------!
! Calcul de la matrice gamma1.

temp1 = ei(1) * e_s1(1) + ei(2) * e_s1(2) + ei(3) * e_s1(3)
norme = dsqrt( 1.0 - temp1 * temp1 )

beta_i(1) = ( ei(2) * e_s1(3) - ei(3) * e_s1(2) ) / norme
beta_i(2) = ( ei(3) * e_s1(1) - ei(1) * e_s1(3) ) / norme
beta_i(3) = ( ei(1) * e_s1(2) - ei(2) * e_s1(1) ) / norme

alpha_i(1) = beta_i(2) * ei(3) - beta_i(3) * ei(2)
alpha_i(2) = beta_i(3) * ei(1) - beta_i(1) * ei(3)
alpha_i(3) = beta_i(1) * ei(2) - beta_i(2) * ei(1) 

temp1 = n1(1) * e_i1(1) + n1(2) * e_i1(2) + n1(3) * e_i1(3)
norme = dsqrt( 1.0 - temp1 * temp1 )

beta_i1(1) = ( n1(2) * e_i1(3) - n1(3) * e_i1(2) ) / norme
beta_i1(2) = ( n1(3) * e_i1(1) - n1(1) * e_i1(3) ) / norme
beta_i1(3) = ( n1(1) * e_i1(2) - n1(2) * e_i1(1) ) / norme

alpha_i1(1) = beta_i1(2) * e_i1(3) - beta_i1(3) * e_i1(2)
alpha_i1(2) = beta_i1(3) * e_i1(1) - beta_i1(1) * e_i1(3)
alpha_i1(3) = beta_i1(1) * e_i1(2) - beta_i1(2) * e_i1(1) 

gamma1(1,1) = alpha_i1(1) * alpha_i(1) + alpha_i1(2) * alpha_i(2) + alpha_i1(3) * alpha_i(3)
gamma1(1,2) = alpha_i1(1) * beta_i(1)  + alpha_i1(2) * beta_i(2)  + alpha_i1(3) * beta_i(3)
gamma1(2,1) = beta_i1(1) * alpha_i(1)  + beta_i1(2) * alpha_i(2)  + beta_i1(3) * alpha_i(3)
gamma1(2,2) = beta_i1(1) * beta_i(1)   + beta_i1(2) * beta_i(2)   + beta_i1(3) * beta_i(3)

!----------------------------------------------------------------------------------------!
! Calcul intermédiaire 1 pour la section efficace d'absorption.

temp1 = ea1(1,1)
temp2 = ea1(1,2)
temp3 = ea1(2,1)
temp4 = ea1(2,2)

ea1(1,1) = gamma1(1,1) * temp1 + gamma1(1,2) * temp2
ea1(1,2) = gamma1(2,1) * temp1 + gamma1(2,2) * temp2
ea1(2,1) = gamma1(1,1) * temp3 + gamma1(1,2) * temp4
ea1(2,2) = gamma1(2,1) * temp3 + gamma1(2,2) * temp4

temp1 = ea2(1,1)
temp2 = ea2(1,2)
temp3 = ea2(2,1)
temp4 = ea2(2,2)

ea2(1,1) = gamma1(1,1) * temp1 + gamma1(1,2) * temp2
ea2(1,2) = gamma1(2,1) * temp1 + gamma1(2,2) * temp2
ea2(2,1) = gamma1(1,1) * temp3 + gamma1(1,2) * temp4
ea2(2,2) = gamma1(2,1) * temp3 + gamma1(2,2) * temp4

!----------------------------------------------------------------------------------------!
! Calcul des coefficients de Fresnel.

temp1 = - ( e_i1(1) * n1(1) + e_i1(2) * n1(2) + e_i1(3) * n1(3) )    
temp2 = - ( e_t1(1) * n1(1) + e_t1(2) * n1(2) + e_t1(3) * n1(3) )

temp3 = inr * temp1 + temp2
temp4 = temp1 + inr * temp2

fresnel_r(1,1,1) = ( inr * temp1 - temp2 ) / temp3 
fresnel_r(1,1,2) = 0.0d0
fresnel_r(1,2,1) = 0.0d0
fresnel_r(1,2,2) = ( temp1 - inr * temp2 ) / temp4

fresnel_t(1,1) = dsqrt( 1.0d0 - fresnel_r(1,1,1) * fresnel_r(1,1,1) )
fresnel_t(1,2) = 0.0d0
fresnel_t(2,1) = 0.0d0
fresnel_t(2,2) = dsqrt( 1.0d0 - fresnel_r(1,2,2) * fresnel_r(1,2,2) )

!----------------------------------------------------------------------------------------!
! Calcul intermédiaire 2 pour la section efficace d'absorption.

temp1 = ea1(1,1)
temp2 = ea1(1,2)
temp3 = ea1(2,1)
temp4 = ea1(2,2)

ea1(1,1) = fresnel_t(1,1) * temp1 + fresnel_t(1,2) * temp2
ea1(1,2) = fresnel_t(2,1) * temp1 + fresnel_t(2,2) * temp2
ea1(2,1) = fresnel_t(1,1) * temp3 + fresnel_t(1,2) * temp4
ea1(2,2) = fresnel_t(2,1) * temp3 + fresnel_t(2,2) * temp4

temp1 = ea2(1,1)
temp2 = ea2(1,2)
temp3 = ea2(2,1)
temp4 = ea2(2,2)

ea2(1,1) = fresnel_t(1,1) * temp1 + fresnel_t(1,2) * temp2
ea2(1,2) = fresnel_t(2,1) * temp1 + fresnel_t(2,2) * temp2
ea2(2,1) = fresnel_t(1,1) * temp3 + fresnel_t(1,2) * temp4
ea2(2,2) = fresnel_t(2,1) * temp3 + fresnel_t(2,2) * temp4

!----------------------------------------------------------------------------------------!
! ray tracing condition d'arret.

rtcond = rtcond * 0.5d0 * ( fresnel_t(1,1) * fresnel_t(1,1) + fresnel_t(2,2) * fresnel_t(2,2) )                   
 
!----------------------------------------------------------------------------------------!
! Calcul intermediaire 1. 

tab1(1,1,1) = fresnel_r(1,1,1) * gamma1(1,1) 
tab1(1,1,2) = fresnel_r(1,1,1) * gamma1(1,2)
tab1(1,2,1) = fresnel_r(1,2,2) * gamma1(2,1)
tab1(1,2,2) = fresnel_r(1,2,2) * gamma1(2,2)

!----------------------------------------------------------------------------------------!
! Calcul de la matrice gammas.

beta_s(1) = beta_i(1)
beta_s(2) = beta_i(2)
beta_s(3) = beta_i(3)

alpha_s(1) = beta_s(2) * e_s1(3) - beta_s(3) * e_s1(2)
alpha_s(2) = beta_s(3) * e_s1(1) - beta_s(1) * e_s1(3)
alpha_s(3) = beta_s(1) * e_s1(2) - beta_s(2) * e_s1(1) 

beta_r1(1) = beta_i1(1)
beta_r1(2) = beta_i1(2)
beta_r1(3) = beta_i1(3)

alpha_r1(1) = beta_r1(2) * e_s1(3) - beta_r1(3) * e_s1(2)
alpha_r1(2) = beta_r1(3) * e_s1(1) - beta_r1(1) * e_s1(3)
alpha_r1(3) = beta_r1(1) * e_s1(2) - beta_r1(2) * e_s1(1) 

gammas(1,1) = alpha_s(1) * alpha_r1(1) + alpha_s(2) * alpha_r1(2) + alpha_s(3) * alpha_r1(3)
gammas(1,2) = alpha_s(1) * beta_r1(1)  + alpha_s(2) * beta_r1(2)  + alpha_s(3) * beta_r1(3)
gammas(2,1) = beta_s(1) * alpha_r1(1)  + beta_s(2) * alpha_r1(2)  + beta_s(3) * alpha_r1(3)
gammas(2,2) = beta_s(1) * beta_r1(1)   + beta_s(2) * beta_r1(2)   + beta_s(3) * beta_r1(3)

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 2. 

tab2(1,1,1) = gammas(1,1) * tab1(1,1,1) + gammas(1,2) * tab1(1,2,1)
tab2(1,1,2) = gammas(1,1) * tab1(1,1,2) + gammas(1,2) * tab1(1,2,2)
tab2(1,2,1) = gammas(2,1) * tab1(1,1,1) + gammas(2,2) * tab1(1,2,1)
tab2(1,2,2) = gammas(2,1) * tab1(1,1,2) + gammas(2,2) * tab1(1,2,2)
   
!----------------------------------------------------------------------------------------!
! Calcul de la boite correspondante a l'angle de diffusion.                    

teta = acos( e_s1(1)*ei(1) + e_s1(2)*ei(2) + e_s1(3)*ei(3) ) / cdr  ! Angle zenithal de diffusion (en degres) 

min = 180.0d0

do i = 1 , nb_teta

temp2 = abs( teta - tab_teta(i) )

if( temp2 < min )then
bin = i
min = temp2
end if

end do 

!----------------------------------------------------------------------------------------!
! Calcul de la matrice de transformation.

ma(1,1) = tab2(1,2,2)
ma(1,2) = tab2(1,1,1)
ma(1,3) = tab2(1,1,2)
ma(1,4) = tab2(1,2,1)

ma(2,1) = 0.0d0
ma(2,2) = 0.0d0
ma(2,3) = 0.0d0
ma(2,4) = 0.0d0

m1 = ma(1,1) * ma(1,1) + ma(2,1) * ma(2,1)
m2 = ma(1,2) * ma(1,2) + ma(2,2) * ma(2,2)
m3 = ma(1,3) * ma(1,3) + ma(2,3) * ma(2,3)
m4 = ma(1,4) * ma(1,4) + ma(2,4) * ma(2,4)

s21 = ma(1,2) * ma(1,1) + ma(2,1) * ma(2,2)
s23 = ma(1,2) * ma(1,3) + ma(2,3) * ma(2,2) 
s24 = ma(1,2) * ma(1,4) + ma(2,4) * ma(2,2)
s31 = ma(1,3) * ma(1,1) + ma(2,1) * ma(2,3) 
s34 = ma(1,3) * ma(1,4) + ma(2,4) * ma(2,3)
s41 = ma(1,4) * ma(1,1) + ma(2,1) * ma(2,4)

d21 = ma(1,2) * ma(2,1) - ma(1,1) * ma(2,2)
d23 = ma(1,2) * ma(2,3) - ma(1,3) * ma(2,2)
d24 = ma(1,2) * ma(2,4) - ma(1,4) * ma(2,2)
d31 = ma(1,3) * ma(2,1) - ma(1,1) * ma(2,3)
d34 = ma(1,3) * ma(2,4) - ma(1,4) * ma(2,3)
d41 = ma(1,4) * ma(2,1) - ma(1,1) * ma(2,4)

F11_RT(bin) = F11_RT(bin) + 0.5d0 * ( m2 + m3 + m4 + m1 )
F12_RT(bin) = F12_RT(bin) + 0.5d0 * ( m2 - m3 + m4 - m1 )
F13_RT(bin) = F13_RT(bin) + s23 + s41
F14_RT(bin) = F14_RT(bin) - d23 - d41
F21_RT(bin) = F21_RT(bin) + 0.5d0 *( m2 + m3 - m4 - m1 )
F22_RT(bin) = F22_RT(bin) + 0.5d0 *( m2 - m3 - m4 + m1 )
F23_RT(bin) = F23_RT(bin) + s23 - s41
F24_RT(bin) = F24_RT(bin) - d23 + d41
F31_RT(bin) = F31_RT(bin) + s24 + s31
F32_RT(bin) = F32_RT(bin) + s24 - s31
F33_RT(bin) = F33_RT(bin) + s21 + s34
F34_RT(bin) = F34_RT(bin) - d21 + d34
F41_RT(bin) = F41_RT(bin) + d24 + d31
F42_RT(bin) = F42_RT(bin) + d24 - d31
F43_RT(bin) = F43_RT(bin) + d21 + d34
F44_RT(bin) = F44_RT(bin) + s21 - s34

!----------------------------------------------------------------------------------------!
! Calcul des reflexions multiples.

beta_t1(1) = beta_i1(1)
beta_t1(2) = beta_i1(2)
beta_t1(3) = beta_i1(3)

alpha_t1(1) = beta_t1(2) * e_t1(3) - beta_t1(3) * e_t1(2)
alpha_t1(2) = beta_t1(3) * e_t1(1) - beta_t1(1) * e_t1(3)
alpha_t1(3) = beta_t1(1) * e_t1(2) - beta_t1(2) * e_t1(1)

beta(1) = beta_t1(1)
beta(2) = beta_t1(2)
beta(3) = beta_t1(3)

alpha(1) = alpha_t1(1)
alpha(2) = alpha_t1(2)
alpha(3) = alpha_t1(3)

tab(1,1,1) = fresnel_t(1,1) 
tab(1,1,2) = fresnel_t(1,2) 
tab(1,2,1) = fresnel_t(2,1) 
tab(1,2,2) = fresnel_t(2,2) 

tab(2,1,1) = 0.0d0
tab(2,1,2) = 0.0d0
tab(2,2,1) = 0.0d0
tab(2,2,2) = 0.0d0

e_ip(1) = e_t1(1)
e_ip(2) = e_t1(2)
e_ip(3) = e_t1(3)

dist_tot = 0.0d0

max_rm = 0

do while( rtcond > eps )           

!----------------------------------------------------------------------------------------!
! Calcul des points d'impacts 2,3,4...

distmin = distmax

do k = 1 , 8

temp1 = e_ip(1) * cn(k,1) + e_ip(2) * cn(k,2) + e_ip(3) * cn(k,3)

if( temp1 < 0.0 ) then

if( k == 7 ) then
distance = length2 + xrt * cn(k,1) + yrt * cn(k,2) + zrt * cn(k,3)
end if

if( k == 8 ) then
distance = length2 + xrt * cn(k,1) + yrt * cn(k,2) + zrt * cn(k,3)
end if

if( k /= 8 .and. k /= 7 ) then
distance = ray * const1 + xrt * cn(k,1) + yrt * cn(k,2) + zrt * cn(k,3) 
end if

distance = - distance / temp1

if( distance  < distmin ) then
distmin = distance
face = k
end if

end if

end do

np(1) = cn(face,1)
np(2) = cn(face,2)
np(3) = cn(face,3)

dist = distmin  

dist_tot = dist_tot + dist

absorp = exp( - ini * wk * dist )

xrt = xrt + dist * e_ip(1)
yrt = yrt + dist * e_ip(2)
zrt = zrt + dist * e_ip(3)

ps = e_ip(1) * np(1) + e_ip(2) * np(2) + e_ip(3) * np(3)

!----------------------------------------------------------------------------------------!
! Contribution à la section efficace d'absorption.

temp1 = ea1(1,1) * ea1(1,1) + ea1(1,2) * ea1(1,2) + ea1(2,1) * ea1(2,1) + ea1(2,2) * ea1(2,2) 
temp2 = ea2(1,1) * ea2(1,1) + ea2(1,2) * ea2(1,2) + ea2(2,1) * ea2(2,1) + ea2(2,2) * ea2(2,2) 

temp3 = dexp( - 2.0d0 * ini * wk * ( dist_tot - dist ) ) 
temp4 = dexp( - 2.0d0 * ini * wk * dist_tot )

abcs = abcs + ( temp3 - temp4 ) * ( temp1 + temp2 )

!----------------------------------------------------------------------------------------!
! Calcul de la matrice gammap.

temp1 = np(1) * e_ip(1) + np(2) * e_ip(2) + np(3) * e_ip(3)

if( temp1 * temp1 .gt. 0.9999999 )then
norme = 1.0
else
norme = dsqrt( 1.0 - temp1 * temp1 )
end if

beta_ip(1) = ( np(2) * e_ip(3) - np(3) * e_ip(2) ) / norme
beta_ip(2) = ( np(3) * e_ip(1) - np(1) * e_ip(3) ) / norme
beta_ip(3) = ( np(1) * e_ip(2) - np(2) * e_ip(1) ) / norme

alpha_ip(1) = beta_ip(2) * e_ip(3) - beta_ip(3) * e_ip(2)
alpha_ip(2) = beta_ip(3) * e_ip(1) - beta_ip(1) * e_ip(3)
alpha_ip(3) = beta_ip(1) * e_ip(2) - beta_ip(2) * e_ip(1) 

gammap(1,1) = alpha_ip(1) * alpha(1) + alpha_ip(2) * alpha(2) + alpha_ip(3) * alpha(3)
gammap(1,2) = alpha_ip(1) * beta(1)  + alpha_ip(2) * beta(2)  + alpha_ip(3) * beta(3)
gammap(2,1) = beta_ip(1) * alpha(1)  + beta_ip(2) * alpha(2)  + beta_ip(3) * alpha(3)
gammap(2,2) = beta_ip(1) * beta(1)   + beta_ip(2) * beta(2)   + beta_ip(3) * beta(3)

!----------------------------------------------------------------------------------------!
! Calcul intermédiaire 1 pour la section efficace d'absorption.

temp1 = ea1(1,1)
temp2 = ea1(1,2)
temp3 = ea1(2,1)
temp4 = ea1(2,2)

ea1(1,1) = gammap(1,1) * temp1 + gammap(1,2) * temp2
ea1(1,2) = gammap(2,1) * temp1 + gammap(2,2) * temp2
ea1(2,1) = gammap(1,1) * temp3 + gammap(1,2) * temp4
ea1(2,2) = gammap(2,1) * temp3 + gammap(2,2) * temp4

temp1 = ea2(1,1)
temp2 = ea2(1,2)
temp3 = ea2(2,1)
temp4 = ea2(2,2)

ea2(1,1) = gammap(1,1) * temp1 + gammap(1,2) * temp2
ea2(1,2) = gammap(2,1) * temp1 + gammap(2,2) * temp2
ea2(2,1) = gammap(1,1) * temp3 + gammap(1,2) * temp4
ea2(2,2) = gammap(2,1) * temp3 + gammap(2,2) * temp4

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 1.

tab1(1,1,1) = gammap(1,1) * tab(1,1,1) + gammap(1,2) * tab(1,2,1)
tab1(1,1,2) = gammap(1,1) * tab(1,1,2) + gammap(1,2) * tab(1,2,2)
tab1(1,2,1) = gammap(2,1) * tab(1,1,1) + gammap(2,2) * tab(1,2,1)
tab1(1,2,2) = gammap(2,1) * tab(1,1,2) + gammap(2,2) * tab(1,2,2)

tab1(2,1,1) = gammap(1,1) * tab(2,1,1) + gammap(1,2) * tab(2,2,1)
tab1(2,1,2) = gammap(1,1) * tab(2,1,2) + gammap(1,2) * tab(2,2,2)
tab1(2,2,1) = gammap(2,1) * tab(2,1,1) + gammap(2,2) * tab(2,2,1)
tab1(2,2,2) = gammap(2,1) * tab(2,1,2) + gammap(2,2) * tab(2,2,2)

!----------------------------------------------------------------------------------------!
! Calcul des coefficients de Fresnel.

temp1 = - ( e_ip(1) * np(1) + e_ip(2) * np(2) + e_ip(3) * np(3) )

if( temp1 * temp1 .gt. 0.9999999 )then
cond2 = 0.0d0
else
cond2 = inr * dsqrt( 1.0 - temp1 * temp1 ) 
end if

if( cond2 < 1.0 ) then 

!----------------------------------------------------------------------------------------!
! Passage glace -> air.

temp1 = dsqrt( 1.0 / ( inr * inr ) - 1.0d0 + ps * ps )      

e_sp(1) = inr * ( e_ip(1) - ps * np(1) - temp1 * np(1) )   
e_sp(2) = inr * ( e_ip(2) - ps * np(2) - temp1 * np(2) )
e_sp(3) = inr * ( e_ip(3) - ps * np(3) - temp1 * np(3) )

!----------------------------------------------------------------------------------------!
! Calcul de la matrice gamma1.

temp1 = ei(1) * e_sp(1) + ei(2) * e_sp(2) + ei(3) * e_sp(3)

if( temp1 * temp1 .gt. 0.9999999 )then

beta_i(1) = t12
beta_i(2) = t22
beta_i(3) = t32

alpha_i(1) = beta_i(2) * ei(3) - beta_i(3) * ei(2)
alpha_i(2) = beta_i(3) * ei(1) - beta_i(1) * ei(3)
alpha_i(3) = beta_i(1) * ei(2) - beta_i(2) * ei(1) 

gamma1(1,1) = alpha_i1(1) * alpha_i(1) + alpha_i1(2) * alpha_i(2) + alpha_i1(3) * alpha_i(3)
gamma1(1,2) = alpha_i1(1) * beta_i(1)  + alpha_i1(2) * beta_i(2)  + alpha_i1(3) * beta_i(3)
gamma1(2,1) = beta_i1(1) * alpha_i(1)  + beta_i1(2) * alpha_i(2)  + beta_i1(3) * alpha_i(3)
gamma1(2,2) = beta_i1(1) * beta_i(1)   + beta_i1(2) * beta_i(2)   + beta_i1(3) * beta_i(3)

else

norme = dsqrt( 1.0d0 - temp1 * temp1 )

beta_i(1) = ( ei(2) * e_sp(3) - ei(3) * e_sp(2) ) / norme
beta_i(2) = ( ei(3) * e_sp(1) - ei(1) * e_sp(3) ) / norme
beta_i(3) = ( ei(1) * e_sp(2) - ei(2) * e_sp(1) ) / norme

alpha_i(1) = beta_i(2) * ei(3) - beta_i(3) * ei(2)
alpha_i(2) = beta_i(3) * ei(1) - beta_i(1) * ei(3)
alpha_i(3) = beta_i(1) * ei(2) - beta_i(2) * ei(1) 

gamma1(1,1) = alpha_i1(1) * alpha_i(1) + alpha_i1(2) * alpha_i(2) + alpha_i1(3) * alpha_i(3)
gamma1(1,2) = alpha_i1(1) * beta_i(1)  + alpha_i1(2) * beta_i(2)  + alpha_i1(3) * beta_i(3)
gamma1(2,1) = beta_i1(1) * alpha_i(1)  + beta_i1(2) * alpha_i(2)  + beta_i1(3) * alpha_i(3)
gamma1(2,2) = beta_i1(1) * beta_i(1)   + beta_i1(2) * beta_i(2)   + beta_i1(3) * beta_i(3)

end if

!----------------------------------------------------------------------------------------!

temp1 = - ( e_ip(1) * np(1) + e_ip(2) * np(2) + e_ip(3) * np(3) )
temp2 = - ( e_sp(1) * np(1) + e_sp(2) * np(2) + e_sp(3) * np(3) )

temp3 = temp1 + inr * temp2
temp4 = inr * temp1 + temp2

fresnel_r(1,1,1) = ( temp1 - inr * temp2 ) / temp3
fresnel_r(1,1,2) = 0.0d0
fresnel_r(1,2,1) = 0.0d0
fresnel_r(1,2,2) = ( inr * temp1 - temp2 ) / temp3

if(fresnel_r(1,1,1)*fresnel_r(1,1,1).gt.1.0d0)then
fresnel_t(1,1) = 0.0d0
else
fresnel_t(1,1) = dsqrt( 1.0d0 - fresnel_r(1,1,1) * fresnel_r(1,1,1) )
end if

fresnel_t(1,2) = 0.0d0
fresnel_t(2,1) = 0.0d0

if(fresnel_r(1,2,2)*fresnel_r(1,2,2).gt.1.0)then
fresnel_t(2,2) = 0.0
else
fresnel_t(2,2) = dsqrt( 1.0d0 - fresnel_r(1,2,2) * fresnel_r(1,2,2) )
end if

!----------------------------------------------------------------------------------------!
! Calcul intermédiaire 2 pour la section efficace d'absorption.

temp1 = ea1(1,1)
temp2 = ea1(1,2)
temp3 = ea1(2,1)
temp4 = ea1(2,2)

ea1(1,1) = fresnel_r(1,1,1) * temp1 + fresnel_r(1,1,2) * temp2
ea1(1,2) = fresnel_r(1,2,1) * temp1 + fresnel_r(1,2,2) * temp2
ea1(2,1) = fresnel_r(1,1,1) * temp3 + fresnel_r(1,1,2) * temp4
ea1(2,2) = fresnel_r(1,2,1) * temp3 + fresnel_r(1,2,2) * temp4

temp1 = ea2(1,1)
temp2 = ea2(1,2)
temp3 = ea2(2,1)
temp4 = ea2(2,2)

ea2(1,1) = fresnel_r(1,1,1) * temp1 + fresnel_r(1,1,2) * temp2
ea2(1,2) = fresnel_r(1,2,1) * temp1 + fresnel_r(1,2,2) * temp2
ea2(2,1) = fresnel_r(1,1,1) * temp3 + fresnel_r(1,1,2) * temp4
ea2(2,2) = fresnel_r(1,2,1) * temp3 + fresnel_r(1,2,2) * temp4

!----------------------------------------------------------------------------------------!
! ray tracing condition d'arret.

rtcond = rtcond * absorp * 0.5 * ( fresnel_r(1,1,1) * fresnel_r(1,1,1) + fresnel_r(1,2,2) * fresnel_r(1,2,2) ) 

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 2. 

tab2(1,1,1) = fresnel_t(1,1) * tab1(1,1,1) + fresnel_t(1,2) * tab1(1,2,1)
tab2(1,1,2) = fresnel_t(1,1) * tab1(1,1,2) + fresnel_t(1,2) * tab1(1,2,2)
tab2(1,2,1) = fresnel_t(2,1) * tab1(1,1,1) + fresnel_t(2,2) * tab1(1,2,1)
tab2(1,2,2) = fresnel_t(2,1) * tab1(1,1,2) + fresnel_t(2,2) * tab1(1,2,2)

tab2(2,1,1) = fresnel_t(1,1) * tab1(2,1,1) + fresnel_t(1,2) * tab1(2,2,1)
tab2(2,1,2) = fresnel_t(1,1) * tab1(2,1,2) + fresnel_t(1,2) * tab1(2,2,2)
tab2(2,2,1) = fresnel_t(2,1) * tab1(2,1,1) + fresnel_t(2,2) * tab1(2,2,1)
tab2(2,2,2) = fresnel_t(2,1) * tab1(2,1,2) + fresnel_t(2,2) * tab1(2,2,2)

tab(1,1,1) = fresnel_r(1,1,1) * tab1(1,1,1) + fresnel_r(1,1,2) * tab1(1,2,1)
tab(1,1,2) = fresnel_r(1,1,1) * tab1(1,1,2) + fresnel_r(1,1,2) * tab1(1,2,2)
tab(1,2,1) = fresnel_r(1,2,1) * tab1(1,1,1) + fresnel_r(1,2,2) * tab1(1,2,1)
tab(1,2,2) = fresnel_r(1,2,1) * tab1(1,1,2) + fresnel_r(1,2,2) * tab1(1,2,2)

tab(2,1,1) = fresnel_r(1,1,1) * tab1(2,1,1) + fresnel_r(1,1,2) * tab1(2,2,1)
tab(2,1,2) = fresnel_r(1,1,1) * tab1(2,1,2) + fresnel_r(1,1,2) * tab1(2,2,2)
tab(2,2,1) = fresnel_r(1,2,1) * tab1(2,1,1) + fresnel_r(1,2,2) * tab1(2,2,1)
tab(2,2,2) = fresnel_r(1,2,1) * tab1(2,1,2) + fresnel_r(1,2,2) * tab1(2,2,2)

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 3.

tab3(1,1,1) = tab2(1,1,1) * gamma1(1,1) + tab2(1,1,2) * gamma1(2,1)
tab3(1,1,2) = tab2(1,1,1) * gamma1(1,2) + tab2(1,1,2) * gamma1(2,2)
tab3(1,2,1) = tab2(1,2,1) * gamma1(1,1) + tab2(1,2,2) * gamma1(2,1)
tab3(1,2,2) = tab2(1,2,1) * gamma1(1,2) + tab2(1,2,2) * gamma1(2,2)

tab3(2,1,1) = tab2(2,1,1) * gamma1(1,1) + tab2(2,1,2) * gamma1(2,1)
tab3(2,1,2) = tab2(2,1,1) * gamma1(1,2) + tab2(2,1,2) * gamma1(2,2)
tab3(2,2,1) = tab2(2,2,1) * gamma1(1,1) + tab2(2,2,2) * gamma1(2,1)
tab3(2,2,2) = tab2(2,2,1) * gamma1(1,2) + tab2(2,2,2) * gamma1(2,2)

!----------------------------------------------------------------------------------------!
! Calcul de la matrice gammas.

beta_s(1) = beta_i(1)
beta_s(2) = beta_i(2)
beta_s(3) = beta_i(3)

alpha_s(1) = beta_s(2) * e_sp(3) - beta_s(3) * e_sp(2)
alpha_s(2) = beta_s(3) * e_sp(1) - beta_s(1) * e_sp(3)
alpha_s(3) = beta_s(1) * e_sp(2) - beta_s(2) * e_sp(1) 

beta_tp(1) = beta_ip(1)
beta_tp(2) = beta_ip(2)
beta_tp(3) = beta_ip(3)

alpha_tp(1) = beta_tp(2) * e_sp(3) - beta_tp(3) * e_sp(2)
alpha_tp(2) = beta_tp(3) * e_sp(1) - beta_tp(1) * e_sp(3)
alpha_tp(3) = beta_tp(1) * e_sp(2) - beta_tp(2) * e_sp(1) 

gammas(1,1) = alpha_s(1) * alpha_tp(1) + alpha_s(2) * alpha_tp(2) + alpha_s(3) * alpha_tp(3)
gammas(1,2) = alpha_s(1) * beta_tp(1)  + alpha_s(2) * beta_tp(2)  + alpha_s(3) * beta_tp(3)
gammas(2,1) = beta_s(1) * alpha_tp(1)  + beta_s(2) * alpha_tp(2)  + beta_s(3) * alpha_tp(3)
gammas(2,2) = beta_s(1) * beta_tp(1)   + beta_s(2) * beta_tp(2)   + beta_s(3) * beta_tp(3)

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 4. 

tab4(1,1,1) = gammas(1,1) * tab3(1,1,1) + gammas(1,2) * tab3(1,2,1)
tab4(1,1,2) = gammas(1,1) * tab3(1,1,2) + gammas(1,2) * tab3(1,2,2)
tab4(1,2,1) = gammas(2,1) * tab3(1,1,1) + gammas(2,2) * tab3(1,2,1)
tab4(1,2,2) = gammas(2,1) * tab3(1,1,2) + gammas(2,2) * tab3(1,2,2)

tab4(2,1,1) = gammas(1,1) * tab3(2,1,1) + gammas(1,2) * tab3(2,2,1)
tab4(2,1,2) = gammas(1,1) * tab3(2,1,2) + gammas(1,2) * tab3(2,2,2)
tab4(2,2,1) = gammas(2,1) * tab3(2,1,1) + gammas(2,2) * tab3(2,2,1)
tab4(2,2,2) = gammas(2,1) * tab3(2,1,2) + gammas(2,2) * tab3(2,2,2)

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 5. 

temp1 = exp( - ini * wk * dist_tot )

tab5(1,1,1) = temp1 * tab4(1,1,1) 
tab5(1,1,2) = temp1 * tab4(1,1,2) 
tab5(1,2,1) = temp1 * tab4(1,2,1)
tab5(1,2,2) = temp1 * tab4(1,2,2) 

tab5(2,1,1) = temp1 * tab4(2,1,1)
tab5(2,1,2) = temp1 * tab4(2,1,2) 
tab5(2,2,1) = temp1 * tab4(2,2,1) 
tab5(2,2,2) = temp1 * tab4(2,2,2) 

!----------------------------------------------------------------------------------------!
! Calcul de la boite correspondante a l'angle de diffusion.

teta = acos( e_sp(1) * ei(1) + e_sp(2) * ei(2) + e_sp(3) * ei(3) ) / cdr  

min = 180.0d0

do i = 1 , nb_teta

temp2 = abs( teta - tab_teta(i) )

if( temp2 < min )then
bin = i
min = temp2
end if

end do 

!----------------------------------------------------------------------------------------!
! Calcul de la matrice de transformation.

ma(1,1) = tab5(1,2,2)
ma(1,2) = tab5(1,1,1)
ma(1,3) = tab5(1,1,2)
ma(1,4) = tab5(1,2,1)

ma(2,1) = tab5(2,2,2)
ma(2,2) = tab5(2,1,1)
ma(2,3) = tab5(2,1,2)
ma(2,4) = tab5(2,2,1)

m1 = ma(1,1) * ma(1,1) + ma(2,1) * ma(2,1)
m2 = ma(1,2) * ma(1,2) + ma(2,2) * ma(2,2)
m3 = ma(1,3) * ma(1,3) + ma(2,3) * ma(2,3)
m4 = ma(1,4) * ma(1,4) + ma(2,4) * ma(2,4)

s21 = ma(1,2) * ma(1,1) + ma(2,1) * ma(2,2)
s23 = ma(1,2) * ma(1,3) + ma(2,3) * ma(2,2) 
s24 = ma(1,2) * ma(1,4) + ma(2,4) * ma(2,2)
s31 = ma(1,3) * ma(1,1) + ma(2,1) * ma(2,3) 
s34 = ma(1,3) * ma(1,4) + ma(2,4) * ma(2,3)
s41 = ma(1,4) * ma(1,1) + ma(2,1) * ma(2,4)

d21 = ma(1,2) * ma(2,1) - ma(1,1) * ma(2,2)
d23 = ma(1,2) * ma(2,3) - ma(1,3) * ma(2,2)
d24 = ma(1,2) * ma(2,4) - ma(1,4) * ma(2,2)
d31 = ma(1,3) * ma(2,1) - ma(1,1) * ma(2,3)
d34 = ma(1,3) * ma(2,4) - ma(1,4) * ma(2,3)
d41 = ma(1,4) * ma(2,1) - ma(1,1) * ma(2,4)

F11_RT(bin) = F11_RT(bin) + 0.5d0 * ( m2 + m3 + m4 + m1 )
F12_RT(bin) = F12_RT(bin) + 0.5d0 * ( m2 - m3 + m4 - m1 )
F13_RT(bin) = F13_RT(bin) + s23 + s41
F14_RT(bin) = F14_RT(bin) - d23 - d41
F21_RT(bin) = F21_RT(bin) + 0.5d0 *( m2 + m3 - m4 - m1 )
F22_RT(bin) = F22_RT(bin) + 0.5d0 *( m2 - m3 - m4 + m1 )
F23_RT(bin) = F23_RT(bin) + s23 - s41
F24_RT(bin) = F24_RT(bin) - d23 + d41
F31_RT(bin) = F31_RT(bin) + s24 + s31
F32_RT(bin) = F32_RT(bin) + s24 - s31
F33_RT(bin) = F33_RT(bin) + s21 + s34
F34_RT(bin) = F34_RT(bin) - d21 + d34
F41_RT(bin) = F41_RT(bin) + d24 + d31
F42_RT(bin) = F42_RT(bin) + d24 - d31
F43_RT(bin) = F43_RT(bin) + d21 + d34
F44_RT(bin) = F44_RT(bin) + s21 - s34

!----------------------------------------------------------------------------------------!
! Reflexion totale.

else  

temp1  = - ( e_ip(1) * np(1) + e_ip(2) * np(2) + e_ip(3) * np(3) )
temp2 =   cond2 * cond2 - 1.0
temp3 =   ( ( 1.0 / inr ) * temp1 )**2
temp4 = - 2.0 * ( 1.0 / inr ) * temp1 * sqrt( temp2 )
temp5 =   1.0 - ( 1.0 / inr )**2

fresnel_r(1,1,1) = ( temp3 - temp2 ) / ( temp3 + temp2 ) 
fresnel_r(1,1,2) = 0.0
fresnel_r(1,2,1) = 0.0
fresnel_r(1,2,2) = ( 2.0 * temp1 * temp1 ) / temp5 - 1.0

fresnel_r(2,1,1) = temp4 / ( temp3 + temp2 )
fresnel_r(2,1,2) = 0.0
fresnel_r(2,2,1) = 0.0
fresnel_r(2,2,2) = temp4 / temp5

fresnel_t(1,1) = 0.0
fresnel_t(1,2) = 0.0
fresnel_t(2,1) = 0.0
fresnel_t(2,2) = 0.0

!----------------------------------------------------------------------------------------!
! Calcul intermédiaire 2 pour la section efficace d'absorption.

temp1 = ea1(1,1)
temp2 = ea1(1,2)
temp3 = ea1(2,1)
temp4 = ea1(2,2)

temp5 = fresnel_r(1,1,1) * temp1 + fresnel_r(1,1,2) * temp2
temp6 = fresnel_r(2,1,1) * temp3 + fresnel_r(2,1,2) * temp4

ea1(1,1) = temp5 - temp6

temp5 = fresnel_r(1,2,1) * temp1 + fresnel_r(1,2,2) * temp2
temp6 = fresnel_r(2,2,1) * temp3 + fresnel_r(2,2,2) * temp4

ea1(1,2) = temp5 - temp6

temp5 = fresnel_r(1,1,1) * temp3 + fresnel_r(1,1,2) * temp4
temp6 = fresnel_r(2,1,1) * temp1 + fresnel_r(2,1,2) * temp2

ea1(2,1) = temp5 + temp6

temp5 = fresnel_r(1,2,1) * temp3 + fresnel_r(1,2,2) * temp4
temp6 = fresnel_r(2,2,1) * temp1 + fresnel_r(2,2,2) * temp2

ea1(2,2) = temp5 + temp6

!----------------------------------------------------------------------------------------!

temp1 = ea2(1,1)
temp2 = ea2(1,2)
temp3 = ea2(2,1)
temp4 = ea2(2,2)

temp5 = fresnel_r(1,1,1) * temp1 + fresnel_r(1,1,2) * temp2
temp6 = fresnel_r(2,1,1) * temp3 + fresnel_r(2,1,2) * temp4

ea2(1,1) = temp5 - temp6

temp5 = fresnel_r(1,2,1) * temp1 + fresnel_r(1,2,2) * temp2
temp6 = fresnel_r(2,2,1) * temp3 + fresnel_r(2,2,2) * temp4

ea2(1,2) = temp5 - temp6

temp5 = fresnel_r(1,1,1) * temp3 + fresnel_r(1,1,2) * temp4
temp6 = fresnel_r(2,1,1) * temp1 + fresnel_r(2,1,2) * temp2

ea2(2,1) = temp5 + temp6

temp5 = fresnel_r(1,2,1) * temp3 + fresnel_r(1,2,2) * temp4
temp6 = fresnel_r(2,2,1) * temp1 + fresnel_r(2,2,2) * temp2

ea2(2,2) = temp5 + temp6

!----------------------------------------------------------------------------------------!
! Calcul intermediaire 5. 

taba(1,1) = fresnel_r(1,1,1) * tab1(1,1,1) 
taba(1,2) = fresnel_r(1,1,1) * tab1(1,1,2) 
taba(2,1) = fresnel_r(1,2,2) * tab1(1,2,1)
taba(2,2) = fresnel_r(1,2,2) * tab1(1,2,2)

tabb(1,1) = fresnel_r(2,1,1) * tab1(2,1,1) 
tabb(1,2) = fresnel_r(2,1,1) * tab1(2,1,2) 
tabb(2,1) = fresnel_r(2,2,2) * tab1(2,2,1)
tabb(2,2) = fresnel_r(2,2,2) * tab1(2,2,2)

tabc(1,1) = fresnel_r(1,1,1) * tab1(2,1,1) 
tabc(1,2) = fresnel_r(1,1,1) * tab1(2,1,2)
tabc(2,1) = fresnel_r(1,2,2) * tab1(2,2,1)
tabc(2,2) = fresnel_r(1,2,2) * tab1(2,2,2)

tabd(1,1) = fresnel_r(2,1,1) * tab1(1,1,1) 
tabd(1,2) = fresnel_r(2,1,1) * tab1(1,1,2) 
tabd(2,1) = fresnel_r(2,2,2) * tab1(1,2,1)
tabd(2,2) = fresnel_r(2,2,2) * tab1(1,2,2)

tab(1,1,1) = taba(1,1) - tabb(1,1)
tab(1,1,2) = taba(1,2) - tabb(1,2)
tab(1,2,1) = taba(2,1) - tabb(2,1)
tab(1,2,2) = taba(2,2) - tabb(2,2)

tab(2,1,1) = tabc(1,1) + tabd(1,1)
tab(2,1,2) = tabc(1,2) + tabd(1,2)
tab(2,2,1) = tabc(2,1) + tabd(2,1)
tab(2,2,2) = tabc(2,2) + tabd(2,2)

end if

!----------------------------------------------------------------------------------------!

e_ip(1) = e_ip(1) - 2.0 * ps * np(1) 
e_ip(2) = e_ip(2) - 2.0 * ps * np(2)
e_ip(3) = e_ip(3) - 2.0 * ps * np(3)

beta_rp(1) = beta_ip(1)
beta_rp(2) = beta_ip(2)
beta_rp(3) = beta_ip(3)

alpha_rp(1) = beta_rp(2) * e_ip(3) - beta_rp(3) * e_ip(2)
alpha_rp(2) = beta_rp(3) * e_ip(1) - beta_rp(1) * e_ip(3)
alpha_rp(3) = beta_rp(1) * e_ip(2) - beta_rp(2) * e_ip(1)

beta(1) = beta_rp(1)
beta(2) = beta_rp(2)
beta(3) = beta_rp(3)

alpha(1) = alpha_rp(1)
alpha(2) = alpha_rp(2)
alpha(3) = alpha_rp(3)

!----------------------------------------------------------------------------------------!
! Fin des reflexions multiples. 

max_rm = max_rm + 1

if( max_rm == 10 )then
rtcond = 0.0d0
end if

end do                                          

!----------------------------------------------------------------------------------------!
! Fin de la boucle sur les rayons.

end do

!write(6,*) ' Nombre de realisations terminees : ' , n , 'sur' , nor

!----------------------
! Fin de la boucle sur les orientations aleatoire.

end do  

!----------------------------------------------------------------------------------------!
! Calcul de la matrice de phase pour l'ensemble de particules.

do i = 1 , nb_teta 

if(i==1)then

tetai = ( tab_teta(i) + 0.0d0 ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

else if(i==nb_teta)then

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + 180.0d0 ) / 2.0d0

else

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

end if

diff_solid_angle = ( cos( tetai * cdr ) - cos( tetaf * cdr ) )

F11_RT(i) = F11_RT(i) / diff_solid_angle
F12_RT(i) = F12_RT(i) / diff_solid_angle
F22_RT(i) = F22_RT(i) / diff_solid_angle
F33_RT(i) = F33_RT(i) / diff_solid_angle
F43_RT(i) = F43_RT(i) / diff_solid_angle
F44_RT(i) = F44_RT(i) / diff_solid_angle

end do    

!----------------------------------------------------------------------------------------!
! Calcul du paramètre norme_rt : 

norme_rt = 0.0

do i = 1 , nb_teta 

if(i==1)then

tetai = ( tab_teta(i) + 0.0d0 ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

else if(i==nb_teta)then

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + 180.0d0 ) / 2.0d0

else

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

end if

norme_rt = norme_rt + 0.5d0 * F11_RT(i) * ( cos( tetai * cdr ) - cos( tetaf * cdr ) )

end do

!----------------------------------------------------------------------------------------!
! Normalized phase matrix is given by :

do i = 1 , nb_teta

P11_RT(i) = F11_RT(i) / norme_rt
P12_RT(i) = F12_RT(i) / norme_rt
P22_RT(i) = F22_RT(i) / norme_rt
P33_RT(i) = F33_RT(i) / norme_rt
P43_RT(i) = F43_RT(i) / norme_rt
P44_RT(i) = F44_RT(i) / norme_rt

end do

!----------------------------------------------------------------------------------------!
! Todo: removed this part:
! Calcul de la fraction de l'énergie de la delta-transmission
! sur l'energie totale diffusée.

!norme_rt = 0.0

!do i = 2 , nb_teta 

!if(i==1)then

!tetai = ( tab_teta(i) + 0.0d0 ) / 2.0d0 
!tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

!else if(i==nb_teta)then

!tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
!tetaf = ( tab_teta(i) + 180.0d0 ) / 2.0d0

!else

!tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
!tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

!end if

!norme_rt = norme_rt + P11_RT(i) * ( cos( tetai * cdr ) - cos( tetaf * cdr ) )

!end do

!f_delta = 1.0d0 - 0.5 * norme_rt 

f_delta = 0.0

!----------------------------------------------------------------------------------------!
! Ecriture des resultats.

open( unit = 1,file = 'hexagomal_pristine_crystal_RT.txt')

do i = 1 , nb_teta

write(1,*) tab_teta(i),P11_RT(i),-P12_RT(i)/P11_RT(i),P22_RT(i)/P11_RT(i)&
&          ,P33_RT(i)/P11_RT(i),-P43_RT(i)/P11_RT(i),P44_RT(i)/P11_RT(i)

end do

close(1)

!----------------------------------------------------------------------------------------!
! Calcul de la matrice de diffusion totale.

abcs = 0.5d0 * rcsi * abcs

scs = 2.0d0 * rcs_tot - abcs
ecs = 2.0d0 * rcs_tot

albedo = scs / ecs

!----------------------------------------------------------------------------------------!

temp1 = 1.0d0 / ( 2.0d0 * albedo )

TAB_temp2(1) = ( ( 1.0d0 - f_delta ) - temp1 ) * P11_RT(1)
TAB_temp2(2) = ( ( 1.0d0 - f_delta ) - temp1 ) * P12_RT(1)
TAB_temp2(3) = ( ( 1.0d0 - f_delta ) - temp1 ) * P22_RT(1)
TAB_temp2(4) = ( ( 1.0d0 - f_delta ) - temp1 ) * P33_RT(1)
TAB_temp2(5) = ( ( 1.0d0 - f_delta ) - temp1 ) * P43_RT(1)
TAB_temp2(6) = ( ( 1.0d0 - f_delta ) - temp1 ) * P44_RT(1)

TAB_temp3(1) = temp1 * P11_MR(1)
TAB_temp3(2) = temp1 * P12_MR(1)
TAB_temp3(3) = temp1 * P22_MR(1)
TAB_temp3(4) = temp1 * P33_MR(1)
TAB_temp3(5) = temp1 * P43_MR(1)
TAB_temp3(6) = temp1 * P44_MR(1)

tetai = ( tab_teta(1) + 0.0d0 ) / 2.0d0 
tetaf = ( tab_teta(1) + tab_teta(2) ) / 2.0d0

num_dirac_function = 1.0d0 / ( dcos( tetai * cdr ) - dcos( tetaf * cdr ) )

temp4 = 2.0d0 * num_dirac_function * f_delta 

P11(1) = TAB_temp2(1) + TAB_temp3(1) + temp4
P12(1) = TAB_temp2(2) + TAB_temp3(2) 
P22(1) = TAB_temp2(3) + TAB_temp3(3) + temp4
P33(1) = TAB_temp2(4) + TAB_temp3(4) + temp4
P43(1) = TAB_temp2(5) + TAB_temp3(5) 
P44(1) = TAB_temp2(6) + TAB_temp3(6) + temp4

do i = 2 , nb_teta

TAB_temp2(1) = ( ( 1.0d0 - f_delta ) - temp1 ) * P11_RT(i)
TAB_temp2(2) = ( ( 1.0d0 - f_delta ) - temp1 ) * P12_RT(i)
TAB_temp2(3) = ( ( 1.0d0 - f_delta ) - temp1 ) * P22_RT(i)
TAB_temp2(4) = ( ( 1.0d0 - f_delta ) - temp1 ) * P33_RT(i)
TAB_temp2(5) = ( ( 1.0d0 - f_delta ) - temp1 ) * P43_RT(i)
TAB_temp2(6) = ( ( 1.0d0 - f_delta ) - temp1 ) * P44_RT(i)

TAB_temp3(1) = temp1 * P11_MR(i)
TAB_temp3(2) = temp1 * P12_MR(i)
TAB_temp3(3) = temp1 * P22_MR(i)
TAB_temp3(4) = temp1 * P33_MR(i)
TAB_temp3(5) = temp1 * P43_MR(i)
TAB_temp3(6) = temp1 * P44_MR(i)

P11(i) = TAB_temp2(1) + TAB_temp3(1) 
P12(i) = TAB_temp2(2) + TAB_temp3(2) 
P22(i) = TAB_temp2(3) + TAB_temp3(3) 
P33(i) = TAB_temp2(4) + TAB_temp3(4) 
P43(i) = TAB_temp2(5) + TAB_temp3(5) 
P44(i) = TAB_temp2(6) + TAB_temp3(6) 

end do

!----------------------------------------------------------------------------------------!

open( unit = 1,file = 'hexagomal_pristine_crystal_GOM.txt')

write(1,*) 'ray, length, length/(2ray) :'
write(1,*) ray,length,length / ( 2.0d0 * ray )
write(1,*) 'wavelength :'
write(1,*) wl
write(1,*) 'refractive index :'
write(1,*) inr,ini
write(1,*) 'extinction cross section :'
write(1,*) ecs
write(1,*) 'scattering cross section :'
write(1,*) scs
write(1,*) 'absorption cross section :'
write(1,*) abcs
write(1,*) 'single scatter albedo :'
write(1,*) albedo
write(1,*) 'delta transmission fraction :'
write(1,*) f_delta
write(1,*) 'phase function (PF) values as fct of scattering angle :'
write(1,*)

do i = 1 , nb_teta

write(1,*) tab_teta(i),P11(i),-P12(i)/P11(i),P22(i)/P11(i)&
&          ,P33(i)/P11(i),-P43(i)/P11(i),P44(i)/P11(i)

end do

close(1)

!----------------------------------------------------------------------------------------!
! Verification de la norme

norme = 0.0

do i = 1 , nb_teta 

if(i==1)then

tetai = ( tab_teta(i) + 0.0d0 ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

else if(i==nb_teta)then

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + 180.0d0 ) / 2.0d0

else

tetai = ( tab_teta(i) + tab_teta(i-1) ) / 2.0d0 
tetaf = ( tab_teta(i) + tab_teta(i+1) ) / 2.0d0

end if

norme = norme + P11(i) * ( cos( tetai * cdr ) - cos( tetaf * cdr ) )

end do

write(6,*) 'norme',norme

close(1)

!----------------------------------------------------------------------------------------!
! Main Code end

contains

!----------------------------------------------------------------------------------------!
! Function sinus cardinal 

real function fsc(x) 

real(8) :: x

if( abs(x) > 0.001d0 ) then
fsc = sin( x ) / x
else
fsc = 1.0d0 - x * x / 6.0d0
end if

return
end function fsc

end program gom_phm

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2017-04-25

Edited / Modifié: 2017-04-25

License / Licence: MIT License

Abstract / Résumé