Changing fsky to sky_area (#314)

Author: philipp128

docs/pipeline/examples/lightcone.yml

@@ -10,4 +10,4 @@ tables:
       phi_star: 0.0017
       alpha: -0.5
       m_lim: 25
-      fsky: 1.0e-6
+      sky_area: 0.05 deg2

examples/gravitational_wave_rates.yml

@@ -1,7 +1,7 @@
 cosmology: !astropy.cosmology.default_cosmology.get
 z_range: !numpy.arange [0, 2.01, 0.1]
 mag_lim: 30
-fsky: 0.0000242407
+sky_area: 1 deg2
 M_star_red: !astropy.modeling.models.Linear1D [-0.70798041, -20.37196157]
 phi_star_red: !astropy.modeling.models.Exponential1D [0.0035097, -1.41649]
 alpha_red: -0.5
@@ -16,7 +16,7 @@ tables:
       phi_star: $phi_star_red
       alpha: $alpha_red
       m_lim: $mag_lim
-      fsky: $fsky
+      sky_area: $sky_area
     magnitude: !skypy.galaxy.luminosity.schechter_lf_magnitude
       redshift: $merger_rates.redshift
       M_star: $M_star_red

examples/mccl_galaxies.yml

@@ -3,7 +3,7 @@ cosmology: !astropy.cosmology.FlatLambdaCDM
   Om0: 0.307
 z_range: !numpy.arange [0, 2.01, 0.1]
 mag_lim: 30
-fsky: 0.0000242407
+sky_area: 1 deg2
 M_star_red: !astropy.modeling.models.Linear1D [-0.70798041, -20.37196157]
 phi_star_red: !astropy.modeling.models.Exponential1D [0.0035097, -1.41649]
 alpha_red: -0.5
@@ -20,7 +20,7 @@ tables:
       phi_star: $phi_star_red
       alpha: $alpha_red
       m_lim: $mag_lim
-      fsky: $fsky
+      sky_area: $sky_area
     spectral_coefficients: !skypy.galaxy.spectrum.dirichlet_coefficients
       alpha0: [2.461, 2.358, 2.568, 2.268, 2.402]
       alpha1: [2.410, 2.340, 2.200, 2.540, 2.464]
@@ -47,7 +47,7 @@ tables:
       phi_star: $phi_star_blue
       alpha: $alpha_blue
       m_lim: $mag_lim
-      fsky: $fsky
+      sky_area: $sky_area
     spectral_coefficients: !skypy.galaxy.spectrum.dirichlet_coefficients
       alpha0: [2.079, 3.524, 1.917, 1.992, 2.536]
       alpha1: [2.265, 3.862, 1.921, 1.685, 2.480]

skypy/galaxy/_schechter.py

@@ -5,20 +5,21 @@
 from ..utils import uses_default_cosmology
 from .redshift import schechter_lf_redshift
 from .luminosity import schechter_lf_magnitude
-
+from astropy import units
 
 __all__ = [
     'schechter_lf',
 ]
 
 
 @uses_default_cosmology
-def schechter_lf(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmology, noise=True):
+@units.quantity_input(sky_area=units.sr)
+def schechter_lf(redshift, M_star, phi_star, alpha, m_lim, sky_area, cosmology, noise=True):
     r'''Sample redshifts and magnitudes from a Schechter luminosity function.
 
     Sample the redshifts and magnitudes of galaxies following a Schechter
     luminosity function with potentially redshift-dependent parameters, limited
-    by an apparent magnitude `m_lim`, for a fraction `fsky` of the sky.
+    by an apparent magnitude `m_lim`, for a sky area `sky_area`.
 
     Parameters
     ----------
@@ -37,8 +38,8 @@ def schechter_lf(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmology, nois
         values for each `redshift`, or a function of redshift.
     m_lim : float
         Limiting apparent magnitude.
-    fsky : array_like
-        Sky fraction over which galaxies are sampled.
+    sky_area : `~astropy.units.Quantity`
+        Sky area over which galaxies are sampled. Must be in units of solid angle.
     cosmology : Cosmology, optional
         Cosmology object to convert apparent to absolute magnitudes. If not
         given, the default cosmology is used.
@@ -61,7 +62,7 @@ def schechter_lf(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmology, nois
     '''
 
     # sample galaxy redshifts
-    z = schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmology, noise)
+    z = schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, sky_area, cosmology, noise)
 
     # if a function is NOT given for M_star, phi_star, alpha, interpolate to z
     if not callable(M_star) and np.ndim(M_star) > 0:

skypy/galaxy/redshift.py

@@ -7,6 +7,7 @@
 import numpy as np
 import scipy.integrate
 import scipy.special
+from astropy import units
 
 from ..utils import uses_default_cosmology, broadcast_arguments, dependent_argument
 
