diff_... and integrate..._complex_medium() now take 1D thetas and phis

vnmanoharan · vnmanoharan · commit 33a0471d8705 · 2025-08-09T21:30:01.000-04:00
diff --git a/pymie/mie.py b/pymie/mie.py
@@ -1148,20 +1148,18 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, phis=None,
     x : complex, array-like
         size parameter, x = ka = 2*pi*n_med/lambda * a (sphere radius a)
     thetas : array-like
-        scattering angles.  Should be 2D, as output from np.meshgrid, if
-        cartesian=True
+        Scattering angles.  Must be in radians.
     kd : float
         k * distance, where k = 2*np.pi*n_matrix/wavelen, and distance is the
         distance away from the center of the particle. The standard far-field
         solutions are obtained when distance >> radius in a non-absorbing
         medium.
     phis : None or ndarray
-        azimuthal angles for which to calculate the diff scat intensity. In the
+        Azimuthal angles for which to calculate the diff scat intensity. In the
         'scattering plane' coordinate system, the scattering matrix does not
         depend on phi, so phi should be set to None. In the 'cartesian'
         coordinate system, the scattering matrix does depend on phi, so an
-        array of values should be provided.  For 'cartesian' both thetas and
-        phis should be 2D, as output from np.meshgrid.
+        array of values should be provided.
     cartesian : boolean (default False)
         If False (default), scattering calculations will be carried out in the
         'scattering plane' coordinate system, defined by basis vectors parallel
