jleinonen / pytmatrix

Python code for T-matrix scattering calculations
MIT License
101 stars 45 forks source link

PSD integrator doesnt work with beta=90 #11

Closed mrozkamil closed 7 years ago

mrozkamil commented 7 years ago

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

jleinonen commented 7 years ago

I think what's missing is that you're not specifying the geometries in the PSDIntegrator. You need to do something like this: scatterer.psd_integrator.geometries = (tmatrix_aux.geom_horiz_back, tmatrix_aux.geom_horiz_forw)

You can replace the predefined geometries with custom ones, but it's a good idea include both backward and forward scattering geometries. Check out the last example in the wiki for an example of using PSDIntegrator.

Let me know if this helped - or if it doesn't, feel free to ask further questions.

mrozkamil commented 7 years ago

Your suggestion works. The script is not crushing any more but what worries me is a fact that I get strange results, i.e. with scatterer.psd_integrator.geometries = ((90.0, 90.0, 0.0, 180.0, 0.0, 90.0),) scatterer.psd = GammaPSD(D0=1.0, Nw=1e3, mu=3) I get ZDR of 4.59653202257 dB (only 16 sizes are calculated but it shouldn't affect the result). I was sure that I must get something around machine precision 0.

The same geometry setting for a single particle: 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,radius=1.0) gives ZDR of 0 dB.

jleinonen commented 7 years ago

I think your problem is that you are enabling orientation averaging in the size distribution integrated version, on these lines: scatterer.or_pdf = orientation.gaussian_pdf(20.0) scatterer.orient = orientation.orient_averaged_fixed These will override the beta parameter for Scatterer so that your beta will be close to zero, which leads to a nonzero ZDR because you have an oblate scatterer.

mrozkamil commented 7 years ago

Please correct me if I am mistaken. When I put beta = 90 it means that the axis of rotation of the particle is aligned along the direction of the propagation, therefore ZDR is equal zero, but then orientation.gaussian_pdf(20.0) overwrites beta, so the question is how can I keep the mean orientation of particles and add some variability. I tried to use orientation.gaussian_pdf(std=10.0, mean=90.0) but it also resulted in nonzero ZDR, for scatterer.psd = GammaPSD(D0=1.0, Nw=1e3, mu=3), ZDR=-2.2926145941.

jleinonen commented 7 years ago

This is actually a new feature that is so far missing from the wiki. If you do scatterer.or_pdf = orientation.gaussian_pdf(std=20.0, mean=90.0) that should do what you want.

mrozkamil commented 6 years ago

Thanks for your reply, I found another "missing functionality", i.e. when we integrate over a PSD then we should also integrate over the alpha parameter, otherwise we are assuming that particles oscillate only around the axis perpendicular to the direction of propagation. In my case I again get ZDR much different than zero.

jleinonen commented 6 years ago

The alpha parameter is always integrated over when orientation averaging is enabled.

mrozkamil commented 6 years ago

I understand that it is integrated to get e.g. the asymmetry parameter but I am talking about integration over the PSD. If we use orientation.gaussian_pdf(std=20.0) then all of our particles have different canting angles oscillating around 0deg but they are tilted along the same axis. Isn't it true?

jleinonen commented 6 years ago

orientation.gaussian_pdf(std=20.0) implies that the beta angle follows a normal distribution with a mean of 0 and a standard deviation of 20 degrees, and the alpha angle is uniformly distributed. If orientation averaging is enabled, integration over the alpha angle is always done.