function [cotresult,x]=cot(x);
cotresult=[];
persistent ainty ainty2 cot cotcs first firstCall ifn nterms pi2rec prodbg sqeps xmax xmin xsml y yrem ; if isempty(firstCall),firstCall=1;end;
if isempty(ainty), ainty=0; end;
if isempty(ainty2), ainty2=0; end;
if isempty(cotresult), cotresult=0; end;
if isempty(cotcs), cotcs=zeros(1,8); end;
if isempty(pi2rec), pi2rec=0; end;
if isempty(prodbg), prodbg=0; end;
if isempty(sqeps), sqeps=0; end;
if isempty(xmax), xmax=0; end;
if isempty(xmin), xmin=0; end;
if isempty(xsml), xsml=0; end;
if isempty(y), y=0; end;
if isempty(yrem), yrem=0; end;
if isempty(ifn), ifn=0; end;
if isempty(nterms), nterms=0; end;
%***BEGIN PROLOGUE COT
%***PURPOSE Compute the cotangent.
%***LIBRARY SLATEC (FNLIB)
%***CATEGORY C4A
%***TYPE SINGLE PRECISION (COT-S, DCOT-D, CCOT-C)
%***KEYWORDS COTANGENT, ELEMENTARY FUNCTIONS, FNLIB, TRIGONOMETRIC
%***AUTHOR Fullerton, W., (LANL)
%***DESCRIPTION
%
% COT(X) calculates the cotangent of the real argument X. X is in
% units of radians.
%
% Series for COT on the interval 0. to 6.25000D-02
% with weighted error 3.76E-17
% log weighted error 16.42
% significant figures required 15.51
% decimal places required 16.88
%
%***REFERENCES (NONE)
%***ROUTINES CALLED CSEVL, INITS, R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 770601 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)
% 920618 Removed space from variable names. (RWC, WRB)
%***end PROLOGUE COT
if isempty(first), first=false; end;
if firstCall, cotcs(1)=[.24025916098295630e0]; end;
if firstCall, cotcs(2)=[-.016533031601500228e0]; end;
if firstCall, cotcs(3)=[-.000042998391931724e0]; end;
if firstCall, cotcs(4)=[-.000000159283223327e0]; end;
if firstCall, cotcs(5)=[-.000000000619109313e0]; end;
if firstCall, cotcs(6)=[-.000000000002430197e0]; end;
if firstCall, cotcs(7)=[-.000000000000009560e0]; end;
if firstCall, cotcs(8)=[-.000000000000000037e0]; end;
if firstCall, pi2rec=[.0116197723675813430e0]; end;
if firstCall, first=[true]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT COT
if( first )
[nterms ,cotcs]=inits(cotcs,8,0.1.*r1mach(3));
xmax = 1.0./r1mach(4);
xsml = sqrt(3.0.*r1mach(3));
xmin = exp(max(log(r1mach(1)),-log(r1mach(2)))+0.01);
sqeps = sqrt(r1mach(4));
end;
first = false;
%
y = abs(x);
if( abs(x)<xmin )
xermsg('SLATEC','COT','ABS(X) IS ZERO OR SO SMALL COT OVERFLOWS',2,2);
end;
if( y>xmax )
xermsg('SLATEC','COT','NO PRECISION BECAUSE ABS(X) IS TOO BIG',3,2);
end;
%
% CAREFULLY COMPUTE Y * (2/PI) = (AINT(Y) + REM(Y)) * (.625 + PI2REC)
% = AINT(.625*Y) + REM(.625*Y) + Y*PI2REC = AINT(.625*Y) + Z
% = AINT(.625*Y) + AINT(Z) + REM(Z)
%
ainty = fix(y);
yrem = y - ainty;
prodbg = 0.625.*ainty;
ainty = fix(prodbg);
y =(prodbg-ainty) + 0.625.*yrem + y.*pi2rec;
ainty2 = fix(y);
ainty = ainty + ainty2;
y = y - ainty2;
%
ifn = fix(rem(ainty,2.));
if( ifn==1 )
y = 1.0 - y;
end;
%
if( abs(x)>0.5 && y<abs(x).*sqeps )
xermsg('SLATEC','COT',['ANSWER LT HALF PRECISION, ABS(X) TOO BIG OR X NEAR N*PI ','(N.NE.0)'],1,1);
end;
%
if( y<=0.25 )
cotresult = 1.0./x;
if( y>xsml )
cotresult =(0.5+csevl(32.0.*y.*y-1.,cotcs,nterms))./y;
end;
%
elseif( y>0.5 ) ;
%
cotresult =(0.5+csevl(2.0.*y.*y-1.,cotcs,nterms))./(0.25.*y);
cotresult =(cotresult.^2-1.0).*0.5./cotresult;
cotresult =(cotresult.^2-1.0).*0.5./cotresult;
else;
cotresult =(0.5+csevl(8.0.*y.*y-1.,cotcs,nterms))./(0.5.*y);
cotresult =(cotresult.^2-1.0).*0.5./cotresult;
end;
%
if( x~=0. )
cotresult = (abs(cotresult).*sign(x));
end;
if( ifn==1 )
cotresult = -cotresult;
end;
%
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
end
%DECK CPADD