一維等離子體FDTD的Matlab源代碼兩種方法

1、一維等離子體FDTD的Matlab源代碼(兩種方法）% 1=V-V<  % *$g-:ILRuZ  %        1D              % oCz/HQoBk  % k9L;!TH1K  % 'D1xh  %初始化 QVgl(;lX  clear; 3yK!-Wp(  % utV_

2、 W&  %系統(tǒng)參數(shù) uwGcxOgg,  TimeT=3000;%迭代次數(shù) qITg"%t  KE=2000;%網(wǎng)格樹木 p4Z(+Aa  kc=450;%源的位置 f3y=Wxk  kpstart=500;%等離子體開始位置 |2A:eI8  kpstop=1000;%等離子體終止位置 'LDQgC*%  % A#E ;lm  %物理參數(shù) V !wj  c0=3e8;%真空中波速 3Jn ;  zdelta=1e-9;%網(wǎng)格大小 #GFro0$  dt=z

3、delta/(2*c0);%時間間隔 ) )Za&S*<  f=900e12;%Gause脈沖的載頻 <C*hokqqP  d=3e-15%脈沖底座寬度 .e-#yET  t0=2.25/f;%脈沖中心時間 uXiNj &Be  u0=57e12%碰撞頻率 m&SNz=  fpe=2000e12;%等離子體頻率 K (|dl:  wpe=2*pi*fpe;%等離子體圓頻率 m4Zk,1m.|  epsz=1/(4*pi*9*109); % 真空介電常數(shù) $/ ,tSm  mu=1

4、/(c02*epsz);%磁常數(shù) ;9#KeA _  ex_low_m1=0; yt2PU_),  ex_low_m2=0; CvdN"k  ex_high_m1=0; 8 FhdN  ex_high_m2=0; v r:=K  a0=2*u0/dt+(2/dt)2; 'hf8ZEW9'  a1=-8/(dt)2; +H2Qk4XFB  a2=-2*u0/dt+(2/dt)2; 2t,zLwBdnJ  b0=wpe2+2*u0/dt+(2/dt)2; *lb<$E="! &

5、#160;b1=2*wpe2-8/(dt)2; F:ELPs4"  b2=wpe2-2*u0/dt+(2/dt)2; Vw"  % &m vSiyKX  %初始化電磁場 ?X;RLpEc|A  Ex=zeros(1,KE); %XTI-B/K  Ex_Pre=zeros(1,KE); rM "l3hP  Hy=zeros(1,KE); i6N',&jFU  Hy_Pre=zeros(1,KE); E .kc(4  Dx=zeros(1,KE); XHh8_ &a

6、mp;  Dx_Pre=zeros(1,KE); K+iP 6B  Sx1=zeros(1,KE); (iGTACoF  Sx2=zeros(1,KE); We z 5N  Sx3=zeros(1,KE); U ;I9 bK8  Sx=zeros(1,KE); -.3wD"l  Dx=Ex; nF/OPd  Dx1=Ex; Qcii)s$js  Dx2=Ex; 3N:D6w-R  Dx3=Ex; :i7;w%B  Ex1=Ex; cS+>JL  Ex2=Ex; WjjB

7、<YKzF  % L.WljNo  %開始計算 MZI  for T=1:TimeT L_s:l9!r      %保存前一時間的電磁場 #o2hibq      Ex_Pre=Ex; o? $.fhD      Hy_Pre=Hy; bYPKh      %中間差分計算Dx 8sCv|cn      for i

8、=2:KE O| hpXkV          Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1); cFWc<55aX6      end Ip;F!s      %加入源 !Rt >xD      Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*(T*dt-t0)/d)2); 1!gbTeVlY

9、      >dGG>      %計算電場Ex w+ LAS      for i=1:kpstart-1 09Cez0          Ex(i)=Dx(i)/epsz; tNX|U:Y*      end DDH:)=;z      for

10、 i=kpstop+1:KE Um54fU          Ex(i)=Dx(i)/epsz; _f:W?$ho      end J9r|gJ(      Dx3=Dx2; C 6AUNRpl  Dx2=Dx1; w*JGUk  Dx1=Dx; > "=>3      for i=kpstart:kpstop F

11、G*r'tCr                     Ex(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex2(i); )TH# 1      end Em&6!      Ex2=Ex1; 'CkIz"Wd

12、 Ex1=Ex; vOpK Np      Sx3=Sx2; kq,ucU%>p      Sx2=Sx1; dKe_Q0      Ex(1)=ex_low_m2; RtIh-Z.9      ex_low_m2=ex_low_m1; 4u5-7TZ      ex_low_m1=Ex(2); HqT#$rv  DG:Z

13、=LuJr      Ex(KE)=ex_high_m2; 76h ,xi      ex_high_m2=ex_high_m1; SmSH2m-      ex_high_m1=Ex(KE-1); X=fYWjH,      %計算磁場 O*)Vhw'pK      for i=1:KE-1 XBu"-(   

