[zr,fnu,kode,n,y,nz,cw,tol,elim,alim]=cwrsk(zr,fnu,kode,n,y,nz,cw,tol,elim,alim);

function [zr,fnu,kode,n,y,nz,cw,tol,elim,alim]=cwrsk(zr,fnu,kode,n,y,nz,cw,tol,elim,alim);
%***BEGIN PROLOGUE  CWRSK
%***SUBSIDIARY
%***PURPOSE  Subsidiary to CBESI and CBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CWRSK-A, ZWRSK-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
%     NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
%
%***SEE ALSO  CBESI, CBESK
%***ROUTINES CALLED  CBKNU, CRATI, R1MACH
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  CWRSK
persistent act acw ascle c1 c2 cinu cscl ct i nw rct s1 s2 st yy ; 

if isempty(cinu), cinu=0; end;
if isempty(cscl), cscl=0; end;
if isempty(ct), ct=0; end;
if isempty(c1), c1=0; end;
if isempty(c2), c2=0; end;
if isempty(rct), rct=0; end;
if isempty(st), st=0; end;
if isempty(act), act=0; end;
if isempty(acw), acw=0; end;
if isempty(ascle), ascle=0; end;
if isempty(s1), s1=0; end;
if isempty(s2), s2=0; end;
if isempty(yy), yy=0; end;
if isempty(i), i=0; end;
if isempty(nw), nw=0; end;
%***FIRST EXECUTABLE STATEMENT  CWRSK
%-----------------------------------------------------------------------
%     I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
%     Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
%     WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
%-----------------------------------------------------------------------
nz = 0;
[zr,fnu,kode,dumvar4,cw,nw,tol,elim,alim]=cbknu(zr,fnu,kode,2,cw,nw,tol,elim,alim);
if( nw~=0 )
nz = -1;
if( nw==(-2) )
nz = -2;
end;
else;
[zr,fnu,n,y,tol]=crati(zr,fnu,n,y,tol);
%-----------------------------------------------------------------------
%     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
%     R(FNU+J-1,Z)=Y(J),  J=1,...,N
%-----------------------------------------------------------------------
cinu = complex(1.0e0,0.0e0);
if( kode~=1 )
yy = imag(zr);
s1 = cos(yy);
s2 = sin(yy);
cinu = complex(s1,s2);
end;
%-----------------------------------------------------------------------
%     ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
%     THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
%     SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
%     THE RESULT IS ON SCALE.
%-----------------------------------------------------------------------
acw = abs(cw(2));
ascle = 1.0e+3.*r1mach(1)./tol;
cscl = complex(1.0e0,0.0e0);
if( acw>ascle )
ascle = 1.0e0./ascle;
if( acw>=ascle )
cscl = complex(tol,0.0e0);
end;
else;
cscl = complex(1.0e0./tol,0.0e0);
end;
c1 = cw(1).*cscl;
c2 = cw(2).*cscl;
st = y(1);
%-----------------------------------------------------------------------
%     CINU=CINU*(CONJG(CT)/ABS(CT))*(1.0E0/ABS(CT) PREVENTS
%     UNDER- OR OVERFLOW PREMATURELY BY SQUARING ABS(CT)
%-----------------------------------------------------------------------
ct = zr.*(c2+st.*c1);
act = abs(ct);
rct = complex(1.0e0./act,0.0e0);
ct = conj(ct).*rct;
cinu = cinu.*rct.*ct;
y(1) = cinu.*cscl;
if( n==1 )
return;
end;
for i = 2 : n;
cinu = st.*cinu;
st = y(i);
y(i) = cinu.*cscl;
end; i = fix(n+1);
return;
end;
end