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);
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')
Catégories : Matlab · Radars
Tagué : Bistatic, RCS, scattering, SER, sphere
#!/usr/bin/python
# -*- coding: UTF-8 -*-
import numpy
from scipy import special
from scipy import *
from pylab import *
def sphjn(n,x):
return sqrt(pi/2./x)*special.jv(n+0.5,x)
def sphyn(n,x):
return sqrt(pi/2./x)*special.yv(n+0.5,x)
def sphh1n(n,x):
return sphjn(n,x)+1j*sphyn(n,x)
def sphh2n(n,x):
return sphjn(n,x)-1j*sphyn(n,x)
def MieCoeffSphPerfCond(a,ordre,freq):
c=3e8;
k0a=a*2./c*pi*freq;
A=numpy.zeros((freq.shape[0],ordre),complex)
B=numpy.zeros((freq.shape[0],ordre),complex)
for n in range(1,ordre+1):
N=float(2*n+1)/float(n)/float(n+1)
pj=pow(0-1j,n)
A[:,n-1] = -pj*N*sphjn(n,k0a)/sphh1n(n,k0a)
B[:,n-1] = (0-1j)*pj*N*(k0a*sphjn(n-1,k0a)-n*sphjn(n,k0a))/(k0a*sphh1n(n-1,k0a)-n*sphh1n(n,k0a))
return A, B
def Stheta(theta,A,B,ordre,freq):
S1=numpy.zeros((freq.shape[0],theta.shape[0]),complex)
S2=numpy.zeros((freq.shape[0],theta.shape[0]),complex)
for nT in range(theta.shape[0]) :
ct=cos(theta[nT])
st=sin(theta[nT])
for n in range(1,ordre+1):
pn=numpy.array(special.lpmn(1,n,ct))
pn1=pn[0][1][n]
dpn1=pn[1][1][n]
pj=pow((0-1j),n+1)
S1[:,nT]=S1[:,nT]+pj*(A[:,n-1]*pn1/st-(0+1j)*B[:,n-1]*st*dpn1)
S2[:,nT]=S2[:,nT]+pj*(-st*A[:,n-1]*dpn1+(0+1j)*B[:,n-1]*pn1/st)
return S1, S2
def bistaticRCSPerfCscipyondSph(a,ordre,theta,phi,freq):
c=3e8
k0=2.*pi*freq/c
A=numpy.zeros((freq.shape[0],ordre),complex)
B=numpy.zeros((freq.shape[0],ordre),complex)
A,B = MieCoeffSphPerfCond(a,ordre,freq)
S1=numpy.zeros((freq.shape[0],theta.shape[0]),complex)
S2=numpy.zeros((freq.shape[0],theta.shape[0]),complex)
S1,S2=Stheta(theta,A,B,ordre,freq)
sigmaTheta=numpy.zeros((freq.shape[0],theta.shape[0],phi.shape[0]))
sigmaPhi=numpy.zeros((freq.shape[0],theta.shape[0],phi.shape[0]))
for nT in range(theta.shape[0]):
for nP in range(phi.shape[0]):
sigmaTheta[:,nT,nP]=4*pi*pow(abs(S1[:,nT]),2)*pow(cos(phi[nP]),2)/pow(k0,2)
sigmaPhi[:,nT,nP]=4*pi*pow(abs(S2[:,nT]),2)*pow(sin(phi[nP]),2)/pow(k0,2)
return sigmaTheta, sigmaPhi
def test():
c=3e8
a=0.01
n=60
nbptsf=500
fmin=0.001e9
fmax=100e9
stepf=(fmax-fmin)/float(nbptsf)
freq=numpy.arange(fmin,fmax,stepf)
nbptst=100
tmin=0.001
tmax=pi
stept=(tmax-tmin)/(nbptst)
theta=numpy.arange(tmin,tmax,stept)
nbptsp=100
pmin=0.001
pmax=pi
stepp=(pmax-pmin)/(nbptsp)
phi= numpy.arange(pmin,pmax,stepp)
k0a=a*2.0*pi*freq/c
sigmaTheta, sigmaPhi=bistaticRCSPerfCscipyondSph(a,n,theta,phi,freq)
Eplane=sigmaTheta[:,0,0]/pi/a/a
Hplane=sigmaPhi[:,0,(nbptsp-1)/2]/pi/a/a
plot(k0a, Eplane, linewidth=1.0)
xlabel('k0a (Sphere circonference in wavelength)')
ylabel('Bistatic RCS/pia2 Normalized echo area')
title("E plane Bistatic Backscattering Cross Section, metallic sphere")
grid(True)
show()
plot(k0a, Hplane, linewidth=1.0)
xlabel('k0a (Sphere circonference in wavelength)')
ylabel('Bistatic RCS/pia2 Normalized echo area')
title("H plane Bistatic Backscattering Cross Section, metallic sphere")
grid(True)
show()
plot(k0a, sigmaTheta[:,(nbptst-1)/2,30]/pi/a/a, linewidth=1.0)
xlabel('k0a (Sphere circonference in wavelength)')
ylabel('Bistatic RCS/pia2 Normalized echo area')
title("H plane Bistatic Backscattering Cross Section, metallic sphere")
grid(True)
show()
test()
Catégories : Python · Radars
Tagué : Bistatic, RCS, scattering, SER, sphere
