#!/usr/bin/env python

# Author: Benjamin Marchant
# Ref: Hong Du's paper: "Mie-Scattering Calculation, Applied Optics, Vol. 43, Issue 9, pp. 1951-1956 (2004)

import matplotlib.pyplot as plt

import numpy as np
import math

#----------------------------------------------------------------------------------------#
# user inputs

ri = complex(1.33,-0.0) # refractive index

x = 5.0 # Mie Parameter (2*pi*radius/wavelength)

n = 100   

angular_resolution = 1.0

#----------------------------------------------------------------------------------------#
# Calculations of the Riccati-Bessel functions by upward recurrence.

PHI = np.zeros((n+1))

a = math.cos(x) 
PHI[0] = math.sin(x) 
b = PHI[0]
for i in range(1,n+1):
	PHI[i] = ( ( 2 * i - 1 ) / x ) * b  - a
	a = b
	b = PHI[i]

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

CHI = np.zeros((n+1))

a = math.cos(x)
CHI[0] = a
CHI[1] = math.cos(x) / x + math.sin(x)
b = CHI[1]
for i in range(2,n+1):
	CHI[i] = ( ( 2 * i - 1 ) / x ) * b  - a
	a = b
	b = CHI[i]

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

ZETA = np.zeros((n+1), dtype=complex)
for i in range(0,n+1):
	ZETA[i] = complex(PHI[i], CHI[i])

#----------------------------------------------------------------------------------------#
# Calculations of the complex functions rn(mx) by upward recurrence.

RN = np.zeros((n+1), dtype=complex)	

a = ri.real * x
b = ri.imag * x	

c1 = math.tan(a) - math.exp(-2.0*b) * math.tan(a)
c2 = 1.0 + math.exp(-2.0*b)
c3 = -1.0 + math.exp(-2.0*b)
c4 = math.tan(a) + math.exp(-2.0*b) * math.tan(a) 

a = c1 * c3 + c2 * c4
b = c2 * c3 - c1 * c4
c = c3 * c3 + c4 * c4

RN[0] = complex(a/c, b/c) 
for i in range(1,n+1):
	a = ( 2.0 * i - 1.0 ) / ( ri * x ) - RN[i-1]
	RN[i] = 1.0 / a	

#----------------------------------------------------------------------------------------#
# Calculations of the Mie coefficients an et bn.

MA = np.zeros((n), dtype=complex)
MB = np.zeros((n), dtype=complex)

for i in range(0,n):
	a = RN[i+1] / ri + (i + 1) * ( 1.0 - 1.0 / ( ri * ri ) ) / x 
	MA[i] = ( a * PHI[i+1] - PHI[i]) / ( a * ZETA[i+1] - ZETA[i] )
	MB[i] = ( RN[i+1] * ri * PHI[i+1] - PHI[i] ) / ( RN[i+1] * ri * ZETA[i+1] - ZETA[i] ) 

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

def angular_functions(x,angle,n,MA,MB):

	#------------------------------------------------------------------------------------#
	# Calculations of the Mie angular functions.

	PI = np.zeros((n))

	a = math.cos(angle) 
	b = 0.0
	c = 1.0

	PI[0] = 1.0
	for j in range(1,n):
		s = a * c
		t = s - b 
		PI[j] = s + t + t / j  
		b = c
		c = PI[j]

	#------------------------------------------------------------------------------------#
	
	TAU = np.zeros((n))
	
	TAU[0] = a
	for j in range(1,n):
		t = a * PI[j] - PI[j-1]
		TAU[j] = ( j + 1 ) * t - PI[j-1]
		
	#------------------------------------------------------------------------------------#
	# scattering amplitudes

	b = complex(0.0,0.0)
	c = complex(0.0,0.0)
	for i in range(0,n):
		a = ( 2.0 * ( i + 1 ) + 1.0 ) / ( ( i + 1 ) * ( i + 2 ) )
		b = b + a * ( MA[i] * PI[i] + MB[i] * TAU[i] )
		c = c + a * ( MB[i] * PI[i] + MA[i] * TAU[i] )
	sa1 = b
	sa2 = c

	return sa1, sa2

