function mjy01b
%This program is a direct conversion of the corresponding Fortran program in
%S. Zhang & J. Jin "Computation of Special Functions" (Wiley, 1996).
%online: http://iris-lee3.ece.uiuc.edu/~jjin/routines/routines.html
%
%Converted by f2matlab open source project:
%online: https://sourceforge.net/projects/f2matlab/
% written by Ben Barrowes (barrowes@alum.mit.edu)
%

%       =========================================================
%       Purpose: This program computes Bessel functions Jn(x)
%                and Yn(x,n=0,1)and their derivatives
%                using subroutine JY01B
%       Input :  x   --- Argument of Jn(x)& Yn(x,x ò 0)
%       Output:  BJ0 --- J0(x)
%                DJ0 --- J0'(x)
%                BJ1 --- J1(x)
%                DJ1 --- J1'(x)
%                BY0 --- Y0(x)
%                DY0 --- Y0'(x)
%                BY1 --- Y1(x)
%                DY1 --- Y1'(x)
%       Example:
%        x       J0(x)J0'(x)J1(x)J1'(x)
%       ---------------------------------------------------------
%        1     .76519769   -.44005059    .44005059    .32514710
%        5    -.17759677    .32757914   -.32757914   -.11208094
%       10    -.24593576   -.04347275    .04347275   -.25028304
%       20     .16702466   -.06683312    .06683312    .16368301
%       30    -.08636798    .11875106   -.11875106   -.08240961
%       40     .00736689   -.12603832    .12603832    .00421593
%       50     .05581233    .09751183   -.09751183    .05776256
%        x       Y0(x)Y0'(x)Y1(x)Y1'(x)
%      ---------------------------------------------------------
%        1     .08825696    .78121282   -.78121282    .86946979
%        5    -.30851763   -.14786314    .14786314   -.33809025
%       10     .05567117   -.24901542    .24901542    .03076962
%       20     .06264060    .16551161   -.16551161    .07091618
%       30    -.11729573   -.08442557    .08442557   -.12010992
%       40     .12593642    .00579351   -.00579351    .12608125
%       50    -.09806500    .05679567   -.05679567   -.09692908
%       =========================================================
x=[];bj0=[];dj0=[];bj1=[];dj1=[];by0=[];dy0=[];by1=[];dy1=[];
fprintf(1,'%s \n','please enter x');
%        READ(*,*)X
x=50;
fprintf(1,[repmat(' ',1,3),'x =','%5.1g' ' \n'],x);
fprintf(1,'%s ','  x          j0(x)j0''(x)j1(x)');fprintf(1,'%s \n','          j1''(x)');
fprintf(1,'%s ','------------------------------------------');fprintf(1,'%s \n','----------------------------');
[x,bj0,dj0,bj1,dj1,by0,dy0,by1,dy1]=jy01b(x,bj0,dj0,bj1,dj1,by0,dy0,by1,dy1);
fprintf(1,[repmat(' ',1,1),'%5.1g',repmat('%16.8g',1,4) ' \n'],x,bj0,dj0,bj1,dj1);
fprintf(1,'%0.15g \n');
fprintf(1,'%s ','  x          y0(x)y0''(x)y1(x)');fprintf(1,'%s \n','          y1''(x)');
fprintf(1,'%s ','------------------------------------------');fprintf(1,'%s \n','----------------------------');
fprintf(1,[repmat(' ',1,1),'%5.1g',repmat('%16.8g',1,4) ' \n'],x,by0,dy0,by1,dy1);
%format(1x,f5.1,4e16.8);
%format(3x,',f5.1);
end
function [x,bj0,dj0,bj1,dj1,by0,dy0,by1,dy1]=jy01b(x,bj0,dj0,bj1,dj1,by0,dy0,by1,dy1,varargin);
%       =======================================================
%       Purpose: Compute Bessel functions J0(x), J1(x), Y0(x),
%                Y1(x), and their derivatives
%       Input :  x   --- Argument of Jn(x)& Yn(x,x ò 0)
%       Output:  BJ0 --- J0(x)
%                DJ0 --- J0'(x)
%                BJ1 --- J1(x)
%                DJ1 --- J1'(x)
%                BY0 --- Y0(x)
%                DY0 --- Y0'(x)
%                BY1 --- Y1(x)
%                DY1 --- Y1'(x)
%       =======================================================
pi=3.141592653589793d0;
if(x == 0.0d0);
bj0=1.0d0;
bj1=0.0d0;
dj0=0.0d0;
dj1=0.5d0;
by0=-1.0d+300;
by1=-1.0d+300;
dy0=1.0d+300;
dy1=1.0d+300;
return;
elseif(x <= 4.0d0);
t=x./4.0d0;
t2=t.*t;
bj0=((((((-.5014415d-3.*t2+.76771853d-2).*t2-.0709253492d0).*t2+.4443584263d0).*t2 -1.7777560599d0).*t2+3.9999973021d0).*t2-3.9999998721d0).*t2+1.0d0;
bj1=t.*(((((((-.1289769d-3.*t2+.22069155d-2).*t2-.0236616773d0).*t2+.1777582922d0).*t2-.8888839649d0).*t2+2.6666660544d0).*t2 -3.9999999710d0).*t2+1.9999999998d0);
by0=(((((((-.567433d-4.*t2+.859977d-3).*t2-.94855882d-2).*t2+.0772975809d0).*t2 -.4261737419d0).*t2+1.4216421221d0).*t2-2.3498519931d0).*t2+1.0766115157).*t2 +.3674669052d0;
by0=2.0d0./pi.*log(x./2.0d0).*bj0+by0;
by1=((((((((.6535773d-3.*t2-.0108175626d0).*t2+.107657606d0).*t2-.7268945577d0).*t2+3.1261399273d0).*t2-7.3980241381d0).*t2+6.8529236342d0).*t2+.3932562018d0).*t2 -.6366197726d0)./x;
by1=2.0d0./pi.*log(x./2.0d0).*bj1+by1;
else;
t=4.0d0./x;
t2=t.*t;
a0=sqrt(2.0d0./(pi.*x));
p0=((((-.9285d-5.*t2+.43506d-4).*t2-.122226d-3).*t2+.434725d-3).*t2-.4394275d-2).*t2+.999999997d0;
q0=t.*(((((.8099d-5.*t2-.35614d-4).*t2+.85844d-4).*t2-.218024d-3).*t2+.1144106d-2).*t2-.031249995d0);
ta0=x-.25d0.*pi;
bj0=a0.*(p0.*cos(ta0)-q0.*sin(ta0));
by0=a0.*(p0.*sin(ta0)+q0.*cos(ta0));
p1=((((.10632d-4.*t2-.50363d-4).*t2+.145575d-3).*t2-.559487d-3).*t2+.7323931d-2).*t2+1.000000004d0;
q1=t.*(((((-.9173d-5.*t2+.40658d-4).*t2-.99941d-4).*t2+.266891d-3).*t2-.1601836d-2).*t2+.093749994d0);
ta1=x-.75d0.*pi;
bj1=a0.*(p1.*cos(ta1)-q1.*sin(ta1));
by1=a0.*(p1.*sin(ta1)+q1.*cos(ta1));
end;
dj0=-bj1;
dj1=bj0-bj1./x;
dy0=-by1;
dy1=by0-by1./x;
return;
end