@@ -77,12 +78,14 @@ def smail(z_median, alpha, beta, size=None):
 @dependent_argument('phi_star', 'redshift')
 @dependent_argument('alpha', 'redshift')
 @broadcast_arguments('redshift', 'M_star', 'phi_star', 'alpha')
-def schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmology, noise=True):
+@units.quantity_input(sky_area=units.sr)
+def schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, sky_area,
+                          cosmology, noise=True):
     r'''Sample redshifts from Schechter luminosity function.
 
     Sample the redshifts of galaxies following a Schechter luminosity function
     with potentially redshift-dependent parameters, limited by an apparent
-    magnitude `m_lim`, for a fraction `fsky` of the sky.
+    magnitude `m_lim`, for a sky area `sky_area`.
 
     Parameters
     ----------
@@ -101,8 +104,8 @@ def schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmol
         values for each `redshift`, or a function of redshift.
     m_lim : float
         Limiting apparent magnitude.
-    fsky : array_like
-        Sky fraction over which galaxies are sampled.
+    sky_area : `~astropy.units.Quantity`
+        Sky area over which galaxies are sampled. Must be in units of solid angle.
     cosmology : Cosmology, optional
         Cosmology object to convert apparent to absolute magnitudes. If not
         given, the default cosmology is used.
@@ -123,11 +126,13 @@ def schechter_lf_redshift(redshift, M_star, phi_star, alpha, m_lim, fsky, cosmol
     the sky.
 
     >>> from skypy.galaxy.redshift import schechter_lf_redshift
+    >>> from astropy import units
     >>> z = [0., 5.]
     >>> M_star = -20.5
     >>> phi_star = 3.5e-3
     >>> alpha = -1.3
-    >>> z_gal = schechter_lf_redshift(z, M_star, phi_star, alpha, 22, 1/41253)
+    >>> sky_area = 1*units.deg**2
+    >>> z_gal = schechter_lf_redshift(z, M_star, phi_star, alpha, 22, sky_area)
 
     '''
 
@@ -150,16 +155,17 @@ def f(lnx, a):
 
     # sample redshifts from the comoving density
     return redshifts_from_comoving_density(redshift=redshift, density=density,
-                                           fsky=fsky, cosmology=cosmology, noise=noise)
+                                           sky_area=sky_area, cosmology=cosmology, noise=noise)
 
 
 @uses_default_cosmology
-def redshifts_from_comoving_density(redshift, density, fsky, cosmology, noise=True):
+@units.quantity_input(sky_area=units.sr)
+def redshifts_from_comoving_density(redshift, density, sky_area, cosmology, noise=True):
     r'''Sample redshifts from a comoving density function.
 
     Sample galaxy redshifts such that the resulting distribution matches a past
     lightcone with comoving galaxy number density `density` at redshifts
-    `redshift`. The comoving volume sampled corresponds to a sky fraction `fsky`
+    `redshift`. The comoving volume sampled corresponds to a sky area `sky_area`
     and transverse comoving distance given by the cosmology `cosmology`.
 
     If the `noise` parameter is set to true, the number of galaxies has Poisson
@@ -171,8 +177,8 @@ def redshifts_from_comoving_density(redshift, density, fsky, cosmology, noise=Tr
         Redshifts at which comoving number densities are provided.
     density : array_like
         Comoving galaxy number density at each redshift in Mpc-3.
-    fsky : array_like
-        Sky fraction over which galaxies are sampled.
+    sky_area : `~astropy.units.Quantity`
+        Sky area over which galaxies are sampled. Must be in units of solid angle.
     cosmology : Cosmology, optional
         Cosmology object for conversion to comoving volume. If not given, the
         default cosmology is used.
@@ -191,15 +197,16 @@ def redshifts_from_comoving_density(redshift, density, fsky, cosmology, noise=Tr
     redshift 1 for a survey of 1 square degree = 1/41253 of the sky.
 
     >>> from skypy.galaxy.redshift import redshifts_from_comoving_density
+    >>> from astropy import units
     >>> z_range = np.arange(0, 1.01, 0.1)
-    >>> z_gal = redshifts_from_comoving_density(z_range, 1e-3, 1/41253)
+    >>> sky_area = 1*units.deg**2
+    >>> z_gal = redshifts_from_comoving_density(z_range, 1e-3, sky_area)
 
     '''
 
     # redshift number density
-    dN_dz = cosmology.differential_comoving_volume(redshift).to_value('Mpc3/sr')
+    dN_dz = (cosmology.differential_comoving_volume(redshift) * sky_area).to_value('Mpc3')
     dN_dz *= density
-    dN_dz *= 4*np.pi*fsky
 
     # integrate density to get expected number of galaxies
     N = np.trapz(dN_dz, redshift)

