#Constants
H0=72.0 #km/s/Mpc
c = 3*10^5 #km/s


#Hubble's law
Hubble=function(z,Omegam) H0*sqrt(Omegam*(1+z)^3+(1.0-Omegam))


#H^(-1)
Hinv=function(zint,Omegam) 1.0/Hubble(zint,Omegam)


#Luminosity distance
dL=function(z,Omegami){
	dLsol=z
	for (i in (1:length(z))){
		zarg=z[i]
		I=integrate(Hinv,0.0,zarg,Omegam=Omegami)
		dLsol[i] = c*(1+zarg)*I$value
		}
	return(dLsol)
	}