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

[post]一维等离子体FDTD的Matlab源代码 mkrVeBp  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ?7J::}R  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9A/bA|
$  
%%%%%%%%%%%%%%%%%%                1D                   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% *z.rOY= 8  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% /)`]p1c1%w  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 9,JWi{lIv  
%%%%%%%%%初始化 FZIC|uz  
clear; a6O <t;&  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [HL>Lp&A?  
%%%%%%%%%%%%%%%%系统参数 b*/Mco 9O  
'+3C2!  
TimeT=200;%迭代次数 p#_5w  
KE=2000;%网格树木 UtQCTNjC{  
kc=450;%源的位置 ).A9>^6?{  
kpstart=500;%等离子体开始位置 )dh`aQ%N "  
kpstop=1000;%等离子体终止位置 e m0 hTxb  
DispE=zeros(1,TimeT); drk BW}_  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% )lz~Rt;1i  
%%%%%%%%%%%%%物理参数 %wI)uJ2  
c0=3e8;%真空中波速 c((bUjS'=Y  
f=150e6; ~qIr'?D  
f1=75e6; U<E]c 4*  
lamda=c0/f; yH>C7M7t  
WL=50; 67 ~pn  
OMIGA=pi/WL; 3s67)n  
zdelta=lamda/WL;%网格大小 Z
`U+a  
dt=zdelta/(2*c0);%时间间隔 FCWk8/  
L u'<4 R  
=;E0PB_w  
u0=1e7;%碰撞频率 !IA\c(c^  
fpe=1e7;%等离子体频率 oZ^,*  
wpe=2*pi*fpe;%等离子体圆频率 <q>d@Foi  
epsz=1/(4*pi*9*10^9); % 真空介电常数 5&O%0`t  
mu=1/(c0^2*epsz);%磁常数 s|`wi}"x  
ex_low_m1=0; "a3?m)  
ex_low_m2=0; 6 jm@`pYbE  
ex_high_m1=0; yz5! >|EB  
ex_high_m2=0; pOh<I{r1  
F^aD#  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% &`m$Zzl;  
%%%%%%%%%%%%%初始化电磁场 +r9neS.l  
Ex=zeros(1,KE); WW>m`RU`  
Ex_Pre=zeros(1,KE); u[oV Jvc  
Hy=zeros(1,KE); 9NNXj^7
  
Hy_Pre=zeros(1,KE); |lZp5MOc  
Dx=zeros(1,KE); w=a$]`  
Dx_Pre=zeros(1,KE); ;NrPMz  
Sx1=zeros(1,KE); P}qpy\/(4  
Sx2=zeros(1,KE); Px9 K  
Sx3=zeros(1,KE); ;(A-  
Sx=zeros(1,KE); y)a)VvU":  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FOsxId[f9  
%%%%%%%%%%%%%%%%%%开始计算 Pai8r%Zfu  
for T=1:TimeT &%;n9K  
    %%%保存前一时间的电磁场 #S
x  
    Ex_Pre=Ex; rHk,OC  
    Hy_Pre=Hy; C"%B>e  
    %%%%中间差分计算Dx >/lB%<$/  
    for i=2:KE u6Wan*I?  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); `@v;QLD"d<  
    end e+D]9wM8  
    %%%%%%%%加入源 %u_dxpx  
    %Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); }`%ks  
       %>6ilGQ+  
    Dx(kc)=sin(OMIGA*T)+sin(0.5*OMIGA*T); 9%"`9j~H>  
    if T<WL |0?v4%g  
         Dx(kc)=Dx(kc)*0.5*(1-cos(OMIGA*T));%开关函数，升余弦函数 CC;^J-h/  
    end 3HW&\:q5'M  
    %%%计算电场Ex SM2
N3"\  
    for i=1:kpstart-1 Bq1}"092  
        Ex(i)=Dx(i)/epsz; |lg jI!iK  
    end <;O^3_'  
    for i=kpstop+1:KE ?