14、;       Hy(i)=Hy(i)-(dt/(mu*zdelta)*(Ex(i+1)-Ex(i); GM f A,>      end *kDCliL      plot(Ex); )g#T9tx2D      grid on; CxOob1      pause(0.01); :hk5 .  end% FD

15、TD Main Function Jobs to Workers %*% 3-D FDTD code with PEC boundaries%*% This MATLAB M-file implements the finite-difference time-domain% solution of Maxwell's curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells.% % To illustrate the algorithm,

16、 an air-filled rectangular cavity % resonator is modeled. The length, width, and height of the % cavity are X cm (x-direction), Y cm (y-direction), and % Z cm (z-direction), respectively.% The computational domain is truncated using PEC boundary % conditions:% ex(i,j,k)=0 on the j=1, j=jb, k=1, and

17、k=kb planes% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossless walls of the cavity.% The cavity is excited by an additive current source oriented% along the z-direction. The source waveform is a differ

18、entiated % Gaussian pulse given by % J(t)=-J0*(t-t0)*exp(-(t-t0)2/tau2), % where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-% content pulse is approximately 7 GHz. The grid resolution % (dx = 2 mm) was chosen to provide at least 10 samples per % wavelength up through 15 GHz.% To execute

19、this M-file, type "fdtd3D" at the MATLAB prompt.% This M-file displays the FDTD-computed Ez fields at every other% time step, and records those frames in a movie matrix, M, which % is played at the end of the simulation using the "movie" command.%*function Ex,Ey,Ez=FDTD3D_Main(ha

20、ndles)global SimRunStop% if isdir('C:MATLAB7workcavityfigures')% mkdir 'C:MATLAB7workcavityfigures'% end%*% Grid Partition%*p.ip = get(handles.XdirPar,'Value');p.jp = get(handles.YdirPar,'Value');p.PN = get(handles.PartNo,'Value');%*% Grid Dimensons%*ie = get(

21、handles.xslider,'Value'); %number of grid cells in x-directionje = get(handles.yslider,'Value'); %number of grid cells in y-directionke = get(handles.zslider,'Value'); %number of grid cells in z-directionib=ie+1; jb=je+1; kb=ke+1;%*% All Domains Fields Ini.%*Ex=zeros(ie,jb,kb

22、);Ey=zeros(ib,je,kb);Ez=zeros(ib,jb,ke);Hx=zeros(ib,je,ke);Hy=zeros(ie,jb,ke);Hz=zeros(ie,je,kb);%*% Fundamental constants%*param.cc=2.99792458e8; %speed of light in free spaceparam.muz=4.0*pi*1.0e-7; %permeability of free spaceparam.epsz = 1.0/(param.cc*param.cc*param.muz); %permittivity of free sp

23、ace%*% Grid parameters%*param.is = get(handles.xsource,'Value'); %location of z-directed current sourceparam.js = get(handles.ysource,'Value'); %location of z-directed current sourceparam.kobs = floor(ke/2); %Surface of observationparam.dx = get(handles.CellSize,'Value'); %sp

24、ace increment of cubic latticeparam.dt = param.dx/(2.0*param.cc); %time stepparam.nmax = get(handles.TimeStep,'Value'); %total number of time steps%*% Differentiated Gaussian pulse excitation%*param.rtau=get(handles.tau,'Value')*100e-12;param.tau=param.rtau/param.dt;param.ndelay=3*pa

25、ram.tau;param.Amp = get(handles.sourceamp,'Value')*10e11;param.srcconst=-param.dt*param.Amp;%*% Material parameters%*param.eps= get(handles.epsilon,'Value');param.sig= get(handles.sigma,'Value'); %*% Updating coefficients%*param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*pa

26、ram.eps)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps);param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps);param.da=1.0;param.db=param.dt/param.muz/param.dx;%*% Calling FDTD Algorithm%*ex=zeros(ib,jb,kb);ey=zeros(ib,jb,kb);ez=zeros(ib,jb,kb);

27、hx=zeros(ib,jb,kb);hy=zeros(ib,jb,kb);hz=zeros(ib,jb,kb);X,Y,Z = meshgrid(1:ib,1:jb,1:kb); % Grid coordinatesPsim = zeros(param.nmax,1);Panl = zeros(param.nmax,1); if (p.ip = 1)&(p.jp = 0) x = ceil(ie/p.PN) param.a = 1:x-1:ie-x; param.b = x+1:x-1:ie; param.c = 1,1; param.d = je,je; m2 = 1; for n

28、=1:1:param.nmax for m1=1:1:p.PN ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); hx,hy,hz = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); end Psim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n); field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);

