均匀介质圆柱对平面波的散射(Mie级数)

2023年12月30日发(作者：重点人员管控措施)

均匀介质圆柱对平面波的散射

1. TMz极化

假设TMz极化均匀平面波垂直入射半径为a的无限长均匀介质圆柱，相对介电常数为εr，相对磁导率为μr，波的传播方向为+x，入射电场和入射磁场用柱面波展开，分别表示为

E1jincˆzE0eajkxˆzE0eajkcosˆzE0anjnJnkejn (1)

HincEincE01n1kE0njnˆˆanjJnkeajJ'nkejnjnjn(2)

散射场朝外传播，因此，散射电场和磁场用柱第二类Hankel函数展开，分别表示如下

EscaˆzE0anHnan2k (3)

HscaE01dankE022ˆˆaHnkaanHn'k

jndjn (4)

透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，因此

EtranˆzE0bnJnk1

an (5)

HtranE01dbnkE0ˆˆaJnk1abnJn'k1

jndjn (6)

根据介质表面的边界条件，切向电场和切向磁场连续，可以得到

2jnJnkejnanHnkbnJnk1

(7)

(8)

1jnJ'nkejn12anHn'k11bnJn'k1

求解方程组，从而得到展开项的系数为

anjn1Jn'kaJnk1aJnkaJn'k1ajne

22Jn'k1aHnka1Jnk1aHn'ka (9)

bnjnJnkaHn'kaJn'kaHn22kaJnk1aH2n2'k1aJn'k1aHnk1a1ejn (1

0)

ˆn令a1Jn'kaJnk1aJnkaJn'k1a，将系数带入展开式，得到散射电22Jn'k1aHnka1Jnk1aHn'ka场和磁场的表达式为

EscaˆzE0jnaˆnHnan2kejn (11)

HscaˆaE01nnjnˆnHna2kejnkE0n2ˆˆnHnaja'kejn(12)

jn2对于远区散射场，kρ → ∞，Hnk2jnjk，则相应的电场和磁场jek为

EscaˆzE0a2jjkˆnejn

eakn (13)

H

scakE02jjk2jjkjnˆˆˆneaˆnejn

aenaeakknnE01(14)

c**********************************************************************

c Compute TMz Scattering from homogenerous lossless dielectric Circular

c Cylinder by Mie Series

c

c a INPUT, real(8)

c On entry, 'a' specifies the radius of the circular cylinder

c epsR INPUT, real(8)

c On entry, 'epsR' specifies the relative permittivity of

c the homogenerous dielectric circular cylinder

c MuR INPUT, real(8)

c On entry, 'muR' specifies the relative permeability of

c the homogenerous dielectric circular cylinder

c f INPUT, real(8)

c On entry, 'f' specifies the incident frequency

c r INPUT, real(8)

c On entry, 'r' specifies the distance between the obrvation

c point and the origin of coordinates

c ph INPUT, real(8)

c On entry, 'ph' specifies the obrvation angle

c Ez OUTPUT, complex(8)

c On exit, 'Ez' specifies the z component of the electric

c scattering field

c Hpho OUTPUT, complex(8)

c On exit, 'Hpho' specifies the pho component of the magnetic

c scattering field

c Hphi OUTPUT, complex(8)

c On exit, 'Hphi' specifies the phi component of the magnetic

c scattering field

c

c Programmed by Panda Brewmaster

c**********************************************************************

subroutine dSca_TM_DIE_Cir_Cyl_Mie(a, epsR, muR, f, r, ph, Ez, Hpho, Hphi)

c**********************************************************************

implicit none

c - Input Parameters

real(8) a, epsR, muR, f, r, ph

complex(8) Ez, Hpho, Hphi

c - Constant Numbers

real(8), parameter :: pi = 3.9793

real(8), parameter :: eps0 = 8.854187817d-12

real(8), parameter :: mu0 = pi * 4.d-7

complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)

c - Temporary Variables

integer k, nmax

real(8) eta0, wavek0, ka0, kr

real(8) eta1, wavek1, ka1

complex(8) coe1, coe2

real(8), allocatable, dimension (:) :: Jnka0, Ynka0, DJnka0, Jnkr, Ynkr

real(8), allocatable, dimension (:) :: Jnka1, Ynka1, DJnka1

