均匀介质圆柱对TE和TM平面波的散射(附程序)

均匀介质圆柱对TE和TM平面波的散射1、求解均匀介质圆柱对TM波的雷达散射截面以及表面电、磁流假设TM极化均匀平面波垂直入射半径为a的无限长均匀介质圆柱，其中介质圆柱沿z轴放置，相对介电常数为εr，相对磁导率为μr，波的传播方向如图所示为+x方向。入射电场用柱面波展开，可表示为 (1)由Maxwell方程,得到(2)其中，为真空中的磁导率,为真空中的波数。当时，介质外散射场朝外传播。因此，散射电场用柱第二类Hankel函数展开，表示如下(3)同理由Maxwell方程,得到（4）当时，介质内为透射场，透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，则自变量为0时，纽曼函数和汉克尔函数趋于无穷，因此电场和磁场只选用贝塞尔函数表示如下 (5) (6)其中，为介质中的磁导率,为介质中的波数。当时，根据介质表面的边界条件，切向电场和切向磁场连续，可以得到电场边界条件则有 (7)磁场边界条件则有 (8)求解方程组，从而得到展开项的系数为 (9) (10)其中，因为，所以为真空中波阻抗的倒数,为介质中波阻抗的倒数。 如下求介质表面的电、磁流。求表面电流由边界条件有所以（13）求表面磁流由边界条件有所以（14）另外对于远区散射场，kρ→∞，则散射电场为（11）又将上式代入二维雷达散射截面的定义式有（12）MATLAB编程求解%**********TM波照射无限长介质柱*******%***********初始化************clearall;closeall;clc;tic;wlen=1.0;k0=2.0*pi/wlen;ur=1;%相对磁导率er=4*10^12;%相对介电常数eta0=120.0*pi;%自由空间波阻抗eta1=eta0*sqrt(ur/er);%介质波阻抗radius=10.0;Npwave=10;NPL=2.0*pi*radius*Npwavepalen=2.0*pi/NPL;ka=k0*radius;ka1=ka*sqrt(ur*er);%介质波数%!*********计算真空和介质中贝塞尔函数和汉克尔函数的值，存储到数组中************ jn0=besselj(0,ka);%真空中h2n0=besselh(0,2,ka);jn(1)=besselj(1,ka); h2n(1)=besselh(1,2,ka);jn(2)=2.0*jn(1)/ka-jn0;h2n(2)=2.0*h2n(1)/ka-h2n0;Djn0=besselj(0,ka1);%介质中Djn(1)=besselj(1,ka1);Djn(2)=2.0*Djn(1)/ka1-Djn0;forn=3:2000000 jn(n)=2.0*(n-1.0)*jn(n-1)/ka-jn(n-2);%真空中贝塞尔函数 h2n(n)=2.0*(n-1.0)*h2n(n-1)/ka-h2n(n-2);%真空中汉克尔函数Djn(n)=2.0*(n-1.0)*Djn(n-1)/ka1-Djn(n-2);%介质中贝塞尔函数if(abs(h2n(n))>1.0*10^10)break;endendntotal=n%!*****计算真空和介质中贝塞尔函数和汉克尔函数的一次导的值，存储到数组中*****Jn0=-jn(1);%真空中贝塞尔函数的一次导H2n0=-h2n(1);%真空中汉克尔函数的一次导DJn0=-Djn(1);%介质中贝塞尔函数的一次导Jn(1)=(jn0-jn(2))/2;H2n(1)=(h2n0-h2n(2))/2;DJn(1)=(Djn0-Djn(2))/2;jn(ntotal+1)=2.0*ntotal*jn(ntotal)/ka-jn(ntotal-1);h2n(ntotal+1)=2.0*ntotal*h2n(ntotal)/ka-h2n(ntotal-1);Djn(ntotal+1)=2.0*ntotal*Djn(ntotal)/ka1-Djn(ntotal-1);forn=2:ntotalJn(n)=(jn(n-1)-jn(n+1))/2;H2n(n)=(h2n(n-1)-h2n(n+1))/2;DJn(n)=(Djn(n-1)-Djn(n+1))/2;end%!*********计算an**********a0=(eta0*DJn0*jn0-eta1*Jn0*Djn0)/(eta1*H2n0*Djn0-eta0*h2n0*DJn0);forn=1:ntotala(n)=(eta0*DJn(n)*jn(n)-eta1*Jn(n)*Djn(n))/...(eta1*H2n(n)*Djn(n)-eta0*h2n(n)*DJn(n));end%!*******************计算每个节点处的解析电流**********form=1:NPL pp=(m-0.5)*palen;%散射角Φ cp1=sin(pp)-i*cos(pp); cp2=-sin(pp)-i*cos(pp);ctemp1=1.0;ctemp2=1.0;sumc0=Jn0+a0*H2n0; sumf0=jn0+a0*h2n0; forn=1:ntotal ctemp1=ctemp1*cp1; ctemp2=ctemp2*cp2; sumc0=sumc0+(Jn(n)+a(n)*H2n(n))*(ctemp1+ctemp2); sumf0=sumf0+(jn(n)+a(n)*h2n(n))*(ctemp1+ctemp2);endJc(m)=sumc0/eta0*(-j);%电流密度Jf(m)=sumf0;%磁流密度endfigure,plot(abs(Jc));grid;title('电流分布');xlabel('散射角Φ(弧度)');ylabel('电流密度(A/m)');figure,plot(abs(Jf));grid;title('磁流分布');xlabel('散射角Φ(弧度)');ylabel('磁流密度(Wb/m)');%*********************计算每个散射方向的RCS*********************************form=1:360 pp=m*2.0*pi/360.0;cp1=cos(pp)+i*sin(pp); cp2=cos(pp)-i*sin(pp); ctemp1=1.0;ctemp2=1.0; sum=a0; forn=1:ntotal ctemp1=ctemp1*cp1; ctemp2=ctemp2*cp2;sum=sum+(ctemp1+ctemp2)*a(n); endRCS_exact(m)=10.0*log10(4.0/k0*abs(sum)^2);endfigure,plot(abs(RCS_exact));grid;title('RCS分布');xlabel('散射角Φ(度)');ylabel('RCS(dbm)');toc;运行结果：NPL=628.3185ntotal=100Elapsedtimeis1.295000seconds.2、求解均匀介质圆柱对TE波的雷达散射截面以及表面电、磁流假设TE极化均匀平面波垂直入射半径为a的无限长均匀介质圆柱，其中介质圆柱沿z轴放置，相对介电常数为εr，相对磁导率为μr，波的传播方向如图所示为+x方向。入射磁场用柱面波展开，可表示为 (1)由Maxwell方程,得到(2)其中，为真空中的介电常数,为真空中的波数。当时，介质外散射场朝外传播。因此，散射磁场用柱第二类Hankel函数展开，表示如下(3)同理由Maxwell方程,得到（4）当时，介质内为透射场，透射场则由柱面基本波函数的线性组合表示，由于透射场在介质内部均为有限大，则自变量为0时，纽曼函数和汉克尔函数趋于无穷，因此电场和磁场只选用贝塞尔函数表示如下 (5) (6)其中，为介质中的介电常数,为介质中的波数。当时，根据介质表面的边界条件，切向电场和切向磁场连续，可以得到磁场边界条件则有 (7)电场边界条件则有 (8)求解方程组，从而得到展开项的系数为 (9) (10)其中，因为，所以为真空中波阻抗,为介质中波阻抗。 如下求介质表面的电、磁流。1）求表面电流由边界条件有所以（13）2)求表面磁流由边界条件有所以（14）另外对于远区散射场，kρ→∞，则散射电场为（11）又将上式代入二维雷达散射截面的定义式有（12）MATLAB编程求解%**********TE波照射无限长介质柱*******%***********初始化************clearall;closeall;clc;tic;wlen=1.0;k0=2.0*pi/wlen;ur=1;%相对磁导率er=4*10^12;%相对介电常数eta0=120.0*pi;%自由空间波阻抗eta1=eta0*sqrt(ur/er);%介质波阻抗radius=10.0;Npwave=10;NPL=2.0*pi*radius*Npwavepalen=2.0*pi/NPL;ka=k0*radius;ka1=ka*sqrt(ur*er);%介质波数%!*********计算真空和介质中贝塞尔函数和汉克尔函数的值，存储到数组中************ jn0=besselj(0,ka);%真空中h2n0=besselh(0,2,ka);jn(1)=besselj(1,ka); h2n(1)=besselh(1,2,ka);jn(2)=2.0*jn(1)/ka-jn0;h2n(2)=2.0*h2n(1)/ka-h2n0;Djn0=besselj(0,ka1);%介质中Djn(1)=besselj(1,ka1);Djn(2)=2.0*Djn(1)/ka1-Djn0;forn=3:2000000 jn(n)=2.0*(n-1.0)*jn(n-1)/ka-jn(n-2);%真空中贝塞尔函数 h2n(n)=2.0*(n-1.0)*h2n(n-1)/ka-h2n(n-2);%真空中汉克尔函数Djn(n)=2.0*(n-1.0)*Djn(n-1)/ka1-Djn(n-2);%介质中贝塞尔函数if(abs(h2n(n))>1.0*10^10)break;endendntotal=n%!*****计算真空和介质中贝塞尔函数和汉克尔函数的一次导的值，存储到数组中*****Jn0=-jn(1);%真空中贝塞尔函数的一次导H2n0=-h2n(1);%真空中汉克尔函数的一次导DJn0=-Djn(1);%介质中贝塞尔函数的一次导Jn(1)=(jn0-jn(2))/2;H2n(1)=(h2n0-h2n(2))/2;DJn(1)=(Djn0-Djn(2))/2;jn(ntotal+1)=2.0*ntotal*jn(ntotal)/ka-jn(ntotal-1);h2n(ntotal+1)=2.0*ntotal*h2n(ntotal)/ka-h2n(ntotal-1);Djn(ntotal+1)=2.0*ntotal*Djn(ntotal)/ka1-Djn(ntotal-1);forn=2:ntotalJn(n)=(jn(n-1)-jn(n+1))/2;H2n(n)=(h2n(n-1)-h2n(n+1))/2;DJn(n)=(Djn(n-1)-Djn(n+1))/2;end%!*********计算an**********a0=(eta0*Djn0*Jn0-eta1*jn0*DJn0)/(eta1*h2n0*DJn0-eta0*H2n0*Djn0);forn=1:ntotala(n)=(eta0*Djn(n)*Jn(n)-eta1*jn(n)*DJn(n))/...(eta1*h2n(n)*DJn(n)-eta0*H2n(n)*Djn(n));end%!*******************计算每个节点处的解析电流**********form=1:NPL pp=(m-0.5)*palen;%散射角Φ cp1=sin(pp)-i*cos(pp); cp2=-sin(pp)-i*cos(pp);ctemp1=1.0;ctemp2=1.0;sumc0=jn0+a0*h2n0; sumf0=Jn0+a0*H2n0; forn=1:ntotal ctemp1=ctemp1*cp1; ctemp2=ctemp2*cp2; sumc0=sumc0+(jn(n)+a(n)*h2n(n))*(ctemp1+ctemp2); sumf0=sumf0+(Jn(n)+a(n)*H2n(n))*(ctemp1+ctemp2);