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