Zsh\^k.g  
        Ex(i)=Dx(i)/epsz; DdUw~n,  
    end *Ms"{+C  
    if T>50 ICr.Gwe3_  
        gw=0; 6}!1a?X  
    end zSU,le  
    i=kpstart; k_ywwkG9lU  
         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); Cb`,N  
         Ex(i)=Dx(i)/epsz-Sx1(i); c))?9H ,e)  
         Ex(i)=Ex(i); \nPf\6;M  
     :FfEjNil  
    for i=kpstart+1:kpstop &[_@f
#  
         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); CV_M |  
         Ex(i)=Dx(i)/epsz-Sx1(i); odC"#Rb  
         %Ex(i)=Ex(i); 2|d^#8)ZC  
    end &q9=0So4\  
    Sx3=Sx2; S'kgpF"bm  
    Sx2=Sx1; }f14# y;  
    Ex(1)=ex_low_m2; 9V[}#(f$  
    ex_low_m2=ex_low_m1; q\|RI;W  
    ex_low_m1=Ex(2); F",TP,X  

0a^bAEP  
    Ex(KE)=ex_high_m2; 4>LaA7)v  
    ex_high_m2=ex_high_m1; N)AlQ'Lwx  
    ex_high_m1=Ex(KE-1); G%;>_E  
    %%%%%%%%%%%%%%%%%%计算磁场 u[Si=)`VPk  
    for i=1:KE-1 < .. (Y8LyY  
&5${k'  

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% f_[dFKoX  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 5X+`aB  
%%%%%%%%%%%%%%%%%%%%        1D              %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% QOY
MT(j  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2|& S2uq  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %_4#WI  
%%%%%%%%%初始化 IF|;;*Z8  
clear; saP%T~  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% [f6BA|   
%%%%%%%%%%%%%%%%系统参数 z,x
)Xx  
TimeT=3000;%迭代次数 GdNhEv  
KE=2000;%网格树木 [J}eNprg  
kc=450;%源的位置 VrP{U-`  
kpstart=500;%等离子体开始位置 Ys}^hy  
kpstop=1000;%等离子体终止位置 WPNw")t!  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% zEtsMU  
%%%%%%%%%%%%%物理参数 aK;OzB)  
c0=3e8;%真空中波速 r(uo-/7z  
zdelta=1e-9;%网格大小 qve'Gm)  
dt=zdelta/(2*c0);%时间间隔 P(b[|QF  
f=900e12;%Gause脉冲的载频 av|T|J/(  
d=3e-15%脉冲底座宽度 wn!=G~nB  
t0=2.25/f;%脉冲中心时间 sL~4~178  
u0=57e12%碰撞频率 e1h7~ j  
fpe=2000e12;%等离子体频率 ^@RvCJ+  
wpe=2*pi*fpe;%等离子体圆频率 nb=mY&q}~  
epsz=1/(4*pi*9*10^9); % 真空介电常数 +~iiy;i(  
mu=1/(c0^2*epsz);%磁常数 H@G$K@L  
ex_low_m1=0; EYj~Xj8_  
ex_low_m2=0; ,3T"fT-(  
ex_high_m1=0; =Q<7
[  
ex_high_m2=0; `vAcCahM  
a0=2*u0/dt+(2/dt)^2; em3
+V  
a1=-8/(dt)^2; rZ3ji(4HS  
a2=-2*u0/dt+(2/dt)^2; {cOx0=  
b0=wpe^2+2*u0/dt+(2/dt)^2; i]? Eq?k  
b1=2*wpe^2-8/(dt)^2; p<L{e~{!7f  
b2=wpe^2-2*u0/dt+(2/dt)^2; yTg|L9  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DGb1_2ZQ  
%%%%%%%%%%%%%初始化电磁场 z{\tn.67  
Ex=zeros(1,KE); TmLCmy!  
Ex_Pre=zeros(1,KE); 0>td[f  
Hy=zeros(1,KE); 6Yodx$  
Hy_Pre=zeros(1,KE); YSt*uOZK  
Dx=zeros(1,KE); R@jMFh;  
Dx_Pre=zeros(1,KE); Z^%a 1>`  
Sx1=zeros(1,KE); 0{z8pNrc  
Sx2=zeros(1,KE); XF)N_}X^  
Sx3=zeros(1,KE); gFHBIN;u  
Sx=zeros(1,KE); S@u46X>  
Dx=Ex; "IzAvKPM  
Dx1=Ex; d[,Rgdd@I  
Dx2=Ex; ~ E6e~  
Dx3=Ex; Q\kWQOB_  
Ex1=Ex;  d9k`  
Ex2=Ex; v/rBjUc+X  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 6JZ>&HA  
%%%%%%%%%%%%%%%%%%开始计算 S
N ?Z7  
for T=1:TimeT 6_QAE6A  
    %%%保存前一时间的电磁场 xCXsyZ2h  
    Ex_Pre=Ex; T%6JVFD  
    Hy_Pre=Hy; Os$E,4,py  
    %%%%中间差分计算Dx kOD
=H-vSi  
    for i=2:KE `b8nz7  
        Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1)); 7AT8QC`u  
    end EF\OM?R  
    %%%%%%%%加入源 aHuMm&  
    Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-t0)/d)^2); C/Z#NP~ *  
     Vlz\n  
    %%%计算电场Ex s) U1U6O  
    for i=1:kpstart-1 GLecB
F+>F  
        Ex(i)=Dx(i)/epsz; n8;p]{  
    end $RY-yKmi  
    for i=kpstop+1:KE :FS5BT$=  
        Ex(i)=Dx(i)/epsz; "J+L]IC?AD  
    end iv+a5  
    Dx3=Dx2; y@I9>}"y  
Dx2=Dx1; m=@xZw<  
Dx1=Dx; erUK;+2g  
    for i=kpstart:kpstop <Qih&P9;>  
                   Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i)); !;*flr`/  
    end 7,p.M)t)  
    Ex2=Ex1; EniV-Uj\D  
Ex1=Ex; f|w;u!U(  
    Sx3=Sx2; PaQ
lQ#  
    Sx2=Sx1; B//*hH >F  
    Ex(1)=ex_low_m2; Ya\:C]   
    ex_low_m2=ex_low_m1; J(d+EjC  
    ex_low_m1=Ex(2); 0>SA90Q  
"%D"h  
    Ex(KE)=ex_high_m2; G-9i   
    ex_high_m2=ex_high_m1; $%DoLpE>  
    ex_high_m1=Ex(KE-1); j]kgdAq>  
    %%%%%%%%%%%%%%%%%%计算磁场 T})q/oUqK  
    for i=1:KE-1 "o`?-bQ:  
        Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i)); NN'pBUR  
    end ;a1DIUm'  
    plot(Ex); \h#aPG<yo  
    grid on; a5t&{ajJ  
    pause(0.01); P7=`P  
end

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

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

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

3楼是？ n(tx'&U"R  
呵呵，不过感兴趣的人好像不多啊。

可不可以说一下Dx和Sx是什么来的呢？ ^Cy=L]  
最好有相关的物理公式，会更好理解点 F0x'^Z}Q;  
谢谢了

要说明的话，可能得整理一份资料出来才行，改天吧。 RA1K$D ?A  
现在天天写code，发现原来写的真是一个烂字啊。 >u+%H vzc  
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); (f;.
`W  
这里面竟然不把exp(-1*u0*dt))那些先提取出来，浪费运行时间啊。 z1nKj\AM2  
现在看到乘法反应都是有没有优化的方法，呵呵。

可以先S1与Ex的关系吗？ vN
ju|=Lo  
主要是 ~G1B}c]  
     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); <G'M/IR a  
         Ex(i)=Dx(i)/epsz-Sx1(i); 45k.U$<|  
         Ex(i)=Ex(i); <}T7;knO  
和 >qBJK)LHOv  
    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); N:;z~`  
         Ex(i)=Dx(i)/epsz-Sx1(i); 7AlL,&+  
h>F"GR?U_(  
这两组方程的来历不是太明白，一般看到的等离子波好像都不涉及碰撞速度的，但我在执行的时候发现碰撞速度对结果的影响和等离子振荡频率是差不多的

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