besjnu.m

很多matlab的源代码

function Jn = besjnu(n1,x1)
% BESJNU Spherical bessel function j(x)
%
%	Y = BESJNU(N,X) returns the spherical Bessel function j(x) 
%	N = order of polynomial. N must be a SCALAR
%
%	NOTE: For the definition, see Abramowitz and Stegun
%
%	BESJNU (with no input arguments) invokes the following example:
%
%	% Plot j(x) of order 3 over (0 12) at intervals of 0.3
%	  >>x = 0:0.3:12;
%	  >>y = besjnu(3,x);
%	  >>plot(x,y),grid


% ADSP Toolbox: Version 2.0 
% For use with "Analog and Digital Signal Processing", 2nd Ed.
% Published by PWS Publishing Co.
%
% Ashok Ambardar, EE Dept. MTU, Houghton, MI 49931, USA
% http://www.ee.mtu/faculty/akambard.html
% e-mail: [e-mail]
% Copyright (c) 1998



if nargin==0,help besjnu,disp('Strike a key to see results of the example')
pause,x=0:0.3:12;y=besjnu(3,x);plot(x,y),axesn,grid,return,end

x=(x1==0);z=x1+x;n=n1+0.5;

if exist('besselj')==2, 
Jn = sqrt(pi/2)*besselj(n,z)./sqrt(z);
else
tol=2e-16;
flag=0;
if (rem(n,1)==0 & n<0),n=abs(n);flag=1;end
Jn=0*z;s=prod(size(z));
i1=find((z<50 | z<n+20) & n>=0);
s1=prod(size(i1));
if s1>0
z2=z(i1);a=exp(n*lnn(z2/2))/gm(n+1);
C=sqrt((pi/2)./z2);y=z2.*z2/4;z1=y.*y/4;
J1=ghgf(0,3,[1 1 1 1],[(1+n)/2 1+n/2 n+1 1],-z1,tol);
In=ghgf(0,1,[1 1 1 1],[n+1 1 1 1],y,tol);
Jn(i1)=C.*J1.*a./In;
end
if s>s1
i2=(1:s);i2(i1)=[];z=z(i2);m=.5*n;
y=z-(m+.25)*pi;z1=z.*z;
A=ghgf(4,1,[.25+m .25-m .75+m .75-m],[.5 1 1 1],(-1)./z1,tol);
B=.5*(.25-n*n)*ghgf(4,1,[.75+m .75-m 5/4+m 5/4-m],[3/2 1 1 1],(-1)./z1,tol)./z;
Jn(i2)=(A.*cos(y)+B.*sin(y))./z;
end
if flag==1,Jn=((-1)^n)*Jn;end

end
Jn=(1-x).*Jn+x*(n1==0);
Jn=real(Jn);