均勻介質(zhì)圓柱對平面波的散射(Mie級數(shù))

2023年12月30日發(fā)(作者：重點人員管控措施)

均勻介質(zhì)圓柱對平面波的散射

1. TMz極化

假設(shè)TMz極化均勻平面波垂直入射半徑為a的無限長均勻介質(zhì)圓柱，相對介電常數(shù)為ε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函數(shù)展開，分別表示如下

Esca??zE0?anHn?an????2??k?? (3)

HscaE01?dankE0??2??2???????aHn?k???aanHn'?k??

??j???n???d?j??n??? (4)

透射場則由柱面基本波函數(shù)的線性組合表示，由于透射場在介質(zhì)內(nèi)部均為有限大，因此

Etran?zE0?bnJn?k1??

?an???? (5)

HtranE01?dbnkE0???????aJn?k1???abnJn'?k1??

??j???n???d?j??n??? (6)

根據(jù)介質(zhì)表面的邊界條件，切向電場和切向磁場連續(xù)，可以得到

?2?j?nJn?k??ejn??anHn?k???bnJn?k1??

(7)

(8)

1?j?nJ'n?k??ejn??1??2?anHn'?k???1?1bnJn'?k1??

求解方程組，從而得到展開項的系數(shù)為

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?，將系數(shù)帶入展開式，得到散射電?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?對于遠區(qū)散射場，kρ → ∞，Hn?k???2jn?jk?，則相應(yīng)的電場和磁場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 均勻介質(zhì)圓柱的雙站散射－TMz極化：(a) εr = 4，μr = 1；(b) εr = 4 ? j，μr =

1；

2. TEz極化

假設(shè)TEz極化均勻平面波垂直入射半徑為a的無限長均勻介質(zhì)圓柱，相對介電常數(shù)為ε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函數(shù)展開，分別表示如下

Hsca?zH0?an????anHn??2??k??

? (17)

EscaH01??ankH0?2?????aHk??a???n?j???n?????j??n????aH??'?k??

2nn (18)

透射場則由柱面基本波函數(shù)的線性組合表示，由于透射場在介質(zhì)內(nèi)部均為有限大，因此

Htran??zH0?an????bJ?k??

nn1 (19)

EtranH01?dbnkH0?????aJn?k1???a?j???n???d?j??n????bJ'?k??

nn1? (20)

根據(jù)介質(zhì)表面的邊界條件，切向電場和切向磁場連續(xù)，可以得到

?2?j?nJn?k??ejn??anHn?k???bnJn?k1??

(21)

(22)

?2??j?nJ'n?k??ejn???anHn'?k????1bnJn'?k1??

求解方程組，從而得到展開項的系數(shù)為

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?，將系數(shù)帶入展開式，得到散射電?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?對于遠區(qū)散射場，kρ → ∞，Hn?k???2jn?jk?，則相應(yīng)的散射磁場為

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均勻介質(zhì)圓柱的雙站散射－TEz極化：(a) εr = 4，μr = 1；(b) εr = 4 ? j，μr = 1；