fdtd.asv

一维中实现等离子体FDTD编程

clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数

TimeT=2000;%迭代次数
KE=200;%网格树木
kc=450;%源的位置
kpstart=100;%等离子体开始位置
kpstop=200;%等离子体终止位置
DispE=zeros(1,TimeT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
f=150e6;
f1=75e6;
lamda=c0/f;
WL=50;
OMIGA=pi/WL;
zdelta=lamda/WL;%网格大小
dt=zdelta/(2*c0);%时间间隔
t0=50.0;
spread=20;



u0=1e7;%碰撞频率
fpe=1e7;%等离子体频率
wpe=2*pi*fpe;%等离子体圆频率
epsz=1/(4*pi*9*10^9); % 真空介电常数
mu=1/(c0^2*epsz);%磁常数
ex_low_m1=0;
ex_low_m2=0;
ex_high_m1=0;
ex_high_m2=0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
Ex=zeros(1,KE);
Ex_Pre=zeros(1,KE);
Hy=zeros(1,KE);
Hy_Pre=zeros(1,KE);
Dx=zeros(1,KE);
Dx_Pre=zeros(1,KE);
Sx1=zeros(1,KE);
Sx2=zeros(1,KE);
Sx3=zeros(1,KE);
Sx=zeros(1,KE);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT

    %%%%中间差分计算Dx
    for i=2:KE
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));
    end
    %%%%%%%%加入源
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2);
      
    Dx(3)=sin(OMIGA*dt*T)*exp(-0.5*((t0-T)^2/spread));

    %%%计算电场Ex
    for i=1:KE
        Ex(i)=Dx(i)/epsz;
    end
    i=kpstart;
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(2*u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         Ex(i)=Ex(i);
    
    for i=kpstart+1:kpstop
         Sx1(i)=(1+exp(-1*u0*dt))*Sx2(i)-exp(-1*u0*dt)*Sx3(i)+(wpe^2*dt/(u0))*(1-exp(-1*u0*dt))*Ex(i);
         Ex(i)=Dx(i)/epsz-Sx1(i);
         %Ex(i)=Ex(i);
    end
    Sx3=Sx2;
    Sx2=Sx1;
    Ex(1)=ex_low_m2;
    ex_low_m2=ex_low_m1;
    ex_low_m1=Ex(2);

    Ex(KE)=ex_high_m2;
    ex_high_m2=ex_high_m1;
    ex_high_m1=Ex(KE-1);
    %%%%%%%%%%%%%%%%%%计算磁场
    for i=1:KE-1
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));
    end
    %if round(T/100)*100==T
      plot(Ex);
      axis([0 200 -1 1]);
      grid on;
      pause(0.00001);
      T
      % end
    
end