skypy/galaxy/tests/test_redshift.py

@@ -6,8 +6,9 @@
 @pytest.mark.flaky
 def test_schechter_lf_redshift():
 
-    from skypy.galaxy.redshift import schechter_lf_redshift, redshifts_from_comoving_density
+    from skypy.galaxy.redshift import schechter_lf_redshift
     from astropy.cosmology import FlatLambdaCDM
+    from astropy import units
     from scipy.special import gamma, gammaincc
 
     # fix this cosmology
@@ -19,10 +20,10 @@ def test_schechter_lf_redshift():
     phi_star = 1e-3
     alpha = -0.5
     m_lim = 30.
-    fsky = 1/41253
+    sky_area = 1.0 * units.deg**2
 
     # sample redshifts
-    z_gal = schechter_lf_redshift(z, M_star, phi_star, alpha, m_lim, fsky, cosmo, noise=False)
+    z_gal = schechter_lf_redshift(z, M_star, phi_star, alpha, m_lim, sky_area, cosmo, noise=False)
 
     # the absolute magnitude limit as function of redshift
     M_lim = m_lim - cosmo.distmod(z).value
@@ -34,7 +35,7 @@ def test_schechter_lf_redshift():
     density = phi_star*gamma(alpha+1)*gammaincc(alpha+1, x_min)
 
     # turn into galaxies/surface area
-    density *= 4*np.pi*fsky*cosmo.differential_comoving_volume(z).to_value('Mpc3/sr')
+    density *= (sky_area * cosmo.differential_comoving_volume(z)).to_value('Mpc3')
 
     # integrate total number
     n_gal = np.trapz(density, z, axis=-1)
@@ -58,6 +59,7 @@ def test_redshifts_from_comoving_density():
 
     from skypy.galaxy.redshift import redshifts_from_comoving_density
     from astropy.cosmology import LambdaCDM
+    from astropy import units
 
     # random cosmology
     H0 = np.random.uniform(50, 100)
@@ -69,10 +71,10 @@ def test_redshifts_from_comoving_density():
     Ngal = 1000
     redshift = np.arange(0.0, 2.001, 0.1)
     density = Ngal/cosmo.comoving_volume(redshift[-1]).to_value('Mpc3')
-    fsky = 1.0
+    sky_area = 4 * np.pi * units.steradian
 
     # sample galaxies over full sky without Poisson noise
-    z_gal = redshifts_from_comoving_density(redshift, density, fsky, cosmo, noise=False)
+    z_gal = redshifts_from_comoving_density(redshift, density, sky_area, cosmo, noise=False)
 
     # make sure number of galaxies is correct (no noise)
     assert np.isclose(len(z_gal), Ngal, atol=1, rtol=0)
@@ -82,6 +84,16 @@ def test_redshifts_from_comoving_density():
     D, p = kstest(z_gal, lambda z: cosmo.comoving_volume(z).to_value('Mpc3')/V_max)
     assert p > 0.01, 'D = {}, p = {}'.format(D, p)
 
+    # Test that an error is raised if sky_area has the wrong units
+    sky_area = 1 * units.degree
+    with pytest.raises(units.UnitsError):
+        redshifts_from_comoving_density(redshift, density, sky_area, cosmo, noise=False)
+
+    # Test that an error is raised if sky_area has no units
+    sky_area = 1
+    with pytest.raises(TypeError):
+        redshifts_from_comoving_density(redshift, density, sky_area, cosmo, noise=False)
+
 
 @pytest.mark.flaky
 def test_smail():

skypy/galaxy/tests/test_schechter.py

@@ -7,6 +7,7 @@ def test_schechter_lf():
 
     from pytest import raises
     from skypy.galaxy import schechter_lf
+    from astropy import units
 
     # redshift and magnitude distributions are tested separately
     # only test that output is consistent here
@@ -17,10 +18,10 @@ def test_schechter_lf():
     phi_star = 1e-3
     alpha = -0.5
     m_lim = 30.
-    fsky = 1/41253
+    sky_area = 1.0 * units.deg**2
 
     # sample redshifts and magnitudes
-    z_gal, M_gal = schechter_lf(z, M_star, phi_star, alpha, m_lim, fsky)
+    z_gal, M_gal = schechter_lf(z, M_star, phi_star, alpha, m_lim, sky_area)
 
     # check length
     assert len(z_gal) == len(M_gal)
@@ -31,4 +32,4 @@ def test_schechter_lf():
     # sample s.t. arrays need to be interpolated
     # alpha array not yet supported
     with raises(NotImplementedError):
-        z_gal, M_gal = schechter_lf(z, M_star, phi_star, alpha, m_lim, fsky)
+        z_gal, M_gal = schechter_lf(z, M_star, phi_star, alpha, m_lim, sky_area)