besjn.m

很多matlab的源代码

function Jn = besjn(n1,x1)
% BESJN Bessel function J(x) for integer and non-integer.
%
%	Y = BESJN(N,X) returns the Bessel function J(X)
%	N = order (integer or non-integer). N must be a scalar. 
%
%	NOTE: For a definition, see Abramowitz and Stegun
%
%       BESJN (with no input arguments) invokes the following example:
%
%       % Plot J(x) of order 3 over (0 12) at intervals 0f 0.3
%         >>x = 0:0.3:12;
%         >>y = besjn(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 besjn,disp('Strike a key to see results of the example')
pause,x=0:0.3:12;y=besjn(3,x);plot(x,y),axesn,grid,return,end

if exist('besselj')==2,Jn=real(besselj(n1,x1));return,end

tol=2e-16;
x=(x1==0);z=x1+x;n=n1;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);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)=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;C=sqrt((2/pi)./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)=C.*(A.*cos(y)+B.*sin(y));
end
if flag==1,Jn=((-1)^n1)*Jn;end
Jn=(1-x).*Jn+x*(n1==0);
Jn=real(Jn);