% 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 [bj1]=besselj1(x);
ixeq0 = find(x == 0);
bj1(ixeq0)=0.0e0;
ixle4 = find(abs(x) > 0 & abs(x) <= 4);
% elseif(x <= 4.0e0);
t=abs(x(ixle4))./4.0e0;
t2=t.*t;
bj1(ixle4)=t.*(((((((-.1289769e-3.*t2+.22069155e-2).*t2-.0236616773e0).*t2+.1777582922e0).*t2-.8888839649e0).*t2+2.6666660544e0).*t2 -3.9999999710e0).*t2+1.9999999998e0);
ixgt4 = find(abs(x) > 4);
% else;
t=4.0e0./abs(x(ixgt4));
t2=t.*t;
a0=sqrt(2.0e0./(pi.*abs(x(ixgt4))));
p1=((((.10632e-4.*t2-.50363e-4).*t2+.145575e-3).*t2-.559487e-3).*t2+.7323931e-2).*t2+1.000000004e0;
q1=t.*(((((-.9173e-5.*t2+.40658e-4).*t2-.99941e-4).*t2+.266891e-3).*t2-.1601836e-2).*t2+.093749994e0);
ta1=abs(x(ixgt4))-.75e0.*pi;
bj1(ixgt4)=a0.*(p1.*cos(ta1)-q1.*sin(ta1));

% J1 is  odd 
ixneg = find(x < 0);
bj1(ixneg) = -bj1(ixneg);