complex(8), allocatable, dimension (:) :: Hnka0, DHnka0, Hnkr, DHnkr

complex(8), allocatable, dimension (:) :: Hnka1, DHnka1

complex(8), allocatable, dimension (:) :: an

eta0 = dsqrt(mu0 / eps0)

wavek0 = 2.d0 * pi * f * dsqrt(mu0 * eps0)

ka0 = wavek0 * a

kr = wavek0 * r

eta1 = dsqrt(muR * mu0 / epsR / eps0)

wavek1 = 2.d0 * pi * f * dsqrt(muR * mu0 * epsR * eps0)

ka1 = wavek1 * a

nmax = ka1 + 10.d0 * ka1 ** (1.d0 / 3.d0) + 1

if(nmax < 5) nmax = 5

allocate(Jnka0(- 1 : nmax + 1), Ynka0(- 1 : nmax + 1), Hnka0(- 1 : nmax + 1), DJnka0(0 :

nmax), DHnka0(0 : nmax))

call dBES(nmax + 1, ka0, Jnka0(0 : nmax + 1), Ynka0(0 : nmax + 1))

Jnka0(- 1) = - Jnka0(1); Ynka0(- 1) = - Ynka0(1)

Hnka0(-1 : nmax + 1) = dcmplx(Jnka0(-1 : nmax + 1), - Ynka0(-1 : nmax + 1))

do k = 0, nmax

DJnka0(k) = (Jnka0(k - 1) - Jnka0(k + 1)) / 2.d0

DHnka0(k) = (Hnka0(k - 1) - Hnka0(k + 1)) / 2.d0

enddo

allocate(Jnka1(- 1 : nmax + 1), Ynka1(- 1 : nmax + 1), Hnka1(- 1 : nmax + 1), DJnka1(0 :

nmax), DHnka1(0 : nmax))

call dBES(nmax + 1, ka1, Jnka1(0 : nmax + 1), Ynka1(0 : nmax + 1))

Jnka1(- 1) = - Jnka1(1); Ynka1(- 1) = - Ynka1(1)

Hnka1(-1 : nmax + 1) = dcmplx(Jnka1(-1 : nmax + 1), - Ynka1(-1 : nmax + 1))

do k = 0, nmax

DJnka1(k) = (Jnka1(k - 1) - Jnka1(k + 1)) / 2.0

DHnka1(k) = (Hnka1(k - 1) - Hnka1(k + 1)) / 2.0

enddo

allocate(an(0 : nmax))

do k = 0, nmax

coe1 = (eta1 * DJnka0(k) * Jnka1(k) - eta0 * Jnka0(k) * DJnka1(k))

coe2 = (eta0 * DJnka1(k) * Hnka0(k) - eta1 * Jnka1(k) * DHnka0(k))

an(k) = coe1 / coe2

enddo

if(r >= 0) then ! Finite Distance

