| [x,n,ktrms,t,ansmlv,ind,ms,gmrn,h,ierr]=dbkias(x,n,ktrms,t,ansmlv,ind,ms,gmrn,h,ierr); |
function [x,n,ktrms,t,ansmlv,ind,ms,gmrn,h,ierr]=dbkias(x,n,ktrms,t,ansmlv,ind,ms,gmrn,h,ierr);
%***BEGIN PROLOGUE DBKIAS
%***SUBSIDIARY
%***PURPOSE Subsidiary to DBSKIN
%***LIBRARY SLATEC
%***TYPE doubleprecision (BKIAS-S, DBKIAS-D)
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% DBKIAS computes repeated integrals of the K0 Bessel function
% by the asymptotic expansion
%
%***SEE ALSO DBSKIN
%***ROUTINES CALLED D1MACH, DBDIFF, DGAMRN, DHKSEQ
%***REVISION HISTORY (YYMMDD)
% 820601 DATE WRITTEN
% 890531 Changed all specific intrinsics to generic. (WRB)
% 890911 Removed unnecessary intrinsics. (WRB)
% 891214 Prologue converted to Version 4.0 format. (BAB)
% 900328 Added TYPE section. (WRB)
% 910722 Updated AUTHOR section. (ALS)
%***end PROLOGUE DBKIAS
persistent b bnd den1 den2 den3 er err firstCall fj fk fln fm1 g1 gs gt hn hrtpi i ii j jmi jn k kk km mm mp rat rg1 rxp rz rzx s ss sumi sumj tol v w xp z ; if isempty(firstCall),firstCall=1;end;
if isempty(i), i=0; end;
if isempty(ii), ii=0; end;
if isempty(j), j=0; end;
if isempty(jmi), jmi=0; end;
if isempty(jn), jn=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(km), km=0; end;
if isempty(mm), mm=0; end;
if isempty(mp), mp=0; end;
if isempty(gt), gt=0; end;
if isempty(b), b=zeros(1,120); end;
if isempty(bnd), bnd=zeros(1,15); end;
if isempty(den1), den1=0; end;
if isempty(den2), den2=0; end;
if isempty(den3), den3=0; end;
if isempty(er), er=0; end;
if isempty(err), err=0; end;
if isempty(fj), fj=0; end;
if isempty(fk), fk=0; end;
if isempty(fln), fln=0; end;
if isempty(fm1), fm1=0; end;
if isempty(g1), g1=0; end;
if isempty(gs), gs=0; end;
if isempty(hn), hn=0; end;
if isempty(hrtpi), hrtpi=0; end;
if isempty(rat), rat=0; end;
if isempty(rg1), rg1=0; end;
if isempty(rxp), rxp=0; end;
if isempty(rz), rz=0; end;
if isempty(rzx), rzx=0; end;
if isempty(s), s=zeros(1,31); end;
if isempty(ss), ss=0; end;
if isempty(sumi), sumi=0; end;
if isempty(sumj), sumj=0; end;
if isempty(tol), tol=0; end;
if isempty(v), v=zeros(1,52); end;
if isempty(w), w=zeros(1,52); end;
if isempty(xp), xp=zeros(1,16); end;
if isempty(z), z=0; end;
h_shape=size(h);h=reshape(h,1,[]);
%-----------------------------------------------------------------------
% COEFFICIENTS OF POLYNOMIAL P(J-1,X), J=1,15
%-----------------------------------------------------------------------
if firstCall, b(1) =[1.00000000000000000d+00]; end;
if firstCall, b(2) =[1.00000000000000000d+00]; end;
if firstCall, b(3) =[-2.00000000000000000d+00]; end;
if firstCall, b(4) =[1.00000000000000000d+00]; end;
if firstCall, b(5) =[-8.00000000000000000d+00]; end;
if firstCall, b(6) =[6.00000000000000000d+00]; end;
if firstCall, b(7) =[1.00000000000000000d+00]; end;
if firstCall, b(8) =[-2.20000000000000000d+01]; end;
if firstCall, b(9) =[5.80000000000000000d+01]; end;
if firstCall, b(10) =[-2.40000000000000000d+01]; end;
if firstCall, b(11) =[1.00000000000000000d+00]; end;
if firstCall, b(12) =[-5.20000000000000000d+01]; end;
if firstCall, b(13) =[3.28000000000000000d+02]; end;
if firstCall, b(14) =[-4.44000000000000000d+02]; end;
if firstCall, b(15) =[1.20000000000000000d+02]; end;
if firstCall, b(16) =[1.00000000000000000d+00]; end;
if firstCall, b(17) =[-1.14000000000000000d+02]; end;
if firstCall, b(18) =[1.45200000000000000d+03]; end;
if firstCall, b(19) =[-4.40000000000000000d+03]; end;
if firstCall, b(20) =[3.70800000000000000d+03]; end;
if firstCall, b(21) =[-7.20000000000000000d+02]; end;
if firstCall, b(22) =[1.00000000000000000d+00]; end;
if firstCall, b(23) =[-2.40000000000000000d+02]; end;
if firstCall, b(24)=[5.61000000000000000d+03]; end;
if firstCall, b(25) =[-3.21200000000000000d+04]; end;
if firstCall, b(26) =[5.81400000000000000d+04]; end;
if firstCall, b(27) =[-3.39840000000000000d+04]; end;
if firstCall, b(28) =[5.04000000000000000d+03]; end;
if firstCall, b(29) =[1.00000000000000000d+00]; end;
if firstCall, b(30) =[-4.94000000000000000d+02]; end;
if firstCall, b(31) =[1.99500000000000000d+04]; end;
if firstCall, b(32) =[-1.95800000000000000d+05]; end;
if firstCall, b(33) =[6.44020000000000000d+05]; end;
if firstCall, b(34) =[-7.85304000000000000d+05]; end;
if firstCall, b(35) =[3.41136000000000000d+05]; end;
if firstCall, b(36) =[-4.03200000000000000d+04]; end;
if firstCall, b(37) =[1.00000000000000000d+00]; end;
if firstCall, b(38) =[-1.00400000000000000d+03]; end;
if firstCall, b(39) =[6.72600000000000000d+04]; end;
if firstCall, b(40) =[-1.06250000000000000d+06]; end;
if firstCall, b(41) =[5.76550000000000000d+06]; end;
if firstCall, b(42) =[-1.24400640000000000d+07]; end;
if firstCall, b(43) =[1.10262960000000000d+07]; end;
if firstCall, b(44) =[-3.73392000000000000d+06]; end;
if firstCall, b(45) =[3.62880000000000000d+05]; end;
if firstCall, b(46) =[1.00000000000000000d+00]; end;
if firstCall, b(47) =[-2.02600000000000000d+03]; end;
if firstCall, b(48)=[2.18848000000000000d+05]; end;
if firstCall, b(49) =[-5.32616000000000000d+06]; end;
if firstCall, b(50) =[4.47650000000000000d+07]; end;
if firstCall, b(51) =[-1.55357384000000000d+08]; end;
if firstCall, b(52) =[2.38904904000000000d+08]; end;
if firstCall, b(53) =[-1.62186912000000000d+08]; end;
if firstCall, b(54) =[4.43390400000000000d+07]; end;
if firstCall, b(55) =[-3.62880000000000000d+06]; end;
if firstCall, b(56) =[1.00000000000000000d+00]; end;
if firstCall, b(57) =[-4.07200000000000000d+03]; end;
if firstCall, b(58) =[6.95038000000000000d+05]; end;
if firstCall, b(59) =[-2.52439040000000000d+07]; end;
if firstCall, b(60) =[3.14369720000000000d+08]; end;
if firstCall, b(61) =[-1.64838430400000000d+09]; end;
if firstCall, b(62) =[4.00269508800000000d+09]; end;
if firstCall, b(63) =[-4.64216395200000000d+09]; end;
if firstCall, b(64) =[2.50748121600000000d+09]; end;
if firstCall, b(65) =[-5.68356480000000000d+08]; end;
if firstCall, b(66) =[3.99168000000000000d+07]; end;
if firstCall, b(67) =[1.00000000000000000d+00]; end;
if firstCall, b(68) =[-8.16600000000000000d+03]; end;
if firstCall, b(69) =[2.17062600000000000d+06]; end;
if firstCall, b(70) =[-1.14876376000000000d+08]; end;
if firstCall, b(71) =[2.05148277600000000d+09]; end;
if firstCall, b(72)=[-1.55489607840000000d+10]; end;
if firstCall, b(73) =[5.60413987840000000d+10]; end;
if firstCall, b(74) =[-1.01180433024000000d+11]; end;
if firstCall, b(75) =[9.21997902240000000d+10]; end;
if firstCall, b(76) =[-4.07883018240000000d+10]; end;
if firstCall, b(77) =[7.82771904000000000d+09]; end;
if firstCall, b(78) =[-4.79001600000000000d+08]; end;
if firstCall, b(79) =[1.00000000000000000d+00]; end;
if firstCall, b(80) =[-1.63560000000000000d+04]; end;
if firstCall, b(81) =[6.69969600000000000d+06]; end;
if firstCall, b(82) =[-5.07259276000000000d+08]; end;
if firstCall, b(83) =[1.26698177760000000d+10]; end;
if firstCall, b(84) =[-1.34323420224000000d+11]; end;
if firstCall, b(85) =[6.87720046384000000d+11]; end;
if firstCall, b(86) =[-1.81818864230400000d+12]; end;
if firstCall, b(87) =[2.54986547342400000d+12]; end;
if firstCall, b(88) =[-1.88307966182400000d+12]; end;
if firstCall, b(89) =[6.97929436800000000d+11]; end;
if firstCall, b(90) =[-1.15336085760000000d+11]; end;
if firstCall, b(91) =[6.22702080000000000d+09]; end;
if firstCall, b(92) =[1.00000000000000000d+00]; end;
if firstCall, b(93) =[-3.27380000000000000d+04]; end;
if firstCall, b(94) =[2.05079880000000000d+07]; end;
if firstCall, b(95) =[-2.18982980800000000d+09]; end;
if firstCall, b(96)=[7.50160522280000000d+10]; end;
if firstCall, b(97) =[-1.08467651241600000d+12]; end;
if firstCall, b(98) =[7.63483214939200000d+12]; end;
if firstCall, b(99) =[-2.82999100661120000d+13]; end;
if firstCall, b(100) =[5.74943734645920000d+13]; end;
if firstCall, b(101) =[-6.47283751398720000d+13]; end;
if firstCall, b(102) =[3.96895780558080000d+13]; end;
if firstCall, b(103) =[-1.25509040179200000d+13]; end;
if firstCall, b(104) =[1.81099255680000000d+12]; end;
if firstCall, b(105) =[-8.71782912000000000d+10]; end;
if firstCall, b(106) =[1.00000000000000000d+00]; end;
if firstCall, b(107) =[-6.55040000000000000d+04]; end;
if firstCall, b(108) =[6.24078900000000000d+07]; end;
if firstCall, b(109) =[-9.29252692000000000d+09]; end;
if firstCall, b(110)=[4.29826006340000000d+11]; end;
if firstCall, b(111) =[-8.30844432796800000d+12]; end;
if firstCall, b(112) =[7.83913848313120000d+13]; end;
if firstCall, b(113) =[-3.94365587815520000d+14]; end;
if firstCall, b(114) =[1.11174747256968000d+15]; end;
if firstCall, b(115) =[-1.79717122069056000d+15]; end;
if firstCall, b(116) =[1.66642448627145600d+15]; end;
if firstCall, b(117) =[-8.65023253219584000d+14]; end;
if firstCall, b(118)=[2.36908271543040000d+14]; end;
if firstCall, b(119) =[-3.01963769856000000d+13]; end;
if firstCall, b(120)=[1.30767436800000000d+12]; end;
%-----------------------------------------------------------------------
% BOUNDS B(M,K) , K=M-3
%-----------------------------------------------------------------------
if firstCall, bnd(1) =[1.0d0]; end;
if firstCall, bnd(2) =[1.0d0]; end;
if firstCall, bnd(3) =[1.0d0]; end;
if firstCall, bnd(4) =[1.0d0]; end;
if firstCall, bnd(5) =[3.10d0]; end;
if firstCall, bnd(6) =[5.18d0]; end;
if firstCall, bnd(7)=[11.7d0]; end;
if firstCall, bnd(8) =[29.8d0]; end;
if firstCall, bnd(9) =[90.4d0]; end;
if firstCall, bnd(10) =[297.0d0]; end;
if firstCall, bnd(11) =[1070.0d0]; end;
if firstCall, bnd(12) =[4290.0d0]; end;
if firstCall, bnd(13) =[18100.0d0]; end;
if firstCall, bnd(14) =[84700.0d0]; end;
if firstCall, bnd(15)=[408000.0d0]; end;
if firstCall, hrtpi=[8.86226925452758014d-01]; end;
firstCall=0;
%
%***FIRST EXECUTABLE STATEMENT DBKIAS
ierr = 0;
tol = max(d1mach(4),1.0d-18);
fln = n;
rz = 1.0d0./(x+fln);
rzx = x.*rz;
z = 0.5d0.*(x+fln);
if( ind<=1 )
[ gmrn ,z]=dgamrn(z);
end;
gs = hrtpi.*gmrn;
g1 = gs + gs;
rg1 = 1.0d0./g1;
gmrn =(rz+rz)./gmrn;
gt=0;
if( ind<=1 )
%-----------------------------------------------------------------------
% EVALUATE ERROR FOR M=MS
%-----------------------------------------------------------------------
hn = 0.5d0.*fln;
den2 = ktrms + ktrms + n;
den3 = den2 - 2.0d0;
den1 = x + den2;
err = rg1.*(x+x)./(den1-1.0d0);
if( n~=0 )
rat = 1.0d0./(fln.*fln);
end;
if( ktrms~=0 )
fj = ktrms;
rat = 0.25d0./(hrtpi.*den3.*sqrt(fj));
end;
err = err.*rat;
fj = -3.0d0;
for j = 1 : 15;
if( j<=5 )
err = err./den1;
end;
fm1 = max(1.0d0,fj);
fj = fj + 1.0d0;
er = bnd(j).*err;
if( ktrms==0 )
er = er.*(1.0d0+hn./fm1);
if( er<tol )
ms = fix(j);
gt=1;
break;
else;
if( j>=5 )
err = err./fln;
end;
end;
else;
er = er./fm1;
if( er<tol )
ms = fix(j);
gt=1;
break;
else;
if( j>=5 )
err = err./den3;
end;
end;
end;
end;
if(gt==0)
ierr = 2;
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end;
end;
mm = fix(ms + ms);
mp = fix(mm + 1);
%-----------------------------------------------------------------------
% H(K)=(-Z)**(K)*(PSI(K-1,Z)-PSI(K-1,Z+0.5))/GAMMA(K) , K=1,2,...,MM
%-----------------------------------------------------------------------
if( ind>1 )
rat = z./(z-0.5d0);
rxp = rat;
for i = 1 : mm;
h(i) = rxp.*(1.0d0-h(i));
rxp = rxp.*rat;
end; i = fix(mm+1);
else;
[z,mm,h,ierr]=dhkseq(z,mm,h,ierr);
end;
%-----------------------------------------------------------------------
% SCALED S SEQUENCE
%-----------------------------------------------------------------------
s(1) = 1.0d0;
fk = 1.0d0;
for k = 2 : mp;
ss = 0.0d0;
km = fix(k - 1);
i = fix(km);
for j = 1 : km;
ss = ss + s(j).*h(i);
i = fix(i - 1);
end; j = fix(km+1);
s(k) = ss./fk;
fk = fk + 1.0d0;
end; k = fix(mp+1);
%-----------------------------------------------------------------------
% SCALED S-TILDA SEQUENCE
%-----------------------------------------------------------------------
if( ktrms~=0 )
fk = 0.0d0;
ss = 0.0d0;
rg1 = rg1./z;
for k = 1 : ktrms;
v(k) = z./(z+fk);
w(k) = t(k).*v(k);
ss = ss + w(k);
fk = fk + 1.0d0;
end; k = fix(ktrms+1);
s(1) = s(1) - ss.*rg1;
for i = 2 : mp;
ss = 0.0d0;
for k = 1 : ktrms;
w(k) = w(k).*v(k);
ss = ss + w(k);
end; k = fix(ktrms+1);
s(i) = s(i) - ss.*rg1;
end; i = fix(mp+1);
end;
%-----------------------------------------------------------------------
% SUM ON J
%-----------------------------------------------------------------------
sumj = 0.0d0;
jn = 1;
rxp = 1.0d0;
xp(1) = 1.0d0;
for j = 1 : ms;
jn = fix(jn + j - 1);
xp(j+1) = xp(j).*rzx;
rxp = rxp.*rz;
%-----------------------------------------------------------------------
% SUM ON I
%-----------------------------------------------------------------------
sumi = 0.0d0;
ii = fix(jn);
for i = 1 : j;
jmi = fix(j - i + 1);
kk = fix(j + i + 1);
for k = 1 : jmi;
v(k) = s(kk).*xp(k);
kk = fix(kk + 1);
end; k = fix(jmi+1);
[jmi,v]=dbdiff(jmi,v);
sumi = sumi + b(ii).*v(jmi).*xp(i+1);
ii = fix(ii + 1);
end; i = fix(j+1);
sumj = sumj + sumi.*rxp;
end; j = fix(ms+1);
ansmlv = gs.*(s(1)-sumj);
h_shape=zeros(h_shape);h_shape(:)=h(1:numel(h_shape));h=h_shape;
return;
end %subroutine dbkias
%DECK DBKISR
|
|
