Special functions

[Eyal Doron]

Hi,
   All right, here they are:

jybess.m - Bessel functions of the first and second kind, integer positive
           order and real arguments. Its highly vectorized both in order
           and argument. However, it fails if the order is much, much larger
           than the argument, so that more than 300 orders of magnitude are
           spanned. It can be fixed, but I haven't the energy. Volunteers?
jybess_d.m - Derivatives of Bessel functions.
spbess.m - Spherical Bessel functions of the first kind. The second kind
           can be easily added.

legendre.m and spharm.m - Associated Legendre functions and spherical 
harmonics, can be picked up from ftp.mathworks.com/pub/contrib/math
(I think they are freely downloadable, but copyrighted, so I won't post
them here).

Have fun!

Eyal Doron

-----------------------jybess.m------------------------------
function [Jn,Yn] = jybess(n,x)
%jybess Integer order Hankel functions.
%       JYBESS(n,X) = Bessel functions, J-sub-n (X) and Y-sub-n(X).
%
%       The function is evaluated for every point in the array x. The results
%       are matrices of (length(x)) rows by (n+1) columns.
%       The maximal order, n, must be a nonnegative integer scalar.

%       Reference: Abramowitz and Stegun, Handbook of Mathematical
%         Functions, National Bureau of Standards, Applied Math.
%         Series #55, sections 9.1.1, 9.1.89 and 9.12.
   if nargout==2 & sum(x==0)
     error('x=0 not permitted when caclulating Y!');
   end
   if n<0
     error('Negative argument!')
   end
   [r c]=size(x); if r>c, x=x'; end
   % Backwards three term recurrence for J-sub-n(x).
   % Temporarily replace x=0 with x=1
   z = x==0;
   x = x + z;

   % Starting index for backwards recurrence.
   c = [17.1032 0.2639 0.6487 -0.0018 0.6457];
   m = round(c(1) + c(2)*n + c(3)*abs(x) + c(4)*n*abs(x) + c(5)*max(n,abs(x)));
   % Make sure rem(m,4)=0.
   m = 4*ceil(m/4);
   mm = max(m(:));

tiny = 16^(-250);    %   if computer=='VAX', tiny = 16^(-30);
                     %      else tiny = 16^(-250); end
   k = mm;
   bkp1 = 0*x;
   bk = tiny*(m==k);
   t = 2*bk;
   y = 2*bk/k;
   sig = 1;
   for k = mm-1:-1:0
     bkp2 = bkp1;
     bkp1 = bk;
     bk = 2*(k+1)*bkp1./x - bkp2 + tiny*(m==k);
     if k<=n, Jn(k+1,:)=bk;end
     if (k > 0 & floor(k/2)*2-k==0)
        t = t + 2*bk;
        sig = -sig;
        y = y + sig*2*bk/k;
     end
     if k == 0 , t = t + bk; end
   end
   if nargout==2
     gamma = 0.57721566490153286;
     y = 2/pi*(log(x/2) + gamma).*bk - 4/pi*y;
   end
   % Normalizing condition, J0(x) + 2*J2(x) + 2*J4(x) + ... = 1.
%   Jn = bn./t;
   nrm=ones(n+1,1)*(1 ./t);
   Jn=Jn.*nrm;
   % Restore results for x = 0; j0(0) = 1, jn(0) = 0 for n > 0.
   if any(z)==1
     Ii=find(z); LI=max(size(Ii));
     Jn(:,Ii)=[ones(1,LI);zeros(n,LI)];
   end
   J0 = bk./t;
   J1 = bkp1./t;

   if nargout==1, return; end

   Y0 = y./t;

   % Forward recurrence for Yn.
   if n>0
      ykm1 = Y0;
      yk = (J1.*Y0 - (2)./(pi*x))./J0;
      Yn=[Y0;yk];
      for k = 2:n
        ykm2 = ykm1;
        ykm1 = yk;
        yk = 2*(k-1)*ykm1./x - ykm2;
        Yn=[Yn;yk];
      end
%      Yn = yk;
   else
      Yn = Y0;
   end

Jn=Jn.';
if nargout==2, Yn=Yn.'; end

-----------------------------------jybess_d.m-------------------------
function [Jnd,Ynd] = jybess_d(n,x)
%jybess Integer order Hankel functions.
%       [Jnd,Ynd]=JYBESS_d(n,X) = Derivatives of Bessel functions,
%                                 J-sub-n (X) and Y-sub-n(X).
%
%       The function is evaluated for every point in the array X.
%       The maximal order, n, must be a nonnegative integer scalar.

%       Reference: Abramowitz and Stegun, Handbook of Mathematical
%         Functions, National Bureau of Standards, Applied Math.
%         Series #55, sections 9.1.1, 9.1.89 and 9.12.

if nargout==1
  J=jybess(n+1,x);
else
  [J,Y]=jybess(n+1,x);
end
Jnd=[-J(:,2) (J(:,1:n)-J(:,3:n+2))/2];
if nargout==2
  Ynd=[-Y(:,2) (Y(:,1:n)-Y(:,3:n+2))/2];
end

-------------------------sphbess.m------------------------------
function Jn = sphbess(n,x)
%SPHBESS  Spherical Bessel functions of the first kind.
%         SPHBESS(n,X) = Spherical Bessel function, j-sub-n (X).
%         SPHBESS returns a matrix of (n+1) orders by (length(X)) arguments.
%         If X is a scalar, it returns a row vector.

   [r c]=size(x); if r>c, x=x'; end
   % Backwards three term recurrence for j-sub-n(x).
   % Treat x==0:
%   if x==0, Jn=[1 zeros(1,n)]; return; end
   % Temporarily replace x=0 with x=1
   z = x==0;
   x = x + z;

   % Starting index for backwards recurrence.
   c = [17.1032 0.2639 0.6487 -0.0018 0.6457];
   m = round(c(1) + c(2)*n + c(3)*abs(x) + c(4)*n*abs(x) + c(5)*max(n,abs(x)));
   % Make sure m is even.
   m = 2*ceil(m/2);
   mm = max(m(:));

   tiny = 16^(-100);
   k = mm;
   bn=zeros(n+1,length(x));
   bkp1 = 0*x;
   bk = tiny*(m==k);
   t = (2*k+1)*bk.^2;
   for k = mm-1:-1:0
     bkp2 = bkp1;
     bkp1 = bk;
     bk = (2*(k+1)+1)*bkp1./x - bkp2 + tiny*(m==k);
     if k <= n , bn(k+1,:) = bk; end
     t = t + (2*k+1)*bk.^2;
   end
   % Normalizing condition: sum[(2n+1)*jn(x)^2] = 1.
   Jn = bn.*(ones(n+1,1)*(1 ./sqrt(t)));

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