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

[post]一维等离子体FDTD的Matlab源代码 77F
I&*q  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% F;4*,Ap  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% o$*aAgS+  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6zo'w Wc3  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% q#{.8H-X'  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% S>-x<'Os  
%%%%%%%%%初始化 8T7ex(w  
clear; DqfWu*  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% {:D8@jb[  
%%%%%%%%%%%%%%%%系统参数 `&\jOve   
S(B$[
)(  
TimeT=200;%迭代次数 4pvT?s>68  
KE=2000;%网格树木 w\"~*(M  
kc=450;%源的位置 S}O>@%  
kpstart=500;%等离子体开始位置 tx,_0[hZi  
kpstop=1000;%等离子体终止位置 BHVC&F*>  
DispE=zeros(1,TimeT); Lro[ |A  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% m*Zq3j  
%%%%%%%%%%%%%物理参数 #.<Uy."z2  
c0=3e8;%真空中波速 (=2-*((&(A  
f=150e6; #ujry.m  
f1=75e6; Mnu8d:$  
lamda=c0/f; :%6OFO$z  
WL=50; N~""Lc&  
OMIGA=pi/WL; WH>=*\  
zdelta=lamda/WL;%网格大小 #%4XZ3j#j;  
dt=zdelta/(2*c0);%时间间隔 m\4V;F  
TY."?` [FK  
QKW\z aG  
u0=1e7;%碰撞频率 jGg,)~)Y  
fpe=1e7;%等离子体频率
F9ys.Bc  
wpe=2*pi*fpe;%等离子体圆频率 -0Q
oVGw  
epsz=1/(4*pi*9*10^9); % 真空介电常数 P3Wnso  
mu=1/(c0^2*epsz);%磁常数 y*F !k{P  
ex_low_m1=0; 0}'xoYv f  
ex_low_m2=0; 04K[U9W3  
ex_high_m1=0; @q98ac*{  
ex_high_m2=0; v qt#JdPp9  
<r9L-4  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Wq5}SM  
%%%%%%%%%%%%%初始化电磁场 S:bY
eD4  
Ex=zeros(1,KE); zBQV2.@  
Ex_Pre=zeros(1,KE); jR1o<]
?  
Hy=zeros(1,KE); 'YKzs;y$  
Hy_Pre=zeros(1,KE); #x?Ku\ts  
Dx=zeros(1,KE); &d%\&fCm(  
Dx_Pre=zeros(1,KE); Oe)d|6=  
Sx1=zeros(1,KE); oZD+AF$R  
Sx2=zeros(1,KE); PBeBI:  
Sx3=zeros(1,KE); l}bAwJ?  
Sx=zeros(1,KE); B(p
xyv)  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% O'm5k l  
%%%%%%%%%%%%%%%%%%开始计算 #r=Jc8J_  
for T=1:TimeT $U_M|Xa  
    %%%保存前一时间的电磁场 2 c 2lK  
    Ex_Pre=Ex; i;flK*HOZ9  
    Hy_Pre=Hy; </hv{<  
    %%%%中间差分计算Dx Ulf'gD4e  
    for i=2:KE t5| }0ID-  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); r}**^"mFy  
    end Wc*jTip  
    %%%%%%%%加入源 m4 k:uk7N  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); y*G3dWb  
       Fb!Ew`;QT  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); b<\2j5  
    if T<WL 5NR@<FE  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 diz=|g=w  
    end o>6c?Xi&  
    %%%计算电场Ex *:5S*E&}V  
    for i=1:kpstart-1 -JEiwi,  
        Ex(i)=Dx(i)/epsz; _Gu;=H,~&  
    end NJNS8\4  
    for i=kpstop+1:KE ~WjK'N4n5  
        Ex(i)=Dx(i)/epsz; 3s]aXz:  
    end @%'1Jd7-Wp  
    if T>50 bu?4$O  
        gw=0; U<YcUmX  
    end %.h&W;  
    i=kpstart; wD\viuq0  
         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); fT2F$U  
         Ex(i)=Dx(i)/epsz-Sx1(i); &<=?O a  
         Ex(i)=Ex(i); F@tfbDO?  
     kqH:H~sgD  
    for i=kpstart+1:kpstop jIL+^{
K<  
         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); CN{xh=2qY[  
         Ex(i)=Dx(i)/epsz-Sx1(i); khc1<BBsT  
         %Ex(i)=Ex(i); %eE0a4^".  
    end W?F Q  
    Sx3=Sx2; .E:3I!dH7  
    Sx2=Sx1; e=cb%  
    Ex(1)=ex_low_m2; *8bj3A]vf  
    ex_low_m2=ex_low_m1; t/wo G9N  
    ex_low_m1=Ex(2); VLfc6:Yg  
hCOCX_  
    Ex(KE)=ex_high_m2; C09rgEB\B  
    ex_high_m2=ex_high_m1; E0*KKo%  
    ex_high_m1=Ex(KE-1); kl~)<,/@  
    %%%%%%%%%%%%%%%%%%计算磁场 [n2+`A  
    for i=1:KE-1 < .. @*9c2\"k  
