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

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()

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()

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()

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()

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