function rf=rf(x,y,z)
% RF   Carlson's elliptic integral of the first kind.
% usage:: Rf = RF(x,y,z)
%   RF(x,y,z) returns the value of Carlson's elliptic integral of
%   the first kind evaluated for each element of x,y,z.  As currently
%   implemented, x, y, and z should be the same size or constants,
%   x,y,z should be non-negitive and no more than one can be zero,
%   Also they cannot be too small or too large (within factor of 5 of
%   realmin/max).  These condition are not checked
%
%   RF is related to the Legendre ellipitic integral of the first
%   kind F(phi,m)=sin(phi)*RF(cos(phi)^2,1-m*sin(phi)^2,1);
%

% revision history:
% 07/18/01 Mark D. Shattuck <mds> rf.m
%          translated from Numerical Recipies in C

tol=0.0025;
%TINY=10*realmin;
%BIG=realmin/10;
C1=1/24;
C2=0.1;
C3=3/44;
C4=1/14;

%	if (min(min(x,y),z) < 0.0 || FMIN(FMIN(x+y,x+z),y+z) < TINY ||
%		FMAX(FMAX(x,y),z) > BIG)
%			nrerror("invalid arguments in rf");

nn=zeros(2,3);

nn(:,1)=size(x)';
nn(:,2)=size(y)';
nn(:,3)=size(z)';
NN=prod(nn);
[N pp]=max(NN);
if (NN(1)==1)
  xt=x*ones(1,N)';
else
  xt=x(:);
end
if (NN(2)==1)
  yt=y*ones(1,N)';
else
  yt=y(:);
end
if (NN(3)==1)
  zt=z*ones(1,N)';
else
  zt=z(:);
end

mm=1;
while mm>tol
  sqrtx=sqrt(xt);
  sqrty=sqrt(yt);
  sqrtz=sqrt(zt);
  alamb=sqrtx.*(sqrty+sqrtz)+sqrty.*sqrtz;
  xt=0.25*(xt+alamb);
  yt=0.25*(yt+alamb);
  zt=0.25*(zt+alamb);
  ave=(xt+yt+zt)/3;
  delx=(ave-xt)./ave;
  dely=(ave-yt)./ave;
  delz=(ave-zt)./ave;
  mm=sqrt(max(flat(([delx(:) dely(:) delz(:)]).^2)));
end
e2=delx.*dely-delz.*delz;
e3=delx.*dely.*delz;
rf=reshape((1+(C1*e2-C2-C3*e3).*e2+C4*e3)./sqrt(ave),nn(:,pp)');
