[资料共享]一维等离子体FDTD的Matlab源代码(两种方法）

[post]一维等离子体FDTD的Matlab源代码 j=up7395  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% gb^<6BYUG  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v;WfcpWq2  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /S #Z.T~~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% r,r"?}Z  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Bo4MoSF}  
%%%%%%%%%初始化 0^25uAD=  
clear; f;`7}7C  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1C5~GI`  
%%%%%%%%%%%%%%%%系统参数 a !yBEpMo  
.3 S9=d?  
TimeT=200;%迭代次数 ~p!=w#/  
KE=2000;%网格树木 ?^by3\,VZ  
kc=450;%源的位置 N0V`xrS  
kpstart=500;%等离子体开始位置 1)BIh~1{p  
kpstop=1000;%等离子体终止位置 N
?{.}-Q  
DispE=zeros(1,TimeT); p[8H!=`K  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PiVp(; rtQ  
%%%%%%%%%%%%%物理参数 MwHxn%  
c0=3e8;%真空中波速 1G}\IK1+  
f=150e6; z=j,-d%9  
f1=75e6; W%-`  
lamda=c0/f; O4g2s8k  
WL=50; ;uhpo  
OMIGA=pi/WL; d1>L&3HKx  
zdelta=lamda/WL;%网格大小 3P|z`}Ka  
dt=zdelta/(2*c0);%时间间隔 ?X'l&k>  
u$Wv*;TT%  
r}4  
u0=1e7;%碰撞频率 |I2~@RfpO:  
fpe=1e7;%等离子体频率 ,{jF)NQaP  
wpe=2*pi*fpe;%等离子体圆频率 rUW/d3y  
epsz=1/(4*pi*9*10^9); % 真空介电常数 +UX~TT:  
mu=1/(c0^2*epsz);%磁常数 GTocN1,Z~a  
ex_low_m1=0; PN"=P2e/ 6  
ex_low_m2=0; ,GYK3+}Z  
ex_high_m1=0; [!S%nYs&8L  
ex_high_m2=0; ($X2SIZh  
m
:W+s4!E  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% r]B`\XWz  
%%%%%%%%%%%%%初始化电磁场 G@4n]c_  
Ex=zeros(1,KE); XE`u  
Ex_Pre=zeros(1,KE); WX&0;Kr  
Hy=zeros(1,KE); ,ieew`  
Hy_Pre=zeros(1,KE); =
KW~k7TaN  
Dx=zeros(1,KE); iI$;%uY3g  
Dx_Pre=zeros(1,KE); l,k.Jo5  
Sx1=zeros(1,KE); @M)"  
Sx2=zeros(1,KE); 6#;u6@+}yy  
Sx3=zeros(1,KE); Blox~=cW  
Sx=zeros(1,KE); WxJV zHtR  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% u|sdQ  
%%%%%%%%%%%%%%%%%%开始计算 }Ryrd!3bY  
for T=1:TimeT 9zNMv-  
    %%%保存前一时间的电磁场 FfM,~s<Efz  
    Ex_Pre=Ex; ``,q[|  
    Hy_Pre=Hy; e% #?B *  
    %%%%中间差分计算Dx  r<1.'F  
    for i=2:KE I~&*8)x
M  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); C,) e7  
    end e8U6D+jY  
    %%%%%%%%加入源 KN:V:8:J  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); o9+fAH`D  
       [N_)V kpr  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); /_V4gwb}|-  
    if T<WL &9khIJIn  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 *XkgwJq  
    end :+v4,=fHy  
    %%%计算电场Ex <y<   
    for i=1:kpstart-1 R,8460e7  
        Ex(i)=Dx(i)/epsz; h/NI5  
    end 3Lm7{s?=Z-  
    for i=kpstop+1:KE 3RP\w~?  
        Ex(i)=Dx(i)/epsz; 0I}c|V'P  
    end ;6q`c!p7  
    if T>50 Id?2(Tg  
        gw=0; LrhQG  
    end /&o<kY  
    i=kpstart; '"9Wt@ .  
         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); 6Zmzo,{  
         Ex(i)=Dx(i)/epsz-Sx1(i); O6b.oS'-  
         Ex(i)=Ex(i); Ih%LKFT  
     Y.#:l<  
    for i=kpstart+1:kpstop 4v#A#5+O E  
         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); }UWi[UgA  
         Ex(i)=Dx(i)/epsz-Sx1(i); L[. )!c8k  
         %Ex(i)=Ex(i); n7,LfO#  
    end qFwAzW
;"  
    Sx3=Sx2; ?1]h5Uh[b  
    Sx2=Sx1;  Wo,fHY  
    Ex(1)=ex_low_m2; 13MB1n  
    ex_low_m2=ex_low_m1; n=%D}W  
    ex_low_m1=Ex(2); qj*77  
$
sb `BS  
    Ex(KE)=ex_high_m2; iz:O]kI  
    ex_high_m2=ex_high_m1; kp8kp
