[catanresult,z]=catan(z);

function [catanresult,z]=catan(z);
catanresult=[];
persistent first firstCall i nterms pi2 r r2 rmax rmin sqeps twoi x xans y yans z2 ; if isempty(firstCall),firstCall=1;end; 

;
if isempty(i), i=0; end;
if isempty(nterms), nterms=0; end;
if isempty(pi2), pi2=0; end;
if isempty(r), r=0; end;
if isempty(r2), r2=0; end;
if isempty(rmax), rmax=0; end;
if isempty(rmin), rmin=0; end;
if isempty(sqeps), sqeps=0; end;
if isempty(twoi), twoi=0; end;
if isempty(x), x=0; end;
if isempty(xans), xans=0; end;
if isempty(y), y=0; end;
if isempty(yans), yans=0; end;
%***BEGIN PROLOGUE  CATAN
%***PURPOSE  Compute the complex arc tangent.
%***LIBRARY   SLATEC (FNLIB)
%***CATEGORY  C4A
%***TYPE      COMPLEX (CATAN-C)
%***KEYWORDS  ARC TANGENT, ELEMENTARY FUNCTIONS, FNLIB, TRIGONOMETRIC
%***AUTHOR  Fullerton, W., (LANL)
%***DESCRIPTION
%
% CATAN(Z) calculates the complex trigonometric arc tangent of Z.
% The result is in units of radians, and the real part is in the first
% or fourth quadrant.
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  R1MACH, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   770801  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%***end PROLOGUE  CATAN
if isempty(z2), z2=0; end;
if isempty(first), first=false; end;
if firstCall,   pi2=[1.57079632679489661923e0];  end;
if firstCall,   first=[true];  end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT  CATAN
if( first )
% NTERMS = LOG(EPS)/LOG(RBND) WHERE RBND = 0.1
nterms = fix(-0.4343.*log(r1mach(3)) + 1.0);
sqeps = sqrt(r1mach(4));
rmin = sqrt(3.0.*r1mach(3));
rmax = 1.0./r1mach(3);
end;
first = false;
%
r = abs(z);
if( r<=0.1 )
%
catanresult = z;
if( r<rmin )
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',z); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
catanresult =complex(0.0,0.0);
z2 = z.*z;
for i = 1 : nterms;
twoi = 2.*(nterms-i) + 1;
catanresult = 1.0./twoi - z2.*catanresult;
end; i = fix(nterms+1);
catanresult = z.*catanresult;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',z); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
%
elseif( r>rmax ) ;
%
catanresult = complex(pi2,0.);
if( real(z)<0.0 )
catanresult = complex(-pi2,0.0);
end;
else;
x = real(z);
y = imag(z);
r2 = r.*r;
if( r2==1.0 && x==0.0 )
xermsg('SLATEC','CATAN','Z IS +I OR -I',2,2);
end;
if( abs(r2-1.0)<=sqeps )
if( abs(complex(1.0,0.0)+z.*z)<sqeps )
xermsg('SLATEC','CATAN','ANSWER LT HALF PRECISION, Z**2 CLOSE TO -1',1,1);
end;
end;
%
xans = 0.5.*atan2(2.0.*x,1.0-r2);
yans = 0.25.*log((r2+2.0.*y+1.0)./(r2-2.0.*y+1.0));
catanresult = complex(xans,yans);
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',z); evalin('caller',[inputname(1),'=FUntemp;']); end
return;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',z); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK CATANH