Weighted random sampling with a reservoir example in fortran 90

					
program WeightRandomSamplingReservoir

implicit none 

integer :: i,j

integer*4 timeArray(3) ! Holds the hour, minute, and second

integer, dimension(100,2) :: V ! Population

real, dimension(10,2) :: R ! Reservoir

integer, dimension(1) :: MinIndex
integer, dimension(2) :: R_dim

real :: ui,ki
real :: ki_min

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

call itime(timeArray) ! Generate sequence of random number 

i = rand ( timeArray(1)+timeArray(2)+timeArray(3) )

R_dim = shape(R)

!----------------------------------------------------------------------------------------!
! Population V: (vi,wi) 

 write(6,*) '!----- Population -----!'

V = 0 
do i = 1, 100
  V(i,1) = i
  if( i < 90 )then
    V(i,2) = 1
  else
    V(i,2) = 2
  end if
  write(6,*) 'Data, Weight: ',V(i,1),V(i,2)
end do

!----------------------------------------------------------------------------------------!
! Weighted Random Sampling

R = 9999.0
ki_min = 9999.0

ui = rand(0)
R(1,1) = V(1,1)
R(1,2) = ui**(1.0/V(1,2))
ki_min = R(1,2)
do i = 2, R_dim(1)
  ui = rand(0)
  R(i,1) = V(i,1)
  R(i,2) = ui**(1.0/V(i,2))
  if( R(i,2) < ki_min ) ki_min = R(i,2)
end do

do i = R_dim(1)+1, 100
  ui = rand(0)
  ki = ui**(1.0/V(i,2))

  !write(6,*) V(i,1),V(i,2),ki
  !write(6,*) (R(j,1), j=1,10)
  !write(6,*) (R(j,2), j=1,10)
  !write(6,*) minloc(R(:,2)),'dimension',shape(minloc(R)),MinIndex(1)
  !write(6,*) "------------"

  MinIndex = minloc(R(:,2))

  if( ki > ki_min )then
    j = MinIndex(1)
    R(j,1) = V(i,1)
    R(j,2) = ki

    MinIndex = minloc(R(:,2))
    j = MinIndex(1)
    ki_min = R(j,2)
  end if

  !if( rand(0) < 0.2 )then
  !write(6,*) D(i,1),D(i,2)
  !end if
end do

!----------------------------------------------------------------------------------------!
! Print results

 write(6,*) '!----- Sampling -----!'

do i=1,R_dim(1)

	write(6,*) R(i,1), R(i,2)

end do

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

end program WeightRandomSamplingReservoir

Cholesky decomposition in fortran 90

					
program choleski

implicit none

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

integer :: i,j,k

real(pr) :: a,b,c,d

real(pr), dimension(4,4) :: MA,ML,MLT,R

real(pr), dimension(4) :: MX,MY,MB

data MA / 1.0, 1.0, 1.0, 1.0, &
        & 1.0, 5.0, 5.0, 5.0, &
        & 1.0, 5.0, 14.0, 14.0, &
        & 1.0, 5.0, 14.0, 15.0  /        

data MB / 9.0, 25.0, 43.0, 44.0/

!----------------------------------------------------------------------------------------!
! Print A MAtrix

write(6,*) '---------- Matrix A ----------'

do i = 1, 4
	write(6,*) MA(i,1),MA(i,2),MA(i,3),MA(i,4)
end do

!----------------------------------------------------------------------------------------!
! Choleski factorisation 

ML = 0.0

ML(1,1) = sqrt( MA(1,1) )
do i = 2 , 4
	ML(i,1) = MA(1,i) / ML(1,1)
enddo

do j = 2 , 4
	do i = 1 , 4
		if( i .lt. j ) then
			ML(i,j) = 0.0
		else
			if( i .eq. j) then
				a = 0.0
				do k = 1 , i -1
					a = a + ML(i,k)**2
				enddo
				ML(i,j) = sqrt( MA(i,i) - a )
			else
				a = 0.0
				do k = 1 , i - 1
					a = a + ML(j,k) * ML(i,k)
				enddo
				ML(i,j) = ( MA(j,i) - a ) / ML(j,j)
			endif
		endif
	enddo