S4_C8  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% D4N(FZ
0~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !' 0PM[  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ckb(+*+l  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% N^M6*,F,J  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 8Vjv #pm  
%%%%%%%%%初始化 )MF 4b][  
clear; ~Zn|(  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% njZJp|y6  
%%%%%%%%%%%%%%%%系统参数 >,QCKZH  
TimeT=3000;%迭代次数 lCgzQZ  
KE=2000;%网格树木 W'~s  
kc=450;%源的位置 BIS.,  
kpstart=500;%等离子体开始位置 $NCR V:J  
kpstop=1000;%等离子体终止位置 (< >Lfn  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% VJP
#  
%%%%%%%%%%%%%物理参数 O')=]6CQ*  
c0=3e8;%真空中波速 2rW9ja  
zdelta=1e-9;%网格大小 2+/r~LwbK  
dt=zdelta/(2*c0);%时间间隔 #].qjOj  
f=900e12;%Gause脉冲的载频 m] yUcj{F  
d=3e-15%脉冲底座宽度 \*&?o51!e  
t0=2.25/f;%脉冲中心时间 Eg&:yF}?(  
u0=57e12%碰撞频率 }d>X
h8:%)  
fpe=2000e12;%等离子体频率 *fs[]q'Q  
wpe=2*pi*fpe;%等离子体圆频率 $@vB<(sk  
epsz=1/(4*pi*9*10^9); % 真空介电常数 ^We}i  
mu=1/(c0^2*epsz);%磁常数 gnkeJ}K  
ex_low_m1=0; tDavp:M1v  
ex_low_m2=0; | TG6-e_  
ex_high_m1=0; ~;z] _`_Va  
ex_high_m2=0; %#o@c  
a0=2*u0/dt+(2/dt)^2; NJSbS<O  
a1=-8/(dt)^2; ]imVIu  
a2=-2*u0/dt+(2/dt)^2; "UreV  
b0=wpe^2+2*u0/dt+(2/dt)^2; &uF~t |!c  
b1=2*wpe^2-8/(dt)^2; ;B;@MD,B  
b2=wpe^2-2*u0/dt+(2/dt)^2; 3d]~e  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +V(^"Z~  
%%%%%%%%%%%%%初始化电磁场 0y t36Du  
Ex=zeros(1,KE); sv&^sARN  
Ex_Pre=zeros(1,KE); Wx8:GBM$2  
Hy=zeros(1,KE); =1k%T{>  
Hy_Pre=zeros(1,KE); E=CAWj\  
Dx=zeros(1,KE); a~J!G
:(  
Dx_Pre=zeros(1,KE); Td|u-9OM  
Sx1=zeros(1,KE); rT=C/SKP  
Sx2=zeros(1,KE); *~YdL7f)J  
Sx3=zeros(1,KE); Caz5q|Oo  
Sx=zeros(1,KE); ~]#-S20  
Dx=Ex; pY[b[ezb  
Dx1=Ex; Uf|uFGb  
Dx2=Ex; $P#x>#+[A  
Dx3=Ex; eEfGH  
Ex1=Ex; eXD~L&s[  
Ex2=Ex; 9XYm8g'X  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% v%c/eAF  
%%%%%%%%%%%%%%%%%%开始计算 V']Z_$_  
for T=1:TimeT *>"NUHq  
    %%%保存前一时间的电磁场 .zl[nx[9"D  
    Ex_Pre=Ex; R?;mu^B  
    Hy_Pre=Hy; _VRxI4q  
    %%%%中间差分计算Dx $)$r  
    for i=2:KE B'Jf&v  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); {* :^K\-  
    end &ds+9A  
    %%%%%%%%加入源 xJAQ'ANr  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); O$;#GpR  
     10h;N[  
    %%%计算电场Ex Rnoz[1y?0  
    for i=1:kpstart-1 M>|ZBEK  
        Ex(i)=Dx(i)/epsz; $U.|  
    end =+UtAf<n  
    for i=kpstop+1:KE G3`9'-2q@c  
        Ex(i)=Dx(i)/epsz; )/vom6y*   
    end ,*dLE   
    Dx3=Dx2; G[JWG  
Dx2=Dx1; k\qFWFR  
Dx1=Dx; 9l"=]7~%  
    for i=kpstart:kpstop 5jso)`IL  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); _WvVF*Q"k  
    end J}[[tl  
    Ex2=Ex1; pW@W-k:u  
Ex1=Ex; /V{UTMSz  
    Sx3=Sx2; Ww60-d}}Q  
    Sx2=Sx1; SD:D8"8  
    Ex(1)=ex_low_m2;
:.-z!  
    ex_low_m2=ex_low_m1; i-w$-2w  
    ex_low_m1=Ex(2); ,`RX~ H=C  
vb?.`B_>&  
    Ex(KE)=ex_high_m2; }:{9!RMO  
    ex_high_m2=ex_high_m1; j{r@>g;3  
    ex_high_m1=Ex(KE-1); ']>9/r#  
    %%%%%%%%%%%%%%%%%%计算磁场 NA/`LaJ  
    for i=1:KE-1 dt@c,McN|Q  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); :hs~;vn)  
    end Hi=</ Wy;  
    plot(Ex); 343d`FRa}  
    grid on; OJ"./*H  
    pause(0.01); .D ^~!A  
end

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

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

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

3楼是？ N,B!D~@  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ n6t@ e^  
最好有相关的物理公式，会更好理解点 }@t"B9D  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 :[
_msd  
现在天天写code，发现原来写的真是一个烂字啊。 dV'EiNpf  
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); }}tbOD)t  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 S$Zi{bU`G  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ XnQR(r)pR2  
主要是 g F*AS(9  
     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); hGz_F/  
         Ex(i)=Dx(i)/epsz-Sx1(i); Kp`{-dUf  
         Ex(i)=Ex(i); "r-P[EKpL  
和 r0XEB,}  
    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); Mqr_w!8d  
         Ex(i)=Dx(i)/epsz-Sx1(i); Cd}^&z  
?q:|vt  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

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