set up an array of z-values, with something like zarr=0.001*(findgen(10000)+1.) Then I did the following (note that this should include multiplying by the factor lH/c for physical length l, Hubble constant H; that amounts to 0.052 for l=1 kpc, H=75 km/s Mpc) plot,zarr,hqfun(zarr,0.5),/xlog,/ylog,xtitle="Redshift z",ytitle="Size, arcseconds",yrange=[1,100] oplot,zarr,hqfun(zarr,0.01) oplot,zarr,hqfun(zarr,1.0) xyouts,5.5,1.8,'0.01' xyouts,6.,4,'0.5' xyouts,3.,10,'q0=1.0' xyouts,0.12,40,'Angular extent of 1 kpc' xyouts,0.12,32,'Ho=75 km/s Mpc' after defining this function to do all the work: function hqfun,z,q top=q^2 * (1+z)^2 bottom=q*z+((q-1.)*(sqrt(1.+2.*q*z)-1.)) return, top/bottom end