RCS-FDTD方柱-球的RCS程序

1、精选优质文档-倾情为你奉上% FDTD_2D_RCSclear;clc;%*初始化*%tic;% v=3*108; % 波速f=15*10(9); % 频率lamda=v/f; % 波长k=2*pi/lamda; % 波数epsz=1/(4*pi*9*109); % 真空介电常数mu=4*pi*10.(-7); % 真空磁导率Z=sqrt(mu/epsz); % 真空波阻抗epsilon=1; % 相对介电常数sigma=0.0; % 电导率%N=120; % 网格数量L=1000; % 迭代次数ddx=lamda/20; % 网格尺寸dt=ddx/(2*v); % 时间间隔lstart=1;

2、 % 空间起始行坐标 lend=N; % 空间终止行坐标rstart=1; % 空间起始列坐标rend=N; % 空间终止列坐标ia=N/4; % 总场区域x左ib=3*N/4; % 总场区域x右ja=ia; % 总场区域x下jb=ib; % 总场区域x上pa=ia-5; % 外推区域x左pb=ib+5; % 外推区域x右qa=pa; % 外推区域x下qb=pb; % 外推区域x上length=pb-pa+1; % 每个边上的总长度npml=N/8; % PML点数% 方柱参数side=2*lamda; % 方柱边长halfgrid=round(side/(2*ddx); % 方柱占据网格大小

3、% 媒质参数for i=lstart:lend; % 控制媒质分部区域 for j=rstart:rend; ga(i,j)=1/(epsilon+sigma*dt/epsz); % 求和参量 gb(i,j)=sigma*dt/epsz; % 求和参量 end;end;% 源参数spread=8; % 脉冲宽度t0=25; % 脉冲高度is=N/2; % 源的X位置js=N/2; % 源的Y位置% 输入平面波ez_inc=zeros(1,N);hx_inc=zeros(1,N);% 平面波吸收条件变量初始化ez_v1=0;ez_v2=0;ez_v3=0;ez_v4=0;% 迭待电磁场参量dz=

4、zeros(N,N); % z方向电荷密度ez=zeros(N,N); % z方向电场iz=zeros(N,N); % z方向电场求和参量hx=zeros(N,N); % x方向磁场hy=zeros(N,N); % y方向磁场ihx=zeros(N,N); % x方向磁场参量ihy=zeros(N,N); % y方向磁场参量% 相位和幅度提取(傅里叶变换)ine=0;inh=0;ez_out=zeros(4,length+1);hx_out=zeros(4,length);hy_out=zeros(4,length);%*PML设置*% PML初始化for i=1:N; gi2(i)=1; g

5、i3(i)=1; fi1(i)=0; fi2(i)=1; fi3(i)=1;end;for j=1:N; gj2(j)=1; gj3(j)=1; fj1(j)=0; fj2(j)=1; fj3(j)=1;end;% 阻抗渐变设置for i=1:npml+1; xnum=npml-i+1; xxn=xnum/npml; xn=0.33*(xxn3); gi2(i)=1/(1+xn); gi2(N-i+1)=1/(1+xn); gi3(i)=(1-xn)/(1+xn); gi3(N-i+1)=(1-xn)/(1+xn); xxn=(xnum-0.5)/npml; xn=0.25*(xxn3); f

6、i1(i)=xn; fi1(N-i)=xn; fi2(i)=1/(1+xn); fi2(N-i)=1/(1+xn); fi3(i)=(1-xn)/(1+xn); fi3(N-i)=(1-xn)/(1+xn);end;for j=1:npml+1; xnum=npml-j+1; xxn=xnum/npml; xn=0.33*(xxn3); gj2(j)=1/(1+xn); gj2(N-j+1)=1/(1+xn); gj3(j)=(1-xn)/(1+xn); gj3(N-j+1)=(1-xn)/(1+xn); xxn=(xnum-0.5)/npml; xn=0.25*(xxn3); fj1(j)=

7、xn; fj1(N-j)=xn; fj2(j)=1/(1+xn); fj2(N-j)=1/(1+xn); fj3(j)=(1-xn)/(1+xn); fj3(N-j)=(1-xn)/(1+xn);end;%*迭代求解电场和磁场*%for T=1:L; disp(T);% %*电场由磁场更新*% TM波Y方向传播 for j=2:N-1; ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j); end; % 平面波吸收条件 ez_inc(1)=ez_v2; ez_v2=ez_v1; ez_v1=ez_inc(2); ez_inc(N)=ez_v3; ez_

8、v3=ez_v4; ez_v4=ez_inc(N-1); % 入射电场傅里叶变换(必须用边界加入) % 注意：当计算频谱用的是其他的位置的入射波时候，会产生很大的误差； ine=ine+exp(-sqrt(-1)*2*pi*f*T*dt)*ez_inc(ja-1); % 电荷密度Dz for i=2:N-1; % 为了使每个电场周围都有磁场进行数组下标处理 for j=2:N; dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*( hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1) ); end; end; % 脉冲或正弦源的加

