-
Notifications
You must be signed in to change notification settings - Fork 44
Home
- Introduction
- Installation
- Testing
-
Usage
- Minimal example
-
Using the
Scattererclass - Orientation averaging: the
orientationmodule - Particle size distribution integration: the
psdmodule - The
tmatrix_auxmodule - Computation of scattering parameters: the
scattermodule - Microwave-specific functionality: the
radarmodule - Refractive indices: the
refractivemodule
- Examples
- Code documentation
- Current limitations and known problems
- Citing the code
- Original code
- References
- See also
This is a Python package for computing the electromagnetic scattering properties of nonspherical particles using the T-matrix method.
The package can be used particles of any size as long as the wavelength is not too short compared to the particle size. However, there is some functionality that is intended specifically for the calculation of scattering properties of hydrometeors at microwave frequencies, such as those used in weather radars.
Make sure you have NumPy, SciPy and a Python environment (version >=2.6) installed. You also need a Fortran 77 compiler (e.g. gfortran) to install the package. For example, if you are on a system with apt, this should get you what you need:
$ sudo apt-get install gfortran python-numpy python-scipy
The code is easiest to install from PyPI. If you have pip installed, install simply with
$ pip install pytmatrix
(using sudo if installing globally).
If you can't use pip, install manually:
- Download the latest release package and unzip
- Go to the root directory of the package and run (with
sudoif you need the permissions):
$ python setup.py install
You can also pull the latest development version (which may or may not be stable) from the repository:
git clone https://github.com/jleinonen/pytmatrix.git
This is to verify that you have a working package after installation.
Open a Python shell. Run the following:
from pytmatrix.test import test_tmatrix
test_tmatrix.run_tests()An output like this should be printed:
test_adaptive_orient (pytmatrix.test.test_tmatrix.TMatrixTests)
Test an adaptive orientation averaging case ... ok
test_against_mie (pytmatrix.test.test_tmatrix.TMatrixTests)
Test scattering parameters against Mie results ... ok
test_asymmetry (pytmatrix.test.test_tmatrix.TMatrixTests)
Test calculation of the asymmetry parameter ... ok
test_fixed_orient (pytmatrix.test.test_tmatrix.TMatrixTests)
Test a fixed-point orientation averaging case ... ok
test_integrated_x_sca (pytmatrix.test.test_tmatrix.TMatrixTests)
Test Rayleigh scattering cross section integrated over sizes. ... ok
test_optical_theorem (pytmatrix.test.test_tmatrix.TMatrixTests)
Optical theorem: test that for a lossless particle, Csca=Cext ... ok
test_psd (pytmatrix.test.test_tmatrix.TMatrixTests)
Test a case that integrates over a particle size distribution ... ok
test_radar (pytmatrix.test.test_tmatrix.TMatrixTests)
Test that the radar properties are computed correctly ... ok
test_rayleigh (pytmatrix.test.test_tmatrix.TMatrixTests)
Test match with Rayleigh scattering for small spheres ... ok
test_single (pytmatrix.test.test_tmatrix.TMatrixTests)
Test a single-orientation case ... ok
----------------------------------------------------------------------
Ran 10 tests in 16.940s
OK
The running time may vary; the most important thing is that "ok" is printed after each test. If some of the tests fail, please notify the author.
Please note that these tests only verify the basic functionality of the package, they do not (as of the latest version) test for convergence in edge cases.
This calculates the amplitude scattering matrix S for a sphere with a volume-equivalent radius of 1.0 and horizontal-to-vertical axis ratio of 0.6 and refractive index 1.5+0.5i, at a wavelength of 10.0:
>>> from pytmatrix import tmatrix
>>> scatterer = tmatrix.Scatterer(radius=1.0, wavelength=10.0, m=complex(1.5, 0.5), axis_ratio=1.0/0.6)
>>> scatterer.get_S()
array([[ 9.66435719e-02 +6.79614074e-02j,
6.16862803e-25 +7.07135826e-25j],
[ -6.16862447e-25 -7.07128493e-25j,
-1.01600111e-01 -1.06748868e-01j]])The Scatterer class, found in the pytmatrix.tmatrix module, is an object-oriented interface to the Fortran 77 T-matrix code. To use it, create an instance of the class, set the properties of the scatterer, and then run one of the commands to retrieve the amplitude and/or phase matrices.
The properties of the scatterer, the radiation and the computations can be set using the attributes of the TMatrix instance. You can also set any of them as keyword arguments to the constructor, as shown in the example above. The available properties are:
| Attribute | Description | Default value |
|---|---|---|
radius |
The radius of the equivalent sphere (see rat). | 1.0 |
radius_type |
If this is set to Scatterer.RADIUS_EQUAL_VOLUME, radius is the equivalent-volume radius. If set to Scatterer.RADIUS_EQUAL_AREA, radius is the equivalent-surface-area radius. If set to Scatterer.RADIUS_MAXIMUM, radius is the maximum radius. |
Scatterer.RADIUS_EQUAL_VOLUME |
wavelength |
The wavelength of the incident radiation. | 1.0 |
m |
The complex refractive index of the scatterer. | complex(2,0) |
axis_ratio |
The horizontal-to-rotational axis ratio. | 1.000001 |
shape |
The type of the particles. Use shape=Scatterer.SHAPE_SPHEROID for spheroids (default) and shape=Scatterer.SHAPE_CYLINDER for cylinders. Note that Chebyshev particles are not currently supported. |
Scatterer.SHAPE_SPHEROID |
ddelt |
The accuracy of the computations. | 1e-3 |
ndgs |
Number of division points used to integrate over the particle surface. Try increasing this if the computations do not converge. | 2 |
alpha |
The Euler angle α (alpha) of the scatterer orientation. | 0.0 |
beta |
The Euler angle β (beta) of the scatterer orientation. | 0.0 |
thet0 |
The zenith angle of the incident beam. | 90.0 |
thet |
The zenith angle of the scattered beam. | 90.0 |
phi0 |
The azimuth angle of the incident beam. | 0.0 |
phi |
The azimuth angle of the scattered beam. | 180.0 |
Kw_sqr |
The reference value of the water dielectric factor |K|² (used for calculating radar reflectivity). | 0.93 |
orient |
The method of performing orientation averaging. | orientation.orient_single |
or_pdf |
Probability distribution function for the particle orientation. | orientation.gaussian_pdf() |
n_alpha |
Number of integration points for averaging over the alpha angle when fixed-point orientation averaging is used. | 5 |
n_beta |
As above, but for the beta angle. |
10 |
psd_integrator |
Set this to a PSDIntegrator instance to enable size distribution integration. If this is None, size distribution integration is not used. See the PSD integration documentation for more information. |
None |
psd |
Set to a callable object giving the PSD value for a given diameter (for example a GammaPSD instance); default None. Has no effect if psd_integrator is None. |
None |
Note on units: radius and wavelength should always have the same unit. Some components of the radar module expect them to be given in millimeters, but otherwise, other units can be used. The results of the scattering computations will be in corresponding units. Angles should always be given in degrees (not radians).
The names of the parameters are defined such as to be compatible with the original T-Matrix code.
The Scatterer class encapsulates the computation of the amplitude (S) and phase (Z) matrices of the scatterer. There are three instance methods that can be used to access the results, and an additional convenience function for setting the geometry.
| Method | Arguments | Description |
|---|---|---|
get_SZ() |
None | Returns the amplitude and phase matrices for the currently defined scatterer, packed in a tuple. |
get_S() |
None | As get_SZ, but only returns the amplitude matrix. |
get_Z() |
None | As get_SZ, but only returns the phase matrix. |
set_geometry(geom) |
geom: A tuple of (thet0, thet, phi0, phi, alpha, beta). |
A convenience function for setting all angles of the scattering geometry at the same time (see attributes). |
One of the advantages of the T-matrix method is that once the T matrix is computed, it is relatively inexpensive to compute the amplitude and phase matrices for any scattering geometry. The Scatterer class implements a lazy evaluation scheme that exploits this property. Calling any of get_SZ, get_S or get_Z will reuse the results of previous calculations if they are still applicable. Changing any of the following parameters will trigger a recalculation of the T matrix : radius, rat, wavelength, m, axis_ratio, np, ddelt, ndgs. The recalculation is not done immediately, but instead on the next call to one of the above-mentioned methods. Changing any of the other attributes will only recalculate the amplitude and phase matrices. If no parameters are changed between calls to those methods, the amplitude and phase matrices will also be reused.
This module contains functions that allow the use of different methods for computing the scattering properties. Most importantly, they allow the use of orientation averaging when computing the scattering properties.
To enable orientation averaging, two things are needed. Firstly, you need to set the orient property of your Scatterer object to a function that is able to perform averaging. The function orientation.orient_single, the default, only computes for a single orientation. For orientation averaging, set orient to either orientation.orient_averaged_fixed or orientation.orient_averaged_adaptive. The former uses a weighted fixed-point quadrature integration for averaging, while the latter, much slower, employs an adaptive integration method. Both orientation averaging methods fully utilize the possibility to reuse the T matrix. If you enable orientation averaging, the Scatterer attributes alpha and beta are ignored.
Secondly, you should set the orientation distribution function. You can either use a function from the orientation module, or define your own. Currently, two different distributions are implemented: a uniform distribution returned by orientation.uniform_pdf() and a Gaussian distribution in orientation.gaussian_pdf(). The latter takes one argument, which specifies the standard deviation of the angle with respect to vertical orientation (this is often called the canting angle).
If you define your own orientation distribution, it should be represented by a callable object (e.g. a function) which takes one argument, the deviation from vertical orientation in degrees, and returns the orientation probability for that angle. Take a look at the functions in the orientation module for an example. The distribution function must be defined in the interval [0,180] and must integrate to 1 over that interval. These criteria are not automatically checked, so you are responsible for the correct behavior of your function.
To enable integration over particle size distributions (PSD), you should set the psd_integrator and psd attributes of the Scatterer instance. When PSDs are used, it is common to use the same particles for several different PSDs. For this reason, the PSDIntegrator class computes and stores a lookup table of the amplitude and phase matrices over a particle size interval; the PSD-weighted integration is then performed over this table.
The psd_integrator attribute should be set to an instance of the PSDIntegrator class from the psd module. The following attributes can be set:
| Attribute | Description | Default value |
|---|---|---|
num_points |
The number of different (equally spaced) particle diameters at which to store the amplitude and phase matrices. | 1024 |
D_max |
The maximum diameter for which to store the amplitude and phase matrices. | psd.Dmax |
m_func |
The refractive index as a function of size (see below). | None |
axis_ratio_func |
The horizontal-to-rotational axis ratio as a function of size (see below). | None |
geometries |
A tuple of the scattering geometries for which to store the amplitude and phase matrices. The format of a geometry specification is the same as that for the set_geometry function. The geometry of the Scatterer instance must be set to one of these geometries. |
(tmatrix_aux.geom_horiz_back,) |
The m_func and axis_ratio_func attributes specify the refractive index and the axis ratio as a function of size. These should be callables that take the diameter as an argument and return the corresponding variable. They must be defined in the interval (0,D_max]. If set to None, the constant values of m and axis_ratio are used, respectively.
The psd attribute specifies the particle size distribution. Currently, several PSDs are pre-defined in the psd module. These are summarized below; please refer to the docstrings of the classes for details.
- The normalized gamma PSD, as given by Eq. (7.62) of Bringi and Chandrasekar (2001), and a binned PSD for computing results for experimental data. Can be constructed as an instance of the
psd.GammaPSDclass, which takes three keyword arguments:D0,Nwandmu. - Exponential PSD
psd.ExponentialPSDwith argumentsLambdaandmu. - Unnormalized gamma PSD
psd.UnnormalizedGammaPSDwith argumentsLambda,N0andmu. - The
psd.BinnedPSDtakes as its arguments the bin edges and the bin values.
If you define your own PSD, it should be a callable object with two properties. Firstly, it should be called with a single argument that specifies the diameter, and it should return the value of the PSD for that diameter. The PSD must be defined in the interval (0, D_max) and should integrate to the number concentration of particles in a unit volume. Secondly, the PSD may optionally define a way to determine if two PSDs are equal. This avoids reintegrating over the PSD if it has not changed.
To use the PSDIntegrator class, set the above attributes as needed, then call the init_scatter_table method of the PSDIntegrator instance using your Scatterer instance as the argument. After this, the PSD-integrated amplitude and phase matrices can be retrieved as usual from the Scatterer class.
Because computing the lookup tables can sometimes take a long time, it is possible to save the lookup tables to a file. Call the save_scatter_table method to do this and load_scatter_table to retrieve the saved tables. After a call to load_scatter_table, you can use the PSDIntegrator object just as if you had called init_scatter_table.
Starting from version 0.3, the PSD integration also supports integrating some quantities that are obtained by integrating over all scattering angles: the scattering cross section, the extinction cross section and the asymmetry parameter. This is useful, for example, for a single scattering model in radiative transfer applications. You can choose to enable the PSD integration of these quantities by setting angular_integration to True when calling init_scatter table().
The methods of PSDIntegrator instances are summarized below:
| Method | Arguments | Description |
|---|---|---|
init_scatter_table(tm, angular_integration=False, verbose=False) |
tm: A scatterer instance. angular_integration: If True, compute the angular-integrated quantities (see above). verbose: Print information of the progress of the computation. |
Initialize the matrix lookup tables. |
save_scatter_table(fn, desription="") |
fn: the name and path of the file to save todescription: A description the table. |
Save the matrix lookup tables to a file. |
load_scatter_table(fn) |
fn: the name and path of the file to load from. |
Load the matrix lookup tables from a file. |
This module has some convenient definitions of constants that may make the code easier to use. These are listed below:
| Constant | Description |
|---|---|
VERSION |
The current version of pytmatrix
|
wl_S, wl_C, wl_X, wl_Ku, wl_Ka, wl_W
|
Typical wavelengths (millimeters) for various radar bands. |
K_w_sqr |
The reference dielectric factors |K|² for calculating the radar reflectivity. This is a dictionary with keys for each of the above mentioned wavelengths. |
geom_horiz_back |
A geometry tuple in the form accepted by TMatrix.set_geometry. Defines the geometry for backscattering of horizontal incident radiation. |
geom_horiz_forw |
As above, but for forward scattering of horizontal incident radiation. |
geom_vert_back |
As above, but for backscattering of vertical incident radiation. |
geom_vert_forw |
As above, but for forward scattering of vertical incident radiation. |
The following raindrop shape relationship functions in tmatrix_aux may also be useful for atmospheric applications: dsr_thurai_2007, dsr_pb, dsr_bc. See the docstrings for each function for a more detailed description.
This module contains functions that are used to compute scattering quantities of general interest in many applications. These all take any Scatterer object as their first argument and compute the results for the selected scatterer parameters. Some functions operate only on specific geometries, which are also described below.
| Method | Arguments | Description | Geometry |
|---|---|---|---|
sca_intensity(scatterer, h_pol=True) |
scatterer: a Scatterer objecth_pol: True for horizontal polarization of incident radiation, False for vertical polarization |
Scattering intensity (phase function) to the direction specified by the thet and phi arguments of scatterer. |
Any |
sca_xsect(scatterer, h_pol=True) |
scatterer: As aboveh_pol: As above. |
Scattering cross section | Any (thet and phi parameters ignored). |
ext_xsect(scatterer, h_pol=True) |
scatterer: As aboveh_pol: As above. |
Extinction cross section | Forward scattering (must be set as this is computed using the optical theorem). |
ssa(scatterer, h_pol=True) |
scatterer: As aboveh_pol: As above. |
Single scattering albedo, or sca_xsect/ext_xsect. |
Forward scattering (must be set as the extinction cross section is computed using the optical theorem). |
asym(scatterer, h_pol=True) |
scatterer: As aboveh_pol: As above. |
Asymmetry parameter <cos Θ> | Any (thet and phi parameters ignored). |
ldr(scatterer, h_pol=True) |
scatterer: As aboveh_pol: If True, return LDR_h. If False, return LDR_v. |
Linear depolarization ratio | Backscattering |
This module contains functions that are mostly interesting for use with radars and other microwave instruments. These all take any Scatterer object as their first argument, although many of them only make sense if the Scatterer object has PSD integration defined by enabling the psd_integrator attribute. Many of the functions also make sense only for a specific scattering geometry; you are responsible for setting an appropriate geometry yourself. Most of these functions expect that you give all sizes in millimeters.
| Method | Arguments | Description | Geometry |
|---|---|---|---|
radar_xsect(scatterer, h_pol=True) |
scatterer: a Scatterer objecth_pol: True for horizontal polarization, False for vertical polarization |
Radar cross section | Any (use non-backscattering geometries for bistatic cross section) |
refl(scatterer, h_pol=True) |
scatterer: As aboveh_pol: As above |
Radar reflectivity | Backscattering |
ldr(scatterer, h_pol=True) |
scatterer: As aboveh_pol: If True, return LDR_h. If False, return LDR_v. |
Linear depolarization ratio | Backscattering |
Zdr(scatterer) |
scatterer: As above |
Differential reflectivity | Backscattering |
delta_hv(scatterer) |
scatterer: As above |
Scattering differential phase | Backscattering |
rho_hv(scatterer) |
scatterer: As above |
Co-polar correlation | Backscattering |
Kdp(scatterer) |
scatterer: As above |
Specific differential phase (deg/km) | Forward scattering |
Ai(scatterer, h_pol=True) |
scatterer: As aboveh_pol: True for horizontal polarization, False for vertical polarization |
Specific attenuation (dB/km) | Forward scattering |
This module has various functions to aid with refractive indices, especially in the context of microwave measurements.
Two methods are available for calculating effective medium approximations (EMA):
| Method | Arguments | Description |
|---|---|---|
mg_refractive(m, mix) |
m: A tuple giving the refractive indices of the componentsmix: A tuple with the volume fractions of the components (if their sum is not 1, they are normalized with the sum). This should be the same size as m. |
The Maxwell-Garnett effective medium approximation. The arguments m and mix should both have at least two items. If len(m)==2, the first element is taken as the matrix and the second as the inclusion. If len(m)>2, the media are mixed recursively so that the last element is used as the inclusion and the second to last as the matrix, then this mixture is used as the last element on the next iteration, and so on. The effective complex refractive index is then returned. |
mg_bruggeman(m, mix) |
m: As abovemix: As above. |
The Bruggeman effective medium approximation. This works as above, but only supports two components. |
There are also definitions of the microwave complex refractive index of water for various wavelengths:
| Constant | Description |
|---|---|
m_w_0C |
The complex refractive indices of water in the microwave at 0 °C temperature. This is a dictionary with keys for each of the wavelengths mentioned in the tmatrix_aux description. |
m_w_10C |
As above, but for 10 °C temperature. |
m_w_20C |
As above, but for 20 °C temperature. |
See also the more detailed example on computing the the specific differential phase Kdp.
The following examples are adopted from the test.test_tmatrix module. They assume that the following imports have been made:
from pytmatrix.tmatrix import Scatterer
from pytmatrix.psd import PSDIntegrator, GammaPSD
from pytmatrix import orientation, radar, tmatrix_aux, refractiveA single-orientation case:
>>> scatterer = Scatterer(radius=2.0, wavelength=6.5, m=complex(1.5,0.5), axis_ratio=1.0/0.6)
>>> scatterer.get_SZ()
(array([[ 3.89338755e-02 -2.43467777e-01j,
-1.11474042e-24 -3.75103868e-24j],
[ 1.11461702e-24 +3.75030914e-24j,
-8.38637654e-02 +3.10409912e-01j]]),
array([[ 8.20899248e-02, -2.12975199e-02, -1.94051304e-24,
2.43057373e-25],
[ -2.12975199e-02, 8.20899248e-02, 2.00801268e-25,
-1.07794906e-24],
[ 1.94055633e-24, -2.01190190e-25, -7.88399525e-02,
8.33266362e-03],
[ 2.43215306e-25, -1.07799010e-24, -8.33266362e-03,
-7.88399525e-02]]))An orientation-averaged case with a 20° standard deviation Gaussian distribution:
>>> scatterer = Scatterer(radius=2.0, wavelength=6.5, m=complex(1.5,0.5), axis_ratio=1.0/0.6)
>>> scatterer.or_pdf = orientation.gaussian_pdf(std=20.0)
>>> scatterer.orient = orientation.orient_averaged_adaptive
>>> scatterer.get_S()
array([[ 6.49005717e-02 -2.42488000e-01j,
-6.13348029e-16 -4.11094415e-15j],
[ -1.50045335e-14 -1.63765222e-15j,
-9.54176591e-02 +2.84758322e-01j]])The same, but using fixed-point averaging instead (note the slightly different results):
>>> scatterer = Scatterer(radius=2.0, wavelength=6.5, m=complex(1.5,0.5), axis_ratio=1.0/0.6)
>>> scatterer.or_pdf = orientation.gaussian_pdf(20.0)
>>> scatterer.orient = orientation.orient_averaged_fixed
>>> scatterer.get_S()
array([[ 6.49006144e-02 -2.42487917e-01j,
1.20257280e-11 -5.23021952e-11j],
[ 6.21754954e-12 +2.95662694e-11j,
-9.54177106e-02 +2.84758152e-01j]])A simple example of using particle size distribution integration:
>>> scatterer = Scatterer(wavelength=6.5, m=complex(1.5,0.5), axis_ratio=1.0/0.6)
>>> scatterer.psd_integrator = PSDIntegrator()
>>> scatterer.psd = GammaPSD(D0=1.0, Nw=1e3, mu=4)
>>> scatterer.psd_integrator.D_max = 10.0
>>> scatterer.psd_integrator.init_scatter_table(scatterer)
>>> scatterer.get_Z()
array([[ 7.20539942e-02, -1.54020511e-02, -9.96222004e-25,
8.34245116e-26],
[ -1.54020511e-02, 7.20539942e-02, 1.23279117e-25,
1.40047875e-25],
[ 9.96224481e-25, -1.23290932e-25, -6.89738802e-02,
1.38873117e-02],
[ 8.34136779e-26, 1.40047656e-25, -1.38873117e-02,
-6.89738802e-02]])A somewhat more complex example of computing radar scattering properties for a C-band radar, showing more features of the size distribution module as well as radar-specific functions:
# For testing variable aspect ratio
# This is from Thurai et al., J. Atmos. Ocean Tech., 24, 2007
def drop_ar(D_eq):
if D_eq < 0.7:
return 1.0;
elif D_eq < 1.5:
return 1.173 - 0.5165*D_eq + 0.4698*D_eq**2 - 0.1317*D_eq**3 - \
8.5e-3*D_eq**4
else:
return 1.065 - 6.25e-2*D_eq - 3.99e-3*D_eq**2 + 7.66e-4*D_eq**3 - \
4.095e-5*D_eq**4
>>> scatterer = Scatterer(wavelength=tmatrix_aux.wl_C, m=refractive.m_w_10C[tmatrix_aux.wl_C])
>>> scatterer.psd_integrator = PSDIntegrator()
>>> scatterer.psd_integrator.axis_ratio_func = lambda D: 1.0/drop_ar(D)
>>> scatterer.psd_integrator.D_max = 10.0
>>> scatterer.psd_integrator.geometries = (tmatrix_aux.geom_horiz_back, tmatrix_aux.geom_horiz_forw)
>>> scatterer.or_pdf = orientation.gaussian_pdf(20.0)
>>> scatterer.orient = orientation.orient_averaged_fixed
>>> scatterer.psd_integrator.init_scatter_table(scatterer)
>>> scatterer.psd = GammaPSD(D0=2.0, Nw=1e3, mu=4)
>>> radar.refl(scatterer)
6382.9718136764586
>>> radar.refl(scatterer, False)
5066.7211159239041
>>> radar.Zdr(scatterer)
1.2598450941838721
>>> radar.ldr(scatterer)
0.0021935578301794474
>>> radar.rho_hv(scatterer)
0.0021935578301794474
>>> scatterer.set_geometry(tmatrix_aux.geom_horiz_forw)
>>> radar.Kdp(scatterer)
0.1933072695008976
>>> radar.Ai(scatterer)
0.01892301354749585The public interface of the code is documented with Python docstrings. If you have IPython, you can type a question mark after an object to get the docstring. For example:
In [1]: from pytmatrix import radar
In [2]: radar.Kdp?
Type: function
Base Class: <type 'function'>
String Form:<function Kdp at 0x40017d0>
Namespace: Interactive
Definition: radar.Kdp(tm)
Docstring:
Specific differential phase (K_dp) for the current setup.
Args:
tm: a TMatrix instance.
Returns:
K_dp [deg/km].
NOTE: This only returns the correct value if the particle diameter and
wavelength are given in [mm]. The tm object should be set to forward
scattering geometry before calling this function.- The main problem with
pytmatrixat the moment is that if an error occurs in the Fortran T-matrix code, this will cause the calculation to exit, also exiting the Python interpreter at the same time. This can happen, e.g., with invalid input values or if the T-matrix code fails to converge. While this does not affect the correctness of the results, it is an inconvenience that one should be aware of. - Chebyshev particles are currently not supported as particle shapes (cylinders and spheroids should work normally, but most of the testing has been done with spheroids).
If you use PyTMatrix in a published work, please cite the following paper:
Leinonen, J., High-level interface to T-matrix scattering calculations: architecture, capabilities and limitations, Opt. Express, vol. 22, issue 2, 1655-1660 (2014), doi: 10.1364/OE.22.001655.
If you want to cite the code directly, the following reference is additionally suggested:
Leinonen, J., Python code for T-matrix scattering calculations. Available at https://github.com/jleinonen/pytmatrix/.
Relevant references should be also made to the Fortran T-matrix code that forms the core of this package. The following should be cited, as appropriate:
- Mishchenko, M. I. and L. D. Travis, T-matrix computations of light scattering by large spheroidal particles, Opt. Commun., vol. 109, 16-21 (1994).
- Mishchenko, M. I., L. D. Travis, and A. Macke, Scattering of light by polydisperse, randomly oriented, finite circular cylinders, Appl. Opt., vol. 35, 4927-4940 (1996).
- Wielaard, D. J., M. I. Mishchenko, A. Macke, and B. E. Carlson, Improved T-matrix computations for large, nonabsorbing and weakly absorbing nonspherical particles and comparison with geometrical optics approximation, Appl. Opt., vol. 36, 4305-4313 (1997).
The original Fortran 77 code is available at: http://www.giss.nasa.gov/staff/mmishchenko/t_matrix.html.
The following books were used extensively as source material to build this code:
-
Mishchenko, M. I. , J. W. Hovenier, and L. D. Travis (ed.), Light Scattering by Nonspherical Particles, Academic Press (2000).
-
Bringi, V. N., and V. Chandrasekar, Polarimetric Doppler weather radar: principles and applications, Cambridge University Press (2001).
A general review of the T-matrix approach can be found in:
- Mishchenko, M. I., L. D. Travis, and D. W. Mackowski, T-matrix computations of light scattering by nonspherical particles: a review, J. Quant. Spectrosc. Radiat. Transfer, vol. 55, 535-575 (1996).
Additional useful information is contained in the paper:
- Mishchenko, M. I. and L. D. Travis, Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers, J. Quant. Spectrosc. Radiat. Transfer, vol. 60, 309-324 (1998).
The definitions and notation used can also be found in:
- Mishchenko, M. I., Calculation of the amplitude matrix for a nonspherical particle in a fixed orientation, Appl. Opt., vol. 39, 1026-1031 (2000).
0.3:
- Added ability to integrate angular-integrated quantities over size distributions
- Numerous minor enhancements
0.2:
- Added calculation of generic scattering parameters (scattering cross section, extinction cross section, single-scattering albedo, asymmetry parameter).
- Changed how PSD integration works (see documentation above).
- Replaced quadrature-point code with a pure Python version.
- Renamed a number of things to be more descriptive. Old names now raise deprecation warnings.
0.1: First release
Python code for calculating Mie scattering from single- and dual-layered spheres: pymiecoated.