@@ -1233,9 +1231,10 @@ def diff_scat_intensity_complex_medium(m, x, thetas, kd, phis=None,
     # ensure that broadcasting will work correctly by adding an axis
     # corresponding to theta
     kd = np.atleast_1d(kd)[..., np.newaxis]
-    if cartesian:
+    if phis is not None:
         # add another axis to correspond to phi
         kd = kd[..., np.newaxis]
+        thetas, phis = np.meshgrid(thetas, phis, indexing="ij")
 
     if near_field:
         if not cartesian:
@@ -1362,7 +1361,7 @@ def integrate_intensity_complex_medium(dscat, thetas, kd,
                           "depend on azimuthal angle, so specified values "
                           "will be ignored")
 
-        # strip units from integrand
+        # include Jacobian
         integrand_par = dsigma_1 * np.abs(np.sin(thetas))
         integrand_perp = dsigma_2 * np.abs(np.sin(thetas))
 
diff --git a/pymie/tests/test_mie.py b/pymie/tests/test_mie.py
@@ -730,9 +730,7 @@ def test_diff_scat_intensity_complex_medium_cartesian():
     n_particle = 1.59 + 1e-4 * 1.0j
     thetas = np.linspace(np.pi/2, np.pi, 4)
     phis = np.linspace(0, 2*np.pi, 3)
-    thetas_2d, phis_2d = np.meshgrid(thetas, phis) # be careful with meshgrid shape.
-                                                   # for integration, theta dimension must always come first,
-                                                   # which is not how it is done here
+    thetas_2d, _ = np.meshgrid(thetas, phis, indexing="ij")
 
     # parameters for calculating scattering
     m = index_ratio(n_particle, n_matrix)
@@ -748,8 +746,8 @@ def test_diff_scat_intensity_complex_medium_cartesian():
     # calculate differential scattered intensity in xy basis
     # if incident vector is unpolarized (1,1), then the resulting differential
     # scattered intensity should be the same as I_par, I_perp
-    I_xy = mie.diff_scat_intensity_complex_medium(m, x, thetas_2d, kd,
-                            cartesian=True, phis = phis_2d,
+    I_xy = mie.diff_scat_intensity_complex_medium(m, x, thetas, kd,
+                            cartesian=True, phis = phis,
                             near_field=False, incident_vector = (1, 1))
 
     # calculate magnitudes
@@ -771,8 +769,7 @@ def test_integrate_intensity_complex_medium_cartesian():
     n_particle = 1.59 + 1e-4 * 1.0j
     thetas = np.linspace(0, np.pi, 500)
     phis = np.linspace(0, 2*np.pi, 550)
-    phis_2d, thetas_2d = np.meshgrid(phis, thetas) # remember, meshgrid shape is (len(thetas), len(phis))
-                                                   # and theta dimension MUST come first in these calculations
+
     # parameters for calculating scattering
     m = index_ratio(n_particle, n_matrix)
     x = size_parameter(wavelen, n_matrix, radius)
@@ -781,7 +778,7 @@ def test_integrate_intensity_complex_medium_cartesian():
     kd = (k*distance).to("").magnitude
 
     # calculate the differential scattered intensities
-    I_xy = mie.calc_ang_scat(m, x, thetas_2d, kd=kd, phis=phis_2d)
+    I_xy = mie.calc_ang_scat(m, x, thetas, kd=kd, phis=phis)
     I_parperp = mie.calc_ang_scat(m, x, thetas, kd=kd)
 
     # integrate the differential scattered intensities
@@ -814,7 +811,6 @@ def test_value_errors():
     n_particle = 1.59 + 1e-4 * 1.0j
     thetas = np.linspace(np.pi/2, np.pi, 4)
     phis = np.linspace(0, 2*np.pi, 3)
-    thetas_2d, phis_2d = np.meshgrid(thetas, phis)
 
     # parameters for calculating scattering
     m = index_ratio(n_particle, n_matrix)
@@ -825,13 +821,13 @@ def test_value_errors():
 
     with pytest.raises(ValueError):
         # try to calculate near field in cartesian
-        _ = mie.diff_scat_intensity_complex_medium(m, x, thetas_2d, kd,
+        _ = mie.diff_scat_intensity_complex_medium(m, x, thetas, kd,
                                                    cartesian=True,
-                                                   phis=phis_2d,
+                                                   phis=phis,
                                                    near_field=True)
     # calculate the differential scattered intensities
-    I_xy = mie.diff_scat_intensity_complex_medium(m, x, thetas_2d, kd,
-                                                  cartesian=True, phis=phis_2d,
+    I_xy = mie.diff_scat_intensity_complex_medium(m, x, thetas, kd,
+                                                  cartesian=True, phis=phis,
                                                   near_field=False)
 
     _ = mie.diff_scat_intensity_complex_medium(m, x, thetas, kd,
@@ -842,7 +838,7 @@ def test_value_errors():
         _ = mie.integrate_intensity_complex_medium(I_xy, thetas, kd,
                                                    cartesian=True)[0]
     with pytest.raises(ValueError):
-        _ = mie.vector_scattering_amplitude(m, x, thetas_2d,
+        _ = mie.vector_scattering_amplitude(m, x, thetas,
                                             cartesian=True)
 
 def test_dwell_time_and_energy():
diff --git a/pymie/tests/test_mie_vectorized.py b/pymie/tests/test_mie_vectorized.py
@@ -599,20 +599,30 @@ def test_vectorized_angular_functions(self, num_wavelen,
             mx(num_wavelen=num_wavelen, num_layer=num_layer, **self.mxargs,
                return_all=True)
 
+        num_phi = 13
+
+        # while diff_scat_intensity_complex_medium() and
+        # integrate_intensity_complex_medium() take 1D arrays of phis and
+        # thetas, vector_scattering_amplitude() and
+        # amplitude_scattering_matrix() do not.  So we have to set up
+        # meshgrid arrays and non-meshgrid and send them to the appropriate
+        # functions
+        thetas = self.thetas
         if not cartesian:
             phis = None
-            thetas = self.thetas
+            phis_2d = None
+            thetas_2d = thetas
         else:
-            phis = np.linspace(0, 2*np.pi, self.num_theta)
-            thetas, phis = np.meshgrid(self.thetas, phis)
+            phis = np.linspace(0, 2*np.pi, num_phi)
+            thetas_2d, phis_2d = np.meshgrid(self.thetas, phis, indexing="ij")
 
-        vsa = mie.vector_scattering_amplitude(m, x, thetas,
+        vsa = mie.vector_scattering_amplitude(m, x, thetas_2d,
                                               cartesian=cartesian,
-                                              phis = phis)
+                                              phis = phis_2d)
 
-        mat = mie.amplitude_scattering_matrix(m, x, thetas,
+        mat = mie.amplitude_scattering_matrix(m, x, thetas_2d,
                                               cartesian=cartesian,
-                                              phis = phis)
+                                              phis = phis_2d)
 
         # choose distance reasonably close to the particle for differential
         # scattering calculations
@@ -632,14 +642,14 @@ def test_vectorized_angular_functions(self, num_wavelen,
 
         # check that shapes of all the computed quantities are correct
         for element in vsa + mat:
-            assert element.shape == (num_wavelen, ) + thetas.shape
-        assert i12.shape == (2, num_wavelen) + thetas.shape
+            assert element.shape == (num_wavelen, ) + thetas_2d.shape
+        assert i12.shape == (2, num_wavelen) + thetas_2d.shape
 
         # check that vectorized calculations match looped calculations over
         # scalars
-        amp0 = np.zeros((num_wavelen, ) + thetas.shape, dtype=complex)
+        amp0 = np.zeros((num_wavelen, ) + thetas_2d.shape, dtype=complex)
         amp1 = np.zeros_like(amp0)
-        S1 = np.zeros((num_wavelen, ) + thetas.shape, dtype=complex)
+        S1 = np.zeros((num_wavelen, ) + thetas_2d.shape, dtype=complex)
         S2 = np.zeros_like(S1)
         S3 = np.zeros_like(S1)
         S4 = np.zeros_like(S1)
@@ -649,21 +659,23 @@ def test_vectorized_angular_functions(self, num_wavelen,
         sigma = np.zeros(num_wavelen)
         sigma_1 = np.zeros_like(sigma)
         sigma_2 = np.zeros_like(sigma)
-        dsigma_1 = np.zeros((num_wavelen, ) + thetas.shape)
+        dsigma_1 = np.zeros((num_wavelen, ) + thetas_2d.shape)
         dsigma_2 = np.zeros_like(dsigma_1)
 
         m = np.atleast_1d(m)
         x = np.atleast_1d(x)
         kd = np.atleast_1d(kd)
         for i in range(num_wavelen):
             # m[[i]] preserves 2D array
-            mat_loop = mie.amplitude_scattering_matrix(m[[i]], x[[i]], thetas,
+            mat_loop = mie.amplitude_scattering_matrix(m[[i]], x[[i]],
+                                                       thetas_2d,
                                                        cartesian=cartesian,
-                                                       phis = phis)
+                                                       phis = phis_2d)
 
-            vsa_loop = mie.vector_scattering_amplitude(m[[i]], x[[i]], thetas,
+            vsa_loop = mie.vector_scattering_amplitude(m[[i]], x[[i]],
+                                                       thetas_2d,
                                                        cartesian=cartesian,
-                                                       phis = phis)
+                                                       phis = phis_2d)
             i_loop = mie.diff_scat_intensity_complex_medium(m[[i]], x[[i]],
                                                             thetas,
                                                             kd[i],
