Bessel functions

[Eyal Doron]

Hi,
   Ok, here they are. The following .m file calculates the Bessel and
Newmann functions for real orders and complex arguments. It returns
all orders from {\nu} to \nu in increments of one, where \nu is the
maximal requested order. Also works for negative \nu, and can return
orders from -\nu to \nu .

This .m file uses various variations on the subject of forward-backward
recurrences + normalization. It seems to work to about 12 digits at least,
usually to 14-15 digits, and is fully vectorized.

I have taken care of various sections in the complex plane, but I'm not
sure that all of them are OK, so use at your own risk. Also, it should
be possible to extend it to complex orders, but I'm too lazy.

Enjoy!

Eyal Doron

-----------------------jybess.m------------------------
function [Jn,Yn] = jybess(n,x,ToDo);
% FUNCTION [Jn,Yn] = jybess(n,x,ToDo):
%
% returns the Bessel and Newmann functions for complex arguments and real
% order. Returns a size [length(x), |n|+1] or [length(x), 2|n|+1] matrix.
%
% ToDo - Optional command string. 'Y' returns only Y_n(x) instead of
%        J_n(x), 'L' returns only the highest computed order.
% n    - Maximal order, real scalar.
% x    - Real or complex argument vector.
% ToDo - Optional string of command characters:
%        'S' - returns the orders -n to n rather than 0 to n. (also ToDo=1).
%        'Y' - if nargout<2, return Y_n(x) instead of J_n(x).
%        'L' - return only the highest computed order.
% -----------------------------------------
% Jn   - Bessel functions of the first kind.
% Yn   - Bessel functions of the second kind (Newmann functions), optional.

% Algorithm:
% ----------
% 1) Calculate J_n by backwards recursion from a calculated starting point.
%    Minimal order is \nu=|n|-floor(|n|).
% 2) Normalization: Normalize J_\nu(x) by:
%    a) for integer n and |\Im(x)|<log(10), use
%                         1 + 2 sum_{k=1}^\infty J_{2k}(x) = 1 
%    b) for integer n and |\Im(x)|>log(10), use
%                         1 + 2 sum_{k=1}^\infty (-)^k J_{2k}(x) = cos(x)
%    c) for half-integer n, use explicit formula for the spherical Bessels.
%    d) for real n and |\Im(x)|<5, use
%       sum_{k=0}^\infty (\nu+2k)\Gamma(\nu+k)/k! J_{\nu+2k}(x) = (x/2)^\nu
%    e) for real n and |\Im(x)|>5, use
%       e^{i\phi\nu} \sum_{k=0}^\infty [r/2*(1-e^{2i\phi})]^k/k! J_{\nu+k}(r)
%                   = J_\nu(r e^i\phi)
%       This requires a recursive call to jybess to obtain the J_{\nu+k}(r).
% 3) If Y_n, or J_n for negative non-integer n, are required, calculate
%    the Y_n(x) as follows:
% 4) Evaluate Y_\nu:
%    a) for nu=0, use
%      2/\pi(\log(x/2)+\gamma)J_0(x) - 4/\pi\sum_{k=1}^\infty (-)^k J_{2k}(x)/k
%                  = Y_0(x)
%    b) for nu=0.5, use explicit formula for the spherical Bessels.
%    c) otherwise, do the following:
%        i) Calculate J_{1-\nu}(x) and J_{2-\nu}(x) using a recursive call.
%       ii) Calculate J_{-\nu}(x) using one backwards recursion step.
%      iii) Use [J_\nu(x)\cos(\nu\pi)-J_{-\nu}(x)]/\sin(\nu\pi) = Y_\nu(x)
% 5) Obtain the rest of the Y_{\nu+k}(x) using the Wronskian relation
%    J_{\nu+1}(x)Y_\nu(x) - J_\nu(x)Y_{\nu+1}(x) = 2/(\pi x)
%    This turns out to be much more stable for complex arguments than the
%    more conventional forward recurrence technique.
% 6) If J for negative orders is required:
%    a) for nu=0, use J_n(x) = (-)^n J_{-n}(x)
%    b) otherwise, use 4.c.iii
% 7) If Y for negative orders is required:
%    a) for nu=0, use Y_n(x) = (-)^n Y_{-n}(x)
%    b) otherwise, use 4.c.iii

if nargin==0
  help jybess
  return
end
if nargin<2
  error('Not enough input parameters!');
end
if max(max(size(n)))>1 | any(imag(n)~=0)
  error('Max order must be a real scalar!')
