Units

(* km in metres *)
km=(10^3);
(* AU in meters (IAU 2012 definition) *)
AU=149597870700;
(* parsec in metres *)
pc=((60 60 180)/π)AU;
(* Megaparsec *)
Mpc=(10^6)pc;
(* year in seconds *)
yr=3600 24 365.25;
(* speed of light in m/s *)
c=299792458;
(* H0 in (s^(-1)) *)
H0=70km/Mpc;

(* 1/H0 in yr *)
(1/(H0 yr))

Distances

Clear[dL,xpar,x];
(* xpar = geodesic comoving distance in units c/H0 *)
xpar[Ωr_,Ωm_,Ωk_,ΩΛ_,z_]:=NIntegrate[(1/((a^2)(Sqrt[Ωr (a^(-4))+Ωm (a^(-3))+Ωk (a^(-2))+ΩΛ]))),{a,(1/(1+z)),1}];
(* x = circumferential comoving distance in units c/H0 *)
x[Ωr_,Ωm_,Ωk_,ΩΛ_,z_]:=(Piecewise[{{(1/(Sqrt[Ωk]))Sinh[(Sqrt[Ωk])xpar[Ωr,Ωm,Ωk,ΩΛ,z]], Ωk>0 (* open *)}, {xpar[Ωr,Ωm,Ωk,ΩΛ,z], Ωk==0 (* flat *)}, {(1/(Sqrt[-Ωk]))Sin[(Sqrt[-Ωk])xpar[Ωr,Ωm,Ωk,ΩΛ,z]], Ωk<0 (* closed *)}}]);
(* dL = luminosity distance in units c/H0 *)
dL[Ωr_,Ωm_,Ωk_,ΩΛ_,z_]:=(1+z)x[Ωr,Ωm,Ωk,ΩΛ,z]

Theoretical Hubble diagram

(* Hubble diagram for various cosmological parameters *)
theoryplot=LogLogPlot[{
(* ΩΛ = 1 *)dL[0,0,0,1,z],
(* Ωk = 1 *)dL[0,0,1,0,z],
(* Ωm = .3, ΩΛ = .7 *)dL[0,.3,0,.7,z],
(* Ωm = .3, Ωk = .7 *)dL[0,.3,.7,0,z],
(* Ωm = 1 *)dL[0,1,0,0,z]
},{z,.01,2}]

(* Hubble diagram for various cosmological parameters, relative to the Ωm = .3, ΩΛ = .7 model *)
theoryplotrel=LogLogPlot[{
(* ΩΛ = 1 *)dL[0,0,0,1,z]/dL[0,.3,0,.7,z],
(* Ωk = 1 *)dL[0,0,1,0,z]/dL[0,.3,0,.7,z],
(* Ωm = .3, ΩΛ = .7 *)dL[0,.3,0,.7,z]/dL[0,.3,0,.7,z],
(* Ωm = .3, Ωk = .7 *)dL[0,.3,.7,0,z]/dL[0,.3,0,.7,z],
(* Ωm = 1 *)dL[0,1,0,0,z]/dL[0,.3,0,.7,z]
},{z,.01,2}]

Supernova data

Betoule et al 2014

(* read data from file *)
Clear[file,sndataBetoule];
file=OpenRead["betoule_sne_14.txt"];
(* skip first 77 lines of file *)
Skip[file,Record,77];
(* read data for 740 supernovae *)
sndataBetoule=ReadList[file,Record,740];
Close[file];

(* data for i'th supernova *)
Clear[sn,zcmb,mB,shape,color,logMstellar,MB,μ,d];
sn[i_]:=ReadList[StringToStream[sndataBetoule[[i]]],{Word,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Number,Word,Word}][[1]];
(* redshift in CMB frame *)
zcmb[i_]:=sn[i][[2]];
(* apparent B magnitude *)
mB[i_]:=sn[i][[5]];
(* shape parameter of supernova lightcurve *)
shape[i_]:=sn[i][[7]];
(* color parameter of supernova lightcurve *)
color[i_]:=sn[i][[9]];
(* log10 of the mass of the host galaxy in solar masses *)
logMstellar[i_]:=sn[i][[11]];
MB[i_]:=(Piecewise[{{-19.05, logMstellar[i]<10}, {-19.12, logMstellar[i]≥10}}]);
(* corrected distance modulus *)
μ[i_]:=mB[i]+α shape[i]-β color[i]-MB[i];
(* luminosity distance in units of c/H0 *) 
d[i_]:=((10 pc)/(c/H0))(10^(μ[i]/5));

(* Hubble diagram of the data *)
dataplot=Block[{α=0.141,β=3.101},ListLogLogPlot[Table[{zcmb[i],d[i]},{i,1,Length[sndataBetoule]}]]]

(* Hubble diagram of the data relative to the standard ΛCDM model *)
dataplotrel=Block[{α=0.141,β=3.101},ListLogLogPlot[Table[{zcmb[i],d[i]/dL[0,.3,0,.7,zcmb[i]]},{i,1,Length[sndataBetoule]}]]]

Scolnic et al 2018

(* read data from file *)
Clear[file,sndataScolnic];
file=OpenRead["scolnic/hlsp_ps1cosmo_panstarrs_gpc1_all_model_v1_lcparam-full.txt"];
(* skip first 3 lines of file *)
Skip[file,Record,3];
(* read data for 1048 supernovae *)
sndataScolnic=ReadList[file,Record,1048];
Close[file];

(* data for i'th supernova *)
Clear[sn,zcmb,mB,dmB,d];
sn[i_]:=ReadList[StringToStream[sndataScolnic[[i]]],{Word,Number,Number,Number,Number,Number}][[1]];
(* redshift in CMB frame *)
zcmb[i_]:=sn[i][[2]];
(* apparent B magnitude *)
mB[i_]:=sn[i][[5]];
(* uncertainty in apparent B magnitude *)
dmB[i_]:=sn[i][[7]];
(* corrected distance modulus *)
μ[i_]:=mB[i]+19.33;
(* luminosity distance in units of c/H0 *) 
d[i_]:=((10 pc)/(c/H0))(10^(μ[i]/5));

(* Hubble diagram of the data *)
dataplot=ListLogLogPlot[Table[{zcmb[i],d[i]},{i,1,Length[sndataScolnic]}]]

(* Hubble diagram of the data relative to the standard ΛCDM model *)
dataplotrel=ListLogLogPlot[Table[{zcmb[i],d[i]/dL[0,.3,0,.7,zcmb[i]]},{i,1,Length[sndataScolnic]}]]

Hubble diagram, theory and data

(* Hubble diagram *)
Show[theoryplot,dataplot,ImageSize->Large]

(* Hubble diagram relative to the Ωm = .3, ΩΛ = .7 model *)
Show[theoryplotrel,dataplotrel,ImageSize->Large]