function rd=rd(x,y,z)
% RD   Carlson's elliptic integral of the second kind.
% usage:: Rd = rd(x,y,z)
%   rd(x,y,z) returns the value of Carlson's elliptic integral of the
%   second kind, RD(x; y; z). x and y must be nonnegative, and at most one
%   can be zero. z must be positive. TINY must be at least twice the
%   negative 2/3 power of the machine overflow limit. BIG must be at most
%   0:1  ERRTOL times the negative 2/3 power of the machine underflow
%   limit.  These condition are not checked.
%
%   RD is related to the Legendre ellipitic integral of the second
%   kind:
%   E(phi,m)=sin(phi)*RF(cos(phi)^2,1-m*sin(phi)^2,1)
%    -1/3*m*sin(phi)^3*RD(cos(phi)^2,1-m*sin(phi)^2,1);


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

tol=0.05;
%TINY=1.0e-25;
%BIG=4.5e21;
C1=3/14;
C2=1/6;
C3=9/22;
C4=3/26;
C5=C3/4;
C6=1.5*C4;

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

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

sum=0;
fac=1;
mm=1;
while (mm>tol)
  sqrtx=sqrt(xt);
  sqrty=sqrt(yt);
  sqrtz=sqrt(zt);
  alamb=sqrtx.*(sqrty+sqrtz)+sqrty.*sqrtz;
  sum=sum+fac./(sqrtz.*(zt+alamb));
  fac=0.25*fac;
  xt=0.25*(xt+alamb);
  yt=0.25*(yt+alamb);
  zt=0.25*(zt+alamb);
  ave=0.2*(xt+yt+3.0*zt);
  delx=(ave-xt)./ave;
  dely=(ave-yt)./ave;
  delz=(ave-zt)./ave;
  mm=sqrt(max(flat(([delx(:) dely(:) delz(:)]).^2)));
end
ea=delx.*dely;
eb=delz.*delz;
ec=ea-eb;
ed=ea-6*eb;
ee=ed+ec+ec;
rd=3*sum+fac.*(1+ed.*(-C1+C5*ed-C6*delz.*ee)...
+delz.*(C2.*ee+delz.*(-C3*ec+C4*delz.*ea)))./(ave.*sqrt(ave));

rd=reshape(rd,nn(:,pp)');