In [1]:
import numpy as np

In [2]:
import snsims
import healpy as hp

In [3]:
from astropy.cosmology import Planck15 as cosmo

In [7]:
help(snsims.PowerLawRates)


Help on class PowerLawRates in module snsims.paramDistribution:

class PowerLawRates(snsims.populationParamSamples.RateDistributions)
 |  This class is a concrete implementation of `RateDistributions` with the
 |  following properties:
 |  - The SN rate : The SN rate is a single power law with numerical
 |      coefficients (alpha, beta)  passed into the instantiation. The rate is
 |      the number of SN at redshift z per comoving volume per unit observer
 |      time over the entire sky expressed in units of numbers/Mpc^3/year 
 |  - A binning in redshift is used to perform the calculation of numbers of SN.
 |      This is assumed
 |  - The expected number of SN in each of these redshift bins is computed using
 |      the rate above, and a cosmology to compute the comoving volume for the
 |      redshift bin
 |  - The numbers of SN are determined by a Poisson Distribution about the
 |      expected number in each redshift bin,  determined with a random state
 |      passed in as an argument. This number must be integral.
 |  - It is assumed that the change of rates and volume within a redshift bin
 |      is negligible enough that samples to the true distribution may be drawn
 |      by obtaining number of SN samples of z from a uniform distribution
 |      within the z bin.
 |  
 |  Method resolution order:
 |      PowerLawRates
 |      snsims.populationParamSamples.RateDistributions
 |      __builtin__.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, rng, alpha=2.6e-05, beta=1.5, zbinEdges=None, zlower=1e-08, zhigher=1.4, numBins=20, surveyDuration=10.0, fieldArea=None, skyFraction=None, cosmo=FlatLambdaCDM(name="Planck15", H0=67.7 km / (Mpc...ff=3.05, m_nu=[ 0.    0.    0.06] eV, Ob0=0.0486))
 |      Parameters
 |      ----------
 |      rng : instance of `np.random.RandomState`
 |      cosmo : Instance of `astropy.cosmology` class, optional, defaults to Planck15
 |          data structure specifying the cosmological parameters 
 |      alpha : float, optional, defaults to 2.6e-5
 |      beta : float, optional, defaults to 1.5
 |          parameter in expression
 |  
 |  numSN(self)
 |      Return the number of expected supernovae in time DeltaT, with a rate snrate
 |      in a redshift range zlower, zhigher divided into numBins equal redshift
 |      bins. The variation of the rate within a bin is ignored.
 |      
 |      Parameters
 |      ----------
 |      zlower : mandatory, float
 |          lower limit on redshift range
 |      zhigher : mandatory, float
 |          upper limit on redshift range
 |      numBins : mandatory, integer
 |          number of bins
 |      cosmo : `astropy.cosmology` instance, mandatory
 |          cosmological parameters
 |      fieldArea : mandatory, units of radian square
 |          sky area considered
 |  
 |  snRate(self, z)
 |      The rate of SN at a redshift z in units of number of SN/ comoving
 |      volume in Mpc^3/yr in earth years according to the commonly used
 |      power-law expression 
 |      
 |      .. math:: rate(z) = lpha (h/0.7)^3 (1.0 + z)eta
 |      
 |      Parameters
 |      ----------
 |      z : array-like, mandatory 
 |          redshifts at which the rate is evaluated
 |      
 |      Examples
 |      --------
 |  
 |  zSampleSize(self)
 |      Parameters
 |      ----------
 |      zbinEdges : `nunpy.ndarray` of edges of zbins, defaults to None
 |          Should be of the form np.array([z0, z1, z2]) which will have
 |          zbins (z0, z1) and (z1, z2)
 |      skyFraction : np.float, optional, 
 |      fieldArea : optional, units of degrees squared
 |          area of sky considered.
 |      zlower : float, optional, defaults to None
 |          lower edge of z range
 |      zhigher : float, optional, defaults to None
 |          higher edge of z range
 |      numBins : int, optional, defaults to None
 |         if not None, overrides zbinEdges
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  randomState
 |  
 |  skyFraction
 |  
 |  zSamples
 |  
 |  zbinEdges
 |  
 |  ----------------------------------------------------------------------
 |  Data and other attributes defined here:
 |  
 |  __abstractmethods__ = frozenset([])
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors inherited from snsims.populationParamSamples.RateDistributions:
 |  
 |  __dict__
 |      dictionary for instance variables (if defined)
 |  
 |  __weakref__
 |      list of weak references to the object (if defined)


In [31]:
zdist = snsims.PowerLawRates(rng=np.random.RandomState(1), 
                             fieldArea=9.6,
                             surveyDuration=10.,
                             zbinEdges=np.arange(0.010001, 1.1, 0.1))

In [32]:
# ten years
zdist.DeltaT


Out[32]:
10.0

In [33]:
# The sky is >~ 40000 sq degrees ~ 4000 * LSST field of view
zdist.skyFraction * 2000 * 2


Out[33]:
0.93084226773030898

In [34]:
zdist.zbinEdges


Out[34]:
array([ 0.010001,  0.110001,  0.210001,  0.310001,  0.410001,  0.510001,
        0.610001,  0.710001,  0.810001,  0.910001,  1.010001])

In [36]:
zdist.zSampleSize().sum()


Out[36]:
11989.600801547122

In [37]:
zdist.zbinEdges


Out[37]:
array([ 0.010001,  0.110001,  0.210001,  0.310001,  0.410001,  0.510001,
        0.610001,  0.710001,  0.810001,  0.910001,  1.010001])

In [47]:
fig, ax = plt.subplots(1, 2)
_ = ax[0].hist(zdist.zSamples, bins=np.arange(0., 1.1, 0.1), histtype='step', lw=2)
_ = ax[0].errorbar(0.5*(zdist.zbinEdges[:-1]+zdist.zbinEdges[1:]), zdist.zSampleSize(),
                   yerr=np.sqrt(zdist.zSampleSize()), fmt='o')
_ = ax[1].plot(0.5*(zdist.zbinEdges[:-1]+zdist.zbinEdges[1:]), 
               18000*zdist.zSampleSize()/ zdist.fieldArea, 'ko')


On the left subplot, the histogram shows the samples of redshift actually simulated as a function of z. The green points are the number of samples that were expected from the rate, time, and area. The error bars show the square root of the number.

On the right subplot, we show the numbers expected in 18000 degrees of a sky from this redshift distribution