function M = radecz2xxyyzz (radecz)

% RADECZ2XYZ converts ra,dec,z to x-y-z positions
%
% 
% M = radecz2M (radecz)
%
% radecz is a N x L matrix representing N data points, 
%   with L values per point.
%
% (Alternatively, radecz is a string e.g. 'blah.dat' which
%  can be loaded into a N x L matrix.)
%
%  Each row represents a data point.
%  Column 1 has ra  (in degrees, between 0 and 360)
%  Column 2 has declination (in degrees, between -90 and 90)
%  Column 3 has redshift
%  radecz may have other columns (4 to L)
%
% M is a N x (L+1) matrix with rows corresponding to those of radecz.
%  Column 1 has x
%  Column 2 has y
%  Column 3 has z
%  Columns 4 to L (iff L > 3) are identical to those of radecz
%  Column L+1 has distance from origin, i.e. sqrt(x*x + y*y + z*z)
%  Column L+2 has redshift (column 3 of original matrix)
%
% Vectorized version - under 10 seconds even for 250K galaxies from Sloan DR2
%
% The ra/dec/r -> x/y/z conversion used is as below:
%
%  x=r*sin(90-dec)*cos(ra);
%  y=r*sin(90-dec)*sin(ra);
%  z=r*cos(90-dec) ;
%
% NOTE : there are other variants of this formula!  
%         (I used it because Mark Subbarao did when he was getting his data
%          from Sloan in an earlier incarnation of this project.)
%
% Dinoj Surendran dinoj@cs.uchicago.edu, 29 Feb 2004, Modified 3 Mar 04


if ischar (radecz)
  A = load (radecz);
else
  A = radecz;
end
N = size(A,1);
L = size(A,2);

% convert degrees -> radians (and dec->90-dec)

ra = A(:,1);
dec = A(:,2);
redshift = A(:,3);
dist = redshift2dist (redshift);

dc = (90-dec)*(3.14159/180);
ra  = ra * (3.14159/180);

M(:,1) = dist.*sin(dc).*cos(ra);
M(:,2) = dist.*sin(dc).*sin(ra);
M(:,3) = dist.*cos(dc) ;
M(:,4:L) = A(:,4:L);
M(:,L+1) = dist;
M(:,L+2) = redshift;