| [zr,fnu,n,y,nz,rz,ascle,tol,elim]=ckscl(zr,fnu,n,y,nz,rz,ascle,tol,elim); |
function [zr,fnu,n,y,nz,rz,ascle,tol,elim]=ckscl(zr,fnu,n,y,nz,rz,ascle,tol,elim);
%***BEGIN PROLOGUE CKSCL
%***SUBSIDIARY
%***PURPOSE Subsidiary to CBKNU, CUNK1 and CUNK2
%***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 CBKNU, CUNK1, CUNK2
%***ROUTINES CALLED CUCHK
%***REVISION HISTORY (YYMMDD)
% ?????? DATE WRITTEN
% 910415 Prologue converted to Version 4.0 format. (BAB)
%***end PROLOGUE CKSCL
persistent aa acs alas as celm ck cs csi csr cy czero elm firstCall fn helim i ic k kk nn nw s1 s2 xx zd zri ; if isempty(firstCall),firstCall=1;end;
if isempty(ck), ck=0; end;
if isempty(cs), cs=0; end;
if isempty(cy), cy=zeros(1,2); end;
if isempty(czero), czero=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(aa), aa=0; end;
if isempty(acs), acs=0; end;
if isempty(as), as=0; end;
if isempty(csi), csi=0; end;
if isempty(csr), csr=0; end;
if isempty(fn), fn=0; end;
if isempty(xx), xx=0; end;
if isempty(zri), zri=0; end;
if isempty(elm), elm=0; end;
if isempty(alas), alas=0; end;
if isempty(helim), helim=0; end;
if isempty(i), i=0; end;
if isempty(ic), ic=0; end;
if isempty(k), k=0; end;
if isempty(kk), kk=0; end;
if isempty(nn), nn=0; end;
if isempty(nw), nw=0; end;
if firstCall, czero=[complex(0.0e0,0.0e0)]; end;
firstCall=0;
%***FIRST EXECUTABLE STATEMENT CUCHK
nz = 0;
ic = 0;
xx = real(zr);
nn = fix(min(2,n));
for i = 1 : nn;
s1 = y(i);
cy(i) = s1;
as = abs(s1);
acs = -xx + log(as);
nz = fix(nz + 1);
y(i) = czero;
if( acs>=(-elim) )
cs = -zr + log(s1);
csr = real(cs);
csi = imag(cs);
aa = exp(csr)./tol;
cs = complex(aa,0.0e0).*complex(cos(csi),sin(csi));
[cs,nw,ascle,tol]=cuchk(cs,nw,ascle,tol);
if( nw==0 )
y(i) = cs;
nz = fix(nz - 1);
ic = fix(i);
end;
end;
end; i = fix(nn+1);
if( n==1 )
return;
end;
if( ic<=1 )
y(1) = czero;
nz = 2;
end;
if( n==2 )
return;
end;
if( nz==0 )
return;
end;
fn = fnu + 1.0e0;
ck = complex(fn,0.0e0).*rz;
s1 = cy(1);
s2 = cy(2);
helim = 0.5e0.*elim;
elm = exp(-elim);
celm = complex(elm,0.0e0);
zri = imag(zr);
zd = zr;
%
% FIND TWO CONSECUTIVE Y VALUES ON SCALE. SCALE RECURRENCE IF
% S2 GETS LARGER THAN EXP(ELIM/2)
%
for i = 3 : n;
kk = fix(i);
cs = s2;
s2 = ck.*s2 + s1;
s1 = cs;
ck = ck + rz;
as = abs(s2);
alas = log(as);
acs = -xx + alas;
nz = fix(nz + 1);
y(i) = czero;
if( acs>=(-elim) )
cs = -zd + log(s2);
csr = real(cs);
csi = imag(cs);
aa = exp(csr)./tol;
cs = complex(aa,0.0e0).*complex(cos(csi),sin(csi));
[cs,nw,ascle,tol]=cuchk(cs,nw,ascle,tol);
if( nw==0 )
y(i) = cs;
nz = fix(nz - 1);
if( ic==(kk-1) )
nz = fix(kk - 2);
for k = 1 : nz;
y(k) = czero;
end; k = fix(nz+1);
return;
else;
ic = fix(kk);
continue;
end;
end;
end;
if( alas>=helim )
xx = xx - elim;
s1 = s1.*celm;
s2 = s2.*celm;
zd = complex(xx,zri);
end;
end; i = fix(n+1);
nz = fix(n);
if( ic==n )
nz = fix(n - 1);
end;
for k = 1 : nz;
y(k) = czero;
end; k = fix(nz+1);
end
%DECK CLBETA
|
|