program mie_code

! reference: Hong Du: Mie-Scattering Calculation, Applied Optics, 
!                     Vol. 43, Issue 9, pp. 1951-1956 (2004)

implicit none

integer :: i,j,k

integer, parameter :: n = 100  ! number of coefficients 

double precision :: x = 5.0  ! Mie Parameter (2*pi*r/wl)

double complex :: ri = (1.33d0,-0.0d0)  ! refractive index 

double complex, dimension(n) :: ma , mb  

double precision, dimension(:,:), allocatable :: pi , tau  

double precision, dimension(0:n) :: phi , chi 
double complex, dimension(0:n) :: zeta  

double complex, dimension(0:n) :: rn  

double complex, dimension(:), allocatable :: sa1 , sa2 

double precision, dimension(:), allocatable :: p11,p12,p33,p34  

double precision :: qe

double precision :: albedo , gfact

integer :: nb_angles

double precision :: teta
double precision :: teta_resol = 1.0 ! scattering angle resolution

double precision :: s , t

double precision :: temp1 , temp2 , temp3
double complex :: ctemp1 , ctemp2 ,ctemp3

double precision :: coef1 , coef2 , coef3 , coef4

double precision :: norm

double precision :: cdr = 0.01745329252

double complex :: ci = (0.0d0,1.0d0) 

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

if (n < int(x)) write(6,*) 'Warning: number of coeffs n is not enough !'

nb_angles = int(180/teta_resol) + 1

allocate( pi(nb_angles,n) )
allocate( tau(nb_angles,n) )

allocate( sa1(nb_angles) )
allocate( sa2(nb_angles) )

allocate( p11(nb_angles) )
allocate( p12(nb_angles) )
allocate( p33(nb_angles) )
allocate( p34(nb_angles) )

do i = 1 , n
	ma(i)  = (0.0d0,0.0d0)
	mb(i)  = (0.0d0,0.0d0)
end do

do i = 1 , nb_angles
do j = 1 , n
	pi(i,j)  = 0.0d0
	tau(i,j) = 0.0d0
end do
end do

do i = 1 ,nb_angles
	sa1(i) = 0.0d0
	sa2(i) = 0.0d0
end do

!----------------------------------------------------------------------------------------!
! Step 1: Calculations of the Mie angular functions.

teta = 0.0d0

do i = 1 , nb_angles

!---- pi function

temp1 = dcos( teta * cdr )

temp2 = 0.0d0
temp3 = 1.0d0
pi(i,1) = temp3
do j = 2 , n
	s = temp1 * temp3
	t = s - temp2
	pi(i,j) = s + t + t / ( j - 1 )
	temp2 = temp3
	temp3 = pi(i,j)
end do

!----- tau function

tau(i,1) = temp1
do j = 2 , n
	t = temp1 * pi(i,j) - pi(i,j-1)
	tau(i,j) = j * t - pi(i,j-1)
end do

teta = teta + teta_resol
end do

!----------------------------------------------------------------------------------------!
! Step 2: Calculations of the Riccati-Bessel functions by upward recurrence.

!----- phi function

temp1 = dcos( x )
phi(0) = dsin( x )
temp2 = phi(0) 
do i = 1 , n
	phi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
	temp1 = temp2
	temp2 = phi(i)
end do

!----- chi function

chi(0) = dcos(x)
temp1  = chi(0) 
chi(1) = dcos(x) / x + dsin(x)
temp2  = chi(1) 
do i = 2 , n
	chi(i) = ( ( 2 * i - 1 ) / x ) * temp2  - temp1
	temp1 = temp2
	temp2 = chi(i)
end do

!----- zeta fcuntion

do i = 0 , n
	zeta(i) = phi(i) + ci * chi(i) 
end do

!----------------------------------------------------------------------------------------!
! Step 3: Calculations of the complex functions rn(mx) by upward recurrence.

temp1 = dble( ri ) * x
temp2 = aimag( ri ) * x

coef1 = dtan( temp1 ) - dexp( - 2.0 * temp2 ) * dtan( temp1 )
coef2 = 1.0d0 + dexp( - 2.0 * temp2 )
coef3 = - 1.0d0 + dexp( - 2.0 * temp2 )
coef4 = dtan( temp1 ) + dexp( - 2.0 * temp2 ) * dtan( temp1 )

temp1 = coef1 * coef3 + coef2 * coef4
temp2 = coef2 * coef3 - coef1 * coef4
temp3 = coef3 * coef3 + coef4 * coef4