#----------------------------------------------------------------------------------------#
# phase matrix

angle_list = []
p11_list = []
p12_list = []
p33_list = []
p34_list = []

for angle in np.arange(0.0,180.0,angular_resolution):

	#print angle
	
	sa1,sa2 = angular_functions(x,math.radians(angle),n,MA,MB)
	
	p11 = ( sa1 * sa1.conjugate() + sa2 * sa2.conjugate() )
	p12 = ( sa2 * sa2.conjugate() - sa1 * sa1.conjugate() ) 
	p33 = ( sa2 * sa1.conjugate() + sa2 * sa1.conjugate() ) 
	p34 = complex(0,1) * ( sa1 * sa2.conjugate() - sa2 * sa1.conjugate() ) 

	angle_list.append(angle)
	p11_list.append(p11.real)
	p12_list.append(p12.real)
	p33_list.append(p33.real)
	p34_list.append(p34.real)		

#----------------------------------------------------------------------------------------#
# normalization

norm = 0 

for i in range(0,n):
	a = abs(MA[i]) * abs(MA[i]) + abs(MB[i]) * abs(MB[i])
	norm = norm +  ( 2.0 * ( i + 1.0 ) + 1.0 ) * a 

p11_list = p11_list / norm 
p12_list = p12_list / norm 
p33_list = p33_list / norm 
p34_list = p34_list / norm 

#----------------------------------------------------------------------------------------#
# extinction efficiency.

a = 0.0
for i in range(0,n):
	a = a + ( 2 * ( i + 1 ) + 1 ) * ( MA[i] + MB[i] ).real 

qe = ( 2.0 / ( x * x ) ) * a

#----------------------------------------------------------------------------------------#
# asymmetry factor.

a = 0.0
for i in range(0,n-1):
	c1 = MA[i] * MA[i+1].conjugate() + MB[i] * MB[i+1].conjugate()
	c2 = MA[i] * MB[i].conjugate()
	b = ( (i+1) * ( (i+1) + 2.0 ) ) / ( (i+1) + 1.0 )
	c = ( 2.0 * (i+1) + 1.0 ) / ( (i+1) * ( (i+1) + 1.0 ) )
	a = a + b * c1.real + c * c2.real

gfact = ( 2.0 / norm ) * a

#----------------------------------------------------------------------------------------#
# single scattering albedo.

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

#----------------------------------------------------------------------------------------#
# Data Visualization

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p11_list,'k-')
plt.yscale('log')

plt.grid(True,which="both")

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

plt.text(0.70,0.85,'x = '+str(x), fontsize=12, transform = ax.transAxes)
plt.text(0.70,0.80,'n=('+str(ri.real)+','+str(ri.imag)+')', \
                    fontsize=12, transform = ax.transAxes)     

plt.text(0.70,0.70,'qe = '+str(round(qe,2)), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.65,'w0 = '+str(round(albedo,2)), fontsize=12, transform = ax.transAxes)     
plt.text(0.70,0.60,'g = '+str(round(gfact,2)), fontsize=12, transform = ax.transAxes)     

#plt.savefig('p11.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p12

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p12_list,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

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

#plt.savefig('p12.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p33

fig = plt.figure()
ax = fig.add_subplot(111)

plt.plot(angle_list,p33_list,'k-')
plt.yscale('log')

plt.grid(True,which="both")

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

#plt.savefig('p33.png')
plt.show()

#----------------------------------------------------------------------------------------#
# plot p34

fig = plt.figure()

ax = fig.add_subplot(111)

plt.plot(angle_list,p34_list,'k-')
#plt.yscale('log')

plt.grid(True,which="both")

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

#plt.savefig('p34.png')
plt.show()

Infos

Main author / Auteur principal:
Ben

Created / Créé: 2017-05-03

Edited / Modifié: 2017-05-03

License / Licence: MIT License

Abstract / Résumé