enddo

!----------------------------------------------------------------------------------------!
! Sanity Check 

MLT = transpose(ML)

R = matmul(ML,MLT)

write(6,*) '---------- Matrix R = L * Transpose(L) ----------'

do i = 1, 4
	write(6,*) R(i,1),R(i,2),R(i,3),R(i,4)
end do

!----------------------------------------------------------------------------------------!
! Solve linear system 1

MY(1) = MB(1) / ML(1,1)
do i = 2 , 4
	a = 0.0
	do k = 1 , i-1
		a = a + ML(i,k) * MY(k)
	enddo
	MY(i) = ( MB(i) - a ) / ML(i,i)
enddo

!----------------------------------------------------------------------------------------!
! Solve linear system 2

MX(4) = MY(4) / ML(4,4)
do i = 4-1, 1, -1
	a = 0.0
	do k = 4, i+1, -1
		a = a + ML(k,i) * MX(k)
	enddo
	MX(i) = ( MY(i) - a ) / ML(i,i)
enddo

write(6,*) '---------- Solution A*X = B ----------'

write(6,*) MX(1),MX(2),MX(3),MX(4)

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

end program choleski

Plot phase function and its Legendre polynomial expansion with python

					
#!/usr/bin/env python

from scipy.special import legendre

import math
import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------------------------#

data  = np.loadtxt("phase_function_legendre_polynomial_expansion_02.txt")

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,1],'r',
      label='P11 phase function reconstructed \n from its Legendre polynomial expansion')
plt.plot(data[:,0],data[:,2],'k--',label='P11 true phase function')

plt.legend()
plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.ylabel(r"$P_{11}$")
plt.title('Phase function and its Legendre polynomial expansion \n (n=120)')

plt.ylim(10**(-2.0),10**6)
plt.yscale('log')
plt.grid(True,which="both")

plt.savefig('phase_function_legendre_polynomial_expansion.png', bbox_inches='tight')
plt.show()

#----------------------------------------------------------------------------------------#

Plot MODIS C6 MYD06 Cloud Fraction with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar
import pprint

#----------------------------------------------------------------------------------------#
# Inputs

file_name = './inputs/myd06_granules_viz/MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Cloud_Fraction')

data = data_selected_id.get()

data_selected_attributes = data_selected_id.attributes()

#pprint.pprint( data_selected_attributes )

scale_factor = data_selected_attributes['scale_factor']
add_offset = data_selected_attributes['add_offset']

#----------------------------------------------------------------------------------------#
# Plot data

data = data * scale_factor

data_shape = data.shape

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0, 0.1, 0.2 , 0.3 , 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')

cbar_bounds = [0.0, 0.1, 0.2 , 0.3 , 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
cbar_ticks = [0.0, 0.1, 0.2 , 0.3 , 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]             
               
cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=10)

l = [int(i) for i in np.linspace(0,data_shape[1],6)]
plt.xticks(l, [i for i in reversed(l)], rotation=0, fontsize=8 )

l = [int(i) for i in np.linspace(0,data_shape[0],9)]
plt.yticks(l, l, rotation=0, fontsize=8 )

title = 'Cloud Fraction \n MYD06 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

plt.savefig("myd06_c6_cloud_fraction_data.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Plot MODIS C6 MYD06 Cloud Water Path with python

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar
import pprint

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD06_L2.A2008015.1435.006.2013342100940.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

data_selected_id = file.select('Cloud_Water_Path')

data = data_selected_id.get()

data_selected_attributes = data_selected_id.attributes()

#pprint.pprint( data_selected_attributes )

scale_factor = data_selected_attributes['scale_factor']
add_offset = data_selected_attributes['add_offset']

#----------------------------------------------------------------------------------------#
# Plot data

data = data * scale_factor

data_shape = data.shape

max_value = 500.0 #np.max( data )
min_value = data[data > 0.0].min( )

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0, 1.0, 100.0, 200.0, 300.0, 400.0, 500.0]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')

cbar_bounds = [ 1.0, 100.0, 200.0, 300.0, 400.0, 500.0]
cbar_ticks = [1.0, 100.0, 200.0, 300.0, 400.0, 500.0]              
cbar_labels = ['1.0', '100.0', '200.0', '300.0', '400.0', '> 500.0']   
               
cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_labels, fontsize=10)

l = [int(i) for i in np.linspace(0,data_shape[1],6)]
plt.xticks(l, [i for i in reversed(l)], rotation=0, fontsize=8 )
plt.yticks([0,250,500,750,1000,1250,1500,1750,2000], 
           ['0','250','500', '750', '1000','1250','1500','1750','2000'], 
           rotation=0, fontsize=8 )

title = 'Cloud Water Path (g/m2) \n MYD06 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

plt.savefig("myd06_c6_cloud_water_path_data.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Plot MODIS C6 MYD04 Optical Depth Land And Ocean with python (+ add a marker)

					
#!/usr/bin/env python

from pyhdf.SD import SD, SDC 

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import matplotlib.cm as cm
import calendar
import pprint

#----------------------------------------------------------------------------------------#
# Inputs

file_name = 'MYD04_L2.A2015190.0010.006.2015190191844.hdf'

#----------------------------------------------------------------------------------------#
# Extract Year, Month, Day, Hour and Minutes from file name

path = file_name.split('/')
name = path[len(path)-1]
list = name.split('.')

hhmm = list[2]
hh = int(hhmm[0:2])
mm = int(hhmm[2:4])

year = int(list[1][1:5])
dayofyear = int(list[1][5:8])

month = 1
day = dayofyear
while day - calendar.monthrange(year,month)[1] > 0 and month <= 12:
	day = day - calendar.monthrange(year,month)[1]
	month = month + 1

#print year, month, day, hh, mm

#----------------------------------------------------------------------------------------#
# Read HDF Files

file = SD(file_name, SDC.READ)

longitude_id = file.select('Longitude')
latitude_id = file.select('Latitude')
data_selected_id = file.select('Optical_Depth_Land_And_Ocean')

longitude_attributes = longitude_id.attributes()
latitude_attributes = latitude_id.attributes()
attributes = data_selected_id.attributes()

longitude = longitude_id.get()
latitude = latitude_id.get()
data = data_selected_id.get()

_FillValue = attributes['_FillValue']

#----------------------------------------------------------------------------------------#
# Find longitude and latitude indexes

long_Halemaumau = -155.28
lat_Halemaumau = 19.41

X = np.sqrt( np.square( longitude - long_Halemaumau ) + \
             np.square( latitude - lat_Halemaumau ) )
LongLat_idx = np.where( X == X.min() )

#----------------------------------------------------------------------------------------#
# Plot data

data[ data == _FillValue ] = 0.0

data = data * attributes['scale_factor']

cmap = [(0.0,0.0,0.0)] + [(cm.jet(i)) for i in xrange(1,256)] 
cmap = mpl.colors.ListedColormap(cmap)

bounds = [0.0,0.0001,0.1,0.2,0.3,0.4,0.5]

norm = mpl.colors.BoundaryNorm(bounds, cmap.N)

img = plt.imshow(np.fliplr(data), cmap=cmap, norm=norm,
                  interpolation='none', origin='lower')                  

plt.scatter(135-LongLat_idx[1],LongLat_idx[0],color='r',s=100)
plt.annotate('Halemaumau Crater',xy=(135-LongLat_idx[1]-75,LongLat_idx[0]+10),\
color='red',bbox=dict(boxstyle='round,pad=0.5', fc='yellow', alpha=0.80))
                         
plt.xticks([0,25,50,75,100,135], 
           ['135','100','50','75','25','0'], rotation=0, fontsize=8 )
plt.yticks([0,50,100,150,203], 
           ['0','50','100','150','203'], 
           rotation=0, fontsize=8 )

title = 'Optical Depth Land And Ocean \n MYD04 C6 ({}-{:02d}-{:02d}; {:02d}h{:02d})'
plt.title(title.format( year, month, day, hh, mm), fontsize=10)

cbar_bounds = [0.0,0.1,0.2,0.3,0.4,0.5]
cbar_ticks = [0.0,0.1,0.2,0.3,0.4,0.5] 

cbar = plt.colorbar(img, cmap=cmap, norm=norm, boundaries=cbar_bounds, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks, fontsize=8)

plt.xlim(0.0,135)
plt.ylim(0.0,203)

plt.savefig("myd04_c6_Optical_Depth_Land_And_Ocean.png", 
            bbox_inches='tight', dpi=200)
plt.show()

Appliquer un filtre de Gabor à une image en fortran 90

					
program apply_gabor_filter

implicit none 

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

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

real(pr) :: pi = dacos(-1.0d0)

real(pr) :: theta,f,sigma_x,sigma_y
integer :: radius  
integer :: x1,x2
real(pr) :: G

real(pr), dimension(512,512) :: IMG,CIMG
real(pr), dimension(:,:), allocatable :: GM

!----------------------------------------------------------------------------------------!
! Read Image Data

open(1,file='lena_gray.txt')
do i = 1, 512
	read(1,*) (IMG(i,j),j=1,512)
end do
close(1)

!----------------------------------------------------------------------------------------!
! User Inputs

theta = 90.0 ! degrees

f = 0.1

sigma_x = 6
sigma_y = 6

radius = 20 

!----------------------------------------------------------------------------------------!
! Create Gabor filter 

theta = ( pi * theta  ) / 180.0 ! converting to radian

allocate( GM(2*radius+1,2*radius+1) )

k = 1
do x1 = - radius, radius
l = 1
do x2 = - radius, radius 
	call Gabor_Filter(x1,x2,theta,f,sigma_x,sigma_y,G)
	GM(k,l) = G
l = l + 1
end do
k = k + 1
end do

open(1,file='gabor_filter_01.txt')
do k = 1, 2*radius
	write(1,*) (GM(k,l),l=1,2*radius)
end do
close(1)

!----------------------------------------------------------------------------------------!
! Apply Gabor Filter

CIMG = 0

do i = 1, 512
do j = 1, 512
	k = 1
	do x1 = - radius, radius
	l = 1
	do x2 = - radius, radius 
		m = i - x1
		n = j - x2
		if( 0 < m .AND. m <= 512)then
		if( 0 < n .AND. n <= 512)then
		CIMG(i,j) = CIMG(i,j) + GM(k,l) * IMG(m,n)
		end if
		end if
	l = l + 1
	end do
	k = k + 1
	end do
end do
end do

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

open(1,file='image_gabor_filter_01.txt')
do i = 1, 512
	write(1,*) (CIMG(i,j),j=1,512)
end do
close(1)

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

end program apply_gabor_filter

!----------------------------------------------------------------------------------------!
! Gabor Filter Rotation 

subroutine Change_Base(x1,x2,theta,x_theta,y_theta)

implicit none 

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

integer :: x1
integer :: x2
real(pr) :: theta
real(pr) :: x_theta
real(pr) :: y_theta

x_theta = x1 * dcos(theta) + x2 * dsin(theta)
y_theta = x2 * dcos(theta) - x1 * dsin(theta)

end subroutine Change_Base

!----------------------------------------------------------------------------------------!
! Gabor Filter 

subroutine Gabor_Filter(x1,x2,theta,f,sigma_x,sigma_y,G)

implicit none 

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

integer :: x1
integer :: x2
real(pr) :: theta
real(pr) :: f
real(pr) :: sigma_x
real(pr) :: sigma_y

real(pr) :: x_theta
real(pr) :: y_theta

real(pr) :: r1,r2,arg

real(pr) :: pi = dacos(-1.0d0)
real(pr) :: G

call Change_Base(x1,x2,theta,x_theta,y_theta)

r1 = x_theta / sigma_x
r2 = y_theta / sigma_y

arg = - 0.5 * ( r1**2 + r2**2 )

G = dexp(arg) * dcos(2.0*pi*f*x_theta)

end subroutine Gabor_Filter

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

Plot scattering phase matrix with matplotlib

					
#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt

#----------------------------------------------------------------------------------------#

data  = np.loadtxt("hexagomal_pristine_crystal_GOM.txt", skiprows=18)

#----------------------------------------------------------------------------------------#

fig = plt.figure()

plt.plot(data[:,0],data[:,1],'k-')
plt.yscale('log')

plt.ylim(10**(-2.0),10**7)
plt.grid(True)

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_1.png', bbox_inches='tight')
plt.close()

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,2],'k-')

