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')

SER bistatique d’une sphère métallique

#!/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()

SER d’une sphère parfaitement conductrice

#!/usr/bin/python
# -*- coding: UTF-8 -*-

import numpy
from scipy import special
from cmath 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 sigBCSPerfCondSph(a,ordre,freq):
    c=3e8
    k0=2.*pi*freq/c
    A=numpy.zeros((freq.shape[0],ordre),complex)
    B=numpy.zeros((freq.shape[0],ordre),complex)
    F0=numpy.zeros(freq.shape,complex)
    A,B = MieCoeffSphPerfCond(a,ordre,freq)
    for  n in range(1,ordre+1):
        pj=pow((0-1j),n-1)
        F0=F0-pj*float(n)*float(n+1)*(A[:,n-1]+(0+1j)*B[:,n-1])/2.0
    sigma=4*pi*abs(F0)*abs(F0)/k0/k0
    return sigma

def test():
    c=3e8
    a=0.1
    n=60
    freq=numpy.arange(0.1e9,9.1e9,0.01e9)
    k0a=a*2.0*pi*freq/c
    sigma=sigBCSPerfCondSph(a,n,freq)/pi/a/a
    plot(k0a, sigma, linewidth=1.0)
    xlabel('k0a (Sphere circonference in wavelength)')
    ylabel('BCS/pia2 Normalized echo area')
    title("Backscattering Cross Section, metallic sphere")
    grid(True)
    show()  

test()