Lorentz-Mie scattering calculation using fortran

Open Science Notebook - Codes

Lorentz-Mie scattering calculation using fortran

added by Ben on May 31, 2016

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
you need to be logged in to post a comment

Abstract


Lorentz-Mie scattering code in fortran90, based on Hong Du's paper: "Mie-Scattering Calculation, Applied Optics, Vol. 43, Issue 9, pp. 1951-1956 (2004)", that I wrote when I was a PHD student

size:100

Licence


MIT License

A short and simple permissive license with conditions only requiring preservation of copyright and license notices. Licensed works, modifications, and larger works may be distributed under different terms and without source code.

Read more