plt.ylim(-0.5,0.5)
plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$-P_{12}/P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_2.png', bbox_inches='tight')

plt.close()

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,3],'k-')

plt.ylim(-1.0,1.0)
plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$P_{22}/P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_3.png', bbox_inches='tight')

plt.close()

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,4],'k-')

plt.ylim(-1.0,1.0)
plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$P_{33}/P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_4.png', bbox_inches='tight')

plt.close()

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,5],'k-')

plt.ylim(-0.5,0.5)
plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$-P_{43}/P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_5.png', bbox_inches='tight')

plt.close()

#----------------------------------------------------------------------------------------#

plt.plot(data[:,0],data[:,6],'k-')

plt.ylim(0.0,1.0)
plt.grid(True,which="both")

plt.xlabel(r"Scattering Angle $\Theta$ ($^\circ$)")
plt.title(r"$P_{44}/P_{11}$")

plt.savefig('geometric_optic_method_hexagonal_pristine_crystals_6.png', bbox_inches='tight')

plt.close()

#----------------------------------------------------------------------------------------#

Geometric optic method (ray-tracing+diffraction) for hexagonal pristine crystals

					
!----------------------------------------------------------------------------------------!
! 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

Write a python script to find corresponding MODIS-CALIOP files over a month

					
#!/usr/bin/env python

