Python module to calculate scattering amplitude and phase matrices of rotationally symmetric particles based on the T-matrix method.
Based on the original Fortran T-matrix code by M.I. Mishchenko. (https://www.giss.nasa.gov/staff/mmishchenko/t_matrix.html). Includes orientation averaging schemes adapted from J. Leinonen (https://github.com/jleinonen/pytmatrix).
This is a pure Python reimplementation of the original routines using NumPy, not a Python wrapper around the Fortran code.
Currently spheroids, Chebyshev shapes and generalized chebyshev shapes are implemented. Also includes several empirical raindrop shapes.
Download and install with pip:
pip install scatterpy
Single generalized Chebyshev particle:
import numpy as np
from scatterpy import *
radius = np.asarray([5., 6.]) # equivalent sphere radius (meter)
wl = 3.1415925 # wavelength of incident light (same unit as radius)
mr = np.asarray([1.5+0.02j, 1.8+0.02j]) # refractive index
alpha = np.asarray([np.deg2rad(145)]) # Drop canting euler angle in rad (rotation around z axis).
beta = np.asarray([np.deg2rad(52)]) # # Drop canting euler angle in rad (rotation around x' axis).
theta0 = np.asarray([np.deg2rad(56), np.deg2rad(60)]) # incomming beam zenith angle
theta = np.asarray([np.deg2rad(65), np.deg2rad(70)]) # outgoing beam zenith angle
phi0 = np.asarray([np.deg2rad(114), np.deg2rad(120)]) # incoming beam azimuth angle
phi = np.asarray([np.deg2rad(128), np.deg2rad(130)]) # outgoing beam azimuth angle
sfunc = lambda x: shapes.gen_chebyshev(np.array([-0.0481, 0.0359, -0.1263, 0.0244,
0.0091, -0.0099, 0.0015, 0.0025,
-0.0016, -0.0002, 0.0010]))
T = calc_T(radius*2, wl, mr, rtol=0.001)
S = calc_S(T, theta0, theta, phi0, phi, alpha, beta)
Z = calc_Z(S)
print(S)Orientation averaged raindrops:
import numpy as np
from scatterpy import *
attrs = {'quantity':'diameter', 'unit':'m'}
D = np.arange(4e-3, 5e-3, 1e-4) # equivalent sphere diameter
wl = np.arange(8e-3, 9.1e-3, 1e-4)
wl = 8E-3 # wavelength of incident light (same unit as radius)
mr = np.asarray([1.5+0.02j]) # refractive index
sfunc = shapes.dropshape_CB90
wfunc = sph_gauss_pdf()
theta0 = np.deg2rad(56) # incomming beam zenith angle
theta = np.deg2rad(65) # outgoing beam zenith angle
phi0 = np.deg2rad(114) # incoming beam azimuth angle
phi = np.deg2rad(128) # outgoing beam azimuth angle
T = calc_T(D, wl, mr, sfunc=sfunc)
S, Z = calc_SZ_oa(T, theta0, theta, phi0, phi, wfunc)
S *= wl / (2 * np.pi)
print(S)To evaluate the empirical raindrop shapes:
D = np.asarray([2e-3, 4e-3])
shape = dropshape_CB90(D)
cheb = shape(np.polynomial.legendre.leggauss(100)[0])
print(cheb)