| [z,id,kode,ai,nz,ierr]=cairy(z,id,kode,ai,nz,ierr); |
function [z,id,kode,ai,nz,ierr]=cairy(z,id,kode,ai,nz,ierr);
%***BEGIN PROLOGUE CAIRY
%***PURPOSE Compute the Airy function Ai(z) or its derivative dAi/dz
% for complex argument z. A scaling option is available
% to help avoid underflow and overflow.
%***LIBRARY SLATEC
%***CATEGORY C10D
%***TYPE COMPLEX (CAIRY-C, ZAIRY-C)
%***KEYWORDS AIRY FUNCTION, BESSEL FUNCTION OF ORDER ONE THIRD,
% BESSEL FUNCTION OF ORDER TWO THIRDS
%***AUTHOR Amos, D. E., (SNL)
%***DESCRIPTION
%
% On KODE=1, CAIRY computes the complex Airy function Ai(z)
% or its derivative dAi/dz on ID=0 or ID=1 respectively. On
% KODE=2, a scaling option exp(zeta)*Ai(z) or exp(zeta)*dAi/dz
% is provided to remove the exponential decay in -pi/3<arg(z)
% <pi/3 and the exponential growth in pi/3<abs(arg(z))<pi where
% zeta=(2/3)*z**(3/2).
%
% While the Airy functions Ai(z) and dAi/dz are analytic in
% the whole z-plane, the corresponding scaled functions defined
% for KODE=2 have a cut along the negative real axis.
%
% Input
% Z - Argument of type COMPLEX
% ID - Order of derivative, ID=0 or ID=1
% KODE - A parameter to indicate the scaling option
% KODE=1 returns
% AI=Ai(z) on ID=0
% AI=dAi/dz on ID=1
% at z=Z
% =2 returns
% AI=exp(zeta)*Ai(z) on ID=0
% AI=exp(zeta)*dAi/dz on ID=1
% at z=Z where zeta=(2/3)*z**(3/2)
%
% Output
% AI - Result of type COMPLEX
% NZ - Underflow indicator
% NZ=0 Normal return
% NZ=1 AI=0 due to underflow in
% -pi/3<arg(Z)<pi/3 on KODE=1
% IERR - Error flag
% IERR=0 Normal return - COMPUTATION COMPLETED
% IERR=1 Input error - NO COMPUTATION
% IERR=2 Overflow - NO COMPUTATION
% (Re(Z) too large with KODE=1)
% IERR=3 Precision warning - COMPUTATION COMPLETED
% (Result has less than half precision)
% IERR=4 Precision error - NO COMPUTATION
% (Result has no precision)
% IERR=5 Algorithmic error - NO COMPUTATION
% (Termination condition not met)
%
% *Long Description:
%
% Ai(z) and dAi/dz are computed from K Bessel functions by
%
% Ai(z) = c*sqrt(z)*K(1/3,zeta)
% dAi/dz = -c* z *K(2/3,zeta)
% c = 1/(pi*sqrt(3))
% zeta = (2/3)*z**(3/2)
%
% when abs(z)>1 and from power series when abs(z)<=1.
%
% In most complex variable computation, one must evaluate ele-
% mentary functions. When the magnitude of Z is large, losses
% of significance by argument reduction occur. Consequently, if
% the magnitude of ZETA=(2/3)*Z**(3/2) exceeds U1=SQRT(0.5/UR),
% then losses exceeding half precision are likely and an error
% flag IERR=3 is triggered where UR=R1MACH(4)=UNIT ROUNDOFF.
% Also, if the magnitude of ZETA is larger than U2=0.5/UR, then
% all significance is lost and IERR=4. In order to use the INT
% function, ZETA must be further restricted not to exceed
% U3=I1MACH(9)=LARGEST INTEGER. Thus, the magnitude of ZETA
% must be restricted by MIN(U2,U3). In IEEE arithmetic, U1,U2,
% and U3 are approximately 2.0E+3, 4.2E+6, 2.1E+9 in single
% precision and 4.7E+7, 2.3E+15, 2.1E+9 in doubleprecision.
% This makes U2 limiting is single precision and U3 limiting
% in doubleprecision. This means that the magnitude of Z
% cannot exceed approximately 3.4E+4 in single precision and
% 2.1E+6 in doubleprecision. This also means that one can
% expect to retain, in the worst cases on 32-bit machines,
% no digits in single precision and only 6 digits in double
% precision.
%
% The approximate relative error in the magnitude of a complex
% Bessel function can be expressed as P*10**S where P=MAX(UNIT
% ROUNDOFF,1.0E-18) is the nominal precision and 10**S repre-
% sents the increase in error due to argument reduction in the
% elementary functions. Here, S=MAX(1,ABS(LOG10(ABS(Z))),
% ABS(LOG10(FNU))) approximately (i.e., S=MAX(1,ABS(EXPONENT OF
% ABS(Z),ABS(EXPONENT OF FNU)) ). However, the phase angle may
% have only absolute accuracy. This is most likely to occur
% when one component (in magnitude) is larger than the other by
% several orders of magnitude. If one component is 10**K larger
% than the other, then one can expect only MAX(ABS(LOG10(P))-K,
% 0) significant digits; or, stated another way, when K exceeds
% the exponent of P, no significant digits remain in the smaller
% component. However, the phase angle retains absolute accuracy
% because, in complex arithmetic with precision P, the smaller
% component will not (as a rule) decrease below P times the
% magnitude of the larger component. In these extreme cases,
% the principal phase angle is on the order of +P, -P, PI/2-P,
% or -PI/2+P.
%
%***REFERENCES 1. M. Abramowitz and I. A. Stegun, Handbook of Mathe-
% matical Functions, National Bureau of Standards
% Applied Mathematics Series 55, U. S. Department
% of Commerce, Tenth Printing (1972) or later.
% 2. D. E. Amos, Computation of Bessel Functions of
% Complex Argument and Large Order, Report SAND83-0643,
% Sandia National Laboratories, Albuquerque, NM, May
% 1983.
% 3. D. E. Amos, A Subroutine Package for Bessel Functions
% of a Complex Argument and Nonnegative Order, Report
% SAND85-1018, Sandia National Laboratory, Albuquerque,
% NM, May 1985.
% 4. D. E. Amos, A portable package for Bessel functions
% of a complex argument and nonnegative order, ACM
% Transactions on Mathematical Software, 12 (September
% 1986), pp. 265-273.
%
%***ROUTINES CALLED CACAI, CBKNU, I1MACH, R1MACH
%***REVISION HISTORY (YYMMDD)
% 830501 DATE WRITTEN
% 890801 REVISION DATE from Version 3.2
% 910415 Prologue converted to Version 4.0 format. (BAB)
% 920128 Category corrected. (WRB)
% 920811 Prologue revised. (DWL)
%***end PROLOGUE CAIRY
persistent aa ad ak alaz alim atrm az az3 bb bk c1 c2 ck coef cone csq cy d1 d2 dig dk elim fid firstCall fnu iflag k k1 k2 mr nn r1m5 rl s1 s2 sfac tol trm1 trm2 tth z3 z3i z3r zi zr zta ; if isempty(firstCall),firstCall=1;end;
if isempty(cone), cone=0; end;
if isempty(csq), csq=0; end;
if isempty(cy), cy=zeros(1,1); end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(trm1), trm1=0; end;
if isempty(trm2), trm2=0; end;
if isempty(zta), zta=0; end;
if isempty(z3), z3=0; end;
if isempty(aa), aa=0; end;
if isempty(ad), ad=0; end;
if isempty(ak), ak=0; end;
if isempty(alim), alim=0; end;
if isempty(atrm), atrm=0; end;
if isempty(az), az=0; end;
if isempty(az3), az3=0; end;
if isempty(bk), bk=0; end;
if isempty(ck), ck=0; end;
if isempty(coef), coef=0; end;
if isempty(c1), c1=0; end;
if isempty(c2), c2=0; end;
if isempty(dig), dig=0; end;
if isempty(dk), dk=0; end;
if isempty(d1), d1=0; end;
if isempty(d2), d2=0; end;
if isempty(elim), elim=0; end;
if isempty(fid), fid=0; end;
if isempty(fnu), fnu=0; end;
if isempty(rl), rl=0; end;
if isempty(r1m5), r1m5=0; end;
if isempty(sfac), sfac=0; end;
if isempty(tol), tol=0; end;
if isempty(tth), tth=0; end;
if isempty(zi), zi=0; end;
if isempty(zr), zr=0; end;
if isempty(z3i), z3i=0; end;
if isempty(z3r), z3r=0; end;
if isempty(bb), bb=0; end;
if isempty(alaz), alaz=0; end;
if isempty(iflag), iflag=0; end;
if isempty(k), k=0; end;
if isempty(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(mr), mr=0; end;
if isempty(nn), nn=0; end;
if firstCall, tth =[6.66666666666666667e-01]; end;
if firstCall, c1 =[3.55028053887817240e-01]; end;
if firstCall, c2 =[2.58819403792806799e-01]; end;
if firstCall, coef=[1.83776298473930683e-01]; end;
if firstCall, cone=[complex(1.0e0,0.0e0)]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT CAIRY
ierr = 0;
nz = 0;
if( id<0 || id>1 )
ierr = 1;
end;
if( kode<1 || kode>2 )
ierr = 1;
end;
if( ierr~=0 )
return;
end;
az = abs(z);
tol = max(r1mach(4),1.0e-18);
fid = id;
if( az>1.0e0 )
%-----------------------------------------------------------------------
% CASE FOR ABS(Z).GT.1.0
%-----------------------------------------------------------------------
fnu =(1.0e0+fid)./3.0e0;
%-----------------------------------------------------------------------
% SET PARAMETERS RELATED TO MACHINE CONSTANTS.
% TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
% ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
% EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
% EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
% UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
% RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
% DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
%-----------------------------------------------------------------------
[k1 ]=i1mach(12);
[k2 ]=i1mach(13);
[r1m5 ]=r1mach(5);
k = fix(min(abs(k1),abs(k2)));
elim = 2.303e0.*(k.*r1m5-3.0e0);
k1 = fix(i1mach(11) - 1);
aa = r1m5.*k1;
dig = min(aa,18.0e0);
aa = aa.*2.303e0;
alim = elim + max(-aa,-41.45e0);
rl = 1.2e0.*dig + 3.0e0;
alaz = log(az);
%-----------------------------------------------------------------------
% TEST FOR RANGE
%-----------------------------------------------------------------------
aa = 0.5e0./tol;
bb = i1mach(9).*0.5e0;
aa = min(aa,bb);
aa = aa.^tth;
if( az>aa )
ierr = 4;
nz = 0;
return;
else;
aa = sqrt(aa);
if( az>aa )
ierr = 3;
end;
csq = sqrt(z);
zta = z.*csq.*complex(tth,0.0e0);
%-----------------------------------------------------------------------
% RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
%-----------------------------------------------------------------------
iflag = 0;
sfac = 1.0e0;
zi = imag(z);
zr = real(z);
ak = imag(zta);
if( zr<0.0e0 )
bk = real(zta);
ck = -abs(bk);
zta = complex(ck,ak);
end;
if( zi==0.0e0 )
if( zr<=0.0e0 )
zta = complex(0.0e0,ak);
end;
end;
aa = real(zta);
if( aa<0.0e0 || zr<=0.0e0 )
if( kode~=2 )
%-----------------------------------------------------------------------
% OVERFLOW TEST
%-----------------------------------------------------------------------
if( aa<=(-alim) )
aa = -aa + 0.25e0.*alaz;
iflag = 1;
sfac = tol;
if( aa>elim )
nz = 0;
ierr = 2;
return;
end;
end;
end;
%-----------------------------------------------------------------------
% CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
%-----------------------------------------------------------------------
mr = 1;
if( zi<0.0e0 )
mr = -1;
end;
[zta,fnu,kode,mr,dumvar5,cy,nn,rl,tol,elim,alim]=cacai(zta,fnu,kode,mr,1,cy,nn,rl,tol,elim,alim);
if( nn<0 )
if( nn==(-1) )
nz = 0;
ierr = 2;
return;
end;
nz = 0;
ierr = 5;
return;
else;
nz = fix(nz + nn);
end;
elseif( kode==2 ) ;
[zta,fnu,kode,dumvar4,cy,nz,tol,elim,alim]=cbknu(zta,fnu,kode,1,cy,nz,tol,elim,alim);
%-----------------------------------------------------------------------
% UNDERFLOW TEST
%-----------------------------------------------------------------------
elseif( aa<alim ) ;
[zta,fnu,kode,dumvar4,cy,nz,tol,elim,alim]=cbknu(zta,fnu,kode,1,cy,nz,tol,elim,alim);
else;
aa = -aa - 0.25e0.*alaz;
iflag = 2;
sfac = 1.0e0./tol;
if( aa<(-elim) )
nz = 1;
ai = complex(0.0e0,0.0e0);
return;
else;
[zta,fnu,kode,dumvar4,cy,nz,tol,elim,alim]=cbknu(zta,fnu,kode,1,cy,nz,tol,elim,alim);
end;
end;
s1 = cy(1).*complex(coef,0.0e0);
if( iflag~=0 )
s1 = s1.*complex(sfac,0.0e0);
if( id==1 )
s1 = -s1.*z;
ai = s1.*complex(1.0e0./sfac,0.0e0);
return;
else;
s1 = s1.*csq;
ai = s1.*complex(1.0e0./sfac,0.0e0);
return;
end;
elseif( id==1 ) ;
ai = -z.*s1;
return;
else;
ai = csq.*s1;
return;
end;
end;
nz = 0;
ierr = 2;
return;
else;
%-----------------------------------------------------------------------
% POWER SERIES FOR ABS(Z).LE.1.
%-----------------------------------------------------------------------
s1 = cone;
s2 = cone;
if( az<tol )
aa = 1.0e+3.*r1mach(1);
s1 = complex(0.0e0,0.0e0);
if( id==1 )
ai = -complex(c2,0.0e0);
aa = sqrt(aa);
if( az>aa )
s1 = z.*z.*complex(0.5e0,0.0e0);
end;
ai = ai + s1.*complex(c1,0.0e0);
return;
else;
if( az>aa )
s1 = complex(c2,0.0e0).*z;
end;
ai = complex(c1,0.0e0) - s1;
return;
end;
else;
aa = az.*az;
if( aa>=tol./az )
trm1 = cone;
trm2 = cone;
atrm = 1.0e0;
z3 = z.*z.*z;
az3 = az.*aa;
ak = 2.0e0 + fid;
bk = 3.0e0 - fid - fid;
ck = 4.0e0 - fid;
dk = 3.0e0 + fid + fid;
d1 = ak.*dk;
d2 = bk.*ck;
ad = min(d1,d2);
ak = 24.0e0 + 9.0e0.*fid;
bk = 30.0e0 - 9.0e0.*fid;
z3r = real(z3);
z3i = imag(z3);
for k = 1 : 25;
trm1 = trm1.*complex(z3r./d1,z3i./d1);
s1 = s1 + trm1;
trm2 = trm2.*complex(z3r./d2,z3i./d2);
s2 = s2 + trm2;
atrm = atrm.*az3./ad;
d1 = d1 + ak;
d2 = d2 + bk;
ad = min(d1,d2);
if( atrm<tol.*ad )
break;
end;
ak = ak + 18.0e0;
bk = bk + 18.0e0;
end;
end;
if( id==1 )
ai = -s2.*complex(c2,0.0e0);
if( az>tol )
ai = ai + z.*z.*s1.*complex(c1./(1.0e0+fid),0.0e0);
end;
if( kode==1 )
return;
end;
zta = z.*sqrt(z).*complex(tth,0.0e0);
ai = ai.*exp(zta);
return;
else;
ai = s1.*complex(c1,0.0e0) - z.*s2.*complex(c2,0.0e0);
if( kode==1 )
return;
end;
zta = z.*sqrt(z).*complex(tth,0.0e0);
ai = ai.*exp(zta);
return;
end;
end;
end;
end %subroutine cairy
%DECK CARG
|
|
