From 95ed4c1b618c5f09b87d3dc36bbe774c22f7fd81 Mon Sep 17 00:00:00 2001 From: Jussi Leinonen Date: Wed, 10 May 2023 12:12:04 +0200 Subject: [PATCH] Fixed calculatiion of angular tables with both polarizations - Issue #30 --- pytmatrix/psd.py | 38 ++++++++++++++++++++++---------------- pytmatrix/scatter.py | 14 +++++++++----- pytmatrix/tmatrix.py | 1 - 3 files changed, 31 insertions(+), 22 deletions(-) diff --git a/pytmatrix/psd.py b/pytmatrix/psd.py index 811d41e..b85ea37 100644 --- a/pytmatrix/psd.py +++ b/pytmatrix/psd.py @@ -334,18 +334,19 @@ def get_SZ(self, psd, geometry): return (self._S_dict[geometry], self._Z_dict[geometry]) - def get_angular_integrated(self, psd, geometry, property_name): + def get_angular_integrated(self, psd, geometry, property_name, h_pol=True): if self._angular_table is None: raise AttributeError( "Initialize or load the table of angular-integrated " + "quantities first." ) + pol = "h_pol" if h_pol else "v_pol" psd_w = psd(self._psd_D) def sca_xsect(geom): return trapz( - self._angular_table["sca_xsect"][geom] * psd_w, + self._angular_table["sca_xsect"][pol][geom] * psd_w, self._psd_D ) @@ -353,15 +354,15 @@ def sca_xsect(geom): sca_prop = sca_xsect(geometry) elif property_name == "ext_xsect": sca_prop = trapz( - self._angular_table["ext_xsect"][geometry] * psd_w, + self._angular_table["ext_xsect"][pol][geometry] * psd_w, self._psd_D ) elif property_name == "asym": sca_xsect_int = sca_xsect(geometry) if sca_xsect_int > 0: sca_prop = trapz( - self._angular_table["asym"][geometry] * \ - self._angular_table["sca_xsect"][geometry] * psd_w, + self._angular_table["asym"][pol][geometry] * \ + self._angular_table["sca_xsect"][pol][geometry] * psd_w, self._psd_D ) sca_prop /= sca_xsect_int @@ -399,8 +400,11 @@ def init_scatter_table(self, tm, angular_integration=False, verbose=False): self._previous_psd = None self._m_table = np.empty(self.num_points, dtype=complex) if angular_integration: - self._angular_table = {"sca_xsect": {}, "ext_xsect": {}, - "asym": {}} + self._angular_table = { + "sca_xsect": {"h_pol": {}, "v_pol": {}}, + "ext_xsect": {"h_pol": {}, "v_pol": {}}, + "asym": {"h_pol": {}, "v_pol": {}} + } else: self._angular_table = None @@ -419,8 +423,9 @@ def init_scatter_table(self, tm, angular_integration=False, verbose=False): if angular_integration: for int_var in ["sca_xsect", "ext_xsect", "asym"]: - self._angular_table[int_var][geom] = \ - np.empty(self.num_points) + for pol in ["h_pol", "v_pol"]: + self._angular_table[int_var][pol][geom] = \ + np.empty(self.num_points) for (i,D) in enumerate(self._psd_D): if verbose: @@ -438,12 +443,14 @@ def init_scatter_table(self, tm, angular_integration=False, verbose=False): self._Z_table[geom][:,:,i] = Z if angular_integration: - self._angular_table["sca_xsect"][geom][i] = \ - scatter.sca_xsect(tm) - self._angular_table["ext_xsect"][geom][i] = \ - scatter.ext_xsect(tm) - self._angular_table["asym"][geom][i] = \ - scatter.asym(tm) + for pol in ["h_pol", "v_pol"]: + h_pol = (pol == "h_pol") + self._angular_table["sca_xsect"][pol][geom][i] = \ + scatter.sca_xsect(tm, h_pol=h_pol) + self._angular_table["ext_xsect"][pol][geom][i] = \ + scatter.ext_xsect(tm, h_pol=h_pol) + self._angular_table["asym"][pol][geom][i] = \ + scatter.asym(tm, h_pol=h_pol) finally: #restore old values (tm.m, tm.axis_ratio, tm.radius, tm.psd_integrator) = \ @@ -451,7 +458,6 @@ def init_scatter_table(self, tm, angular_integration=False, verbose=False): tm.set_geometry(old_geom) - def save_scatter_table(self, fn, description=""): """Save the scattering lookup tables. diff --git a/pytmatrix/scatter.py b/pytmatrix/scatter.py index 945ddb5..1010cb3 100644 --- a/pytmatrix/scatter.py +++ b/pytmatrix/scatter.py @@ -78,7 +78,9 @@ def sca_xsect(scatterer, h_pol=True): if scatterer.psd_integrator is not None: return scatterer.psd_integrator.get_angular_integrated( - scatterer.psd, scatterer.get_geometry(), "sca_xsect") + scatterer.psd, scatterer.get_geometry(), "sca_xsect", + h_pol=h_pol + ) old_geom = scatterer.get_geometry() @@ -112,7 +114,9 @@ def ext_xsect(scatterer, h_pol=True): if scatterer.psd_integrator is not None: try: return scatterer.psd_integrator.get_angular_integrated( - scatterer.psd, scatterer.get_geometry(), "ext_xsect") + scatterer.psd, scatterer.get_geometry(), "ext_xsect", + h_pol=h_pol + ) except AttributeError: # Fall back to the usual method of computing this from S pass @@ -125,8 +129,6 @@ def ext_xsect(scatterer, h_pol=True): finally: scatterer.set_geometry(old_geom) - - if h_pol: return 2 * scatterer.wavelength * S[1,1].imag else: @@ -163,7 +165,9 @@ def asym(scatterer, h_pol=True): if scatterer.psd_integrator is not None: return scatterer.psd_integrator.get_angular_integrated( - scatterer.psd, scatterer.get_geometry(), "asym") + scatterer.psd, scatterer.get_geometry(), "asym", + h_pol=h_pol + ) old_geom = scatterer.get_geometry() diff --git a/pytmatrix/tmatrix.py b/pytmatrix/tmatrix.py index ef920f6..74726fc 100644 --- a/pytmatrix/tmatrix.py +++ b/pytmatrix/tmatrix.py @@ -27,7 +27,6 @@ import pytmatrix.orientation as orientation - class Scatterer(object): """T-Matrix scattering from nonspherical particles.