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