Fixed calculatiion of angular tables with both polarizations - Issue #30

pytmatrix/psd.py

@@ -334,34 +334,35 @@ 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
             )
 
         if property_name == "sca_xsect":
             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,20 +443,21 @@ 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) = \
                 (old_m, old_axis_ratio, old_radius, old_psd_integrator) 
             tm.set_geometry(old_geom)
 
 
-
     def save_scatter_table(self, fn, description=""):
         """Save the scattering lookup tables.
         

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()
 

pytmatrix/tmatrix.py

@@ -27,7 +27,6 @@
 import pytmatrix.orientation as orientation
 
 
-
 class Scatterer(object):
     """T-Matrix scattering from nonspherical particles.