end
sym=0; Return_Y=(nargout==2); Return_Last=0;
if nargin>2
  sym=any(ToDo==1) | any(ToDo=='s') | any(ToDo=='S');
  Return_Y=Return_Y | any(ToDo=='y' | ToDo=='Y');
  Return_Last=any(ToDo=='l' | ToDo=='L');
end
orig_n=n;
n_is_negative=(n<0); n=abs(n);
nu=n-floor(n); n=floor(n);
if nargin<3, sym=0; end
sym=(sym~=0);
x=x(:).'; xlen=length(x);
calc_Y=(Return_Y | sym | n_is_negative);
Acc_Jnorm=(nu~=0.5);
Acc_Ynorm=calc_Y & nu==0;

z = x==0; x = x + z;            % Temporarily replace x=0 with x=1

tiny = 16^(-250);

%                             Starting index for backwards recurrence.
c = [ 0.9507    1.4208   14.1850
      0.9507    1.4208   14.1850
      0.9507    1.4208   14.1850
      0.7629    1.4222   13.9554
      0.7369    1.4289   13.1756
      0.7674    1.4311   12.4523
      0.8216    1.4354   11.2121
      0.8624    1.4397    9.9718
      0.8798    1.4410    8.9217
      0.9129    1.4360    7.8214
      0.9438    1.5387    6.5014
      0.9609    1.5216    5.7256
      0.9693    1.5377    5.3565
      0.9823    1.5220    4.5659
      0.9934    1.5049    3.7902
      0.9985    1.4831    3.2100
      1.0006    1.4474    3.0239
      0.9989    1.4137    2.8604
      0.9959    1.3777    2.7760
      1.0005    1.3500    2.3099]';
j = 1+min(n,19);
m = c(1,j).*max(3,j) + c(2,j).*(max(1,abs(x))-1) + ...
    c(3,j)./(1-log(min(1,abs(x))));
m=max([m; n+10+0*m]);
m = 4*ceil(m/4);

%                                       Prevent underflow
logtiny=log(tiny); logx=log(abs(x)/2)+1;
II=1:xlen;
while ~isempty(II)
  est=-log(2*pi*m(II))/2+m(II).*(logx(II)-log(m(II)))<logtiny;
  II=II(find(est)); m(II)=m(II)-4;
end;
mm = max(m(:));

% Normalization summation coefficients

if nu==0
  Nrm=[1; 2*ones(ceil(mm/2),1)]; 
  UseCos=find(abs(imag(x))>log(10));
  if ~isempty(UseCos)                  % Use norm. to cos(z) instead of 1.
    Nrm=Nrm(:,ones(xlen,1));
    [r,c]=size(Nrm);
    Nrm(2:2:r,UseCos)=-Nrm(2:2:r,UseCos);
  end
elseif nu~=0.5               % nu=0.5 doesn't use normalizations
  k=(1:ceil(mm/2)).';
  Nrm=cumprod([nu*gamma(nu); (nu+2*k).*(nu+k-1)./k./(nu+2*k-2)]);
  Nrm=Nrm(:,ones(xlen,1));
end

%                           Backwards recursion for the Jn
Jn=zeros(n+1,xlen);
k = mm;
bkp1 = 0*x;
bk = tiny*(m==k);
if Acc_Jnorm, t = Nrm(k/2+1,:).*bk; end
if Acc_Ynorm, y = 2*bk/k; end
sig = 1;
for k = mm-1:-1:0
  bkp2 = bkp1;
  bkp1 = bk;
  bk = 2*(k+1+nu)*bkp1./x - bkp2 + tiny*(m==k);
  if k<=n
    Jn(k+1,:)=bk;
  end
  if (floor(k/2)*2-k==0)
     if Acc_Jnorm, t = t + Nrm(k/2+1,:).*bk; end
     if Acc_Ynorm & k>0
       sig = -sig;
       y = y + sig*2*bk/k;
     end
  end
end
clear Nrm bkp2 sig

if nu==0.5               % Use explicit formulae for the spherical Bessels
  si=sin(x);
  II=(abs(si)>0.1);
  F1=find(II); F2=find(1-II);
  nrm=ones(1,xlen);
  SpFact=sqrt(2*x/pi);
  if ~isempty(F1)
    nrm(F1)=SpFact(F1).*(si(F1)./x(F1))./Jn(1,F1);
  end
  if ~isempty(F2)
    nrm(F2)=SpFact(F2).*(si(F2)./x(F2).^2-cos(x(F2))./x(F2))./bkp1(F2);
  end