import os
import glob
import numpy as np
import calendar

year = '2008'
month = '07'
caliop_version = 4

if int(year) % 400 == 0:
	day_month = [31,29,31,30,31,30,31,31,30,31,30,31]
elif int(year) % 4 == 0 and int(year) % 100 != 0:
	day_month = [31,29,31,30,31,30,31,31,30,31,30,31]
else:	
	day_month = [31,28,31,30,31,30,31,31,30,31,30,31]

'''---------- step 1 ----------'''	

myd021km_path = '/Volumes/Data/MYD021KM/'
myd03_path = '/Volumes/Data/Myd03/'
myd35_path = '/Volumes/DATA_03/modis/Myd35/'
myd06_c5_path = '/Volumes/DATA_03/modis/Myd06_c5/' 
myd06_c6_path = '/Volumes/DATA_03/modis/Myd06_c6/'
myd06_baseline_c6_path = '/Volumes/DATA_03/modis/Myd06_baseline_c6/'
myd06_research_c6_path = '/Volumes/DATA_03/modis/Myd06_research_c6/'

if caliop_version == 3:
	caliop_01km_path = '/Volumes/Data/Caliop/01km_CLay/'
	caliop_05km_path = '/Volumes/Data/Caliop/05km_CLay/'

if caliop_version == 4:
	caliop_01km_path = '/Volumes/Data/Caliop/v4/01km_CLay/'
	caliop_05km_path = '/Volumes/Data/Caliop/v4/05km_CLay/'

#dest_path = '/Users/marchant/Research/Colocated_data/Caliop_modis/Scripts/Caliop_modis/01km_CLay/'
dest_path = '/Volumes/Data/Colocated_data/Caliop_modis/01kmCLay/'