`S7  
    ex_high_m1=Ex(KE-1); >]/aG!  
    %%%%%%%%%%%%%%%%%%计算磁场 xX\A&9m  
    for i=1:KE-1 < .. N3&n"w _d  
r~;TId} #  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b}5hqIy  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% H2D j`0  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 7EukrE<b'  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% xN]88L}Tn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 'X1fb:8m8  
%%%%%%%%%初始化 zFGZ;?i  
clear; eXsFPM  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% h]#bPb  
%%%%%%%%%%%%%%%%系统参数 %.n 7+  
TimeT=3000;%迭代次数 09C[B+>h  
KE=2000;%网格树木 :Y>M//0  
kc=450;%源的位置 <aEY=IF4  
kpstart=500;%等离子体开始位置 nLv"ON~  
kpstop=1000;%等离子体终止位置 Z|dng6ck  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% _9Y7.5  
%%%%%%%%%%%%%物理参数 SPU_@ Pk  
c0=3e8;%真空中波速 c!}f\ ]D  
zdelta=1e-9;%网格大小 0G"I}Jp{  
dt=zdelta/(2*c0);%时间间隔 X4}Lg2ts  
f=900e12;%Gause脉冲的载频 Y]!{ nW  
d=3e-15%脉冲底座宽度 c=A)_ZFg  
t0=2.25/f;%脉冲中心时间 K/+w6d  
u0=57e12%碰撞频率 /?Fa<{  
fpe=2000e12;%等离子体频率 =_Y#uE$  
wpe=2*pi*fpe;%等离子体圆频率 @Zd/>'  
epsz=1/(4*pi*9*10^9); % 真空介电常数 }Qo:;&"3  
mu=1/(c0^2*epsz);%磁常数 KgMW  
ex_low_m1=0; ~}F$1;t0  
ex_low_m2=0; vTF_`X  
ex_high_m1=0; [MVG\6Up(  
ex_high_m2=0; PgF7ug%,@C  
a0=2*u0/dt+(2/dt)^2; Uq}-<q  
a1=-8/(dt)^2; 6Y(Vs>  
a2=-2*u0/dt+(2/dt)^2; ^9PB+mz  
b0=wpe^2+2*u0/dt+(2/dt)^2; |UZhMF4/-L  
b1=2*wpe^2-8/(dt)^2; Kv26rY8Q  
b2=wpe^2-2*u0/dt+(2/dt)^2; nkvkHh  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% e.:SBXZ  
%%%%%%%%%%%%%初始化电磁场 <xWBS/K  
Ex=zeros(1,KE); @fwk  
Ex_Pre=zeros(1,KE); !O~5<tA[#1  
Hy=zeros(1,KE); }n91aE3v  
Hy_Pre=zeros(1,KE); +r 2\v  
Dx=zeros(1,KE); k4ijWo{:0  
Dx_Pre=zeros(1,KE); /\&Wk;u3  
Sx1=zeros(1,KE); GMO|A.bzzN  
Sx2=zeros(1,KE); =2yg:D  
Sx3=zeros(1,KE); _N-JRM m<  
Sx=zeros(1,KE); #`9D,+2iB%  
Dx=Ex; 9e:}qO5)  
Dx1=Ex; (e$/@3*  
Dx2=Ex; 4 9zOhG |  
Dx3=Ex; dJ0qg_ U&  
Ex1=Ex; yAt,XG3  
Ex2=Ex; \.7O0Q{  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0Zv<]xO  
%%%%%%%%%%%%%%%%%%开始计算 ;\5^yDv[e  
for T=1:TimeT NFQR  
    %%%保存前一时间的电磁场 Lp?JSMe  
    Ex_Pre=Ex; ,wvzY7%  
    Hy_Pre=Hy; v<qiu>sbz}  
    %%%%中间差分计算Dx )FfJ%oT}  
    for i=2:KE  EL[N%M3  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); .g6DKjy>  
    end p&%M=SzN  
    %%%%%%%%加入源 ix=H=U]Q{  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); :U7m@3czU  
     'F665  
    %%%计算电场Ex {} 11U0  
    for i=1:kpstart-1 )~X*&(7RR}  
        Ex(i)=Dx(i)/epsz; =_/,C  
    end 
