| [rcresult,x,y,ier]=rc(x,y,ier); |
function [rcresult,x,y,ier]=rc(x,y,ier);
rcresult=[];
persistent c1 c2 errtol first firstCall lamda lolim mu s sn uplim xern3 xern4 xern5 xn yn ; if isempty(firstCall),firstCall=1;end;
;
%***BEGIN PROLOGUE RC
%***PURPOSE Calculate an approximation to
% RC(X,Y) = Integral from zero to infinity of
% -1/2 -1
% (1/2,t+X) (t+Y) dt,
% where X is nonnegative and Y is positive.
%***LIBRARY SLATEC
%***CATEGORY C14
%***TYPE SINGLE PRECISION (RC-S, DRC-D)
%***KEYWORDS DUPLICATION THEOREM, ELEMENTARY FUNCTIONS,
% ELLIPTIC INTEGRAL, TAYLOR SERIES
%***AUTHOR Carlson, B. C.
% Ames Laboratory-DOE
% Iowa State University
% Ames, IA 50011
% Notis, E. M.
% Ames Laboratory-DOE
% Iowa State University
% Ames, IA 50011
% Pexton, R. L.
% Lawrence Livermore National Laboratory
% Livermore, CA 94550
%***DESCRIPTION
%
% 1. RC
% Standard FORTRAN function routine
% Single precision version
% The routine calculates an approximation to
% RC(X,Y) = Integral from zero to infinity of
%
% -1/2 -1
% (1/2,t+X) (t+Y) dt,
%
% where X is nonnegative and Y is positive. The duplication
% theorem is iterated until the variables are nearly equal,
% and the function is then expanded in Taylor series to fifth
% order. Logarithmic, inverse circular, and inverse hyper-
% bolic functions can be expressed in terms of RC.
%
%
% 2. Calling Sequence
% RC( X, Y, IER )
%
% Parameters on Entry
% Values assigned by the calling routine
%
% X - Single precision, nonnegative variable
%
% Y - Single precision, positive variable
%
%
%
% On Return (values assigned by the RC routine)
%
% RC - Single precision approximation to the integral
%
% IER - Integer to indicate normal or abnormal termination.
%
% IER = 0 Normal and reliable termination of the
% routine. It is assumed that the requested
% accuracy has been achieved.
%
% IER > 0 Abnormal termination of the routine
%
% X and Y are unaltered.
%
%
% 3. Error Messages
%
% Value of IER assigned by the RC routine
%
% Value Assigned Error Message Printed
% IER = 1 X.LT.0.0E0.OR.Y.LE.0.0E0
% = 2 X+Y.LT.LOLIM
% = 3 MAX(X,Y) .GT. UPLIM
%
%
% 4. Control Parameters
%
% Values of LOLIM, UPLIM, and ERRTOL are set by the
% routine.
%
% LOLIM and UPLIM determine the valid range of X and Y
%
% LOLIM - Lower limit of valid arguments
%
% Not less than 5 * (machine minimum) .
%
% UPLIM - Upper limit of valid arguments
%
% Not greater than (machine maximum) / 5 .
%
%
% Acceptable values for: LOLIM UPLIM
% IBM 360/370 SERIES : 3.0E-78 1.0E+75
% CDC 6000/7000 SERIES : 1.0E-292 1.0E+321
% UNIVAC 1100 SERIES : 1.0E-37 1.0E+37
% CRAY : 2.3E-2466 1.09E+2465
% VAX 11 SERIES : 1.5E-38 3.0E+37
%
% ERRTOL determines the accuracy of the answer
%
% The value assigned by the routine will result
% in solution precision within 1-2 decimals of
% 'machine precision'.
%
%
% ERRTOL - Relative error due to truncation is less than
% 16 * ERRTOL ** 6 / (1 - 2 * ERRTOL).
%
%
% The accuracy of the computed approximation to the inte-
% gral can be controlled by choosing the value of ERRTOL.
% Truncation of a Taylor series after terms of fifth order
% introduces an error less than the amount shown in the
% second column of the following table for each value of
% ERRTOL in the first column. In addition to the trunca-
% tion error there will be round-off error, but in prac-
% tice the total error from both sources is usually less
% than the amount given in the table.
%
%
%
% Sample Choices: ERRTOL Relative Truncation
% error less than
% 1.0E-3 2.0E-17
% 3.0E-3 2.0E-14
% 1.0E-2 2.0E-11
% 3.0E-2 2.0E-8
% 1.0E-1 2.0E-5
%
%
% Decreasing ERRTOL by a factor of 10 yields six more
% decimal digits of accuracy at the expense of one or
% two more iterations of the duplication theorem.
%
% *Long Description:
%
% RC Special Comments
%
%
%
%
% Check: RC(X,X+Z) + RC(Y,Y+Z) = RC(0,Z)
%
% where X, Y, and Z are positive and X * Y = Z * Z
%
%
% On Input:
%
% X and Y are the variables in the integral RC(X,Y).
%
% On Output:
%
% X and Y are unaltered.
%
%
%
% RC(0,1/4)=RC(1/16,1/8)=PI=3.14159...
%
% RC(9/4,2)=LN(2)
%
%
%
% ********************************************************
%
% Warning: Changes in the program may improve speed at the
% expense of robustness.
%
%
% --------------------------------------------------------------------
%
% Special Functions via RC
%
%
%
% LN X X .GT. 0
%
% 2
% LN(X) = (X-1) RC(((1+X)/2) , X )
%
%
% --------------------------------------------------------------------
%
% ARCSIN X -1 .LE. X .LE. 1
%
% 2
% ARCSIN X = X RC (1-X ,1 )
%
% --------------------------------------------------------------------
%
% ARCCOS X 0 .LE. X .LE. 1
%
%
% 2 2
% ARCCOS X = SQRT(1-X ) RC(X ,1 )
%
% --------------------------------------------------------------------
%
% ARCTAN X -INF .LT. X .LT. +INF
%
% 2
% ARCTAN X = X RC(1,1+X )
%
% --------------------------------------------------------------------
%
% ARCCOT X 0 .LE. X .LT. INF
%
% 2 2
% ARCCOT X = RC(X ,X +1 )
%
% --------------------------------------------------------------------
%
% ARCSINH X -INF .LT. X .LT. +INF
%
% 2
% ARCSINH X = X RC(1+X ,1 )
%
% --------------------------------------------------------------------
%
% ARCCOSH X X .GE. 1
%
% 2 2
% ARCCOSH X = SQRT(X -1) RC(X ,1 )
%
% --------------------------------------------------------------------
%
% ARCTANH X -1 .LT. X .LT. 1
%
% 2
% ARCTANH X = X RC(1,1-X )
%
% --------------------------------------------------------------------
%
% ARCCOTH X X .GT. 1
%
% 2 2
% ARCCOTH X = RC(X ,X -1 )
%
% --------------------------------------------------------------------
%
%***REFERENCES B. C. Carlson and E. M. Notis, Algorithms for incomplete
% elliptic integrals, ACM Transactions on Mathematical
% Software 7, 3 (September 1981), pp. 398-403.
% B. C. Carlson, Computing elliptic integrals by
% duplication, Numerische Mathematik 33, (1979),
% pp. 1-16.
% B. C. Carlson, Elliptic integrals of the first kind,
% SIAM Journal of Mathematical Analysis 8, (1977),
% pp. 231-242.
%***ROUTINES CALLED R1MACH, XERMSG
%***REVISION HISTORY (YYMMDD)
% 790801 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 891009 Removed unreferenced statement labels. (WRB)
% 891009 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)
% 900510 Changed calls to XERMSG to standard form, and some
% editorial changes. (RWC))
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE RC
if isempty(xern3), xern3=repmat(' ',1,16); end;
if isempty(xern4), xern4=repmat(' ',1,16); end;
if isempty(xern5), xern5=repmat(' ',1,16); end;
if isempty(c1), c1=0; end;
if isempty(c2), c2=0; end;
if isempty(errtol), errtol=0; end;
if isempty(lamda), lamda=0; end;
if isempty(lolim), lolim=0; end;
if isempty(mu), mu=0; end;
if isempty(s), s=0; end;
if isempty(sn), sn=0; end;
if isempty(uplim), uplim=0; end;
if isempty(xn), xn=0; end;
if isempty(yn), yn=0; end;
if isempty(first), first=false; end;
if firstCall, first=[true]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT RC
if( first )
errtol =(r1mach(3)./16.0e0).^(1.0e0./6.0e0);
lolim = 5.0e0.*r1mach(1);
uplim = r1mach(2)./5.0e0;
%
c1 = 1.0e0./7.0e0;
c2 = 9.0e0./22.0e0;
end;
first = false;
%
% CALL ERROR HANDLER IF NECESSARY.
%
rcresult = 0.0e0;
if( x<0.0e0 || y<=0.0e0 )
ier = 1;
xern3=sprintf([repmat('%15.6f',1,1)], x);
xern4=sprintf([repmat('%15.6f',1,1)], y);
xermsg('SLATEC','RC',['X.LT.0 .OR. Y.LE.0 WHERE X = ',[xern3,[' AND Y = ',xern4]]],1,1);
%
elseif( max(x,y)>uplim ) ;
ier = 3;
xern3=sprintf([repmat('%15.6f',1,1)], x);
xern4=sprintf([repmat('%15.6f',1,1)], y);
xern5=sprintf([repmat('%15.6f',1,1)], uplim);
xermsg('SLATEC','RC',['MAX(X,Y).GT.UPLIM WHERE X = ',[xern3,[' Y = ',[xern4,[' AND UPLIM = ',xern5]]]]],3,1);
%
elseif( x+y<lolim ) ;
ier = 2;
xern3=sprintf([repmat('%15.6f',1,1)], x);
xern4=sprintf([repmat('%15.6f',1,1)], y);
xern5=sprintf([repmat('%15.6f',1,1)], lolim);
xermsg('SLATEC','RC',['X+Y.LT.LOLIM WHERE X = ',[xern3,[' Y = ',[xern4,[' AND LOLIM = ',xern5]]]]],2,1);
else;
%
ier = 0;
xn = x;
yn = y;
%
while( true );
mu =(xn+yn+yn)./3.0e0;
sn =(yn+mu)./mu - 2.0e0;
if( abs(sn)<errtol )
break;
end;
lamda = 2.0e0.*sqrt(xn).*sqrt(yn) + yn;
xn =(xn+lamda).*0.250e0;
yn =(yn+lamda).*0.250e0;
end;
%
s = sn.*sn.*(0.30e0+sn.*(c1+sn.*(0.3750e0+sn.*c2)));
rcresult =(1.0e0+s)./sqrt(mu);
end;
csnil=dbstack(1); csnil=csnil(1).name(1)~='@';
if csnil&&~isempty(inputname(2)), assignin('caller','FUntemp',y); evalin('caller',[inputname(2),'=FUntemp;']); end
if csnil&&~isempty(inputname(1)), assignin('caller','FUntemp',x); evalin('caller',[inputname(1),'=FUntemp;']); end
if csnil&&~isempty(inputname(3)), assignin('caller','FUntemp',ier); evalin('caller',[inputname(3),'=FUntemp;']); end
end %function rc
%DECK RD
|
|