allocate(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 = an(0) * Hnkr(0)

do k = 1, nmax

Ez = Ez + cj ** (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Ez = Ez + cj ** k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hpho = 0.d0

do k = 1, nmax

Hpho = Hpho + cj ** (- k) * k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hpho = Hpho + cj ** k * (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hpho = - Hpho / wavek0 / eta0 / r

Hphi = 0.d0

do k = 1, nmax

Hphi = Hphi + cj ** (- k) * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hphi = Hphi + cj ** k * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hphi = - Hphi / eta0 / cj

deallocate(Jnkr, Ynkr, Hnkr, DHnkr)

el ! Infinite Distance

Ez = an(0)

do k = 1, nmax

Ez = Ez + 2.0 * an(k) * dcosd(k * ph)

enddo

Ez = Ez * dsqrt(2.d0 / pi / wavek0)

Hpho = 0.d0

Hphi = Ez / eta0

endif

deallocate(an)

deallocate(Jnka0, Ynka0, Hnka0, DJnka0, DHnka0)

deallocate(Jnka1, Ynka1, Hnka1, DJnka1, DHnka1)

end subroutine dSca_TM_DIE_Cir_Cyl_Mie

c######################################################################

c**********************************************************************

c Compute TMz Scattering from homogenerous lossy dielectric Circular

c Cylinder by Mie Series

c

c a INPUT, real(8)

c On entry, 'a' specifies the radius of the circular cylinder

c epsR INPUT, complex(8)

c On entry, 'epsR' specifies the relative permittivity of

c the homogenerous dielectric circular cylinder

c MuR INPUT, complex(8)

c On entry, 'muR' specifies the relative permeability of

c the homogenerous dielectric circular cylinder

c f INPUT, real(8)

c On entry, 'f' specifies the incident frequency

c r INPUT, real(8)

c On entry, 'r' specifies the distance between the obrvation

c point and the origin of coordinates

c ph INPUT, real(8)

c On entry, 'ph' specifies the obrvation angle

c Ez OUTPUT, complex(8)

c On exit, 'Ez' specifies the z component of the electric

c scattering field

c Hpho OUTPUT, complex(8)

c On exit, 'Hpho' specifies the pho component of the magnetic

c scattering field

c Hphi OUTPUT, complex(8)

c On exit, 'Hphi' specifies the phi component of the magnetic

c scattering field

c

c Programmed by Panda Brewmaster

c**********************************************************************

subroutine dSca_TM_LDIE_Cir_Cyl_Mie(a, epsR, muR, f, r, ph, Ez, Hpho, Hphi)

c**********************************************************************

implicit none

c - Input Parameters

real(8) a, f, r, ph

complex(8) epsR, muR, Ez, Hpho, Hphi

c - Constant Numbers

real(8), parameter :: pi = 3.9793d0

real(8), parameter :: eps0 = 8.854187817d-12

real(8), parameter :: mu0 = pi * 4.d-7

complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)

c - Temporary Variables

integer k, nmax

real(8) eta0, wavek0, ka0, kr

complex(8) eta1, wavek1, ka1, coe1, coe2

real(8), allocatable, dimension (:) :: Jnka0, Ynka0, DJnka0, Jnkr, Ynkr

complex(8), allocatable, dimension (:) :: Jnka1, Ynka1, DJnka1

complex(8), allocatable, dimension (:) :: Hnka0, DHnka0, Hnkr, DHnkr

complex(8), allocatable, dimension (:) :: Hnka1, DHnka1

complex(8), allocatable, dimension (:) :: an

eta0 = dsqrt(mu0 / eps0)

wavek0 = 2.d0 * pi * f * dsqrt(mu0 * eps0)

ka0 = wavek0 * a

kr = wavek0 * r

eta1 = cdsqrt(muR * mu0 / epsR / eps0)

wavek1 = 2.d0 * pi * f * cdsqrt(muR * mu0 * epsR * eps0)

ka1 = wavek1 * a

nmax = ka0 + 10.d0 * ka0 ** (1.d0 / 3.d0) + 1

if(nmax < 5) nmax = 5

allocate(Jnka0(- 1 : nmax + 1), Ynka0(- 1 : nmax + 1), Hnka0(- 1 : nmax + 1), DJnka0(0 :

nmax), DHnka0(0 : nmax))

call dBES(nmax + 1, ka0, Jnka0(0 : nmax + 1), Ynka0(0 : nmax + 1))

Jnka0(- 1) = - Jnka0(1); Ynka0(- 1) = - Ynka0(1)

Hnka0(-1 : nmax + 1) = dcmplx(Jnka0(-1 : nmax + 1), - Ynka0(-1 : nmax + 1))

do k = 0, nmax

DJnka0(k) = (Jnka0(k - 1) - Jnka0(k + 1)) / 2.d0

DHnka0(k) = (Hnka0(k - 1) - Hnka0(k + 1)) / 2.d0

enddo

allocate(Jnka1(- 1 : nmax + 1), Ynka1(- 1 : nmax + 1), Hnka1(- 1 : nmax + 1), DJnka1(0 :

nmax), DHnka1(0 : nmax))

call cdBES(nmax + 1, ka1, Jnka1(0 : nmax + 1), Ynka1(0 : nmax + 1))

Jnka1(- 1) = - Jnka1(1); Ynka1(- 1) = - Ynka1(1)

do k = -1, nmax + 1

Hnka1(k) = dcmplx(dreal(Jnka1(k)) + dimag(Ynka1(k)), dimag(Jnka1(k)) -

dreal(Ynka1(k)))

enddo

do k = 0, nmax

DJnka1(k) = (Jnka1(k - 1) - Jnka1(k + 1)) / 2.d0

DHnka1(k) = (Hnka1(k - 1) - Hnka1(k + 1)) / 2.d0

enddo

allocate(an(0 : nmax))

do k = 0, nmax

coe1 = (eta1 * DJnka0(k) * Jnka1(k) - eta0 * Jnka0(k) * DJnka1(k))

coe2 = (eta0 * DJnka1(k) * Hnka0(k) - eta1 * Jnka1(k) * DHnka0(k))

an(k) = coe1 / coe2

enddo

if(r >= 0) then ! Finite Distance

allocate(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 = an(0) * Hnkr(0)

do k = 1, nmax

Ez = Ez + cj ** (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Ez = Ez + cj ** k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hpho = 0.d0

do k = 1, nmax

Hpho = Hpho + cj ** (- k) * k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hpho = Hpho + cj ** k * (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hpho = - Hpho / wavek0 / eta0 / r

Hphi = 0.d0

do k = 1, nmax

Hphi = Hphi + cj ** (- k) * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hphi = Hphi + cj ** k * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Hphi = - Hphi / eta0 / cj

deallocate(Jnkr, Ynkr, Hnkr, DHnkr)

el ! Infinite Distance

Ez = an(0)

do k = 1, nmax

Ez = Ez + 2.d0 * an(k) * dcosd(k * ph)

enddo

Ez = Ez * dsqrt(2.d0 / pi / wavek0)

Hpho = 0.d0

Hphi = Ez / eta0

endif

deallocate(an)

deallocate(Jnka0, Ynka0, Hnka0, DJnka0, DHnka0)

deallocate(Jnka1, Ynka1, Hnka1, DJnka1, DHnka1)

end subroutine dSca_TM_LDIE_Cir_Cyl_Mie

c######################################################################

2520151050-5-10-150306090Angle (deg)Scattering

Width

(dB)f = 1 GHzf = 3 GHz120150180

(a)

25Scattering

Width

(dB)20151Angle (deg)f = 1 GHzf = 3 GHz120150180

(b)

图1 均匀介质圆柱的双站散射－TMz极化：(a) εr = 4，μr = 1；(b) εr = 4 − j，μr =

1；

2. TEz极化

假设TEz极化均匀平面波垂直入射半径为a的无限长均匀介质圆柱，相对介电常数为εr，相对磁导率为μr，波的传播方向为+x，入射电场和入射磁场用柱面波展开，分别表示为

HincˆzH0eajkxˆzH0eajkcosˆzH0anjnJnkejn (15)

EincH01n1kH0jnˆˆanjJkeanjnjnjnJ'nkejn (16)

散射场朝外传播，因此，散射电场和磁场用柱第二类Hankel函数展开，分别表示如下

HscaˆzH0ananHn2k

 (17)

EscaH01ankH02ˆˆaHkanjnjnaH'k

2nn (18)

透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，因此

HtranˆzH0anbJk

nn1 (19)

EtranH01dbnkH0ˆˆaJnk1ajndjnbJ'k

nn1 (20)

根据介质表面的边界条件，切向电场和切向磁场连续，可以得到

2jnJnkejnanHnkbnJnk1

(21)

(22)

2jnJ'nkejnanHn'k1bnJn'k1

求解方程组，从而得到展开项的系数为

anjnJn'kaJnk1a1JnkaJn'k1ajne

221Jn'k1aHnkaJnk1aHn'ka22JnkaHn'kaJn'kaHnka (23)

bnjnJnk1aHn22'k1a1Jn'k1aHnk1aejn (24)

ˆn令aJn'kaJnk1a1JnkaJn'k1a，将系数带入展开式，得到散射电221Jn'k1aHnkaJnk1aHn'ka场和磁场的表达式为

HH01scaˆzH0a2nˆnHnjna2kejn (25)

EscaˆanˆnHnjnankejnkH0ˆajn2ˆnHnjna'kejn (26)

2对于远区散射场，kρ → ∞，Hnk2jnjk，则相应的散射磁场为

jek (27)

HH01scaˆzH0a2jjkˆnejn

eaknE

scakH02jjk2jjkjnˆˆˆneaˆnejn

aenaeakknn(28)

c**********************************************************************

c Compute TEz Scattering from homogenerous lossless dielectric Circular

c cylinder by Mie Series

c

c a INPUT, real(8)

c On entry, 'a' specifies the radius of the circular cylinder

c epsR INPUT, real(8)

c On entry, 'epsR' specifies the relative permittivity of

c the homogenerous dielectric circular cylinder

c MuR INPUT, real(8)

c On entry, 'muR' specifies the relative permeability of

c the homogenerous dielectric circular cylinder

c f INPUT, real(8)

c On entry, 'f' specifies the incident frequency

c r INPUT, real(8)

c On entry, 'r' specifies the distance between the obrvation

c point and the origin of coordinates

c ph INPUT, real(8)

c On entry, 'ph' specifies the obrvation angle

c Hz OUTPUT, complex(8)

c On exit, 'Hz' specifies the z component of the magnetic

c scattering field

c Epho OUTPUT, complex(8)

c On exit, 'Epho' specifies the pho component of the electric

c scattering field

c Ephi OUTPUT, complex(8)

c On exit, 'Ephi' specifies the phi component of the electric

c scattering field

c

c Programmed by Panda Brewmaster

c**********************************************************************

subroutine dSca_TE_Die_Cir_Cyl_Mie(a, EpsR, MuR, f, r, ph, Hz, Epho, Ephi)

c**********************************************************************

implicit none

c - Input Parameters

real(8) a, EpsR, MuR, f, r, ph

complex(8) Hz, Epho, Ephi

c - Constant Numbers

real(8), parameter :: pi = 3.9793d0

real(8), parameter :: eps0 = 8.8549d-12

real(8), parameter :: mu0 = pi * 4.d-7

complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)

c - Temporary Variables

integer k, nmax

real(8) eta0, wavek0, ka0, kr

real(8) eta1, wavek1, ka1

complex(8) coe1, coe2

real(8), allocatable, dimension (:) :: Jnka0, Ynka0, DJnka0, Jnkr, Ynkr

real(8), allocatable, dimension (:) :: Jnka1, Ynka1, DJnka1

complex(8), allocatable, dimension (:) :: Hnka0, DHnka0, Hnkr, DHnkr

complex(8), allocatable, dimension (:) :: Hnka1, DHnka1

complex(8), allocatable, dimension (:) :: an

eta0 = dsqrt(mu0 / eps0)

wavek0 = 2.d0 * pi * f * dsqrt(mu0 * eps0)

ka0 = wavek0 * a

kr = wavek0 * r

eta1 = dsqrt(MuR * mu0 / EpsR / eps0)

wavek1 = 2.0 * pi * f * dsqrt(MuR * mu0 * EpsR * eps0)

ka1 = wavek1 * a

nmax = ka1 + 10.0 * ka1 ** (1.0 / 3.0) + 1

if(nmax < 5) nmax = 5

allocate(Jnka0(- 1 : nmax + 1), Ynka0(- 1 : nmax + 1), Hnka0(- 1 : nmax + 1), DJnka0(0 :

nmax), DHnka0(0 : nmax))

call dBES(nmax + 1, ka0, Jnka0(0 : nmax + 1), Ynka0(0 : nmax + 1))

Jnka0(- 1) = - Jnka0(1); Ynka0(- 1) = - Ynka0(1)

Hnka0(-1 : nmax + 1) = dcmplx(Jnka0(-1 : nmax + 1), - Ynka0(-1 : nmax + 1))

do k = 0, nmax

DJnka0(k) = (Jnka0(k - 1) - Jnka0(k + 1)) / 2.d0

DHnka0(k) = (Hnka0(k - 1) - Hnka0(k + 1)) / 2.d0

enddo

allocate(Jnka1(- 1 : nmax + 1), Ynka1(- 1 : nmax + 1), Hnka1(- 1 : nmax + 1), DJnka1(0 :

nmax), DHnka1(0 : nmax))

call dBES(nmax + 1, ka1, Jnka1(0 : nmax + 1), Ynka1(0 : nmax + 1))

Jnka1(- 1) = - Jnka1(1); Ynka1(- 1) = - Ynka1(1)

Hnka1(-1 : nmax + 1) = dcmplx(Jnka1(-1 : nmax + 1), - Ynka1(-1 : nmax + 1))

do k = 0, nmax

DJnka1(k) = (Jnka1(k - 1) - Jnka1(k + 1)) / 2.0

DHnka1(k) = (Hnka1(k - 1) - Hnka1(k + 1)) / 2.0

enddo

allocate(an(0 : nmax))

do k = 0, nmax

coe1 = eta0 * DJnka0(k) * Jnka1(k) - eta1 * Jnka0(k) * DJnka1(k)

coe2 = eta1 * DJnka1(k) * Hnka0(k) - eta0 * Jnka1(k) * DHnka0(k)

an(k) = coe1 / coe2

enddo

if(r >= 0) then ! Finite Distance

allocate(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 = an(0) * Hnkr(0)

do k = 1, nmax

Hz = Hz + cj ** (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hz = Hz + cj ** k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Epho = 0.d0

do k = 1, nmax

Epho = Epho + cj ** (- k) * k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Epho = Epho + cj ** k * (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Epho = - Epho * eta0 / wavek0 / r

Ephi = an(0) * DHnkr(0)

do k = 1, nmax

Ephi = Ephi + cj ** (- k) * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Ephi = Ephi + cj ** k * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Ephi = Ephi * eta0 / cj

deallocate(Jnkr, Ynkr, Hnkr, DHnkr)

el ! Infinite Distance

Hz = an(0)

do k = 1, nmax

Hz = Hz + 2.0 * an(k) * dcosd(k * ph)

enddo

Hz = Hz * dsqrt(2.0 / pi / wavek0)

Epho = 0.d0

Ephi = Hz * eta0

endif

deallocate(an)

deallocate(Jnka0, Ynka0, Hnka0, DJnka0, DHnka0)

deallocate(Jnka1, Ynka1, Hnka1, DJnka1, DHnka1)

end subroutine dSca_TE_Die_Cir_Cyl_Mie

c######################################################################

c**********************************************************************

c Compute TEz Scattering from homogenerous lossy dielectric Circular

c Cylinder by Mie Series

c

c a INPUT, real(8)

c On entry, 'a' specifies the radius of the circular cylinder

c epsR INPUT, complex(8)

c On entry, 'epsR' specifies the relative permittivity of

c the homogenerous dielectric circular cylinder

c MuR INPUT, complex(8)

c On entry, 'muR' specifies the relative permeability of

c the homogenerous dielectric circular cylinder

c f INPUT, real(8)

c On entry, 'f' specifies the incident frequency

c r INPUT, real(8)

c On entry, 'r' specifies the distance between the obrvation

c point and the origin of coordinates

c ph INPUT, real(8)

c On entry, 'ph' specifies the obrvation angle

c Hz OUTPUT, complex(8)

c On exit, 'Hz' specifies the z component of the magnetic

c scattering field

c Epho OUTPUT, complex(8)

c On exit, 'Epho' specifies the pho component of the electric

c scattering field

c Ephi OUTPUT, complex(8)

c On exit, 'Ephi' specifies the phi component of the electric

c scattering field

c

c Programmed by Panda Brewmaster

c**********************************************************************

subroutine dSca_TE_LDie_Cir_Cyl_Mie(a, EpsR, MuR, f, r, ph, Hz, Epho, Ephi)

c**********************************************************************

implicit none

c - Input Parameters

real(8) a, f, r, ph

complex(8) EpsR, MuR, Hz, Epho, Ephi

c - Constant Numbers

real(8), parameter :: pi = 3.9793d0

real(8), parameter :: eps0 = 8.854187817d-12

real(8), parameter :: mu0 = pi * 4.d-7

complex(8), parameter :: cj = dcmplx(0.d0, 1.d0)

c - Temporary Variables

integer k, nmax

real(8) eta0, wavek0, ka0, kr

complex(8) eta1, wavek1, ka1, coe1, coe2

real(8), allocatable, dimension (:) :: Jnka0, Ynka0, DJnka0, Jnkr, Ynkr

complex(8), allocatable, dimension (:) :: Jnka1, Ynka1, DJnka1

complex(8), allocatable, dimension (:) :: Hnka0, DHnka0, Hnkr, DHnkr

complex(8), allocatable, dimension (:) :: Hnka1, DHnka1

complex(8), allocatable, dimension (:) :: an

eta0 = dsqrt(mu0 / eps0)

wavek0 = 2.d0 * pi * f * dsqrt(mu0 * eps0)

ka0 = wavek0 * a

kr = wavek0 * r

eta1 = cdsqrt(MuR * mu0 / EpsR / eps0)

wavek1 = 2.d0 * pi * f * cdsqrt(MuR * mu0 * EpsR * eps0)

ka1 = wavek1 * a

nmax = ka0 + 10.d0 * ka0 ** (1.d0 / 3.d0) + 1

if(nmax < 5) nmax = 5

allocate(Jnka0(- 1 : nmax + 1), Ynka0(- 1 : nmax + 1), Hnka0(- 1 : nmax + 1), DJnka0(0 :

nmax), DHnka0(0 : nmax))

call dBES(nmax + 1, ka0, Jnka0(0 : nmax + 1), Ynka0(0 : nmax + 1))

Jnka0(- 1) = - Jnka0(1); Ynka0(- 1) = - Ynka0(1)

Hnka0(-1 : nmax + 1) = dcmplx(Jnka0(-1 : nmax + 1), - Ynka0(-1 : nmax + 1))

do k = 0, nmax

DJnka0(k) = (Jnka0(k - 1) - Jnka0(k + 1)) / 2.d0

DHnka0(k) = (Hnka0(k - 1) - Hnka0(k + 1)) / 2.d0

enddo

allocate(Jnka1(- 1 : nmax + 1), Ynka1(- 1 : nmax + 1), Hnka1(- 1 : nmax + 1), DJnka1(0 :

nmax), DHnka1(0 : nmax))

call cdBES(nmax + 1, ka1, Jnka1(0 : nmax + 1), Ynka1(0 : nmax + 1))

Jnka1(- 1) = - Jnka1(1); Ynka1(- 1) = - Ynka1(1)

do k = - 1, nmax + 1

Hnka1(k) = dcmplx(dreal(Jnka1(k)) + dimag(Ynka1(k)), dimag(Jnka1(k)) -

dreal(Ynka1(k)))

enddo

do k = 0, nmax

DJnka1(k) = (Jnka1(k - 1) - Jnka1(k + 1)) / 2.d0

DHnka1(k) = (Hnka1(k - 1) - Hnka1(k + 1)) / 2.d0

enddo

allocate(an(0 : nmax))

do k = 0, nmax

coe1 = eta0 * DJnka0(k) * Jnka1(k) - eta1 * Jnka0(k) * DJnka1(k)

coe2 = eta1 * DJnka1(k) * Hnka0(k) - eta0 * Jnka1(k) * DHnka0(k)

an(k) = coe1 / coe2

enddo

if(r >= 0) then ! Finite Distance

allocate(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 = an(0) * Hnkr(0)

do k = 1, nmax

Hz = Hz + cj ** (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Hz = Hz + cj ** k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Epho = 0.d0

do k = 1, nmax

Epho = Epho + cj ** (- k) * k * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Epho = Epho + cj ** k * (- k) * an(k) * Hnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Epho = - Epho * eta0 / wavek0 / r

Ephi = an(0) * DHnkr(0)

do k = 1, nmax

Ephi = Ephi + cj ** (- k) * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, k * ph))

Ephi = Ephi + cj ** k * an(k) * DHnkr(k) * cdexp(dcmplx(0.d0, - k * ph))

enddo

Ephi = Ephi * eta0 / cj

deallocate(Jnkr, Ynkr, Hnkr, DHnkr)

el ! Infinite Distance

Hz = an(0)

do k = 1, nmax

Hz = Hz + 2.d0 * an(k) * dcosd(k * ph)

enddo

Hz = Hz * dsqrt(2.d0 / pi / wavek0)

endif

deallocate(an)

deallocate(Jnka0, Ynka0, Hnka0, DJnka0, DHnka0)

deallocate(Jnka1, Ynka1, Hnka1, DJnka1, DHnka1)

end subroutine dSca_TE_LDie_Cir_Cyl_Mie

c######################################################################

2520151050-5-10-150306090Angle (deg)Scattering

Width

(dB)f = 1 GHzf = 3 GHz120150180

(a)

Scattering

Width

(dB)20100-10-20-30-400306090Angle (deg)f = 1 GHzf = 3 GHz120150180

(b)

图2均匀介质圆柱的双站散射－TEz极化：(a) εr = 4，μr = 1；(b) εr = 4 − j，μr = 1；