pytmatrix icon indicating copy to clipboard operation
pytmatrix copied to clipboard

Something wrong with the computation of Av?

Open jfigui opened this issue 1 year ago • 12 comments

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,

@jfigui

jfigui avatar May 09 '23 14:05 jfigui

Looking at the code I don't see why it would be ignored. radar.Ai calls scatter.ext_xsect which does behave differently (uses a different element of the matrix S) based on h_pol. Could you post a working example of that replicates the problem?

jleinonen avatar May 09 '23 15:05 jleinonen

diam_min = 0.1 diam_max = 7. step = 0.1 diam = np.arange(diam_min, diam_max+step, step) num_points = diam.size canting_angle = 10.

geom_back = (90.0, 90.0, 0.0, 180.0, 0.0, 0.0) geom_forw = (90.0, 90.0, 0.0, 0.0, 0.0, 0.0)

wavelength = tmatrix_aux.wl_C m = refractive.m_w_20C[wavelength] scatterer = Scatterer(wavelength=wavelength, m=m) scatterer.orient = orientation.orient_averaged_fixed scatterer.or_pdf = orientation.gaussian_pdf(canting_angle) scatterer.psd_integrator = PSDIntegrator() scatterer.psd_integrator.axis_ratio_func = (lambda diam: 1.0/tmatrix_aux.dsr_thurai_2007(diam)) scatterer.psd_integrator.D_max = diam_max scatterer.psd_integrator.num_points = num_points scatterer.psd_integrator.geometries = (geom_back, geom_forw) scatterer.psd_integrator.init_scatter_table(scatterer, angular_integration=True, verbose=True)

psd = GammaPSD(D0=3, Nw=10e3, mu=0, D_max=diam_max)
scatterer.psd = psd scatterer.set_geometry(geom_forw) Ah = radar.Ai(scatterer) Av = radar.Ai(scatterer, False)

jfigui avatar May 10 '23 06:05 jfigui

@jleinonen,

radar.Ai(scatterer, h_pol=True) calls scatter.ext_xsect(scatterer, h_pol=h_pol)

If the PSD integrator is not None ext_xsect returns:

scatterer.psd_integrator.get_angular_integrated( scatterer.psd, scatterer.get_geometry(), "sca_xsect")

but in this call the information of which polarization you want is lost. By default you get the horizontal polarization regardless of the intended polarization.

Am I right or I am missing something?

jfigui avatar May 10 '23 06:05 jfigui

Looks like there's indeed a bug. I'm working on a fix.

jleinonen avatar May 10 '23 10:05 jleinonen

Thanks a lot!

jfigui avatar May 10 '23 11:05 jfigui

Hi @jfigui, could you check if it works for you with the latest commit?

jleinonen avatar May 10 '23 12:05 jleinonen

Hi @jleinonen ,

It works for me now. See the output: psd_rain_C_2000_ele00000_refl_h-A psd_rain_C_2000_ele00000_refl_h-Adp

Could you make a new release soon? I work with the pypi package usually

jfigui avatar May 10 '23 13:05 jfigui

I don't manage the conda package (someone created it without my knowledge) so I have little control over it. As I understand it gets updated when I create a new release.

jleinonen avatar May 10 '23 13:05 jleinonen

I think if you create a new Pypi package you should get a new conda package at some point. At least this is how it works for us using conda-forge

jfigui avatar May 10 '23 13:05 jfigui

Thanks for being so reactive. Let me know when the new release is ready

jfigui avatar May 10 '23 13:05 jfigui

Version 0.3.3 is now on PyPI. The conda package should eventually update automatically, as mentioned.

jleinonen avatar May 10 '23 16:05 jleinonen

Great!

Thanks a lot for being so reactive. For info I am using pytmatrix in a software package to compute the scattering properties of different hydrometeors. The code is now public and you can find it here: https://github.com/openradar/hydroscatt

I would be grateful if you have a look and report any bug.

jfigui avatar May 11 '23 07:05 jfigui