Return values of calc_vol

If z is a scalar, result is a scalar double precision $dV/(dz d\Omega)$. If z is an array, result is a double precision array with $dV/(dz d\Omega)$ for each z value.
If H₀ is not specified, it is internally set to 50.
The function calc_vol returns -1000 if the passed parameters describe a bouncing universe (see section V below).