9、入 source=exp(-0.5)*( (t0-T)/spread ).2); % 脉冲源 ez_inc(5)=source; % TM平面波源 % 电荷密度Dz加入射场 for i=ia:ib; dz(i,ja)=dz(i,ja)+0.5*hx_inc(ja-1); dz(i,jb)=dz(i,jb)-0.5*hx_inc(jb); end;% % 电场Ez for i=1:N; for j=1:N; ez(i,j)=ga(i,j)*( dz(i,j)-iz(i,j) ); iz(i,j)=iz(i,j)+gb(i,j)*ez(i,j) ; end; end; % 边界电场置零,作为辅助P

10、ML for j=1:N; ez(1,j)=0; ez(N,j)=0; end; for i=1:N; ez(i,1)=0; ez(i,N)=0; end; % 金属边界条件 for i=N/2-halfgrid:N/2+halfgrid-1; for j=N/2-halfgrid:N/2+halfgrid-1; ez(i,j)=0; end; end; % 傅里叶变换Ez; s=0; for i=pa:pb+1; s=s+1; ez_out(1,s)=ez_out(1,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(i,qa); % 下 ez_out(2,s)=ez_out

11、(2,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(i,qb); % 上 end; s=0; for j=qa:qb+1; s=s+1; ez_out(3,s)=ez_out(3,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(pa,j); % 左 ez_out(4,s)=ez_out(4,s)+exp(-sqrt(-1)*2*pi*f*T*dt)*ez(pb,j); % 右 end;% %*磁场由电场更新*% 更新平面波磁场 for j=1:N-1; hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1); end;

12、 % 入射磁场傅里叶变换 inh=inh+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx_inc(ja-1);% % 磁场Hx for i=1:N; for j=1:N-1; curl_e=ez(i,j)-ez(i,j+1); ihx(i,j)=ihx(i,j)+fi1(i)*curl_e; hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j); end; end; % 磁场Hx加入射场 for i=ia:ib; hx(i,ja-1)=hx(i,ja-1)+0.5*ez_inc(ja); hx(i,jb)=hx(i,jb)

13、-0.5*ez_inc(jb); end; % 傅里叶变换Hx; s=0; % 注意：电场和磁场相差半个时间步 for i=pa:pb; s=s+1; hx_out(1,s)=hx_out(1,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx(i,qa); % 下 hx_out(2,s)=hx_out(2,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx(i,qb); % 上 end; s=0; for j=qa:qb; s=s+1; hx_out(3,s)=hx_out(3,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)

14、*dt)*hx(pa,j); % 左 hx_out(4,s)=hx_out(4,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hx(pb,j); % 右 end;% 磁场Hy for i=1:N-1; for j=1:N; curl_e=ez(i+1,j)-ez(i,j); ihy(i,j)=ihy(i,j)+fj1(j)*curl_e; hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j); end; end; % 磁场Hy加入射场 for j=ja:jb; hy(ia-1,j)=hy(ia-1,j)-0.5*ez_i

15、nc(j); hy(ib,j)=hy(ib,j)+0.5*ez_inc(j); end; % 傅里叶变换Hy; s=0; for i=pa:pb; s=s+1; hy_out(1,s)=hy_out(1,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hy(i,qa); % 下 hy_out(2,s)=hy_out(2,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hy(i,qb); % 上 end; s=0; for j=qa:qb; s=s+1; hy_out(3,s)=hy_out(3,s)+exp(-sqrt(-1)*2*pi*f*(T

16、+1/2)*dt)*hy(pa,j); % 左 hy_out(4,s)=hy_out(4,s)+exp(-sqrt(-1)*2*pi*f*(T+1/2)*dt)*hy(pb,j); % 右 end;% % 显示 plot(hx(40,:); axis(1 120 -0.5 0.5) mesh(ez) drawnow;end;%*近远场外推*% 电场平均for i=1:length ez_outave(:,i)=(ez_out(:,i)+ez_out(:,i+1)/2;end%外推回路上的电磁流分量current=zeros(4,length); % 电流z方向magnetic_x=zeros(

17、4,length); % 磁流x方向magnetic_y=zeros(4,length); % 磁流y方向% 电流zcurrent(1,:)=hx_out(1,:); % 下current(2,:)=-hx_out(2,:); % 上current(3,:)=-hy_out(3,:); % 左current(4,:)=hy_out(4,:); % 右% 磁流xmagnetic_x(1,:)=ez_outave(1,:); % 下magnetic_x(2,:)=-ez_outave(2,:); % 上magnetic_x(3,:)=0; % 左magnetic_x(4,:)=0; % 右% 磁流

18、ymagnetic_y(1,:)=0; % 下magnetic_y(2,:)=0; % 上magnetic_y(3,:)=-ez_outave(3,:); % 左magnetic_y(4,:)=ez_outave(4,:); % 右% 外推边界坐标X,Yposx=zeros(4,length);posy=zeros(4,length);posx(1,:)=linspace(pa*ddx,pb*ddx,length)+ddx/2; % 下posx(2,:)=linspace(pa*ddx,pb*ddx,length)+ddx/2; % 上posx(3,:)=pa*ddx; % 左posx(4,:

19、)=pb*ddx; % 右posy(1,:)=qa*ddx; % 下posy(2,:)=qb*ddx; % 上posy(3,:)=linspace(qa*ddx,qb*ddx,length)+ddx/2; % 左posy(4,:)=linspace(qa*ddx,qb*ddx,length)+ddx/2; % 右% 归一化fz,fmx,fmy求解fz=ddx*current/abs(ine); fmx=ddx*magnetic_x/abs(inh);fmy=ddx*magnetic_y/abs(inh);% RCS求解;co=0;for phi=pi/2:pi/180:pi+pi/2; co=

20、co+1; ss=0; for t=1:4; for s=1:length kr=posx(t,s)*cos(phi)+posy(t,s)*sin(phi); ss=ss+(-fz(t,s)-fmx(t,s)*sin(phi)+fmy(t,s)*cos(phi)*exp(sqrt(-1)*k*kr); % 外推公式 end end; rcs(co)=k/4*(abs(ss).2;end;rcs=10*log10(rcs/(lamda);save fdtdupml rcs;plot(rcs,'g');rcs_upml=rcs;save rcs_upml.mat rcs_upml;

21、toc; % 该程序用于计算球的RCS，包含了精确解以及物理光学法近似解clear alleps = 0.00001;index = 0;% kr 0.05 - 15 => 300 点for kr = 0.05:0.05:15 index = index + 1; sphere_rcs = 0. + 0.*i; f1 = 0. + 1.*i; f2 = 1. + 0.*i; m = 1.; n = 0.; q = -1.; % 初始化叠加项 del =+*i; while(abs(del) > eps) q = -q; n = n + 1; m = m + 2; del = (2.

22、*n-1) * f2 / kr-f1; f1 = f2; f2 = del; del = q * m /(f2 * (kr * f1 - n * f2); sphere_rcs = sphere_rcs + del; end rcs1(index) = (abs(1-1/(2*i*kr)2; rcs(index) = abs(sphere_rcs); sphere_rcsdb1(index)= 10. * log10(rcs1(index); sphere_rcsdb(index) = 10. * log10(rcs(index);end % 画图程序figure(1);n=0.05:.05:

23、15;plot (n,rcs1,'-r');hold on;plot (n,rcs,'k');set (gca,'xtick',1 2 3 4 5 6 7 8 9 10 11 12 13 14 15);xlabel ('Sphere circumference in wavelengths');ylabel ('Normalized sphere RCS');grid;figure (2);plot (n,sphere_rcsdb1,'-r');hold on;plot (n,sphere_rcsd

24、b,'k');set (gca,'xtick',1 2 3 4 5 6 7 8 9 10 11 12 13 14 15);xlabel ('Sphere circumference in wavelengths');ylabel ('Normalized sphere RCS - dB');grid;figure (3);semilogx (n,sphere_rcsdb1,'-r');hold on;semilogx (n,sphere_rcsdb,'k');xlabel ('Sphere

25、circumference in wavelengths');ylabel ('Normalized sphere RCS - dB');grid;function rcs = rcs_cylinder(r1, r2, h, freq, phi, CylinderType)% rcs_cylinder.m% 本程序计算椭圆柱和圆柱RCS% 采用的几何光学和物理光学法两种公式 r = r1; eps =0.00001;dtr = pi/180;phir = phi*dtr;freqGH = num2str(freq*1.e-9);lambda = 3.0e+8 /freq

26、; % wavelengthk=2*pi/lambda;% CylinderType= 'Elliptic' % 'Elliptic' or 'Circular' switch CylinderTypecase 'Circular' % 计算 RCS 0 - (90-.5) 90 °附近由于算法问题不连续，所以不计算 index = 0; for theta = 0.0:.1:90-.5 index = index +1; thetar = theta * dtr; rcs(index) = (lambda * r *

27、sin(thetar) / . (8. * pi * (cos(thetar)2) + eps; rcs1(index) = r*sin(thetar)*sin(k*h*cos(thetar)/(k*cos(thetar)2); rcs2(index) = 2*pi*h2*r/lambda; end % 单独计算90°RCS thetar = pi/2; index = index +1; rcs(index) = (2. * pi * h2 * r / lambda )+ eps; rcs1(index) = (2. * pi * h2 * r / lambda )+ eps; %

28、 90.5180 for theta = 90+.5:.1:180. index = index + 1; thetar = theta * dtr; rcs(index) = ( lambda * r * sin(thetar) / . (8. * pi * (cos(thetar)2) + eps; rcs1(index) = r*sin(thetar)*sin(k*h*cos(thetar)/(k*cos(thetar)2); rcs2(index) = 2*pi*h2*r/lambda; endcase 'Elliptic' r12 = r1*r1; r22 = r2*r2; h2 = h*h; % 计算 RCS 0 - (90-.5) 90 °附近由于算法问题不连续，所以不计算 index = 0; for theta = 0.0:.1:90-.5 i