rn(0) = ( temp1 / temp3 ) + ci * ( temp2 / temp3 )

do i = 1 , n
	ctemp1 = ( 2 * i - 1 ) / ( ri * x ) - rn(i-1)
	rn(i) = 1.0d0 / ctemp1
end do

!----------------------------------------------------------------------------------------!
! Step 4: Calculations of the Mie coefficients an et bn.

do i = 1 , n
	ctemp1 = rn(i) / ri + i * ( 1.0d0 - 1.0d0 / ( ri * ri ) ) / x
	ma(i) = ( ctemp1 * phi(i) - phi(i-1) ) / ( ctemp1 * zeta(i) - zeta(i-1) )
	mb(i) = ( rn(i) * ri  * phi(i) - phi(i-1) ) / ( rn(i) * ri  * zeta(i) - zeta(i-1) )
end do

!----------------------------------------------------------------------------------------!
! Step 5: scattering amplitudes

do i = 1 , nb_angles
	ctemp2 = (0.0d0,0.0d0)
	ctemp3 = (0.0d0,0.0d0)
	do j = 1 , n
		temp1 = DBLE( 2 * j + 1 ) / DBLE( j * ( j + 1 ) )
		ctemp2 = ctemp2 + temp1 * ( ma(j) * pi(i,j) + mb(j) * tau(i,j) )
		ctemp3 = ctemp3 + temp1 * ( mb(j) * pi(i,j) + ma(j) * tau(i,j) )
	end do
	sa1(i) = ctemp2
	sa2(i) = ctemp3
end do

!----------------------------------------------------------------------------------------!
! Step 6: phase function.

!----- norm (need to be improved <- integration too simple)

norm = 0.0
do i = 1 , n
	temp1 = abs( ma(i) ) * abs( ma(i) ) + abs( mb(i) ) * abs( mb(i) )
	norm = norm + ( 2.0d0 * i + 1.0d0 ) * temp1
end do

do i = 1 , nb_angles
	p11(i) = ( sa1(i) * conjg( sa1(i) ) + sa2(i) * conjg( sa2(i) ) ) / norm
	p12(i) = ( sa2(i) * conjg( sa2(i) ) - sa1(i) * conjg( sa1(i) ) ) / norm 
	p33(i) = ( sa2(i) * conjg( sa1(i) ) + sa1(i) * conjg( sa2(i) ) ) / norm 
	p34(i) = ci * ( sa1(i) * conjg( sa2(i) ) - sa2(i) * conjg( sa1(i) ) ) / norm 
end do

!----------------------------------------------------------------------------------------!
! extinction efficiency.

temp1 = 0.0d0
do i = 1 , n
	temp1 = temp1 + ( 2 * i + 1 ) * dble( ma(i) + mb(i) ) 
end do

qe = ( 2.0d0 / ( x * x ) ) * temp1

!----------------------------------------------------------------------------------------!
! asymmetry factor.

temp1 = 0.0d0
do i = 1 , n - 1
	ctemp1 = ma(i) * conjg( ma(i+1) ) + mb(i) * conjg( mb(i+1) ) 
	ctemp2 = ma(i) * conjg( mb(i) )
	temp2 = dble( i * ( i + 2.0d0 ) ) / dble( i + 1.0d0 )
	temp3 = dble( 2.0d0 * i + 1.0d0 ) / dble( i * ( i + 1.0d0 ) )
	temp1 = temp1 + ( temp2 * dble( ctemp1 ) + temp3 * dble( ctemp2 ) )
end do

gfact = ( 2.0d0 / norm ) * temp1

!----------------------------------------------------------------------------------------!
! single scattering albedo.

albedo = ( 2.0d0 / ( x * x ) ) * norm / qe

!----------------------------------------------------------------------------------------!
! save data in file

teta = 0.0d0
open(1,file='mie_code_output.txt')
	write(1,*) x, real(real(ri)), real(AIMAG(ri)) !, wl
	write(1,*) qe,albedo,gfact
	do i = 1 , nb_angles
		write(1,*) teta,p11(i),p12(i),p33(i),p34(i)
		teta = teta + teta_resol
	end do
close(1)

end program mie_code

About / Au sujet de

Ben
Created:
31 mai 2016 21:07:57

Updated:
31 mai 2016 21:07:57

License / Licence

MIT License

Abstract / Résumé