for i in np.arange(calendar.monthrange(int(year),int(month))[1]):
	'''---------- ----------'''	
	if i+1 < 10:
		day = '0' + str(i+1)
	else:
		day = str(i+1)
	'''---------- ----------'''	
	sum_days = 0
	for d in np.arange(int(month)-1):
		sum_days = sum_days + day_month[d]
	if int(month)-1 == 0:
		Julian_day = i+1 
	else:
		Julian_day = sum_days + i+1 	
	if Julian_day < 10:	
		Julian_day_str = '00' + str(Julian_day)
	elif Julian_day < 100:	
		Julian_day_str = '0' + str(Julian_day)
	else:
		Julian_day_str = str(Julian_day)
	'''---------- step 2 ----------'''
	myd021km_full_path = myd021km_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd03_full_path = myd03_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd35_full_path = myd35_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd06_c5_full_path = myd06_c5_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd06_c6_full_path = myd06_c6_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd06_baseline_c6_full_path = myd06_baseline_c6_path + year + '/' + year + '_' + month + '_' + day + '/'
	myd06_research_c6_full_path = myd06_research_c6_path + year + '/' + year + '_' + month + '_' + day + '/'
	caliop_01km_full_path = caliop_01km_path + year + '/' + year + '_' + month + '_' + day + '/'
	caliop_05km_full_path = caliop_05km_path + year + '/' + year + '_' + month + '_' + day + '/'
	'''---------- step 3 ----------'''	
	myd021km_day_list = glob.glob(myd021km_full_path + '/*hdf' )
	myd03_day_list = glob.glob(myd03_full_path + '/*hdf' )
	myd35_day_list = glob.glob(myd35_full_path + '/*hdf' )
	myd06_c5_day_list = glob.glob(myd06_c5_full_path + '/*hdf' )
	myd06_c6_day_list = glob.glob(myd06_c6_full_path + '/*hdf' )
	myd06_baseline_c6_day_list = glob.glob(myd06_baseline_c6_full_path + '/*hdf' )
	myd06_research_c6_day_list = glob.glob(myd06_research_c6_full_path + '/*hdf' )
	caliop_01km_day_list = glob.glob(caliop_01km_full_path + '/*hdf' )
	caliop_05km_day_list = glob.glob(caliop_05km_full_path + '/*hdf' )
	'''---------- step 4 ----------'''	
	for modis_granule_name in myd03_day_list:
		offset = len(myd03_full_path)
		myd03_name = modis_granule_name[offset:]
		#print 'offset', offset
		#print modis_granule_name[offset+15:offset+17]
		modis_granule_hour_str = modis_granule_name[offset+15:offset+17] 
		modis_granule_hour = int( modis_granule_hour_str )
		#print modis_granule_name[offset+17:offset+19]
		modis_granule_minute_str = modis_granule_name[offset+17:offset+19] 
		modis_granule_minute = int( modis_granule_minute_str )
		#print modis_granule_hour, modis_granule_minute
		missing_file = 0 
		'''---------- MYD021KM ----------'''	
		file_name_prefix = 'MYD021KM.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0
		for myd021km_granule_name in myd021km_day_list:
			offset = len(myd021km_full_path)
			if file_name_prefix == myd021km_granule_name[offset:offset+22]:
				found_file = 1
				myd021km_name = myd021km_granule_name[offset:]
				#print myd35_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1
			myd021km_name = 'NoFile'
		'''---------- MYD35 ----------'''	
		file_name_prefix = 'MYD35_L2.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0
		for myd35_granule_name in myd35_day_list:
			offset = len(myd35_full_path)
			if file_name_prefix == myd35_granule_name[offset:offset+22]:
				found_file = 1
				myd35_name = myd35_granule_name[offset:]
				#print myd35_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1
			myd35_name = 'NoFile'
		'''---------- MYD06_c5 ----------'''				
		file_name_prefix = 'MYD06_L2.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0		
		for myd06_c5_granule_name in myd06_c5_day_list:
			offset = len(myd06_c5_full_path)
			if file_name_prefix == myd06_c5_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c5_name = myd06_c5_granule_name[offset:]
				#print myd06_c5_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1
			#print "MISSING MYD06_c5", len(myd06_c5_day_list)
			#print myd06_c5_full_path
			myd06_c5_name = 'NoFile'		
		'''---------- MYD06_c6 ----------'''				
		file_name_prefix = 'MYD06_L2.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0		
		for myd06_c6_granule_name in myd06_c6_day_list:
			offset = len(myd06_c6_full_path)
			if file_name_prefix == myd06_c6_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_c6_name = myd06_c6_granule_name[offset:]
				#print myd06_c6_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1	
			myd06_c6_name = 'NoFile'		
		'''---------- MYD06_baseline_c6 ----------'''				
		file_name_prefix = 'MYD06_L2.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0		
		for myd06_baseline_c6_granule_name in myd06_baseline_c6_day_list:
			offset = len(myd06_baseline_c6_full_path)
			if file_name_prefix == myd06_baseline_c6_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_baseline_c6_name = myd06_baseline_c6_granule_name[offset:]
				#print myd06_baseline_c6_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1		
			myd06_baseline_c6_name = 'NoFile'		
		'''---------- MYD06_research_c6 ----------'''				
		file_name_prefix = 'MYD06_L2.A'+ year + Julian_day_str + '.' \
		                   + modis_granule_hour_str + modis_granule_minute_str
		found_file = 0		
		for myd06_research_c6_granule_name in myd06_research_c6_day_list:
			offset = len(myd06_research_c6_full_path)
			if file_name_prefix == myd06_research_c6_granule_name[offset:offset+22]:
				found_file = 1			
				myd06_research_c6_name = myd06_research_c6_granule_name[offset:]
				#print myd06_research_c6_granule_name[offset:]
		if found_file == 0:
			#missing_file = missing_file + 1			
		    myd06_research_c6_name = 'NoFile'
		'''---------- CAL 01km ----------'''					
		minimum = 1000
		find_one = 0
		for caliop_01km_granule_name in caliop_01km_day_list:
			offset = len(caliop_01km_full_path)
			
			if caliop_version == 3:
				caliop_granule_hour = int( caliop_01km_granule_name[offset+47:offset+49] )
				caliop_granule_minute = int( caliop_01km_granule_name[offset+50:offset+52]   )

			if caliop_version == 4:
				caliop_granule_hour = int( caliop_01km_granule_name[offset+46:offset+48] )
				caliop_granule_minute = int( caliop_01km_granule_name[offset+49:offset+51]   )

			#print caliop_granule_hour, caliop_granule_minute
			difference = modis_granule_hour * 60 + modis_granule_minute \
			             - ( caliop_granule_hour * 60 + caliop_granule_minute )
			if (difference < minimum ) and ( difference > 0 ):
				minimum = difference
				tmp_02 = caliop_01km_granule_name
				find_one = 1
		if find_one == 1:
			caliop_01km_name = tmp_02[offset:]
			#print tmp_02[offset:]
		else:
			missing_file = missing_file + 1	
			#print 'CAL 01 km missing'			
		'''---------- CAL 05km ----------'''				
		minimum = 1000
		find_one = 0
		for caliop_05km_granule_name in caliop_05km_day_list:
			offset = len(caliop_05km_full_path)
			
			if caliop_version == 3:
				caliop_granule_hour = int( caliop_05km_granule_name[offset+42:offset+44] )
				caliop_granule_minute = int( caliop_05km_granule_name[offset+45:offset+47] )

			if caliop_version == 4:
				caliop_granule_hour = int( caliop_05km_granule_name[offset+46:offset+48] )
				caliop_granule_minute = int( caliop_05km_granule_name[offset+49:offset+51] )

			#print caliop_05km_granule_name, caliop_granule_hour, caliop_granule_minute
			difference = modis_granule_hour * 60 + modis_granule_minute \
			             - ( caliop_granule_hour * 60 + caliop_granule_minute )
			if (difference < minimum ) and ( difference > 0 ):
				minimum = difference
				tmp_03 = caliop_05km_granule_name
				find_one = 1
		if find_one == 1:
			caliop_05km_name = tmp_03[offset:]		
			#print tmp_03[offset:]
		else:
			missing_file = missing_file + 1	
			#print 'CAL 05 km missing'				
		'''---------- step 7 ----------'''
		#print "missing_file", missing_file
		if ( missing_file == 0 ):
			#print 'NO MISSING FILES FOUND'
			filename = 'calmod_file_name_' + modis_granule_hour_str + modis_granule_minute_str
			f = open(filename, 'w')
			f.write(myd03_name + ' \n')	
			f.write(myd35_name + ' \n')				
			f.write(myd06_c5_name + ' \n')			
			f.write(myd06_c6_name + ' \n')			
			f.write(myd06_baseline_c6_name + ' \n')
			f.write(myd06_research_c6_name + ' \n')			
			f.write(caliop_01km_name + ' \n')
			f.write(caliop_05km_name + ' \n')
			f.write(myd021km_name + ' \n')	
			f.closed
			command = 'mv ' + filename + ' ' + dest_path + year + '/' + year + '_' + month + '_' + day + '/'
			os.system(command)			
		'''----------  ----------'''			
		#print 'modis granule done', modis_granule_hour, modis_granule_minute
	print day + ' done'