29、 Dyn_FFT pause(0.01); end elseif (p.ip = 0)&(p.jp = 1) y = ceil(je/p.PN); param.c = 1:y-1:je-y; param.d = y+1:y-1:je; param.a = 1,1; param.b = ie,ie; m1 = 1; for n=1:1:param.nmax for m2=1:1:p.PN ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); hx,hy,hz = Hfields(par

30、am,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); end Psim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n); field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); pause(0.01); endelseif (p.ip = 1)&(p.jp = 1) x = ceil(ie/p.PN); param.a = 1:x-1:ie-x; param.b = x+1:x-1:ie; y = ceil(je/p.PN); p

31、aram.c = 1:y-1:je-y; param.d = y+1:y-1:je; for n=1:1:param.nmax for m2=1:1:p.PN for m1=1:1:p.PN ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); hx,hy,hz = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); end end Psim(n),Panl(n) = Cavity_Power(param,handles

32、,ex,ey,ez,n); field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); pause(0.01); end else param.a = 1; param.b = ie; param.c = 1; param.d = je; m1 = 1;m2=1; for n=1:1:param.nmax ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); hx,hy,hz = Hfields(param,hx,hy,hz,ex,ey,ez

33、,ie,je,ke,ib,jb,kb,n,m1,m2,p); SimRunStop = get(handles.Stop_sim,'value'); if SimRunStop = 1 h = warndlg('Simulation Run is Stopped by User !','! Warning !'); waitfor(h); break; end Psim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n); if n>=2 Panl(n)=Panl(n) + Panl(n-

34、1); end field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); pause(0.01); end end%文件2% Cavity Field Viz. %function = field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p)%*% Cross-Section initialization%*tview = squeeze(ez(:,:,param.kobs);sview = squeeze(ez(:,param.js,:);ax1 = handles.axes1;ax2 =

35、handles.axes2;ax3 = handles.axes3;%*% Visualize fields%*timestep=int2str(n);ezview=squeeze(ez(:,:,param.kobs);exview=squeeze(ex(:,:,param.kobs);eyview=squeeze(ey(:,:,param.kobs);xmin = min(X(:);xmax = max(X(:);ymin = min(Y(:);ymax = max(Y(:);zmin = min(Z(:);daspect(2,2,1)xrange = linspace(xmin,xmax,

36、8);yrange = linspace(ymin,ymax,8);zrange = 3:4:15;cx cy cz = meshgrid(xrange,yrange,zrange);% sview=squeeze(ez(:,param.js,:);rect = -50 -35 360 210;rectp = -50 -40 350 260;axes(ax1);axis tightset(gca,'nextplot','replacechildren');E_total = sqrt(ex.2 + ey.2 + ez.2);etview=squeeze(E_to

37、tal(:,:,param.kobs);sview = squeeze(E_total(:,param.js,:);surf(etview');shading interp;caxis(-1.0 1.0);colorbar;axis image;title('Total E-Field, time step = ',timestep,'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD');xlabel('i co

38、ordinate');ylabel('j coordinate');set(gca,'fontname','Times New Roman','fontsize',10);% F1 = getframe(gca,rect);% M1 = frame2im(F1);% filename = fullfile('C:MATLAB7workcavityfigures',strcat('XY_ETotal',num2str(n),'.jpeg');% imwrite(M1,filen

39、ame,'jpeg')axes(ax2);axis tightset(gca,'nextplot','replacechildren');surf(sview');shading interp;caxis(-1.0 1.0);colorbar;axis image;title('Ez(i,j=13,k), time step = ',timestep,'fontname','Times New Roman','fontsize',12,'FontWeight'

40、,'BOLD');xlabel('i coordinate');ylabel('k coordinate');set(gca,'fontname','Times New Roman','fontsize',10);% F2 = getframe(gca,rect);% M2 = frame2im(F2);% filename = fullfile('C:MATLAB7workcavityfigures',strcat('XZ_ETotal',num2str(n),&#

41、39;.jpeg');% imwrite(M2,filename,'jpeg')%*% Cavity Power - Analitic expression %*axes(ax3);% axis tight% set(gca,'nextplot','replacechildren');t = 1:1:param.nmax;Psim = 1e3*Psim./max(Psim);Panl = 1e3*Panl./max(Panl);semilogy(t,Psim,'b.-',t,Panl,'rx-');str(

42、1) = 'Normalized Analytic Vs. Simulated Power'str(2) = 'as function of time-steps'title(str,'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD' );xlabel('Time Step');ylabel('Log(Power)');legend('Simulation',&#

43、39;Analytic','location','SouthEast');set(gca,'fontname','Times New Roman','fontsize',10);% F3 = getframe(gca,rectp);% M3 = frame2im(F3);% filename = fullfile('C:MATLAB7workcavityfigures',strcat('Power',num2str(n),'.jpeg');% imwrite(

44、M3,filename,'jpeg')%文件3% Source waveform functionfunction source=waveform(param,handles,n)ax1 = handles.axes1;type = get(handles.sourcetype,'value');amp = get(handles.sourceamp,'value')*1e4;f = get(handles.Frequency,'value')*1e9;switch type case 2 source = param.srcco

45、nst*(n-param.ndelay)*exp(-(n-param.ndelay)2/param.tau2); case 1 source = param.srcconst*sin(param.dt*n*2*pi*f); case 3 source = param.srcconst*exp(-(n-param.ndelay)2/param.tau2)*sin(2*pi*f*(n-param.ndelay)*param.dt); otherwise source = param.srcconst*(n-param.ndelay)*exp(-(n-param.ndelay)2/param.tau

46、2);end%文件4function hx,hy,hz = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,x,y,p)a = param.a(x);b = param.b(x);c = param.c(y);d = param.d(y);hx(a+1:b,c:d,1:ke)=. hx(a+1:b,c:d,1:ke)+. param.db*(ey(a+1:b,c:d,2:kb)-. ey(a+1:b,c:d,1:ke)+. ez(a+1:b,c:d,1:ke)-. ez(a+1:b,c+1:d+1,1:ke);hy(a:b,c+1:d,1

47、:ke)=. hy(a:b,c+1:d,1:ke)+. param.db*(ex(a:b,c+1:d,1:ke)-. ex(a:b,c+1:d,2:kb)+. ez(a+1:b+1,c+1:d,1:ke)-. ez(a:b,c+1:d,1:ke);hz(a:b,c:d,2:ke)=. hz(a:b,c:d,2:ke)+. param.db*(ex(a:b,c+1:d+1,2:ke)-. ex(a:b,c:d,2:ke)+. ey(a:b,c:d,2:ke)-. ey(a+1:b+1,c:d,2:ke);%文件5function ex,ey,ez=Efields(param,handles,ex

48、,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,x,y,p)a = param.a(x);b = param.b(x);c = param.c(y);d = param.d(y);if (param.is>=a)&(param.is<=b)|(param.js>=c)&(param.js<=d) source = waveform(param,handles,n);else source = 0;endex(a:b,c+1:d,2:ke)=. param.ca*ex(a:b,c+1:d,2:ke)+. param.cb*(hz(a

49、:b,c+1:d,2:ke)-. hz(a:b,c:d-1,2:ke)+. hy(a:b,c+1:d,1:ke-1)-. hy(a:b,c+1:d,2:ke);ey(a+1:b,c:d,2:ke)=. param.ca*ey(a+1:b,c:d,2:ke)+. param.cb*(hx(a+1:b,c:d,2:ke)-. hx(a+1:b,c:d,1:ke-1)+. hz(a:b-1,c:d,2:ke)-. hz(a+1:b,c:d,2:ke);ez(a+1:b,c+1:d,1:ke)=. param.ca*ez(a+1:b,c+1:d,1:ke)+. param.cb*(hx(a+1:b,c

50、:d-1,1:ke)-. hx(a+1:b,c+1:d,1:ke)+. hy(a+1:b,c+1:d,1:ke)-. hy(a:b-1,c+1:d,1:ke);ez(param.is,param.js,1:ke)=ez(param.is,param.js,1:ke)+source;function FDTDonedimensionpipei(L,d,T)%version1.0  終端匹配%FDTDonedimensionpipei(6,0.18,0.5e-9)t0=3*T;c=3e8;u=4*pi*1e-7;e=8.8541878e-12;dz=T*c/10;Nz=L/dz;Nz=f

51、ix(Nz);dt=dz/2/c;Ex=zeros(1,Nz+1);B=zeros(1,Nz+1);Hy=zeros(1,Nz);Nt=2*Nz;for n=0:Nt    t=n*dt;    F=exp(-(t-t0).2./T2);    Ex(1)=F;    for k=1:Nz        Hy(k)=Hy(k)+dt./u.*(Ex(k)-Ex(k+1)./dz;  &#

52、160; end    for k=1:Nz-1        Ex(k+1)=Ex(k+1)+dt./e.*(Hy(k)-Hy(k+1)./dz;    end    Ex(1)=B(2)+(c*dt-dz)./(c*dt+dz).*(Ex(2)-B(1);    Ex(Nz+1)=B(Nz)+(c*dt-dz)./(c*dt+dz).*(Ex(Nz)-B(Nz+1);   

53、 Vref1=d.*Ex(Nz-300);    Vref2=d.*Ex(Nz-100);    plot(t,Vref1,'s');    hold on;    plot(t,Vref2,'rx');    hold on;    B=Ex;endfunction FDTD_debug%constantsc_0 = 3E8;                  % Speed of light in free spaceNz = 100;                    % Number of cells in z-directi