Better deal with division by zero in mean opacity calculations

.travis.yml

@@ -12,7 +12,6 @@ env:
         - SETUP_CMD='test'
         - MPLBACKEND='Agg'
     matrix:
-        - PYTHON_VERSION=2.7
         - PYTHON_VERSION=3.6
         - PYTHON_VERSION=3.7
 

hyperion/dust/mean_opacities.py

@@ -64,9 +64,6 @@ def compute(self, optical_properties, n_temp=1200, temp_min=0.1, temp_max=100000
         kappa_nu = interp1d_fast_loglog(optical_properties.nu,
                                         optical_properties.kappa, nu)
 
-        # Set mean opacity variable
-        self.var_name = 'specific_energy'
-
         # Initialize mean opacity arrays
         self.chi_planck = np.zeros(n_temp)
         self.kappa_planck = np.zeros(n_temp)
@@ -82,29 +79,31 @@ def compute(self, optical_properties, n_temp=1200, temp_min=0.1, temp_max=100000
             b_nu = B_nu(nu, T)
             db_nu_dt = dB_nu_dT(nu, T)
 
+            # We use np.divide here to do the right thing with division by zero
+
             # Compute planck mean opacity
-            self.chi_planck[it] = (integrate_loglog(nu, b_nu * chi_nu) /
-                                   integrate_loglog(nu, b_nu))
+            self.chi_planck[it] = np.divide(integrate_loglog(nu, b_nu * chi_nu),
+                                            integrate_loglog(nu, b_nu))
 
             # Compute planck mean absoptive opacity
-            self.kappa_planck[it] = (integrate_loglog(nu, b_nu * kappa_nu) /
-                                     integrate_loglog(nu, b_nu))
+            self.kappa_planck[it] = np.divide(integrate_loglog(nu, b_nu * kappa_nu),
+                                              integrate_loglog(nu, b_nu))
 
             # Compute reciprocal planck mean opacity
-            self.chi_inv_planck[it] = (integrate_loglog(nu, b_nu) /
-                                       integrate_loglog(nu, b_nu / chi_nu))
+            self.chi_inv_planck[it] = np.divide(integrate_loglog(nu, b_nu),
+                                                integrate_loglog(nu, b_nu / chi_nu))
 
             # Compute reciprocal planck mean aborptive opacity
-            self.kappa_inv_planck[it] = (integrate_loglog(nu, b_nu) /
-                                         integrate_loglog(nu, b_nu / kappa_nu))
+            self.kappa_inv_planck[it] = np.divide(integrate_loglog(nu, b_nu),
+                                                  integrate_loglog(nu, b_nu / kappa_nu))
 
             # Compute rosseland mean opacity
-            self.chi_rosseland[it] = (integrate_loglog(nu, db_nu_dt) /
-                                      integrate_loglog(nu, db_nu_dt / chi_nu))
+            self.chi_rosseland[it] = np.divide(integrate_loglog(nu, db_nu_dt),
+                                               integrate_loglog(nu, db_nu_dt / chi_nu))
 
             # Compute rosseland mean aborptive opacity
-            self.kappa_rosseland[it] = (integrate_loglog(nu, db_nu_dt) /
-                                        integrate_loglog(nu, db_nu_dt / kappa_nu))
+            self.kappa_rosseland[it] = np.divide(integrate_loglog(nu, db_nu_dt),
+                                                 integrate_loglog(nu, db_nu_dt / kappa_nu))
 
         self.temperature = temperatures
         self.specific_energy = 4. * sigma * temperatures ** 4. * self.kappa_planck

hyperion/dust/tests/test_mean_opacities.py

@@ -9,7 +9,7 @@
 from ..optical_properties import OpticalProperties
 from ..mean_opacities import MeanOpacities
 from ...util.functions import virtual_file, B_nu
-from ...util.constants import sigma
+from ...util.constants import sigma, c
 
 
 class TestMeanOpacities(object):
@@ -65,3 +65,17 @@ def test_plot(self):
         m.plot(fig, 111)
 
         plt.close(fig)
+
+
+def test_zero_emissivities():
+
+    wav = np.array([2.0, 0.5])
+    nu = c / (wav * 1.e-4)
+
+    o = OpticalProperties()
+    o.nu = nu
+    o.chi = np.array([1.e-2, 1])
+    o.albedo = np.array([0., 0.5])
+
+    m = MeanOpacities()
+    m.compute(o)