help-octave
[Top][All Lists]

## Special functions

 From: Eyal Doron Subject: Special functions Date: Fri, 11 Aug 1995 11:23:57 +0100 (MET DST)

```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

legendre.m and spharm.m - Associated Legendre functions and spherical
harmonics, can be picked up from ftp.mathworks.com/pub/contrib/math
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
```