function ds = redshift2dist (zs)

% converts redshift to megaparsecs, in comoving radial distance
%
% Based on the formula below (emailed to me by Andrey Zhiglo)
%
% r(z) = 299792.458/72.0*\int_0^z \frac{dw}{\sqrt{0.27*(1+w)^3+0.73}}
%
% zs is a Nx1 vector of redshift values.
% ds is a Nx1 vector of corresponding  distance values.
%
%
% Overestimates answer, e.g. gives 14076  instead of 13997   for z=1100.
%                            gives  9664.4 instead of 9659.7 for z=10
%                            gives  7934.9 instead of 7932.3 for z=5
%                            gives  5245.2 instead of 5244.3 for z=2
%                            gives  3317.4 instead of 3317.1 for z=1
% That is in comparison with http://www.astro.ucla.edu/~wright/CosmoCalc.html
% 
% This is ok for outreach purposes, visualization, and many scientific tasks. 
%
% Dinoj Surendran dinoj@cs.uchicago.edu
%
% The code is vectorized if all redshifts in zs are under 20 ...
%  the furthest quasar is still under z<10.

HUBBLE = 71;           % km/s/Mpc
LIGHT  = 299792.458;   % km/s
DARKENERGY = 0.73;     % fraction of universe that is dark energy
MATTER = 1 - DARKENERGY;

if size(zs,2) ~= 1
  zs = zs';
  if size(zs,2) ~= 1
    error ('input should be a vector');
  end
end

% make it work as z changes by 0.0001 

if length (find(zs>0)) == 0
  ds = zeros(size(zs));
  return;
end

if max(zs) < 20
  NUMINTEGPOINTS = max( 50000, ceil (2/min(zs(zs>0))));
  DEL = 1 / NUMINTEGPOINTS;                    % old - ignore -> DEL = (z / NUMINTEGPOINTS);
  xs = [0 : DEL : max(zs)+2*DEL];
  os = ones(size(xs));
  ys = (MATTER * (os+xs).^3 + os*DARKENERGY).^-0.5;
 cys = cumsum(ys);
 cys = cys';

  is = floor (zs * NUMINTEGPOINTS);    % z(i) occurs between xs(is(i)) and xs(is(i)+1) 
  ds = cys (is) + (cys(is+1)-cys(is)).*mod(is,1);
  ds = ds * DEL * LIGHT / HUBBLE;

else

  NUMINTEGPOINTS = 10000;
  for j = 1 : length(zs)
    DEL = (zs(j) / NUMINTEGPOINTS);
    xs = [0:DEL:zs(j)];
    os = ones(size(xs));
    ys = (MATTER * (os+xs).^3 + os*DARKENERGY).^-0.5;
    ds(j) = trapz (ys) * DEL * LIGHT / HUBBLE;
  end
  ds = ds';
end