| [zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim]=zkscl(zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim); |
function [zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim]=zkscl(zrr,zri,fnu,n,yr,yi,nz,rzr,rzi,ascle,tol,elim);
%***BEGIN PROLOGUE ZKSCL
%***SUBSIDIARY
%***PURPOSE Subsidiary to ZBESK
%***LIBRARY SLATEC
%***TYPE ALL (CKSCL-A, ZKSCL-A)
%***AUTHOR Amos, D. E., (SNL)
%***DESCRIPTION
%
% SET K FUNCTIONS TO ZERO ON UNDERFLOW, CONTINUE RECURRENCE
% ON SCALED FUNCTIONS UNTIL TWO MEMBERS COME ON SCALE, THEN
% RETURN WITH MIN(NZ+2,N) VALUES SCALED BY 1/TOL.
%
%***SEE ALSO ZBESK
%***ROUTINES CALLED ZABS, ZLOG, ZUCHK
%***REVISION HISTORY (YYMMDD)
% 830501 DATE WRITTEN
% 910415 Prologue converted to Version 4.0 format. (BAB)
% 930122 Added ZLOG to EXTERNAL statement. (RWC)
%***end PROLOGUE ZKSCL
% COMPLEX CK,CS,CY,CZERO,RZ,S1,S2,Y,ZR,ZD,CELM
persistent acs alas as celmr cki ckr csi csr cyi cyr elm firstCall fn helim i ic idum igo kk nn nw s1i s1r s2i s2r str zdi zdr zeroi zeror ; if isempty(firstCall),firstCall=1;end;
if isempty(acs), acs=0; end;
if isempty(as), as=0; end;
if isempty(cki), cki=0; end;
if isempty(ckr), ckr=0; end;
if isempty(csi), csi=0; end;
if isempty(csr), csr=0; end;
if isempty(cyi), cyi=zeros(1,2); end;
if isempty(cyr), cyr=zeros(1,2); end;
if isempty(fn), fn=0; end;
if isempty(str), str=0; end;
if isempty(s1i), s1i=0; end;
if isempty(s1r), s1r=0; end;
if isempty(s2i), s2i=0; end;
if isempty(s2r), s2r=0; end;
if isempty(zeroi), zeroi=0; end;
if isempty(zeror), zeror=0; end;
if isempty(zdr), zdr=0; end;
if isempty(zdi), zdi=0; end;
if isempty(celmr), celmr=0; end;
if isempty(elm), elm=0; end;
if isempty(helim), helim=0; end;
if isempty(alas), alas=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(idum), idum=0; end;
if isempty(kk), kk=0; end;
if isempty(nn), nn=0; end;
if isempty(nw), nw=0; end;
if isempty(igo), igo=0; end;
if firstCall, zeror =[0.0d0]; end;
if firstCall, zeroi=[0.0d0]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT ZKSCL
nz = 0;
ic = 0;
nn = fix(min(2,n));
for i = 1 : nn;
s1r = yr(i);
s1i = yi(i);
cyr(i) = s1r;
cyi(i) = s1i;
[as ,s1r,s1i]=zabs(s1r,s1i);
acs = -zrr + log(as);
nz = fix(nz + 1);
yr(i) = zeror;
yi(i) = zeroi;
if( acs>=(-elim) )
[s1r,s1i,csr,csi,idum]=zlog(s1r,s1i,csr,csi,idum);
csr = csr - zrr;
csi = csi - zri;
str = exp(csr)./tol;
csr = str.*cos(csi);
csi = str.*sin(csi);
[csr,csi,nw,ascle,tol]=zuchk(csr,csi,nw,ascle,tol);
if( nw==0 )
yr(i) = csr;
yi(i) = csi;
ic = fix(i);
nz = fix(nz - 1);
end;
end;
end; i = fix(nn+1);
if( n==1 )
return;
end;
if( ic<=1 )
yr(1) = zeror;
yi(1) = zeroi;
nz = 2;
end;
if( n==2 )
return;
end;
if( nz==0 )
return;
end;
fn = fnu + 1.0d0;
ckr = fn.*rzr;
cki = fn.*rzi;
s1r = cyr(1);
s1i = cyi(1);
s2r = cyr(2);
s2i = cyi(2);
helim = 0.5d0.*elim;
elm = exp(-elim);
celmr = elm;
zdr = zrr;
zdi = zri;
%
% FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
% S2 GETS LARGER THAN EXP(ELIM/2)
%
igo=0;
for i = 3 : n;
kk = fix(i);
csr = s2r;
csi = s2i;
s2r = ckr.*csr - cki.*csi + s1r;
s2i = cki.*csr + ckr.*csi + s1i;
s1r = csr;
s1i = csi;
ckr = ckr + rzr;
cki = cki + rzi;
[as ,s2r,s2i]=zabs(s2r,s2i);
alas = log(as);
acs = -zdr + alas;
nz = fix(nz + 1);
yr(i) = zeror;
yi(i) = zeroi;
if( acs>=(-elim) )
[s2r,s2i,csr,csi,idum]=zlog(s2r,s2i,csr,csi,idum);
csr = csr - zdr;
csi = csi - zdi;
str = exp(csr)./tol;
csr = str.*cos(csi);
csi = str.*sin(csi);
[csr,csi,nw,ascle,tol]=zuchk(csr,csi,nw,ascle,tol);
if( nw==0 )
yr(i) = csr;
yi(i) = csi;
nz = fix(nz - 1);
if( ic==kk-1 )
nz = fix(kk - 2);
igo=1;
break;
else;
ic = fix(kk);
continue;
end;
end;
end;
if( alas>=helim )
zdr = zdr - elim;
s1r = s1r.*celmr;
s1i = s1i.*celmr;
s2r = s2r.*celmr;
s2i = s2i.*celmr;
end;
end;
if(igo==0)
nz = fix(n);
if( ic==n )
nz = fix(n - 1);
end;
end;
for i = 1 : nz;
yr(i) = zeror;
yi(i) = zeroi;
end; i = fix(nz+1);
end %subroutine zkscl
%DECK ZLOG
|
|