'Comoving volume in cubic Mpc at redshift ``z``.
This is the volume of the universe encompassed by redshifts less
than ``z``. For the case of omega_k = 0 it is a sphere of radius
`comoving_distance` but it is less intuitive
if omega_k is not 0.
Parameters
z : array-like
Input redshifts.  Must be 1D or scalar.
Returns
V : `~astropy.units.Quantity`
Comoving volume in :math:`Mpc^3` at each input redshift.'

def comoving_volume(self, z):

    Ok0 = self._Ok0
    if (Ok0 == 0):
        return (((4.0 / 3.0) * pi) * (self.comoving_distance(z) ** 3))
    dh = self._hubble_distance.value
    dm = self.comoving_transverse_distance(z).value
    term1 = ((((4.0 * pi) * (dh ** 3)) / (2.0 * Ok0)) * (u.Mpc ** 3))
    term2 = ((dm / dh) * np.sqrt((1 + (Ok0 * ((dm / dh) ** 2)))))
    term3 = ((sqrt(abs(Ok0)) * dm) / dh)
    if (Ok0 > 0):
        return (term1 * (term2 - ((1.0 / sqrt(abs(Ok0))) * np.arcsinh(term3))))
    else:
        return (term1 * (term2 - ((1.0 / sqrt(abs(Ok0))) * np.arcsin(term3))))