elseif nu==0
  nrm=ones(size(x));
  if ~isempty(UseCos)
    nrm(UseCos)=cos(x(UseCos));
  end
  nrm=nrm./t;
else
  II=(abs(imag(x))<5);
  F1=find(II); F2=find(1-II);
  nrm=ones(1,xlen);
  if ~isempty(F1)
    nrm(F1)=(x(F1)/2).^nu ./t(F1);
  end
  if ~isempty(F2)               % Use mult. Theorem to calc J_0
    r=abs(x(F2)); ph=angle(x(F2)); rn=max(ceil(r));
    MaxOrder=ceil(18+1.25*rn);
    Jtmp=jybess(nu+MaxOrder,r);
    k=(0:MaxOrder).';
    Kg=cumsum([0;log(1:MaxOrder).']);          % log(k!)
    v=i*ph-i*pi/2+log(imag(x(F2)));
    
nrm(F2)=exp(i*nu*ph).*sum(exp(k*v-Kg(:,ones(length(F2),1))).*Jtmp.')./Jn(1,F2);
  end
end
Jn=Jn.*nrm(ones(n+1,1),:);           % Normalizing condition.

%                Restore results for x = 0; j0(0) = 1, jn(0) = 0 for n > 0.
if any(z)
  Ii=find(z); LI=max(size(Ii));
  Jn(:,Ii)=[ones(1,LI);zeros(n,LI)];
end
Jn=Jn.';

if calc_Y                                    % Also Yn
%                         First, get Y_0(x)
  if n==0
    J1 = bkp1.*nrm;
  else
    J1=Jn(:,2).';
  end
  if nu==0                          % Use dependence on sums of J
    gamma = 0.57721566490153286;
    lx=log(x/2);
    II=find(imag(x)==0 & real(x)<0);
    if ~isempty(II)
      lx(II)=log(abs(x(II))/2);
    end
    y = 2/pi*(lx + gamma).*bk - 4/pi*y;
    Y0 = y.*nrm;
  elseif nu==0.5                    % explicit formula
    Y0=-cos(x)./x.*SpFact;
  else                              % Use linear relations between J and Y
    Jtmp=jybess(2-nu,x);
    Jminus=2*(1-nu)./x.'.*Jtmp(:,1)-Jtmp(:,2);
    Y0=(Jn(:,1).'*cos(pi*nu)-Jminus.')/sin(pi*nu);
  end
  
  Yn=[Y0.', zeros(xlen,n)];
  for l=1:n                      % Iterate using Wronskian relation
    Yn(:,l+1)=(Jn(:,l+1).*Yn(:,l)-2/pi./x)./Jn(:,l);
  end
  
  if any(z), Yn(find(z),:)=-inf*ones(size(Yn(find(z),:))); end
end

if nu==0
  if n_is_negative                                 % negative order
    Jn(:,2:2:n+1)=-Jn(:,2:2:n+1);
    if Return_Y
      Yn(:,2:2:n+1)=-Yn(:,2:2:n+1);
    end
  end
  if sym & n>0                                     % symmetric output
    Rev=n:-1:1; Rev=2*(floor(Rev/2)*2==Rev)-1;
    Jn=[Jn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Jn];
    if Return_Y
      Yn=[Yn(:,n+1:-1:2).*Rev(ones(xlen,1),:) Yn];
    end
  end
else
  if n_is_negative | sym
    ord=nu+(0:n); co=cos(pi*ord); so=sin(pi*ord);
    Jminus=Jn.*co(ones(xlen,1),:)-Yn.*so(ones(xlen,1),:);
    if Return_Y
      ord=nu+(0:n);co=cos(pi*ord); so=sin(-pi*ord);
      Yminus=(Jminus.*co(ones(xlen,1),:)-Jn)./so(ones(xlen,1),:);
    end
  end
  if sym
    Jn=[fliplr(Jminus) Jn];
    if n_is_negative, Jn=fliplr(Jn); end
    if Return_Y
      Yn=[fliplr(Yminus) Yn];
      if n_is_negative, Yn=fliplr(Yn); end
    end
  elseif n_is_negative
    Jn=Jminus;
    if Return_Y, Yn=Yminus; end
  end
end
if Return_Last
  [r,c]=size(Jn);
  Jn=Jn(:,c);
  if calc_Y, Yn=Yn(:,c); end
end
if nargout<2 & Return_Y, Jn=Yn; end