一维等离子体fdtd的matlab源代码(两种方法）.txt

一维等离子的FDTD的MATLAB程序

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%初始化
clear;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%系统参数
TimeT=3000;%迭代次数
KE=2000;%网格数目
kc=450;%源的位置
kpstart=500;%等离子体开始位置
kpstop=1000;%等离子体终止位置
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%物理参数
c0=3e8;%真空中波速
zdelta=1e-9;%网格大小
dt=zdelta/(2*c0);%时间间隔
f=900e12;%Gause脉冲的载频
d=3e-15%脉冲底座宽度
t0=2.25/f;%脉冲中心时间
u0=57e12%碰撞频率
fpe=2000e12;%等离子体频率
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;
a0=2*u0/dt+(2/dt)^2;
a1=-8/(dt)^2;
a2=-2*u0/dt+(2/dt)^2;
b0=wpe^2+2*u0/dt+(2/dt)^2;
b1=2*wpe^2-8/(dt)^2;
b2=wpe^2-2*u0/dt+(2/dt)^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%初始化电磁场
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);
Dx=Ex;
Dx1=Ex;
Dx2=Ex;
Dx3=Ex;
Ex1=Ex;
Ex2=Ex;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%开始计算
for T=1:TimeT
    %%%保存前一时间的电磁场
    Ex_Pre=Ex;
    Hy_Pre=Hy;
    %%%%中间差分计算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);
    
    %%%计算电场Ex
    for i=1:kpstart-1
        Ex(i)=Dx(i)/epsz;
    end
    for i=kpstop+1:KE
        Ex(i)=Dx(i)/epsz;
    end
    Dx3=Dx2;
Dx2=Dx1;
Dx1=Dx;
    for i=kpstart:kpstop
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i));
    end
    Ex2=Ex1;
Ex1=Ex;
    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
    plot(Ex);
    grid on;
    pause(0.01);
end