均勻介質圓柱對平面波的散射1.TM極化假設TM極化均勻平面波垂直入射半徑為a的無限長均勻介質圓柱，相對介電常數為ε，相對磁導率為μ，波的傳播方向為+x，入射電場和入射磁場用柱面波展開，分別表示為(1)(2)散射場朝外傳播，因此，散射電場和磁場用柱第二類Hankel函數展開，分別表示如下(3)(4)透射場則由柱面基本波函數的線性組合表示，由于透射場在介質內部均為有限大，因此(5)(6)根據介質表面的邊界條件，切向電場和切向磁場連續，可以得到(7)(8)求解方程組，從而得到展開項的系數為(9)(10)令，將系數帶入展開式，得到散射電場和磁場的表達式為(11)(12)對于遠區散射場，kρ→∞，，則相應的電場和磁場為(13)(14)c**********************************************************************cComputeTMzScatteringfromhomogenerouslosslessdielectricCircularcCylinderbyMieSeriesccaINPUT,real(8)cOnentry,'a'specifiestheradiusofthecircularcylindercepsRINPUT,real(8)cOnentry,'epsR'specifiestherelativepermittivityofcthehomogenerousdielectriccircularcylindercMuRINPUT,real(8)cOnentry,'muR'specifiestherelativepermeabilityofcthehomogenerousdielectriccircularcylindercfINPUT,real(8)cOnentry,'f'specifiestheincidentfrequencycrINPUT,real(8)cOnentry,'r'specifiesthedistancebetweentheobservationcpointandtheoriginofcoordinatescphINPUT,real(8)cOnentry,'ph'specifiestheobservationanglecEzOUTPUT,complex(8)cOnexit,'Ez'specifiesthezcomponentoftheelectriccscatteringfieldcHphoOUTPUT,complex(8)cOnexit,'Hpho'specifiesthephocomponentofthemagneticcscatteringfieldcHphiOUTPUT,complex(8)cOnexit,'Hphi'specifiesthephicomponentofthemagneticcscatteringfieldccProgrammedbyPandaBrewmasterc**********************************************************************subroutinedSca_TM_DIE_Cir_Cyl_Mie(a,epsR,muR,f,r,ph,Ez,Hpho,Hphi)c**********************************************************************implicitnonec-InputParametersreal(8)a,epsR,muR,f,r,phcomplex(8)Ez,Hpho,Hphic-ConstantNumbersreal(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-TemporaryVariablesintegerk,nmaxreal(8)eta0,wavek0,ka0,krreal(8)eta1,wavek1,ka1complex(8)coe1,coe2real(8),allocatable,dimension(:)::Jnka0,Ynka0,DJnka0,Jnkr,Ynkrreal(8),allocatable,dimension(:)::Jnka1,Ynka1,DJnka1complex(8),allocatable,dimension(:)::Hnka0,DHnka0,Hnkr,DHnkrcomplex(8),allocatable,dimension(:)::Hnka1,DHnka1complex(8),allocatable,dimension(:)::aneta0=dsqrt(mu0/eps0)wavek0=2.d0*pi*f*dsqrt(mu0*eps0)ka0=wavek0*akr=wavek0*reta1=dsqrt(muR*mu0/epsR/eps0)wavek1=2.d0*pi*f*dsqrt(muR*mu0*epsR*eps0)ka1=wavek1*anmax=ka1+10.d0*ka1**(1.d0/3.d0)+1if(nmax<5)nmax=5allocate(Jnka0(-1:nmax+1),Ynka0(-1:nmax+1),Hnka0(-1:nmax+1),DJnka0(0:nmax),DHnka0(0:nmax))calldBES(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))dok=0,nmaxDJnka0(k)=(Jnka0(k-1)-Jnka0(k+1))/2.d0DHnka0(k)=(Hnka0(k-1)-Hnka0(k+1))/2.d0enddoallocate(Jnka1(-1:nmax+1),Ynka1(-1:nmax+1),Hnka1(-1:nmax+1),DJnka1(0:nmax),DHnka1(0:nmax))calldBES(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))dok=0,nmaxDJnka1(k)=(Jnka1(k-1)-Jnka1(k+1))/2.0DHnka1(k)=(Hnka1(k-1)-Hnka1(k+1))/2.0enddoallocate(an(0:nmax))dok=0,nmaxcoe1=(eta1*DJnka0(k)*Jnka1(k)-eta0*Jnka0(k)*DJnka1(k))coe2=(eta0*DJnka1(k)*Hnka0(k)-eta1*Jnka1(k)*DHnka0(k))an(k)=coe1/coe2enddoif(r>=0)then!FiniteDistanceallocate(Jnkr(-1:nmax+1),Ynkr(-1:nmax+1),Hnkr(-1:nmax+1),DHnkr(0:nmax))calldBES(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),-Yn