Code covered by the BSD License  

Highlights from
slatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[cotresult,x]=cot(x);
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

Contact us at files@mathworks.com