Merge pull request #42 from mirochaj/surface_density

added more accurate integrator for galaxy surface densities

ares/populations/GalaxyEnsemble.py

@@ -2779,10 +2779,14 @@ def get_lf(self, z, bins=None, use_mags=True, wave=1600.,
         else:
             x = np.arange(6.5, 14, 0.25)
 
+        if yok.sum() == 0:
+            return x, np.zeros_like(x)
+
         # Make sure binning range covers the range of luminosities/magnitudes
         if use_mags:
-            assert y[yok==1].min() > x.min()
-            assert y[yok==1].max() < x.max()
+            mi, ma = y[yok==1].min(), y[yok==1].max()
+            assert mi > x.min(), "{} NOT > {}".format(mi, x.min())
+            assert ma < x.max(), "{} NOT < {}".format(ma, x.max())
         else:
             assert y[yok==1].min() < x.min()
             assert y[yok==1].max() > x.max()
@@ -3431,7 +3435,8 @@ def SurfaceDensity(self, z, bins, dz=1., dtheta=1., wave=1600.,
 
     def get_surface_density(self, z, bins=None, dz=1., dtheta=1., wave=1600.,
         cam=None, filters=None, filter_set=None, dlam=20., method='closest',
-        window=1, load=True, presets=None, absolute=False, use_mags=True):
+        window=1, load=True, presets=None, absolute=False, use_mags=True,
+        use_central_z=True, zstep=0.1, return_evol=False, use_volume=True):
         """
         Compute surface density of galaxies [number / deg^2 / dz]
 
@@ -3442,18 +3447,47 @@ def get_surface_density(self, z, bins=None, dz=1., dtheta=1., wave=1600.,
         square degree.
         """
 
-        # First, compute the luminosity function.
-        x, phi = self.get_lf(z, bins=bins, wave=wave, cam=cam,
-            filters=filters, filter_set=filter_set, dlam=dlam, method=method,
-            window=window, load=load, presets=presets, absolute=absolute,
-            use_mags=use_mags)
+        # Simplest thing: take central redshift, assume same UVLF throughout
+        # dz interval along LOS.
+        if use_central_z:
+            # First, compute the luminosity function.
+            x, phi = self.get_lf(z, bins=bins, wave=wave, cam=cam,
+                filters=filters, filter_set=filter_set, dlam=dlam, method=method,
+                window=window, load=load, presets=presets, absolute=absolute,
+                use_mags=use_mags)
+
+            # Compute the volume of the shell we're looking at [cMpc^3]
+            if use_volume:
+                vol = self.cosm.ProjectedVolume(z, angle=dtheta, dz=dz)
+            else:
+                vol = 1
 
-        # Compute the volume of the shell we're looking at [cMpc^3]
-        vol = self.cosm.ProjectedVolume(z, angle=dtheta, dz=dz)
+            # Get total number of galaxies in volume
+            Ngal = phi * vol
+        else:
+            # Sub-sample redshift interval
+            zbin_e = np.arange(z - 0.5 * dz, z + 0.5 * dz, zstep)
+
+            phi = np.zeros((zbin_e.size, bins.size))
+            vol = np.zeros_like(zbin_e)
+            for i, ze in enumerate(zbin_e):
+                zmid = ze + 0.5 * zstep
+
+                # Compute LF at midpoint of this bin.
+                x, phi[i] = self.get_lf(zmid, bins=bins, wave=wave, cam=cam,
+                    filters=filters, filter_set=filter_set, dlam=dlam, method=method,
+                    window=window, load=load, presets=presets, absolute=absolute,
+                    use_mags=use_mags)
+
+                # Compute the volume of the shell we're looking at [cMpc^3]
+                if use_volume:
+                    vol[i] = self.cosm.ProjectedVolume(zmid, angle=dtheta,
+                        dz=zstep)
+                else:
+                    vol[i] = 1
 
-        # Number of galaxies per mag bin in survey area.
-        # Currently neglects evolution of LF along LoS.
-        Ngal = phi * vol
+            # Integrate over the redshift interval
+            Ngal = np.sum(phi * vol[:,None], axis=0)
 
         # Faint to bright
         Ngal_asc = Ngal[-1::-1]
@@ -3468,7 +3502,30 @@ def get_surface_density(self, z, bins=None, dz=1., dtheta=1., wave=1600.,
         ntot = np.trapz(Ngal, x=x)
         nltm = cumtrapz(Ngal, x=x, initial=Ngal[0])
 
-        return x, nltm
+        if return_evol and (not use_central_z):
+            return x, nltm, zbin_e, phi, vol
+        else:
+            return x, nltm
+
+    def get_volume_density(self, z, bins=None, wave=1600.,
+        cam=None, filters=None, filter_set=None, dlam=20., method='closest',
+        window=1, load=True, presets=None, absolute=False, use_mags=True,
+        use_central_z=True, zstep=0.1, return_evol=False):
+        """
+        Return volume density of galaxies in given `dz` chunk.
+
+        .. note :: Just a wrapper around `get_surface_density`, with
+            hack parameter `use_volume` set to False and `use_central_z` to
+            True.
+
+
+        """
+
+        return self.get_surface_density(z, bins=bins, wave=wave,
+            cam=cam, filters=filters, filter_set=filter_set, dlam=dlam,
+            method=method, window=window, load=load, presets=presets,
+            absolute=absolute, use_mags=use_mags, use_central_z=True,
+            zstep=zstep, return_evol=return_evol, use_volume=False)
 
     def load(self):
         """

ares/static/SpectralSynthesis.py

@@ -138,10 +138,10 @@ def OpticalDepth(self, z, owaves):
 
         # Scorched earth option: null all flux at < 912 Angstrom
         if self.pf['tau_clumpy'] == 1:
-            tau[rwaves <= lam_LL] = np.inf
+            tau[rwaves < lam_LL] = np.inf
         # Or all wavelengths < 1216 A (rest)
         elif self.pf['tau_clumpy'] == 2:
-            tau[rwaves <= lam_LyA] = np.inf
+            tau[rwaves < lam_LyA] = np.inf
         else:
             tau = self.madau1995(z, owaves)
 

input/litdata/blue_tides.py

@@ -120,7 +120,6 @@
  13: {'alpha': -2.13, 'log_phi_star': -11.71, 'M_star': -19.11},
 }
 
-
 def schechter_jaacks(M, log_phi_star=-9, M_star=-20, alpha=-2):
     """
     .. note :: There's a typo in Eq. 5 of Y. Feng et al. (2016). They had

tests/test_populations_ensemble.py

@@ -156,6 +156,17 @@ def test():
     x, Sigma = pop.get_surface_density(6, bins=amag_bins)
     assert 1e3 <= Sigma[np.argmin(np.abs(amag_bins - 27))] <= 1e4
 
+    # Test surface density integral sub-sampling. Should be a small effect.
+    x1, Sigma1 = pop.get_surface_density(6, dz=0.1, bins=amag_bins)
+    x2, Sigma2 = pop.get_surface_density(6, dz=0.1, bins=amag_bins,
+        use_central_z=False, zstep=0.05)
+
+    rdiff = np.abs(Sigma1 - Sigma2) / Sigma2
+    assert rdiff[Sigma2 > 0].mean() < 0.1
+
+    # Try volume density
+    x, n = pop.get_volume_density(6, bins=amag_bins)
+
     # Test nebular line emission stuff
     pars['pop_nebular'] = 2
     pop_neb = ares.populations.GalaxyPopulation(**pars)