% Converted to a routine more appropriate for use with Matlab, Octave, and Lyme by
% Robert J. McGough
% Michigan State University
% 18 Dec 2005
% It is unlikely that the vectorized operations are optimal, but they
% appear to work.
% Updated: 14 Dec 2008
% Handles negative arguments correctly (Matlab version doesn't).

%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)
%

function [bj0]=besselj0(x);
ixeq0 = find(x == 0);
bj0(ixeq0)=1.0e0;
ixle4 = find(abs(x) > 0 & abs(x) <= 4);
% elseif(x <= 4.0e0);
t=x(ixle4)./4.0e0;
t2=t.*t;
bj0(ixle4)=((((((-.5014415e-3.*t2+.76771853e-2).*t2-.0709253492e0).*t2+.4443584263e0).*t2 -1.7777560599e0).*t2+3.9999973021e0).*t2-3.9999998721e0).*t2+1.0e0;

ixgt4 = find(abs(x) > 4);
% else;
t=4.0e0./abs(x(ixgt4));
t2=t.*t;
a0=sqrt(2.0e0./(pi.*abs(x(ixgt4))));
p0=((((-.9285e-5.*t2+.43506e-4).*t2-.122226e-3).*t2+.434725e-3).*t2-.4394275e-2).*t2+.999999997e0;
q0=t.*(((((.8099e-5.*t2-.35614e-4).*t2+.85844e-4).*t2-.218024e-3).*t2+.1144106e-2).*t2-.031249995e0);
ta0=abs(x(ixgt4))-.25e0.*pi;
bj0(ixgt4)=a0.*(p0.*cos(ta0)-q0.*sin(ta0));