一维等离子体FDTD的Matlab源代码(两种方法).doc

一维等离子体FDTD的Matlab源代码(两种方法）% 1=V-V % *$g-:ILRuZ %1D% oCz/HQoBk % k9L;!TH1K % D1xh %初始化 QVgl(;lX clear; 3yK!-Wp( % utV_ W& %系统参数 uwGcxOgg, TimeT=3000;%迭代次数 qITg%t KE=2000;%网格树木 p4Z(+Aa kc=450;%源的位置 f3y=Wxk kpstart=500;%等离子体开始位置 |2A:eI8 kpstop=1000;%等离子体终止位置 LDQgC*% % A#E ;lm %物理参数 V !wj c0=3e8;%真空中波速 3Jn ; zdelta=1e-9;%网格大小 #GFro0$ dt=zdelta/(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); % 真空介电常数 $/ ,tSm mu=1/(c02*epsz);%磁常数 ;9#KeA _ ex_low_m1=0; yt2PU_), ex_low_m2=0; CvdNk 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; *lbJL Ex2=Ex; WjjBYKzF % 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=2:KE O| hpXkV Dx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1); cFWcxD Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*(T*dt-t0)/d)2); 1!gbTeVlY dGG %计算电场Ex w+ LAS for i=1:kpstart-1 09Cez0 Ex(i)=Dx(i)/epsz; tNX|U:Y* end DDH:)=;z for 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 FG*rtCr 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; CkIzWd 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=LuJr Ex(KE)=ex_high_m2; 76h ,xi ex_high_m2=ex_high_m1; SmSH2m- ex_high_m1=Ex(KE-1); X=fYWjH, %计算磁场 O*)VhwpK for i=1:KE-1 XBu-( 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% FDTD Main Function Jobs to Workers %*% 3-D FDTD code with PEC boundaries%*% This MATLAB M-file implements the finite-difference time-domain% solution of Maxwells curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells.% % To illustrate the algorithm, 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 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 differentiated % 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 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(handles)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(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);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 space%*% 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); %space 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*param.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*param.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);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=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); 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(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); 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); param.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,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,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-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 = 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,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_total(:,:,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 coordinate);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,filename,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,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),.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(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,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(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.srcconst*(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.tau2);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: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,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=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: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: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=fix(Nz);dt=dz/2/c;Ex=zeros(1,Nz+1);B=zeros(1,Nz+1);Hy=zeros