`xpU  
    for i=kpstop+1:KE 5c~OG6COx  
        Ex(i)=Dx(i)/epsz; TdAHw @(  
    end C/!c?$J  
    Dx3=Dx2; !Y3 *\  
Dx2=Dx1; ~!w()v n  
Dx1=Dx; %)K)h&m  
    for i=kpstart:kpstop -gP4| r8&  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); ]EzX$T  
    end 75pz' Cb  
    Ex2=Ex1; Q*+_%n1 /  
Ex1=Ex; (y xrK  
    Sx3=Sx2; ,^_aqH  
    Sx2=Sx1; A%G \ AT  
    Ex(1)=ex_low_m2; A S#D9o  
    ex_low_m2=ex_low_m1; {,|
J?>{  
    ex_low_m1=Ex(2); *
qLOr6  
kM`!'0kt  
    Ex(KE)=ex_high_m2; !y>MchNv  
    ex_high_m2=ex_high_m1; S}rW=hO  
    ex_high_m1=Ex(KE-1); |# z
znT"  
    %%%%%%%%%%%%%%%%%%计算磁场 !PfIe94{`  
    for i=1:KE-1 ktr l|  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); I=,u7w`m  
    end ,DT=(  
    plot(Ex); EQ j2:9f  
    grid on; v&>TU(x\H  
    pause(0.01); ]qJ6#sAw75  
end

我也发到         求助【48小时速问速答区】 » 色散性的FDTD仿真求助  这个里面了，重复发，不算违规吧。

原来是赵博士啊 刚访问了你的科学网博客啊 以后要多多请教

欢迎到这里发帖，好东西是要奖励的

3楼是？ #t:]a<3Y2  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ LO]D XW 9  
最好有相关的物理公式，会更好理解点 Nv "R'Pps  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 "Pi\I9M3  
现在天天写code，发现原来写的真是一个烂字啊。 u\=gps/Z  
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); jC+>^=J(  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 ^_6%dKLK  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ GbUcNROr  
主要是 k/AcXU%O+  
     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); AbUU#C7  
         Ex(i)=Dx(i)/epsz-Sx1(i); _Ptf^+  
         Ex(i)=Ex(i); EA7]o.Nm*{  
和 qyl~*r*  
    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); fC xN!  
         Ex(i)=Dx(i)/epsz-Sx1(i); KIn^,d0H  
V\ch0i 1  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

还有就是步长的设置为什么是dt=dx/2c而不是dt=dx/c呢。但dt=dx/c在这里好像是不稳定的