理想导体圆柱对平面波的散射(Mie级数)

理想导体圆柱对平面波的散射1. TM z 极化假设TM z 极化均匀平面波垂直入射半径为a 的无限长PEC 圆柱，波的传播方向为+x ，入射电场和入射磁场用柱面波展开，分别表示为()cos 000ˆˆˆincjkxjk n jn z z z n n E e E e E j J k e ρϕϕρ∞---=-∞===∑E aaa(1)()()10011ˆˆ'incincn jn n jn n n n n E kE njJ k ej J k e j j j ϕϕρϕρρωμωμρωμ∞∞-+-=-∞=-∞=-∇⨯=-+∑∑HEaa (2)散射场朝外传播，因此，散射电场和磁场用柱第二类Hankel 函数展开，分别表示如下()()20ˆscaz nnn E a H k ρ∞=-∞=∑Ea(3)()()()()22001ˆˆ'scan nn n n n E da kE Hk a H k j d j ρϕρρωμρϕωμ∞∞=-∞=-∞=-+∑∑Haa (4)根据PEC 的边界条件，切向电场为0，可以得到()()()200njn n n n n E jJ k e a H k ϕρρ∞-=-∞⎡⎤+=⎣⎦∑(5)从而得到展开项的系数为()()()2n n jn n n J ka a j e H ka ϕ-=-，将系数带入展开式，得到散射电场和磁场的表达式为()()()()()202n scanjn znn nJ ka EE j H k e H ka ϕρ∞-=-∞=-∑(6)()()()()()2021n scanjn n n nJ ka E H j nH k e H ka ϕρρωμρ∞-=-∞=∑(7)()()()()()202'n scan jn n n n J ka kE H j H k e j H ka ϕϕρωμ∞-=-∞=-∑ (8)对于远区散射场，kρ → ∞，()()2n jk n H k eρρ-≈，则相应的电场和磁场为()()()2nsca jk jnzn nJ kaE E eH kaρϕ∞-=-∞=-∑(9)()()()2nsca jk jnn nJ kaH n eH kaρϕρ∞-=-∞=≈∑(10)()()()2nsca jk jnn nJ kaH eH kaρϕϕ∞-=-∞=∑(11) 利用远区散射电场计算散射宽度(Scattering Width)为()()()22222lim2scan jnincn nJ kaeH kaϕρλσπρπ∞→∞=-∞==∑EE(12)c**********************************************************************c Compute TMz Scattering from PEC Circular Cylinder by Mie Seriescc a INPUT, real(8)c On entry, 'a' specifies the radius of the circular cylinderc f INPUT, real(8)c On entry, 'f' specifies the incident frequencyc r INPUT, real(8)c On entry, 'r' specifies the distance between the observationc point and the origin of coordinatesc ph INPUT, real(8)c On entry, 'ph' specifies the observation anglec Ez OUTPUT, complex(8)c On exit, 'Ez' specifies the z component of the electricc scattering fieldc Hpho OUTPUT, complex(8)c On exit, 'Hpho' specifies the pho component of the magneticc scattering fieldc Hphi OUTPUT, complex(8)c On exit, 'Hphi' specifies the phi component of the magneticc scattering fieldcc Programmed by Panda Brewmasterc********************************************************************** subroutine dSca_TM_PEC_Cir_Cyl_Mie(a, f, r, ph, Ez, Hpho, Hphi)c********************************************************************** implicit nonec - Input Parametersreal(8) a, f, r, phcomplex(8) Ez, Hpho, Hphic - Constant Numbersreal(8), parameter :: pi = 3.141592653589793real(8), parameter :: eps0 = 8.854187817d-12real(8), parameter :: mu0 = pi * 4.d-7complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)c - Temporary Variablesinteger k, nmaxreal(8) eta0, wavek, ka, krreal(8), allocatable, dimension (:) :: Jnka, Ynka, Jnkr, Ynkrcomplex(8), allocatable, dimension (:) :: Hnka, Hnkrcomplex(8), allocatable, dimension (:) :: DHnkreta0 = dsqrt(mu0 / eps0)wavek = 2.d0 * pi * f * dsqrt(mu0 * eps0)ka = wavek * akr = wavek * rnmax = ka + 10.d0 * ka ** (1.d0 / 3.d0) + 1if(nmax < 5) nmax = 5allocate(Jnka(0 : nmax), Ynka(0 : nmax), Hnka(0 : nmax))call dBES(nmax + 1, ka, Jnka(0 : nmax), Ynka(0 : nmax))Hnka(0 : nmax) = dcmplx(Jnka(0 : nmax), - Ynka(0 : nmax))if(kr >= 0) then ! Finite Distanceallocate(Jnkr(- 1 : nmax + 1), Ynkr(- 1 : nmax + 1), Hnkr(- 1 : nmax + 1), DHnkr(0 : nmax))call dBES(nmax + 2, kr, Jnkr(0 : nmax + 1), Ynkr(0 : nmax + 1))Jnkr(- 1) = - Jnkr(1)Ynkr(- 1) = - Ynkr(1)Hnkr(- 1 : nmax + 1) = dcmplx(Jnkr(- 1 : nmax + 1), - Ynkr(- 1 : nmax + 1))Ez = - Jnka(0) / Hnka(0) * Hnkr(0)do k = 1, nmaxEz = Ez - cj ** (- k) * Jnka(k) / Hnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))Ez = Ez - cj ** k * Jnka(k) / Hnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph)) enddoHpho = 0.d0do k = 1, nmaxHpho = Hpho + cj ** (- k) * k * Jnka(k) / Hnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))Hpho = Hpho + cj ** k * (- k) * Jnka(k) / Hnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0,- k * ph)) enddo Hpho = Hpho / wavek / eta0 / r Hphi = - Jnka(0) / Hnka(0) * DHnkr(0) do k = 1, nmax DHnkr(k) = (Hnkr(k - 1) - Hnkr(k + 1)) / 2.d0 Hphi = Hphi - cj ** (- k) * Jnka(k) / Hnka(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph)) Hphi = Hphi - cj ** k * Jnka(k) / Hnka(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph)) enddo Hphi = Hphi / cj / eta0 deallocate(Jnkr, Ynkr, Hnkr, DHnkr) else ! Infinite Distance Ez = - Jnka(0) / Hnka(0) do k = 1, nmax Ez = Ez - 2.d0 * Jnka(k) / Hnka(k) * dcosd(k * ph) enddo Ez = Ez * dsqrt(2.d0 / pi / wavek) Hpho = 0.d0 Hphi = - Ez / eta0 endif deallocate(Jnka, Ynka, Hnka) end subroutine dSca_TM_PEC_Cir_Cyl_Miec######################################################################-50510152025306090120150180Angle (deg)S c a t t e r i n g W i d t h (d B )图1 PEC 圆柱的双站散射－TM z 极化2. TE z 极化假设TE z 极化均匀平面波垂直入射半径为a 的无限长PEC 圆柱，波的传播方向为+x ，入射电场和入射磁场用柱面波展开，分别表示为()cos 000ˆˆˆincjkxjk n jn z z z n n H e H e H j J k e ρϕϕρ∞---=-∞===∑Haaa(13)()()1001ˆˆ'incn jn n jn n n n n H kH njJ k ej J k e j j ϕϕρϕρρωερωε∞∞-+-=-∞=-∞=-∑∑Eaa(14)散射场朝外传播，因此，散射电场和磁场用柱第二类Hankel 函数展开，分别表示如下()()20ˆscaz n n n H a H k ρ∞=-∞=∑Ha(15)()()()()22001ˆˆ'scan nnnn n H a kH Hk a H k j j ρϕρρωερϕωε∞∞=-∞=-∞∂=-∂∑∑Eaa(16)根据PEC 的边界条件，切向电场为0，可以得到()()()20''0incscanjn n n n n kHE E jJ k e a H k j ϕϕϕρρωε∞-=-∞⎡⎤+=-+=⎣⎦∑(17)从而得到展开项的系数为()()()2''n n jn n n J ka a j e H ka ϕ-=-，将系数带入展开式，得到散射电场和磁场的表达式为()()()()()202''n scanjn znn n J ka HH j H k e H ka ϕρ∞-=-∞=-∑(18)()()()()()202'1'n scanjn n n n J ka H E nj H k e H ka ϕρρωερ∞-=-∞=-∑(19)()()()()()202'''n scanjn n n n J ka kH E jH k e j H ka ϕϕρωε∞-=-∞=∑(20)对于远区散射场，kρ → ∞，()()2n jk n H k eρρ-≈，则相应的散射磁场为()()()2''n sca jk jn zn n J ka HH e H ka ρϕ∞-=-∞=-∑(21)()()()2'0'n scajk jn n n J ka E n e H ka ρϕρ∞-=-∞=≈∑ (22)()()()2''nsca jk jnn nJ kaE H eH kaρϕϕη∞-=-∞=-∑(23) 则散射宽度(Scattering Width)为()()()2222'2lim2'scan jnincn nJ kaeH kaϕρλσπρπ∞→∞=-∞==∑HH(24)c**********************************************************************c Compute TEz Scattering from PEC Circular Cylinder by Mie Seriescc a INPUT, real(8)c On entry, 'a' specifies the radius of the circular cylinderc f INPUT, real(8)c On entry, 'f' specifies the incident frequencyc r INPUT, real(8)c On entry, 'r' specifies the distance between the observationc point and the origin of coordinatesc ph INPUT, real(8)c On entry, 'ph' specifies the observation anglec Hz OUTPUT, complex(8)c On exit, 'Hz' specifies the z component of the magneticc scattering fieldc Epho OUTPUT, complex(8)c On exit, 'Epho' specifies the pho component of the electricc scattering fieldc Ephi OUTPUT, complex(8)c On exit, 'Ephi' specifies the phi component of the electricc scattering fieldcc Programmed by Panda Brewmasterc********************************************************************** subroutine dSca_TE_PEC_Cir_Cyl_Mie(a, f, r, ph, Hz, Epho, Ephi)c********************************************************************** implicit nonec - Input Parametersreal(8) a, f, r, phcomplex(8) Hz, Epho, Ephic - Constant Numbersreal(8), parameter :: pi = 3.141592653589793d0real(8), parameter :: eps0 = 8.854187817620389d-12real(8), parameter :: mu0 = pi * 4.d-7complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)c - Temporary Variablesinteger k, nmaxreal(8) eta0, wavek, ka, krreal(8), allocatable, dimension (:) :: Jnka, Ynka, Jnkr, Ynkrreal(8), allocatable, dimension (:) :: DJnka, DJnkrcomplex(8), allocatable, dimension (:) :: Hnka, Hnkrcomplex(8), allocatable, dimension (:) :: DHnka, DHnkreta0 = dsqrt(mu0 / eps0)wavek = 2.d0 * pi * f * dsqrt(mu0 * eps0)ka = wavek * akr = wavek * rnmax = ka + 10.d0 * ka ** (1.d0 / 3.d0) + 1if(nmax < 5) nmax = 5allocate(Jnka(- 1 : nmax + 1), Ynka(- 1 : nmax + 1), Hnka(- 1 : nmax + 1))allocate(DJnka(0 : nmax), DHnka(0 : nmax))call dBES(nmax + 2, ka, Jnka(0 : nmax + 1), Ynka(0 : nmax + 1))Jnka(- 1) = - Jnka(1); Ynka(- 1) = - Ynka(1)Hnka(- 1 : nmax + 1) = dcmplx(Jnka(- 1 : nmax + 1), - Ynka(- 1 : nmax + 1))do k = 0, nmaxDJnka(k) = (Jnka(k - 1) - Jnka(k + 1)) / 2.d0DHnka(k) = (Hnka(k - 1) - Hnka(k + 1)) / 2.d0enddoif(kr >= 0) then ! Finite Distanceallocate(Jnkr(- 1 : nmax + 1), Ynkr(- 1 : nmax + 1), Hnkr(- 1 : nmax + 1))allocate(DHnkr(0 : nmax))call dBES(nmax + 2, kr, Jnkr(0 : nmax + 1), Ynkr(0 : nmax + 1))Jnkr(- 1) = - Jnkr(1); Ynkr(- 1) = - Ynkr(1)Hnkr(- 1 : nmax + 1) = dcmplx(Jnkr(- 1 : nmax + 1), - Ynkr(- 1 : nmax + 1))Hz = - DJnka(0) / DHnka(0) * Hnkr(0)do k = 1, nmaxHz = Hz - cj ** (- k) * DJnka(k) / DHnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))Hz = Hz - cj ** k * DJnka(k) / DHnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))enddoEpho = 0.d0do k = 1, nmaxEpho = Epho - cj ** (- k) * k * DJnka(k) / DHnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))Epho = Epho - cj ** k * (- k) * DJnka(k) / DHnka(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph)) enddo Epho = Epho * eta0 / wavek / r Ephi = DJnka(0) / DHnka(0) * DHnkr(0) do k = 1, nmax Ephi = Ephi + cj ** (- k) * DJnka(k) / DHnka(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph)) Ephi = Ephi + cj ** k * DJnka(k) / DHnka(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph)) enddo Ephi = Ephi * eta0 / cj deallocate(Jnkr, Ynkr, Hnkr, DHnkr) else ! Infinite Distance Hz = - DJnka(0) / DHnka(0) do k = 1, nmax Hz = Hz - 2.d0 * DJnka(k) / DHnka(k) * dcosd(k * ph) enddo Hz = Hz * dsqrt(2.d0 / pi / wavek) Epho = 0.d0 Ephi = Hz * eta0 endif deallocate(Jnka, Ynka, Hnka) deallocate(DJnka, DHnka) end subroutine dSca_TE_PEC_Cir_Cyl_Miec######################################################################-15-10-50510152025Angle (deg)S c a t t e r i n g W i d t h (d B )图2 PEC 圆柱的双站散射－TE z 极化。