[zrr,zri,fnu,kode,n,yr,yi,nz,cwr,cwi,tol,elim,alim]=zwrsk(zrr,zri,fnu,kode,n,yr,yi,nz,cwr,cwi,tol,elim,alim);

function [zrr,zri,fnu,kode,n,yr,yi,nz,cwr,cwi,tol,elim,alim]=zwrsk(zrr,zri,fnu,kode,n,yr,yi,nz,cwr,cwi,tol,elim,alim);
%***BEGIN PROLOGUE  ZWRSK
%***SUBSIDIARY
%***PURPOSE  Subsidiary to ZBESI and ZBESK
%***LIBRARY   SLATEC
%***TYPE      ALL (CWRSK-A, ZWRSK-A)
%***AUTHOR  Amos, D. E., (SNL)
%***DESCRIPTION
%
%     ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
%     NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN
%
%***SEE ALSO  ZBESI, ZBESK
%***ROUTINES CALLED  D1MACH, ZABS, ZBKNU, ZRATI
%***REVISION HISTORY  (YYMMDD)
%   830501  DATE WRITTEN
%   910415  Prologue converted to Version 4.0 format.  (BAB)
%***end PROLOGUE  ZWRSK
%     COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR
persistent act acw ascle c1i c1r c2i c2r cinui cinur csclr cti ctr i nw pti ptr ract sti str ; 

if isempty(act), act=0; end;
if isempty(acw), acw=0; end;
if isempty(ascle), ascle=0; end;
if isempty(cinui), cinui=0; end;
if isempty(cinur), cinur=0; end;
if isempty(csclr), csclr=0; end;
if isempty(cti), cti=0; end;
if isempty(ctr), ctr=0; end;
if isempty(c1i), c1i=0; end;
if isempty(c1r), c1r=0; end;
if isempty(c2i), c2i=0; end;
if isempty(c2r), c2r=0; end;
if isempty(pti), pti=0; end;
if isempty(ptr), ptr=0; end;
if isempty(ract), ract=0; end;
if isempty(sti), sti=0; end;
if isempty(str), str=0; end;
if isempty(i), i=0; end;
if isempty(nw), nw=0; end;
%***FIRST EXECUTABLE STATEMENT  ZWRSK
%-----------------------------------------------------------------------
%     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;
[zrr,zri,fnu,kode,dumvar5,cwr,cwi,nw,tol,elim,alim]=zbknu(zrr,zri,fnu,kode,2,cwr,cwi,nw,tol,elim,alim);
if( nw~=0 )
nz = -1;
if( nw==(-2) )
nz = -2;
end;
else;
[zrr,zri,fnu,n,yr,yi,tol]=zrati(zrr,zri,fnu,n,yr,yi,tol);
%-----------------------------------------------------------------------
%     RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
%     R(FNU+J-1,Z)=Y(J),  J=1,...,N
%-----------------------------------------------------------------------
cinur = 1.0d0;
cinui = 0.0d0;
if( kode~=1 )
cinur = cos(zri);
cinui = sin(zri);
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 ,cwr(2),cwi(2)]=zabs(cwr(2),cwi(2));
ascle = 1.0d+3.*d1mach(1)./tol;
csclr = 1.0d0;
if( acw>ascle )
ascle = 1.0d0./ascle;
if( acw>=ascle )
csclr = tol;
end;
else;
csclr = 1.0d0./tol;
end;
c1r = cwr(1).*csclr;
c1i = cwi(1).*csclr;
c2r = cwr(2).*csclr;
c2i = cwi(2).*csclr;
str = yr(1);
sti = yi(1);
%-----------------------------------------------------------------------
%     CINU=CINU*(CONJG(CT)/ABS(CT))*(1.0D0/ABS(CT) PREVENTS
%     UNDER- OR OVERFLOW PREMATURELY BY SQUARING ABS(CT)
%-----------------------------------------------------------------------
ptr = str.*c1r - sti.*c1i;
pti = str.*c1i + sti.*c1r;
ptr = ptr + c2r;
pti = pti + c2i;
ctr = zrr.*ptr - zri.*pti;
cti = zrr.*pti + zri.*ptr;
[act ,ctr,cti]=zabs(ctr,cti);
ract = 1.0d0./act;
ctr = ctr.*ract;
cti = -cti.*ract;
ptr = cinur.*ract;
pti = cinui.*ract;
cinur = ptr.*ctr - pti.*cti;
cinui = ptr.*cti + pti.*ctr;
yr(1) = cinur.*csclr;
yi(1) = cinui.*csclr;
if( n==1 )
return;
end;
for i = 2 : n;
ptr = str.*cinur - sti.*cinui;
cinui = str.*cinui + sti.*cinur;
cinur = ptr;
str = yr(i);
sti = yi(i);
yr(i) = cinur.*csclr;
yi(i) = cinui.*csclr;
end; i = fix(n+1);
return;
end;
end