[lun,kprint,ipass]=cqcby(lun,kprint,ipass);

function [lun,kprint,ipass]=cqcby(lun,kprint,ipass);
persistent aa ab aer alim arg atol av cc ci cip coe1 coe2 csgn cspn ct cw cwrk dig elim eps er ertol ffnu film firstCall fnu fnul hpi i i4 icase ierr ifnu il ir irb it itl k k1 k2 kdo keps kk kode lflg mflg mqc n nl nu nul nz nz1 nz2 pi r r1m4 r1m5 r2 rhpi rl rm slak st t tol ts v w xn xnu xx y yn yy z zn ; if isempty(firstCall),firstCall=1;end; 

if isempty(mqc), mqc=1 ; end;
if isempty(ci), ci=0; end;
if isempty(cip), cip=zeros(1,4); end;
if isempty(coe1), coe1=0; end;
if isempty(coe2), coe2=0; end;
if isempty(csgn), csgn=0; end;
if isempty(cspn), cspn=0; end;
if isempty(cw), cw=0; end;
if isempty(cwrk), cwrk=zeros(1,20); end;
if isempty(v), v=zeros(1,20); end;
if isempty(w), w=zeros(1,20); end;
if isempty(y), y=zeros(1,20); end;
if isempty(z), z=0; end;
if isempty(zn), zn=0; end;
if isempty(aa), aa=0; end;
if isempty(ab), ab=0; end;
if isempty(aer), aer=zeros(1,20); end;
if isempty(alim), alim=0; end;
if isempty(arg), arg=0; end;
if isempty(atol), atol=0; end;
if isempty(av), av=0; end;
if isempty(cc), cc=0; end;
if isempty(ct), ct=0; end;
if isempty(dig), dig=0; end;
if isempty(elim), elim=0; end;
if isempty(eps), eps=0; end;
if isempty(er), er=0; end;
if isempty(ertol), ertol=0; end;
if isempty(ffnu), ffnu=0; end;
if isempty(film), film=0; end;
if isempty(fnu), fnu=0; end;
if isempty(fnul), fnul=0; end;
if isempty(hpi), hpi=0; end;
if isempty(pi), pi=0; end;
if isempty(r), r=0; end;
if isempty(rhpi), rhpi=0; end;
if isempty(rl), rl=0; end;
if isempty(rm), rm=0; end;
if isempty(r1m4), r1m4=0; end;
if isempty(r1m5), r1m5=0; end;
if isempty(r2), r2=0; end;
if isempty(slak), slak=0; end;
if isempty(st), st=0; end;
if isempty(t), t=zeros(1,20); end;
if isempty(tol), tol=0; end;
if isempty(ts), ts=0; end;
if isempty(xn), xn=0; end;
if isempty(xnu), xnu=zeros(1,20); end;
if isempty(xx), xx=0; end;
if isempty(yn), yn=0; end;
if isempty(yy), yy=0; end;
if isempty(i), i=0; end;
if isempty(icase), icase=0; end;
if isempty(ierr), ierr=0; end;
if isempty(ifnu), ifnu=0; end;
if isempty(il), il=0; end;
if isempty(ir), ir=0; end;
if isempty(irb), irb=0; end;
if isempty(it), it=0; end;
if isempty(itl), itl=0; end;
if isempty(i4), i4=0; end;
if isempty(k), k=0; end;
if isempty(kdo), kdo=zeros(1,20); end;
if isempty(keps), keps=zeros(1,20); end;
if isempty(kk), kk=0; end;
if isempty(kode), kode=0; end;
if isempty(k1), k1=0; end;
if isempty(k2), k2=0; end;
if isempty(lflg), lflg=0; end;
if isempty(mflg), mflg=0; end;
if isempty(n), n=0; end;
if isempty(nl), nl=0; end;
if isempty(nu), nu=0; end;
if isempty(nul), nul=0; end;
if isempty(nz), nz=0; end;
if isempty(nz1), nz1=0; end;
if isempty(nz2), nz2=0; end;
if firstCall,   cip(1) =[complex(1.0e0,0.0e0)];  end;
if firstCall,  cip(2) =[complex(0.0e0,1.0e0)];  end;
if firstCall,  cip(3) =[complex(-1.0e0,0.0e0)];  end;
if firstCall,  cip(4)=[complex(0.0e0,-1.0e0)];  end;
firstCall=0;
if( kprint>=2 )
writef(lun,[' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ','CBESY', '\n '  ' \n']);
%format (' QUICK CHECK ROUTINE FOR THE Y BESSEL FUNCTION FROM ','CBESY'];
end;
[r1m4 ]=r1mach(4);
tol = max(r1m4,1.0e-18);
atol = 100.0e0.*tol;
aa = -log10(r1m4);
[k1 ]=i1mach(12);
[k2 ]=i1mach(13);
[r1m5 ]=r1mach(5);
k = fix(min(abs(k1),abs(k2)));
elim = 2.303e0.*(k.*r1m5-3.0e0);
ab = aa.*2.303e0;
alim = elim + max(-ab,-41.45e0);
dig = min(aa,18.0e0);
fnul = 10.0e0 + 6.0e0.*(dig-3.0e0);
rl = 1.2e0.*dig + 3.0e0;
slak = 3.0e0 + 4.0e0.*(-log10(tol)-7.0e0)./11.0e0;
slak = max(slak,3.0e0);
ertol = tol.*10.0e0.^slak;
rm = 0.5e0.*(alim+elim);
rm = min(rm,200.0e0);
rm = max(rm,rl+10.0e0);
r2 = min(rm,fnul);
if( kprint>=2 )
writef(lun,[' PARAMETERS', '\n ' ,repmat(' ',1,5),'TOL ',repmat(' ',1,8),'ELIM',repmat(' ',1,8),'ALIM',repmat(' ',1,8),'RL  ',repmat(' ',1,8),'FNUL',repmat(' ',1,8),'DIG' ' \n']);
%format (' PARAMETERS'/5X,'TOL ',8X,'ELIM',8X,'ALIM',8X,'RL  ',8X,'FNUL',8X,'DIG');
writef(lun,[repmat(' ',1,1),repmat('%12.4f',1,6), '\n '  ' \n'], tol , elim , alim , rl , fnul , dig);
%format(1x,6e12.4];
end;
ci =complex(0.0e0,1.0e0);
hpi = 2.0e0.*atan(1.0e0);
rhpi = 1.0e0./hpi;
pi = hpi + hpi;
if( mqc~=2 )
nl = 2;
il = 5;
for i = 1 : il;
keps(i) = 0;
kdo(i) = 0;
end; i = fix(il+1);
kdo(5) = 1;
nul = 5;
xnu(1) = 0.0e0;
xnu(2) = 1.0e0;
xnu(3) = 2.0e0;
xnu(4) = 0.5e0.*fnul;
xnu(5) = fnul + 1.2e0;
else;
nl = 4;
il = 13;
for i = 1 : il;
kdo(i) = 0;
keps(i) = 0;
end; i = fix(il+1);
kdo(2) = 1;
kdo(6) = 1;
kdo(8) = 1;
kdo(11) = 1;
kdo(12) = 1;
kdo(13) = 1;
keps(3) = 1;
keps(4) = 1;
keps(5) = 1;
keps(9) = 1;
nul = 6;
xnu(1) = 0.0e0;
xnu(2) = 0.6e0;
xnu(3) = 1.3e0;
xnu(4) = 2.0e0;
xnu(5) = 0.5e0.*fnul;
xnu(6) = fnul + 1.2e0;
end;
i = 2;
eps = 0.01e0;
film = il - 1;
t(1) = -pi + eps;
for k = 2 : il;
if( kdo(k)==0 )
t(i) = pi.*(-il+2.*k-1)./film;
if( keps(k)~=0 )
ts = t(i);
t(i) = ts - eps;
i = fix(i + 1);
t(i) = ts + eps;
else;
i = fix(i + 1);
end;
end;
end; k = fix(il+1);
itl = fix(i - 1);
if( kprint>=2 )
writef(lun,[' CHECKS IN THE (Z,FNU) SPACE' ' \n']);
%format (' CHECKS IN THE (Z,FNU) SPACE');
end;
lflg = 0;
for kode = 1 : 2;
for n = 1 : nl;
for nu = 1 : nul;
fnu = xnu(nu);
ifnu = fix(fix(fnu));
ffnu = fnu - ifnu;
arg = hpi.*ffnu;
csgn = complex(cos(arg),sin(arg));
i4 = fix(rem(ifnu,4) + 1);
csgn = csgn.*cip(i4);
cspn = conj(csgn).*complex(rhpi,0.0e0);
csgn = csgn.*ci;
for icase = 1 : 3;
irb = fix(min(2,icase));
for ir = irb : 4;
if( icase==2 )
r =(2.0e0.*(4-ir)+r2.*(ir-1))./3.0e0;
elseif( icase==3 ) ;
if( r2>=rm )
break;
end;
r =(r2.*(4-ir)+rm.*(ir-1))./3.0e0;
else;
r =(eps.*(4-ir)+2.0e0.*(ir-1))./3.0e0;
end;
for it = 1 : itl;
ct = cos(t(it));
st = sin(t(it));
if( abs(ct)<atol )
ct = 0.0e0;
end;
if( abs(st)<atol )
st = 0.0e0;
end;
z = complex(r.*ct,r.*st);
xx = real(z);
yy = imag(z);
[z,fnu,kode,n,w,nz2,ierr]=cbesi(z,fnu,kode,n,w,nz2,ierr);
if( nz2==0 )
[z,fnu,kode,n,y,nz1,ierr]=cbesk(z,fnu,kode,n,y,nz1,ierr);
if( nz1==0 )
zn = z.*ci;
xn = -yy;
yn = xx;
[zn,fnu,kode,n,v,nz,cwrk,ierr]=cbesy(zn,fnu,kode,n,v,nz,cwrk,ierr);
if( nz==0 )
coe1 = csgn;
coe2 = cspn;
if( kode==2 )
cc = -xx - abs(xx);
if( cc>(-alim) )
cw = complex(cc,-yy);
coe2 = coe2.*exp(cw);
else;
coe2 = complex(0.0e0,0.0e0);
continue;
end;
end;
for kk = 1 : n;
y(kk) = y(kk).*coe2;
w(kk) = w(kk).*coe1;
coe1 = coe1.*ci;
coe2 = -coe2.*ci;
end; kk = fix(n+1);
mflg = 0;
for i = 1 : n;
ab = fnu + i - 1;
aa = max(0.5e0,ab);
cw = w(i) - y(i);
av = abs(v(i));
er = abs(cw-v(i));
if( av~=0.0e0 )
if( yn~=0.0e0 )
er = er./av;
elseif( xn>0.0e0 ) ;
if( abs(xn)<aa )
er = er./av;
end;
elseif( abs(ffnu-0.5e0)<0.125e0 ) ;
if( abs(xn)<aa )
er = er./av;
end;
else;
er = er./av;
end;
end;
aer(i) = er;
if( er>ertol )
mflg = 1;
end;
end; i = fix(n+1);
if( mflg~=0 )
if( lflg==0 )
if( kprint>=2 )
writef(lun,[ '\n ' ,' CASES WHICH VIOLATE THE RELATIVE ','ERROR TEST WITH ERTOL = ','%12.4f', '\n '  ' \n'], ertol);
%format [' CASES WHICH VIOLATE THE RELATIVE ','ERROR TEST WITH ERTOL = ',e12.4];
writef(lun,[' INPUT TO CBESY   ZN, FNU, KODE, N' ' \n']);
%format (' INPUT TO CBESY   ZN, FNU, KODE, N');
end;
if( kprint>=3 )
writef(lun,[' COMPARE Y(ZN,FNU) WITH COE1*I(Z,FNU)','-COE2*K(Z,FNU)' ' \n']);
%format (' COMPARE Y(ZN,FNU) WITH COE1*I(Z,FNU)','-COE2*K(Z,FNU)');
writef(lun,[' Z = ZN*EXP(-i*PI/2)', '\n ' ,' COE1 = EXP(i*(FNU+1)*PI/2)  ',' COE2 = (2/PI)*EXP(-i*FNU*PI/2)' ' \n']);
%format (' Z = ZN*EXP(-i*PI/2)'/' COE1 = EXP(i*(FNU+1)*PI/2)  ',' COE2 = (2/PI)*EXP(-i*FNU*PI/2)');
writef(lun,[' RESULTS FROM CBESY   V(KK)', '\n ' ,repmat(' ',1,9),'FROM CBESI   W(KK)', '\n ' ,repmat(' ',1,9),'FROM CBESK   Y(KK)' ' \n']);
%format (' RESULTS FROM CBESY   V(KK)'/9X,'FROM CBESI   W(KK)'/9X,'FROM CBESK   Y(KK)');
writef(lun,[' TEST CASE INDICES   IR, IT, ICASE', '\n '  ' \n']);
%format (' TEST CASE INDICES   IR, IT, ICASE'];
end;
end;
lflg = fix(lflg + 1);
if( kprint>=2 )
writef(lun,['   INPUT:   ZN=',repmat('%12.4f',1,2),repmat(' ',1,3),'FNU=','%12.4f',repmat(' ',1,3),'KODE=','%3i',repmat(' ',1,3),'N=','%3i' ' \n'], zn , fnu , kode , n);
%format ('   INPUT:   ZN=',2E12.4,3X,'FNU=',e12.4,3X,'KODE=',i3,3X,'N=',i3);
end;
if( kprint>=3 )
for k=(1):(n), writef(lun,['   ERROR:   AER(K)=',repmat('%12.4f',1,4) ' \n'],aer(k)); end;
%format ('   ERROR:   AER(K)=',4E12.4);
writef(lun,[repmat(' ',1,12),'Z=',repmat('%12.4f',1,2), '\n ' ,repmat(' ',1,12),'COE1=',repmat('%12.4f',1,2),repmat(' ',1,3),'COE2=',repmat('%12.4f',1,2) ' \n'], z , coe1 , coe2);
%format (12X,'Z=',2E12.4/12X,'COE1=',2E12.4,3X,'COE2=',2E12.4);
kk = fix(max(nz1,nz2) + 1);
kk = fix(min(n,kk));
writef(lun,[' RESULTS:   V(KK)=',repmat('%12.4f',1,2), '\n ' ,repmat(' ',1,12),'W(KK)=',repmat('%12.4f',1,2), '\n ' ,repmat(' ',1,12),'Y(KK)=',repmat('%12.4f',1,2) ' \n'], v(kk) , w(kk) , y(kk));
%format (' RESULTS:   V(KK)=',2E12.4/12X,'W(KK)=',2E12.4/12X,'Y(KK)=',2E12.4);
writef(lun,['    CASE:   IR=','%3i',repmat(' ',1,3),'IT=','%3i',repmat(' ',1,3),'ICASE=','%3i', '\n '  ' \n'], ir , it , icase);
%format ('    CASE:   IR=',i3,3X,'IT=',i3,3X,'ICASE=',i3];
end;
end;
end;
end;
end;
end; it = fix(itl+1);
end;
end;
end;
end;
end;
if( kprint>=2 )
if( lflg==0 )
writef(lun,[' QUICK CHECKS OK' ' \n']);
%format (' QUICK CHECKS OK');
else;
writef(lun,[' ***','%5i',' FAILURE(S) FOR CBESY IN THE (Z,FNU) ','PLANE' ' \n'], lflg);
%format (' ***',i5,' FAILURE(S) FOR CBESY IN THE (Z,FNU) ','PLANE');
end;
end;
ipass = 0;
if( lflg==0 )
ipass = 1;
end;
if( ipass==1 && kprint>=2 )
writef(lun,[ '\n ' ,' * CBESY  PASSED ALL TESTS  *', '\n '  ' \n']);
%format [' * CBESY  PASSED ALL TESTS  *'];
end;
if( ipass==0 && kprint>=1 )
writef(lun,[ '\n ' ,' * CBESY  FAILED SOME TESTS *', '\n '  ' \n']);
%format [' * CBESY  FAILED SOME TESTS *'];
end;
end %subroutine cqcby