| [x,n,kode,m,tol,en,nz,ierr]=exint(x,n,kode,m,tol,en,nz,ierr); |
function [x,n,kode,m,tol,en,nz,ierr]=exint(x,n,kode,m,tol,en,nz,ierr);
%***BEGIN PROLOGUE EXINT
%***PURPOSE Compute an M member sequence of exponential integrals
% E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0.
%***LIBRARY SLATEC
%***CATEGORY C5
%***TYPE SINGLE PRECISION (EXINT-S, DEXINT-D)
%***KEYWORDS EXPONENTIAL INTEGRAL, SPECIAL FUNCTIONS
%***AUTHOR Amos, D. E., (SNLA)
%***DESCRIPTION
%
% EXINT computes M member sequences of exponential integrals
% E(N+K,X), K=0,1,...,M-1 for N .GE. 1 and X .GE. 0. The
% exponential integral is defined by
%
% E(N,X)=integral on (1,infinity) of EXP(-XT)/T**N
%
% where X=0.0 and N=1 cannot occur simultaneously. Formulas
% and notation are found in the NBS Handbook of Mathematical
% Functions (ref. 1).
%
% The power series is implemented for X .LE. XCUT and the
% confluent hypergeometric representation
%
% E(A,X) = EXP(-X)*(X**(A-1))*U(A,A,X)
%
% is computed for X .GT. XCUT. Since sequences are computed in
% a stable fashion by recurring away from X, A is selected as
% the integer closest to X within the constraint N .LE. A .LE.
% N+M-1. For the U computation, A is further modified to be the
% nearest even integer. Indices are carried forward or
% backward by the two term recursion relation
%
% K*E(K+1,X) + X*E(K,X) = EXP(-X)
%
% once E(A,X) is computed. The U function is computed by means
% of the backward recursive Miller algorithm applied to the
% three term contiguous relation for U(A+K,A,X), K=0,1,...
% This produces accurate ratios and determines U(A+K,A,X), and
% hence E(A,X), to within a multiplicative constant C.
% Another contiguous relation applied to C*U(A,A,X) and
% C*U(A+1,A,X) gets C*U(A+1,A+1,X), a quantity proportional to
% E(A+1,X). The normalizing constant C is obtained from the
% two term recursion relation above with K=A.
%
% Description of Arguments
%
% Input
% X X .GT. 0.0 for N=1 and X .GE. 0.0 for N .GE. 2
% N order of the first member of the sequence, N .GE. 1
% (X=0.0 and N=1 is an error)
% KODE a selection parameter for scaled values
% KODE=1 returns E(N+K,X), K=0,1,...,M-1.
% =2 returns EXP(X)*E(N+K,X), K=0,1,...,M-1.
% M number of exponential integrals in the sequence,
% M .GE. 1
% TOL relative accuracy wanted, ETOL .LE. TOL .LE. 0.1
% ETOL = single precision unit roundoff = R1MACH(4)
%
% Output
% EN a vector of dimension at least M containing values
% EN(K) = E(N+K-1,X) or EXP(X)*E(N+K-1,X), K=1,M
% depending on KODE
% NZ underflow indicator
% NZ=0 a normal return
% NZ=M X exceeds XLIM and an underflow occurs.
% EN(K)=0.0E0 , K=1,M returned on KODE=1
% IERR error flag
% IERR=0, normal return, computation completed
% IERR=1, input error, no computation
% IERR=2, error, no computation
% algorithm termination condition not met
%
%***REFERENCES M. Abramowitz and I. A. Stegun, Handbook of
% Mathematical Functions, NBS AMS Series 55, U.S. Dept.
% of Commerce, 1955.
% D. E. Amos, Computation of exponential integrals, ACM
% Transactions on Mathematical Software 6, (1980),
% pp. 365-377 and pp. 420-428.
%***ROUTINES CALLED I1MACH, PSIXN, R1MACH
%***REVISION HISTORY (YYMMDD)
% 800501 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)
% 900326 Removed duplicate information from DESCRIPTION section.
% (WRB)
% 910408 Updated the REFERENCES section. (WRB)
% 920207 Updated with code with a revision date of 880811 from
% D. Amos. Included correction of argument list. (WRB)
% 920501 Reformatted the REFERENCES section. (WRB)
%***end PROLOGUE EXINT
persistent a aa aams ah ak at b bk bt cc cnorm ct em emx etol fnm fx gt gt2 i i1m ic icase ict ik ind ix jset k kk kn ks ml mu nd nm p1 p2 pt s tx xcut xlim xtol y y1 y2 yt ;
if isempty(a), a=zeros(1,99); end;
if isempty(aa), aa=0; end;
if isempty(aams), aams=0; end;
if isempty(ah), ah=0; end;
if isempty(ak), ak=0; end;
if isempty(at), at=0; end;
if isempty(b), b=zeros(1,99); end;
if isempty(bk), bk=0; end;
if isempty(bt), bt=0; end;
if isempty(cc), cc=0; end;
if isempty(cnorm), cnorm=0; end;
if isempty(ct), ct=0; end;
if isempty(em), em=0; end;
if isempty(emx), emx=0; end;
if isempty(etol), etol=0; end;
if isempty(fnm), fnm=0; end;
if isempty(fx), fx=0; end;
if isempty(pt), pt=0; end;
if isempty(p1), p1=0; end;
if isempty(p2), p2=0; end;
if isempty(s), s=0; end;
if isempty(tx), tx=0; end;
if isempty(xcut), xcut=0; end;
if isempty(xlim), xlim=0; end;
if isempty(xtol), xtol=0; end;
if isempty(y), y=zeros(1,2); end;
if isempty(yt), yt=0; end;
if isempty(y1), y1=0; end;
if isempty(y2), y2=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(icase), icase=0; end;
if isempty(ict), ict=0; end;
if isempty(ik), ik=0; end;
if isempty(ind), ind=0; end;
if isempty(ix), ix=0; end;
if isempty(i1m), i1m=0; end;
if isempty(jset), jset=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(kn), kn=0; end;
if isempty(ks), ks=0; end;
if isempty(ml), ml=0; end;
if isempty(mu), mu=0; end;
if isempty(nd), nd=0; end;
if isempty(nm), nm=0; end;
if isempty(gt), gt=0; end;
if isempty(gt2), gt2=zeros(1,2); end;
en_shape=size(en);en=reshape(en,1,[]);
%***FIRST EXECUTABLE STATEMENT EXINT
ierr = 0;
nz = 0;
etol = max(r1mach(4),0.5e-18);
if( x<0.0e0 )
ierr = 1;
end;
if( n<1 )
ierr = 1;
end;
if( kode<1 || kode>2 )
ierr = 1;
end;
if( m<1 )
ierr = 1;
end;
if( tol<etol || tol>0.1e0 )
ierr = 1;
end;
if( x==0.0e0 && n==1 )
ierr = 1;
end;
if( ierr==0 )
i1m = fix(-i1mach(12));
pt = 2.3026e0.*r1mach(5).*i1m;
xlim = pt - 6.907755e0;
bt = pt +(n+m-1);
if( bt>1000.0e0 )
xlim = pt - log(bt);
end;
%
xcut = 2.0e0;
if( etol>2.0e-7 )
xcut = 1.0e0;
end;
while (1);
gt2(:)=0;
if( x>xcut )
%-----------------------------------------------------------------------
% BACKWARD RECURSIVE MILLER ALGORITHM FOR
% E(N,X)=EXP(-X)*(X**(N-1))*U(N,N,X)
% WITH RECURSION AWAY FROM N=INTEGER CLOSEST TO X.
% U(A,B,X) IS THE SECOND CONFLUENT HYPERGEOMETRIC FUNCTION
%-----------------------------------------------------------------------
emx = 1.0e0;
if( kode~=2 )
if( x<=xlim )
emx = exp(-x);
else;
nz = fix(m);
for i = 1 : m;
en(i) = 0.0e0;
end; i = fix(m+1);
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
end;
ix = fix(x + 0.5e0);
kn = fix(n + m - 1);
gt=0;
while (1);
if( kn<=ix )
icase = 1;
ks = fix(kn);
ml = fix(m - 1);
mu = -1;
ind = fix(m);
if( kn>1 )
gt=1;
break;
end;
else;
if( n>=ix || ix>=kn )
if( n>=ix )
icase = 2;
ind = 1;
ks = fix(n);
mu = fix(m - 1);
if( n>1 )
gt=1;
break;
end;
if( kn==1 )
break;
end;
ix = 2;
else;
ierr = 2;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
end;
icase = 1;
ks = fix(ix);
ml = fix(ix - n);
ind = fix(ml + 1);
mu = fix(kn - ix);
gt=1;
end;
break;
end;
if(gt==0)
ks = 2;
icase = 3;
end;
ik = fix(fix(ks./2));
ah = ik;
jset = fix(1 + ks -(ik+ik));
%-----------------------------------------------------------------------
% START COMPUTATION FOR
% EN(IND) = C*U( A , A ,X) JSET=1
% EN(IND) = C*U(A+1,A+1,X) JSET=2
% FOR AN EVEN INTEGER A.
%-----------------------------------------------------------------------
ic = 0;
aa = ah + ah;
aams = aa - 1.0e0;
aams = aams.*aams;
tx = x + x;
fx = tx + tx;
ak = ah;
xtol = tol;
if( tol<=1.0e-3 )
xtol = 20.0e0.*tol;
end;
ct = aams + fx.*ah;
em =(ah+1.0e0)./((x+aa).*xtol.*sqrt(ct));
bk = aa;
cc = ah.*ah;
%-----------------------------------------------------------------------
% FORWARD RECURSION FOR P(IC),P(IC+1) AND INDEX IC FOR BACKWARD
% RECURSION
%-----------------------------------------------------------------------
p1 = 0.0e0;
p2 = 1.0e0;
while( ic~=99 );
ic = fix(ic + 1);
ak = ak + 1.0e0;
at = bk./(bk+ak+cc+ic);
bk = bk + ak + ak;
a(ic) = at;
bt =(ak+ak+x)./(ak+1.0e0);
b(ic) = bt;
pt = p2;
p2 = bt.*p2 - at.*p1;
p1 = pt;
ct = ct + fx;
em = em.*at.*(1.0e0-tx./ct);
if( em.*(ak+1.0e0)<=p1.*p1 )
ict = fix(ic);
kk = fix(ic + 1);
bt = tx./(ct+fx);
y2 =(bk./(bk+cc+kk)).*(p1./p2).*(1.0e0-bt+0.375e0.*bt.*bt);
y1 = 1.0e0;
%-----------------------------------------------------------------------
% BACKWARD RECURRENCE FOR
% Y1= C*U( A ,A,X)
% Y2= C*(A/(1+A/2))*U(A+1,A,X)
%-----------------------------------------------------------------------
for k = 1 : ict;
kk = fix(kk - 1);
yt = y1;
y1 =(b(kk).*y1-y2)./a(kk);
y2 = yt;
end; k = fix(ict+1);
%-----------------------------------------------------------------------
% THE CONTIGUOUS RELATION
% X*U(B,C+1,X)=(C-B)*U(B,C,X)+U(B-1,C,X)
% WITH B=A+1 , C=A IS USED FOR
% Y(2) = C * U(A+1,A+1,X)
% X IS INCORPORATED INTO THE NORMALIZING RELATION
%-----------------------------------------------------------------------
pt = y2./y1;
cnorm = 1.0e0 - pt.*(ah+1.0e0)./aa;
y(1) = 1.0e0./(cnorm.*aa+x);
y(2) = cnorm.*y(1);
if( icase~=3 )
en(ind) = emx.*y(jset);
if( m==1 )
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
aa = ks;
if( icase==1 )
gt2(1)=1;
break;
end;
if( icase==2 )
gt2(2)=1;
break;
end;
end;
%-----------------------------------------------------------------------
% RECURSION SECTION N*E(N+1,X) + X*E(N,X)=EMX
%-----------------------------------------------------------------------
en(1) = emx.*(1.0e0-y(1))./x;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
end;
if(any(gt2==1))
break;
end;
ierr = 2;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
elseif( x==0.0e0 && n>1 ) ;
for i = 1 : m;
en(i) = 1.0e0./(n+i-2);
end; i = fix(m+1);
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
else;
%-----------------------------------------------------------------------
% SERIES FOR E(N,X) FOR X.LE.XCUT
%-----------------------------------------------------------------------
tx = x + 0.5e0;
ix = fix(tx);
%-----------------------------------------------------------------------
% ICASE=1 MEANS INTEGER CLOSEST TO X IS 2 AND N=1
% ICASE=2 MEANS INTEGER CLOSEST TO X IS 0,1, OR 2 AND N.GE.2
%-----------------------------------------------------------------------
icase = 2;
if( ix>n )
icase = 1;
end;
nm = fix(n - icase + 1);
nd = fix(nm + 1);
ind = fix(3 - icase);
mu = fix(m - ind);
ml = 1;
ks = fix(nd);
fnm = nm;
s = 0.0e0;
xtol = 3.0e0.*tol;
if( nd~=1 )
xtol = 0.3333e0.*tol;
s = 1.0e0./fnm;
end;
aa = 1.0e0;
ak = 1.0e0;
ic = 35;
if( x<etol )
ic = 1;
end;
gt=0;
for i = 1 : ic;
aa = -aa.*x./ak;
if( i==nm )
s = s + aa.*(-log(x)+psixn(nd));
xtol = 3.0e0.*tol;
else;
s = s - aa./(ak-fnm);
if( abs(aa)>xtol.*abs(s) )
ak = ak + 1.0e0;
continue;
elseif( i>=2 ) ;
if( nd-2>i || i>nd-1 )
gt=1;
break;
end;
ak = ak + 1.0e0;
continue;
end;
end;
ak = ak + 1.0e0;
end;
if( ic~=1 && gt==0)
ierr = 2;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
if( nd==1 )
s = s +(-log(x)+psixn(1));
end;
if( kode==2 )
s = s.*exp(x);
end;
en(1) = s;
emx = 1.0e0;
if( m~=1 )
en(ind) = s;
aa = ks;
if( kode==1 )
emx = exp(-x);
end;
%to100
if( icase==1 )
break;
end;
if( icase==2 )
gt2(2)=1;
break;
end;
end;
if( icase~=2 )
if( kode==1 )
emx = exp(-x);
end;
en(1) =(emx-s)./x;
end;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
break;
end;
if(gt2(2)==0)
k = fix(ind - 1);
for i = 1 : ml;
aa = aa - 1.0e0;
en(k) =(emx-aa.*en(k+1))./x;
k = fix(k - 1);
end; i = fix(ml+1);
if( mu<=0 )
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
return;
end;
aa = ks;
end;
k = fix(ind);
for i = 1 : mu;
en(k+1) =(emx-x.*en(k))./aa;
aa = aa + 1.0e0;
k = fix(k + 1);
end; i = fix(mu+1);
end;
en_shape=zeros(en_shape);en_shape(:)=en(1:numel(en_shape));en=en_shape;
end %subroutine exint
%DECK EXPREL
|
|
