RCS-FDTD方柱-球的RCS程序

1、精選優(yōu)質(zhì)文檔-傾情為你奉上% FDTD_2D_RCSclear;clc;%*初始化*%tic;% v=3*108; % 波速f=15*10(9); % 頻率lamda=v/f; % 波長k=2*pi/lamda; % 波數(shù)epsz=1/(4*pi*9*109); % 真空介電常數(shù)mu=4*pi*10.(-7); % 真空磁導(dǎo)率Z=sqrt(mu/epsz); % 真空波阻抗epsilon=1; % 相對介電常數(shù)sigma=0.0; % 電導(dǎo)率%N=120; % 網(wǎng)格數(shù)量L=1000; % 迭代次數(shù)ddx=lamda/20; % 網(wǎng)格尺寸dt=ddx/(2*v); % 時(shí)間間隔lstart=1;

2、 % 空間起始行坐標(biāo) lend=N; % 空間終止行坐標(biāo)rstart=1; % 空間起始列坐標(biāo)rend=N; % 空間終止列坐標(biāo)ia=N/4; % 總場區(qū)域x左ib=3*N/4; % 總場區(qū)域x右ja=ia; % 總場區(qū)域x下jb=ib; % 總場區(qū)域x上pa=ia-5; % 外推區(qū)域x左pb=ib+5; % 外推區(qū)域x右qa=pa; % 外推區(qū)域x下qb=pb; % 外推區(qū)域x上length=pb-pa+1; % 每個(gè)邊上的總長度npml=N/8; % PML點(diǎn)數(shù)% 方柱參數(shù)side=2*lamda; % 方柱邊長halfgrid=round(side/(2*ddx); % 方柱占據(jù)網(wǎng)格大小

3、% 媒質(zhì)參數(shù)for i=lstart:lend; % 控制媒質(zhì)分部區(qū)域 for j=rstart:rend; ga(i,j)=1/(epsilon+sigma*dt/epsz); % 求和參量 gb(i,j)=sigma*dt/epsz; % 求和參量 end;end;% 源參數(shù)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設(shè)置*% 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;% 阻抗?jié)u變設(shè)置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); % 入射電場傅里葉變換(必須用邊界加入) % 注意：當(dāng)計(jì)算頻譜用的是其他的位置的入射波時(shí)候，會(huì)產(chǎn)生很大的誤差； ine=ine+exp(-sqrt(-1)*2*pi*f*T*dt)*ez_inc(ja-1); % 電荷密度Dz for i=2:N-1; % 為了使每個(gè)電場周圍都有磁場進(jìn)行數(shù)組下標(biāo)處理 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; % 注意：電場和磁場相差半個(gè)時(shí)間步 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;%*近遠(yuǎn)場外推*% 電場平均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,:); % 右% 外推邊界坐標(biāo)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; % 該程序用于計(jì)算球的RCS，包含了精確解以及物理光學(xué)法近似解clear alleps = 0.00001;index = 0;% kr 0.05 - 15 => 300 點(diǎn)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.; % 初始化疊加項(xiàng) 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% 本程序計(jì)算橢圓柱和圓柱RCS% 采用的幾何光學(xué)和物理光學(xué)法兩種公式 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' % 計(jì)算 RCS 0 - (90-.5) 90 °附近由于算法問題不連續(xù)，所以不計(jì)算 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 % 單獨(dú)計(jì)算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; % 計(jì)算 RCS 0 - (90-.5) 90 °附近由于算法問題不連續(xù)，所以不計(jì)算 index = 0; for theta = 0.0:.1:90-.5 i