| [z,fnu,kode,n,y,nz,tol,elim,alim]=cbknu(z,fnu,kode,n,y,nz,tol,elim,alim); |
function [z,fnu,kode,n,y,nz,tol,elim,alim]=cbknu(z,fnu,kode,n,y,nz,tol,elim,alim);
%***BEGIN PROLOGUE CBKNU
%***SUBSIDIARY
%***PURPOSE Subsidiary to CAIRY, CBESH, CBESI and CBESK
%***LIBRARY SLATEC
%***TYPE ALL (CBKNU-A, ZBKNU-A)
%***AUTHOR Amos, D. E., (SNL)
%***DESCRIPTION
%
% CBKNU COMPUTES THE K BESSEL FUNCTION IN THE RIGHT HALF Z PLANE
%
%***SEE ALSO CAIRY, CBESH, CBESI, CBESK
%***ROUTINES CALLED CKSCL, CSHCH, CUCHK, GAMLN, I1MACH, R1MACH
%***REVISION HISTORY (YYMMDD)
% 830501 DATE WRITTEN
% 910415 Prologue converted to Version 4.0 format. (BAB)
%***end PROLOGUE CBKNU
%
persistent a1 a2 aa ak alas as ascle bb bk bry caz cc cch celm ck coef cone crsc cs cscl csh csr css ctwo cy cz czero dnu dnu2 elm etest f fc fhs firstCall fk fks fmu fpi g1 g2 gt250 gt300 gt500 helim hpi i ic idum iflag igo inu inub j k kflag kk kmax koded nw p p1 p2 p2i p2m p2r pi pt q r1 rk rthpi rz s s1 s2 smu spi st t1 t2 tm tth xd xx yd yy zd ; if isempty(firstCall),firstCall=1;end;
if isempty(cch), cch=0; end;
if isempty(ck), ck=0; end;
if isempty(coef), coef=0; end;
if isempty(cone), cone=0; end;
if isempty(crsc), crsc=0; end;
if isempty(cs), cs=0; end;
if isempty(cscl), cscl=0; end;
if isempty(csh), csh=0; end;
if isempty(csr), csr=zeros(1,3); end;
if isempty(css), css=zeros(1,3); end;
if isempty(ctwo), ctwo=0; end;
if isempty(cz), cz=0; end;
if isempty(czero), czero=0; end;
if isempty(f), f=0; end;
if isempty(fmu), fmu=0; end;
if isempty(p), p=0; end;
if isempty(pt), pt=0; end;
if isempty(p1), p1=0; end;
if isempty(p2), p2=0; end;
if isempty(q), q=0; end;
if isempty(rz), rz=0; end;
if isempty(smu), smu=0; end;
if isempty(st), st=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(zd), zd=0; end;
if isempty(celm), celm=0; end;
if isempty(cy), cy=zeros(1,2); end;
if isempty(aa), aa=0; end;
if isempty(ak), ak=0; end;
if isempty(ascle), ascle=0; end;
if isempty(a1), a1=0; end;
if isempty(a2), a2=0; end;
if isempty(bb), bb=0; end;
if isempty(bk), bk=0; end;
if isempty(bry), bry=zeros(1,3); end;
if isempty(caz), caz=0; end;
if isempty(cc), cc=zeros(1,8); end;
if isempty(dnu), dnu=0; end;
if isempty(dnu2), dnu2=0; end;
if isempty(etest), etest=0; end;
if isempty(fc), fc=0; end;
if isempty(fhs), fhs=0; end;
if isempty(fk), fk=0; end;
if isempty(fks), fks=0; end;
if isempty(fpi), fpi=0; end;
if isempty(g1), g1=0; end;
if isempty(g2), g2=0; end;
if isempty(hpi), hpi=0; end;
if isempty(pi), pi=0; end;
if isempty(p2i), p2i=0; end;
if isempty(p2m), p2m=0; end;
if isempty(p2r), p2r=0; end;
if isempty(rk), rk=0; end;
if isempty(rthpi), rthpi=0; end;
if isempty(r1), r1=0; end;
if isempty(s), s=0; end;
if isempty(spi), spi=0; end;
if isempty(tm), tm=0; end;
if isempty(tth), tth=0; end;
if isempty(t1), t1=0; end;
if isempty(t2), t2=0; end;
if isempty(xx), xx=0; end;
if isempty(yy), yy=0; end;
if isempty(helim), helim=0; end;
if isempty(elm), elm=0; end;
if isempty(xd), xd=0; end;
if isempty(yd), yd=0; end;
if isempty(alas), alas=0; end;
if isempty(as), as=0; end;
if isempty(i), i=0; end;
if isempty(idum), idum=0; end;
if isempty(iflag), iflag=0; end;
if isempty(inu), inu=0; end;
if isempty(k), k=0; end;
if isempty(kflag), kflag=0; end;
if isempty(kk), kk=0; end;
if isempty(kmax), kmax=0; end;
if isempty(koded), koded=0; end;
if isempty(nw), nw=0; end;
if isempty(j), j=0; end;
if isempty(ic), ic=0; end;
if isempty(inub), inub=0; end;
if isempty(igo), igo=0; end;
if isempty(gt250), gt250=0; end;
if isempty(gt300), gt300=0; end;
if isempty(gt500), gt500=0; end;
%
if firstCall, kmax=[30]; end;
if firstCall, r1=[2.0e0]; end;
if firstCall, czero =[complex(0.0e0,0.0e0)]; end;
if firstCall, cone =[complex(1.0e0,0.0e0)]; end;
if firstCall, ctwo=[complex(2.0e0,0.0e0)]; end;
%
if firstCall, pi =[3.14159265358979324e0]; end;
if firstCall, rthpi =[1.25331413731550025e0]; end;
if firstCall, spi =[1.90985931710274403e0]; end;
if firstCall, hpi =[1.57079632679489662e0]; end;
if firstCall, fpi =[1.89769999331517738e0]; end;
if firstCall, tth=[6.66666666666666666e-01]; end;
%
if firstCall, cc(1) =[5.77215664901532861e-01]; end;
if firstCall, cc(2) =[-4.20026350340952355e-02]; end;
if firstCall, cc(3) =[-4.21977345555443367e-02]; end;
if firstCall, cc(4) =[7.21894324666309954e-03]; end;
if firstCall, cc(5) =[-2.15241674114950973e-04]; end;
if firstCall, cc(6) =[-2.01348547807882387e-05]; end;
if firstCall, cc(7) =[1.13302723198169588e-06]; end;
if firstCall, cc(8)=[6.11609510448141582e-09]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT CBKNU
xx = real(z);
yy = imag(z);
caz = abs(z);
cscl = complex(1.0e0./tol,0.0e0);
crsc = complex(tol,0.0e0);
css(1) = cscl;
css(2) = cone;
css(3) = crsc;
csr(1) = crsc;
csr(2) = cone;
csr(3) = cscl;
bry(1) = 1.0e+3.*r1mach(1)./tol;
bry(2) = 1.0e0./bry(1);
[bry(3) ]=r1mach(2);
nz = 0;
iflag = 0;
koded = fix(kode);
rz = ctwo./z;
inu = fix(fnu + 0.5e0);
dnu = fnu - inu;
gt300=0;
gt500=0;
while (1);
if( abs(dnu)~=0.5e0 )
dnu2 = 0.0e0;
if( abs(dnu)>tol )
dnu2 = dnu.*dnu;
end;
if( caz<=r1 )
%-----------------------------------------------------------------------
% SERIES FOR ABS(Z).LE.R1
%-----------------------------------------------------------------------
fc = 1.0e0;
smu = log(rz);
fmu = smu.*complex(dnu,0.0e0);
[fmu,csh,cch]=cshch(fmu,csh,cch);
if( dnu~=0.0e0 )
fc = dnu.*pi;
fc = fc./sin(fc);
smu = csh.*complex(1.0e0./dnu,0.0e0);
end;
a2 = 1.0e0 + dnu;
%-----------------------------------------------------------------------
% GAM(1-Z)*GAM(1+Z)=PI*Z/SIN(PI*Z), T1=1/GAM(1-DNU), T2=1/GAM(1+DNU)
%-----------------------------------------------------------------------
t2 = exp(-gamln(a2,idum));
t1 = 1.0e0./(t2.*fc);
if( abs(dnu)>0.1e0 )
g1 =(t1-t2)./(dnu+dnu);
else;
%-----------------------------------------------------------------------
% SERIES FOR F0 TO RESOLVE INDETERMINACY FOR SMALL ABS(DNU)
%-----------------------------------------------------------------------
ak = 1.0e0;
s = cc(1);
for k = 2 : 8;
ak = ak.*dnu2;
tm = cc(k).*ak;
s = s + tm;
if( abs(tm)<tol )
break;
end;
end;
g1 = -s;
end;
g2 = 0.5e0.*(t1+t2).*fc;
g1 = g1.*fc;
f = complex(g1,0.0e0).*cch + smu.*complex(g2,0.0e0);
pt = exp(fmu);
p = complex(0.5e0./t2,0.0e0).*pt;
q = complex(0.5e0./t1,0.0e0)./pt;
s1 = f;
s2 = p;
ak = 1.0e0;
a1 = 1.0e0;
ck = cone;
bk = 1.0e0 - dnu2;
if( inu>0 || n>1 )
%-----------------------------------------------------------------------
% GENERATE K(DNU,Z) AND K(DNU+1,Z) FOR FORWARD RECURRENCE
%-----------------------------------------------------------------------
if( caz>=tol )
cz = z.*z.*complex(0.25e0,0.0e0);
t1 = 0.25e0.*caz.*caz;
while (1);
f =(f.*complex(ak,0.0e0)+p+q).*complex(1.0e0./bk,0.0e0);
p = p.*complex(1.0e0./(ak-dnu),0.0e0);
q = q.*complex(1.0e0./(ak+dnu),0.0e0);
rk = 1.0e0./ak;
ck = ck.*cz.*complex(rk,0.0e0);
s1 = s1 + ck.*f;
s2 = s2 + ck.*(p-f.*complex(ak,0.0e0));
a1 = a1.*t1.*rk;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
if( a1<=tol )
break;
end;
end;
end;
kflag = 2;
bk = real(smu);
a1 = fnu + 1.0e0;
ak = a1.*abs(bk);
if( ak>alim )
kflag = 3;
end;
p2 = s2.*css(kflag);
s2 = p2.*rz;
s1 = s1.*css(kflag);
if( koded~=1 )
f = exp(z);
s1 = s1.*f;
s2 = s2.*f;
end;
%to 200
break;
else;
%-----------------------------------------------------------------------
% GENERATE K(FNU,Z), 0.0D0 .LE. FNU .LT. 0.5D0 AND N=1
%-----------------------------------------------------------------------
if( caz>=tol )
cz = z.*z.*complex(0.25e0,0.0e0);
t1 = 0.25e0.*caz.*caz;
while (1);
f =(f.*complex(ak,0.0e0)+p+q).*complex(1.0e0./bk,0.0e0);
p = p.*complex(1.0e0./(ak-dnu),0.0e0);
q = q.*complex(1.0e0./(ak+dnu),0.0e0);
rk = 1.0e0./ak;
ck = ck.*cz.*complex(rk,0.0);
s1 = s1 + ck.*f;
a1 = a1.*t1.*rk;
bk = bk + ak + ak + 1.0e0;
ak = ak + 1.0e0;
if( a1<=tol )
break;
end;
end;
end;
y(1) = s1;
if( koded==1 )
return;
end;
y(1) = s1.*exp(z);
return;
end;
end;
end;
%-----------------------------------------------------------------------
% IFLAG=0 MEANS NO UNDERFLOW OCCURRED
% IFLAG=1 MEANS AN UNDERFLOW OCCURRED- COMPUTATION PROCEEDS WITH
% KODED=2 AND A TEST FOR ON SCALE VALUES IS MADE DURING FORWARD
% RECURSION
%-----------------------------------------------------------------------
coef = complex(rthpi,0.0e0)./sqrt(z);
kflag = 2;
if( koded~=2 )
if( xx>alim )
%-----------------------------------------------------------------------
% SCALE BY EXP(Z), IFLAG = 1 CASES
%-----------------------------------------------------------------------
koded = 2;
iflag = 1;
kflag = 2;
else;
% BLANK LINE
a1 = exp(-xx).*real(css(kflag));
pt = complex(a1,0.0e0).*complex(cos(yy),-sin(yy));
coef = coef.*pt;
end;
end;
if( abs(dnu)==0.5e0 )
s1 = coef;
s2 = coef;
break;
end;
%-----------------------------------------------------------------------
% MILLER ALGORITHM FOR ABS(Z).GT.R1
%-----------------------------------------------------------------------
ak = cos(pi.*dnu);
ak = abs(ak);
if( ak==0.0e0 )
s1 = coef;
s2 = coef;
break;
end;
fhs = abs(0.25e0-dnu2);
if( fhs==0.0e0 )
s1 = coef;
s2 = coef;
break;
end;
%-----------------------------------------------------------------------
% COMPUTE R2=F(E). IF ABS(Z).GE.R2, USE FORWARD RECURRENCE TO
% DETERMINE THE BACKWARD INDEX K. R2=F(E) IS A STRAIGHT LINE ON
% 12.LE.E.LE.60. E IS COMPUTED FROM 2**(-E)=B**(1-I1MACH(11))=
% TOL WHERE B IS THE BASE OF THE ARITHMETIC.
%-----------------------------------------------------------------------
t1 =(i1mach(11)-1).*r1mach(5).*3.321928094e0;
t1 = max(t1,12.0e0);
t1 = min(t1,60.0e0);
t2 = tth.*t1 - 6.0e0;
if( xx~=0.0e0 )
t1 = atan(yy./xx);
t1 = abs(t1);
else;
t1 = hpi;
end;
while (1);
igo=0;
if( t2>caz )
%-----------------------------------------------------------------------
% COMPUTE BACKWARD INDEX K FOR ABS(Z).LT.R2
%-----------------------------------------------------------------------
a2 = sqrt(caz);
ak = fpi.*ak./(tol.*sqrt(a2));
aa = 3.0e0.*t1./(1.0e0+caz);
bb = 14.7e0.*t1./(28.0e0+caz);
ak =(log(ak)+caz.*cos(aa)./(1.0e0+0.008e0.*caz))./cos(bb);
fk = 0.12125e0.*ak.*ak./caz + 1.5e0;
else;
%-----------------------------------------------------------------------
% FORWARD RECURRENCE LOOP WHEN ABS(Z).GE.R2
%-----------------------------------------------------------------------
etest = ak./(pi.*caz.*tol);
fk = 1.0e0;
if( etest>=1.0e0 )
fks = 2.0e0;
rk = caz + caz + 2.0e0;
a1 = 0.0e0;
a2 = 1.0e0;
for i = 1 : kmax;
ak = fhs./fks;
bk = rk./(fk+1.0e0);
tm = a2;
a2 = bk.*a2 - ak.*a1;
a1 = tm;
rk = rk + 2.0e0;
fks = fks + fk + fk + 2.0e0;
fhs = fhs + fk + fk;
fk = fk + 1.0e0;
tm = abs(a2).*fk;
if( etest<tm )
fk = fk + spi.*t1.*sqrt(t2./caz);
fhs = abs(0.25e0-dnu2);
igo=1;
break;
end;
end;
if(igo==1)
break;
end;
nz = -2;
return;
end;
end;
break;
end;
k = fix(fk);
%-----------------------------------------------------------------------
% BACKWARD RECURRENCE LOOP FOR MILLER ALGORITHM
%-----------------------------------------------------------------------
fk = k;
fks = fk.*fk;
p1 = czero;
p2 = complex(tol,0.0e0);
cs = p2;
for i = 1 : k;
a1 = fks - fk;
a2 =(fks+fk)./(a1+fhs);
rk = 2.0e0./(fk+1.0e0);
t1 =(fk+xx).*rk;
t2 = yy.*rk;
pt = p2;
p2 =(p2.*complex(t1,t2)-p1).*complex(a2,0.0e0);
p1 = pt;
cs = cs + p2;
fks = a1 - fk + 1.0e0;
fk = fk - 1.0e0;
end; i = fix(k+1);
%-----------------------------------------------------------------------
% COMPUTE (P2/CS)=(P2/ABS(CS))*(CONJG(CS)/ABS(CS)) FOR BETTER
% SCALING
%-----------------------------------------------------------------------
tm = abs(cs);
pt = complex(1.0e0./tm,0.0e0);
s1 = pt.*p2;
cs = conj(cs).*pt;
s1 = coef.*s1.*cs;
if( inu>0 || n>1 )
%-----------------------------------------------------------------------
% COMPUTE P1/P2=(P1/ABS(P2)*CONJG(P2)/ABS(P2) FOR SCALING
%-----------------------------------------------------------------------
tm = abs(p2);
pt = complex(1.0e0./tm,0.0e0);
p1 = pt.*p1;
p2 = conj(p2).*pt;
pt = p1.*p2;
s2 = s1.*(cone+(complex(dnu+0.5e0,0.0e0)-pt)./z);
else;
zd = z;
if( iflag~=1 )
gt300=1;
break;
end;
gt500=1;
break;
end;
%-----------------------------------------------------------------------
% FORWARD RECURSION ON THE THREE TERM RECURSION RELATION WITH
% SCALING NEAR EXPONENT EXTREMES ON KFLAG=1 OR KFLAG=3
%-----------------------------------------------------------------------
break;
%pre200
end;
%200
while (1);
while(gt500==0);
if(gt300==0)
ck = complex(dnu+1.0e0,0.0e0).*rz;
if( n==1 )
inu = fix(inu - 1);
end;
if( inu>0 )
inub = 1;
%250
while (1);
if( iflag==1 )
%-----------------------------------------------------------------------
% IFLAG=1 CASES, FORWARD RECURRENCE ON SCALED VALUES ON UNDERFLOW
%-----------------------------------------------------------------------
helim = 0.5e0.*elim;
elm = exp(-elim);
celm = complex(elm,0.0);
ascle = bry(1);
zd = z;
xd = xx;
yd = yy;
ic = -1;
j = 2;
for i = 1 : inu;
gt250=0;
st = s2;
s2 = ck.*s2 + s1;
s1 = st;
ck = ck + rz;
as = abs(s2);
alas = log(as);
p2r = -xd + alas;
if( p2r>=(-elim) )
p2 = -zd + log(s2);
p2r = real(p2);
p2i = imag(p2);
p2m = exp(p2r)./tol;
p1 = complex(p2m,0.0e0).*complex(cos(p2i),sin(p2i));
[p1,nw,ascle,tol]=cuchk(p1,nw,ascle,tol);
if( nw==0 )
j = fix(3 - j);
cy(j) = p1;
if( ic==(i-1) )
kflag = 1;
inub = fix(i + 1);
s2 = cy(j);
j = fix(3 - j);
s1 = cy(j);
if( inub<=inu )
gt250=1;
break;
end;
if( n==1 )
s1 = s2;
end;
gt300=1;
break;
else;
ic = fix(i);
continue;
end;
end;
end;
if( alas>=helim )
xd = xd - elim;
s1 = s1.*celm;
s2 = s2.*celm;
zd = complex(xd,yd);
end;
end;
if(gt300==1)
break;
end;
if(gt250==1)
break;
end;
if( n==1 )
s1 = s2;
end;
gt500=1;
%500
break;
end;
break;
end;
%500
if(gt500==1)
break;
end;
if(gt300==0)
p1 = csr(kflag);
ascle = bry(kflag);
for i = inub : inu;
st = s2;
s2 = ck.*s2 + s1;
s1 = st;
ck = ck + rz;
if( kflag<3 )
p2 = s2.*p1;
p2r = real(p2);
p2i = imag(p2);
p2r = abs(p2r);
p2i = abs(p2i);
p2m = max(p2r,p2i);
if( p2m>ascle )
kflag = fix(kflag + 1);
ascle = bry(kflag);
s1 = s1.*p1;
s2 = p2;
s1 = s1.*css(kflag);
s2 = s2.*css(kflag);
p1 = csr(kflag);
end;
end;
end; i = fix(inu+1);
if( n==1 )
s1 = s2;
end;
%gt300
end;
else;
if( n==1 )
s1 = s2;
end;
zd = z;
%500
if( iflag==1 )
break;
end;
end;
%gt300
end;
gt300=0;
y(1) = s1.*csr(kflag);
if( n==1 )
return;
end;
y(2) = s2.*csr(kflag);
if( n==2 )
return;
end;
kk = 2;
kk = fix(kk + 1);
if( kk>n )
return;
end;
p1 = csr(kflag);
ascle = bry(kflag);
for i = kk : n;
p2 = s2;
s2 = ck.*s2 + s1;
s1 = p2;
ck = ck + rz;
p2 = s2.*p1;
y(i) = p2;
if( kflag<3 )
p2r = real(p2);
p2i = imag(p2);
p2r = abs(p2r);
p2i = abs(p2i);
p2m = max(p2r,p2i);
if( p2m>ascle )
kflag = fix(kflag + 1);
ascle = bry(kflag);
s1 = s1.*p1;
s2 = p2;
s1 = s1.*css(kflag);
s2 = s2.*css(kflag);
p1 = csr(kflag);
end;
end;
end; i = fix(n+1);
return;
%gt500
end;
gt500=0;
y(1) = s1;
if( n~=1 )
y(2) = s2;
end;
ascle = bry(1);
[zd,fnu,n,y,nz,rz,ascle,tol,elim]=ckscl(zd,fnu,n,y,nz,rz,ascle,tol,elim);
inu = fix(n - nz);
if( inu<=0 )
return;
end;
kk = fix(nz + 1);
s1 = y(kk);
y(kk) = s1.*csr(1);
if( inu==1 )
return;
end;
kk = fix(nz + 2);
s2 = y(kk);
y(kk) = s2.*csr(1);
if( inu==2 )
return;
end;
t2 = fnu +(kk-1);
ck = complex(t2,0.0e0).*rz;
kflag = 1;
kk = fix(kk + 1);
if( kk>n )
return;
end;
p1 = csr(kflag);
ascle = bry(kflag);
for i = kk : n;
p2 = s2;
s2 = ck.*s2 + s1;
s1 = p2;
ck = ck + rz;
p2 = s2.*p1;
y(i) = p2;
if( kflag<3 )
p2r = real(p2);
p2i = imag(p2);
p2r = abs(p2r);
p2i = abs(p2i);
p2m = max(p2r,p2i);
if( p2m>ascle )
kflag = fix(kflag + 1);
ascle = bry(kflag);
s1 = s1.*p1;
s2 = p2;
s1 = s1.*css(kflag);
s2 = s2.*css(kflag);
p1 = csr(kflag);
end;
end;
end; i = fix(n+1);
return;
%-----------------------------------------------------------------------
% FNU=HALF ODD INTEGER CASE, DNU=-0.5
%-----------------------------------------------------------------------
s1 = coef;
s2 = coef;
%200
end;
end %subroutine cbknu
%DECK CBLKT1
|
|
