Fonctions de Bessel sphériques pour matlab

Il est utile de redéfinir les fonctions de Bessel sphériques en matlab pour les calculs de SER.
Ces fonctions sont utilisées pour calculer les coefficients de Mie par exemple.

function jn=sphjn(n,x)
jn=sqrt(pi/2./x).*besselj(n+0.5,x);

function yn=sphyn(n,x)
yn=sqrt(pi/2./x).*bessely(n+0.5,x);

function h1n=sphh1n(n,x)
h1n=sqrt(0.5*pi./x).*besselh(n+0.5,x);

function h2n=sphh2n(n,x);
h2n=sqrt(0.5*pi./x).*besselh(n+0.5,2,x);

SER bistatique d’une sphère métallique, version matlab ou presque

Tout d’abord, le résultat dans le plan E :

Puis, dans le plan H :

clear all;
close all;

a=0.08;
ordre=20;
c=3e8;
freq=[0.001e9:0.1e9:9e9]';
k0=2/c*pi.*freq;
k0a=k0.*a;

nbptst=80;
theta_min=0.0001;
theta_max=pi*0.999;
nbptsp=90;
phi_min=0.0001;
phi_max=pi*0.999;

stept=(theta_max-theta_min)/(nbptst-1);
stepp=(phi_max-phi_min)/(nbptsp-1);
theta=theta_min:stept:theta_max;
phi=phi_min:stepp:phi_max;

% Coefficients de Mie pour la sphère
A=zeros(length(freq),ordre);
B=zeros(length(freq),ordre);
for n=1:ordre
N=(2*n+1)/n/(n+1);
pj=(-j)^n;
A(:,n) = -pj*N.*sphjn(n,k0a)./sphh1n(n,k0a);
B(:,n) = -j*pj*N.*(k0a.*sphjn(n-1,k0a)-n.*sphjn(n,k0a))...
./(k0a.*sphh1n(n-1,k0a)-n.*sphh1n(n,k0a));
end

% Coefficient S1 et S2 de Theta
S1=zeros(length(freq),length(theta));
S2=zeros(length(freq),length(theta));
pm=zeros(100+1,100+1);
pd=zeros(100+1,100+1);
for nT=1:length(theta)
ct=cos(theta(nT));
st=sin(theta(nT));
for n=1:ordre
[dumvar1,m,n,x,pm,pd]=lpmn(100,1,n,ct,pm,pd);
pn1=pm(m+1,n+1);
dpn1=pd(m+1,n+1);
pj=(0-1j)^(n+1);
S1(:,nT)=S1(:,nT)+pj.*(A(:,n).*pn1./st-(0+1j).*B(:,n).*st.*dpn1);
S2(:,nT)=S2(:,nT)+pj.*(-st.*A(:,n).*dpn1+(0+1j).*B(:,n).*pn1./st);
end
end

% Caclul de la SER selon theta et phi
sigmaTheta=zeros(length(freq),length(theta),length(phi));
sigmaPhi=zeros(length(freq),length(theta),length(phi));
for nT=1:length(theta)
for nP=1:length(phi)
sigmaTheta(:,nT,nP)=4*pi.*abs(S1(:,nT)).^2.*cos(phi(nP)).^2./k0.^2;
sigmaPhi(:,nT,nP)=4*pi*abs(S2(:,nT)).^2.*sin(phi(nP)).^2./k0.^2;
end
end

%Tracé des SER dans les plans E et H
figure(1);
[Xt,Yt] = meshgrid(k0a,theta_min*180/pi:(theta_max-theta_min)/(nbptst-1)*180/pi:theta_max*180/pi);
surf(Xt,Yt,sigmaPhi(:,:,nbptsp/2)'./pi./a^2);
zlabel('Normalized RCS Hplane');
ylabel('Thetas, scaterring angle in degrees');
xlabel('k0a');
title('Metallic sphere RCS function of thetas')

figure(2);
[Xt,Yt] = meshgrid(k0a,theta_min*180/pi:(theta_max-theta_min)/(nbptst-1)*180/pi:theta_max*180/pi);
surf(Xt,Yt,sigmaTheta(:,:,1)'./pi./a^2);
zlabel('Normalized RCS Eplane');
ylabel('Thetas, scaterring angle in degrees');
xlabel('k0a');
title('Metallic sphere RCS function of thetas')