A Python code for computing the scattering properties of homogeneous nonspherical scatterers with the T-Matrix method.
- See the installation and usage instructions
- Download the code.
Uses the T-Matrix code by M. I. Mishchenko and L. D. Travis.
Python code for T-matrix scattering calculations
License: MIT License
A Python code for computing the scattering properties of homogeneous nonspherical scatterers with the T-Matrix method.
Uses the T-Matrix code by M. I. Mishchenko and L. D. Travis.
I'm very thankful for your work on the T-matrix method.
Because study python is mush more easy than Fortran. It spend me less time to do more works.
My work is use WRF data to simulate Airborne polarimetric weather radar.
Now I'm preparing a paper.I find you add "Chapter 8 Citing the code", I will add your paper in mine.
Thanks again~
Best wishes.
Happy New year 2016
Why is the reflectivity factor in meters smaller than the normal one 300db, and the one in millimeters larger than the normal one about 60db
The modules written in fortran does not build. After I installed it with pip
and with python setup.py install
the there are no compiled files and there are no fortran .f
files. Can you please provide me the already compiled files? In C one can provide the .c
or .h
files. I guess that should be possible with fortran as well.
Thank you very much.
Hello,
I'm wondering if this code is suitable for both lidar and radar to compute the scattering cross section and extinction cross section in bad weather (such as rain, fog, and dust storms), in other words, whether this code is self-adaption for the Rayleigh scattering and Mie scattering?
Thanks a lot!
I am using radar data in vertical pointing mode so how can i perform scattering simulations in case of vertical incidence ?
scatterer.set_geometry(tmatrix_aux.geom_vert_forw)?
According to the code of 'radar.py', when the unit of the wavelength is 'm', by using 'Ah = 8.686 x wavelength x Imag[Shh]' , we can get the specific attenuation in 'dB/km', so I guess, the unit of the scattering amplitude the code defines is 'km', is that right?
But if so, then, the unit of RCS is 'km^2/m^3', which seems unreliable.
Thanks for your help and patience!
I tried installing this package through the pip install pytmatrix
method. The .so
file failed to compile. However, installing from source after cloning the repo worked.
Running MacOSx 10.11.6
pip 8.1.2
It appears that BinnedDSD and scattering is not playing nicely in the newest version of pytmatrix. As an example the following:
from pytmatrix.tmatrix import Scatterer
from pytmatrix.psd import PSDIntegrator
from pytmatrix import orientation, radar, tmatrix_aux, refractive
scatterer = Scatterer(wavelength=pytmatrix.tmatrix_aux.wl_C,
m=pytmatrix.refractive.m_w_10C[pytmatrix.tmatrix_aux.wl_C])
BinnedDSD = pytmatrix.psd.BinnedPSD([0.1,0.2,0.3,0.4,0.5], [10,10,10,10])
scatterer.psd = BinnedDSD
scatterer.psd_integrator = PSDIntegrator()
scatterer.psd_integrator.axis_ratio_func = lambda D: 1.0 / \
pytmatrix.tmatrix_aux.dsr_pb(D)
dsr_func = pytmatrix.tmatrix_aux.dsr_pb
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)
radar.refl(scatterer)
gives the error
AttributeError: 'NoneType' object has no attribute 'bin_edges'
> /Users/hard505/anaconda/envs/py36/lib/python3.6/site-packages/pytmatrix/psd.py(243)__eq__()
241
242 def __eq__(self, other):
--> 243 return len(self.bin_edges) == len(other.bin_edges) and \
244 (self.bin_edges == other.bin_edges).all() and \
245 (self.bin_psd == other.bin_psd).all()
The part of "radar.py" to calculate Specific attenuation (A) means "Ah=8.686wavelengthImag(Shh)", but according to formula(4.7b) in a book named "Polarimetric Doppler Weather Radar Principles and Applications" by V.N.Bringi(which is also the conference of this code), that should be "Ah=-8.686wavelengthImag(Shh)", there is a minus sign difference. So is this part of code wrong?
Thanks for your help!
When I'm running the example from github, everything works fine, but when I'm running the code with my parameters I'm getting the line: "Process finished with exit code 0" and can't get the results of the calculated array.
The example from github:
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]])
my code:
scatterer = tmatrix.Scatterer(radius=50, wavelength=5, m=complex(0.1903, 0), axis_ratio=180/100)
scatterer.get_S()
I have got the amplitude scattering matrix, then how to get volume scattering function or phase function with an incident radiant of none polarization by using this code.? thanks a lot.
It would seem that going forward this package might be limited to Python < 3.12, unless the package is updated to account for the new numpy changes. From this page: https://numpy.org/doc/stable/reference/distutils.html
numpy.distutils is deprecated, and will be removed for Python >= 3.12.
Not sure how much work it is to parse this over or if its worth it.
I set m=refractive.ice_refractive('ice_refr.dat')(8.43,0.2) and mu=0(GammaDSD),the other settings are the same as the example given by Wiki which is to compute radar scattering properties for a C-band radar. The result is too small(ref=9.55dBz WhenD0=1.1,Nw=2500)
It's shown that the RCS calculation of particles is based on the Muller Matrix, so the result of calculation should be RCS with unit of 'm^2/m^3', i.e., the RCS in a unit volume. So I'm wondering how to get the single particle's RCS with unit of 'm^2'.
I'm looking forward to your reply. Thank you for your kindness!
Do we have any function to calculate Nw if we have Median drop diameter (Do) and reflectivity Z ?
Hi, I am running into a problem when i try and create a scatterer with a diameter that is ~40x larger than the wavelength. Is this a bug?
Hello,
I am using this code to calculate scattering coefficient and extinction coefficient. But I find that the execution of the code stops when I try to calculate the T-matrix of a scatterer with a size parameter which is about 80 while the axis_ratio = 1, and the size parameter is reduced to about 18 while the axis_ratio = 3. Is this a bug? Or it's a existing limit of the fortran numerical code?
According to table 1 in the paper named "Capabilities and limitations of a current FORTRAN implementation of the T-matrix method for randomly oriented, rotationally symmetric scatterers" by M. I. Mishchenko, the maximal convergent size parameters seem to be much larger than the limit I find.
By the way, I have tried to increase value of the attribute ndgs in the Scatterer class, but still I can't calculate the T-matrix when size parameter is large.
So, I want to know how large size parameters can this code calculate? And what should I do to calculate T-matrix for larger size parameter?
I'm looking forward to your reply. Thanks for your kindness!
Best,
Lee
Sorry to bother you, but I want to know the effect of scatterer.psd_integrator.geometries = (geom_back, geom_forw)
As far as I know, the cross sections are in units of area. But in the code the calculation of extinction cross section is
2 * scatterer.wavelength * S[1,1].imag
which would give a unit of a length? Can you please provide the reference of this equation? This would be really nice.
Thank you.
It looks like the current conda package is out of date on conda-forge. It still appears to have the bug that the most recent commit fixed (scattering table IO). I'm not sure who is maintaining, but it would be good to re-issue these packages based on newest code.
LICENSE.txt
in source folder claims to be BSD-3-clause.As LICENSE.txt
is not distributed with PyPI, people who install from there might assume this package is MIT, whereas people who are installing from GITHUB might assume this is BSD-3-Clause.
Now, what license should be used when redistributing this package?
I try to create a LUT for rain using your code. The code is dedicated for the nadir looking radar so I use the following settings: beta=90,thet0=90,thet=90,angular_integration=True (correct me if it is wrong)
The code works fine with a single scatterer
scatterer = Scatterer(radius=3.0, wavelength=8, m=complex(1.5,0.5), axis_ratio=1.0/0.6, beta=90,thet0=90,thet=90,angular_integration=True)
scatter.asym(scatterer)
0.6983014838870412
np.log10(radar.refl(scatterer))*10
21.354628070158022
np.log10(radar.Zdr(scatterer))*10
0.0
moreover when I define the PSD everything works just fine
scatterer = Scatterer(wavelength=tmatrix_aux.wl_Ka, m=refractive.m_w_10C[tmatrix_aux.wl_C],axis_ratio=1.0/0.6,beta=90.0)
scatterer.psd_integrator = PSDIntegrator()
scatterer.psd_integrator.D_max = 4.0
scatterer.or_pdf = orientation.gaussian_pdf(20.0)
scatterer.orient = orientation.orient_averaged_fixed
scatterer.psd_integrator.num_points=2**3
scatterer.psd_integrator.init_scatter_table(scatterer, verbose=True,
angular_integration=True)
but when I try to get the reflectivity of the ensambe of partiles then I get the following:
10*np.log10(radar.refl(scatterer))
Traceback (most recent call last):File "", line 1, in 10*np.log10(radar.refl(scatterer))
File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/radar.py", line 62, in refl
radar_xsect(scatterer, h_pol)File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/radar.py", line 38, in radar_xsect
Z = scatterer.get_Z()File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/tmatrix.py", line 341, in get_Z
return self.get_SZ()[1]File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/tmatrix.py", line 329, in get_SZ
self.get_geometry())File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/psd.py", line 307, in call
return self.get_SZ(psd, geometry)File "/opt/miniconda3/lib/python3.6/site-packages/pytmatrix/psd.py", line 334, in get_SZ
return (self._S_dict[geometry], self._Z_dict[geometry])KeyError: (90.0, 90.0, 0.0, 180.0, 0.0, 90.0)
which suggests that the geometry setting is causing problems
Hello,
I would like to inform you that while running "test_tmatrix.run_tests()", I have got the following warning:
I reported it here, as I am not sure whether this waring is expected.
Best,
BA
Hi @jleinonen ,
When I compute Ah and Av of rain with a prescribed canting angle using the psd_integrator I get exactly the same results for Ah and for Av. In fact the parameter h_pol in radar.Ai(scatterer, h_pol=False) seems to be ignored when using the integrator.
I am missing something or is this indeed a bug in the code?
Cheers,
Hi,
there is an error, probably a typo, in the function dsr_bc() in tmatrix_aux: The factor D_eq is missing for the second summand (see for example Bringi and Chandrasekar 2001, equation 7.3).
Regards and thank you for this python interface for the T-matrix method!
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.