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

1、v1.0可编辑可修改维等离子体 FDTD的Matlab源代码（两种方法）% 1D%。%clear;%W!%TimeT=3000;%迭代次数KE=2000;%网格树木kc=450;%源的位置kpstart=500;% 等离子体开始位置kpstop=1000;%等离子体终止位置%c0=3e8;%真空中波速v1.0可编辑可修改zdelta=1e-9;% 网格大小dt=zdelta/(2*c0);% 时间间隔f=900e12;%Gause脉冲的载频d=3e-15%脉冲底座宽度t0=f;%脉冲中心时间u0=57e12% 撞频率fpe=2000e12;%等离子体频率wpe=2*pi*fpe;%等离子体圆频

2、率epsz=1/(4*pi*9*10A9); % 真空介电常数mu=1/(c0A2*epsz);% 磁常数ex_low_m1=0;ex_low_m2=0;ex_high_m1=0;ex_high_m2=0;a0=2*u0/dt+(2/dt)A2;a1=-8/(dt)A2;a2=-2*u0/dt+(2/dt)A2;b0=wpeA2+2*u0/dt+(2/dt)A2;b1=2*wpeA2-8/(dt)A2;b2=wpeA2-2*u0/dt+(2/dt)A2;% 磁场Ex=zeros(1,KE);Ex_Pre=zeros(1,KE);Hy=zeros(1,KE);Hy_Pre=zeros(1,KE)

3、;v1.0可编辑可修改Dx=zeros(1,KE);Dx_Pre=zeros(1,KE);Sx1=zeros(1,KE);Sx2=zeros(1,KE);Sx3=zeros(1,KE);Sx=zeros(1,KE);Dx=Ex;Dx1=Ex;Dx2=Ex;Dx3=Ex;Ex1=Ex;Ex2=Ex;% %for T=1:TimeT%:存前一时间的电磁场Ex_Pre=Ex;Hy_Pre=Hy;%间差分计算 Dxfor i=2:KEDx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1);end%Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*(T*dt-t0

4、)/dF2);%算电场Exv1.0可编辑可修改for i=1:kpstart-1Ex(i)=Dx(i)/epsz;endfor i=kpstop+1:KEEx(i)=Dx(i)/epsz;endDx3=Dx2;Dx2=Dx1;Dx1=Dx;for i=kpstart:kpstopEx(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)- b1*Ex1(i)-b2*Ex2(i);endEx2=Ex1;Ex1=Ex;Sx3=Sx2;Sx2=Sx1;Ex(1)=ex_low_m2;ex_low_m2=ex_low_m1;ex_low_m1=Ex(2);Ex(KE)=ex

5、_high_m2;ex_high_m2=ex_high_m1;ex_high_m1=Ex(KE-1);%for i=1:KE-1v1.0可编辑可修改Hy(i)=Hy(i)-(dt/(mu*zdelta)*(Ex(i+1)-Ex(i);end plot(Ex);grid on;pause;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

6、 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 (

7、z-direction), respectively.%v1.0可编辑可修改% The computational domain is truncated using PEC boundary% conditions:%ex(i,j,k)=0onthe j=1, j=jb, k=1, and k=kb planes%ey(i,j,k)=0onthe i=1, i=ib, k=1, and k=kb planes%ez(i,j,k)=0onthe i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossle

8、ss 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)A2/tauA2),% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-% content pulse is approx

9、imately 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 mat

10、rix, M, which% is played at the end of the simulation using the movie command.%o片 *function Ex,Ey,Ez=FDTD3D_Main(handles)global SimRunStop% if isdir(C:MATLAB7workcavityfigures)v1.0可编辑可修改% mkdir C:MATLAB7workcavityfigures% end%*% Grid Partition%*=get,Value);=get,Value);=get,Value);%*% Grid Dimensons%

11、*ie = get,Value);%number of grid cells in x-directionje = get,Value);%number of grid cells in y-directionke = get,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);v1.0可编辑可修改Ez=zeros(ib,jb,ke);Hx=zeros(ib,je,ke);Hy=

12、zeros(ie,jb,ke);Hz=zeros(ie,je,kb);0片 *% Fundamental constantso片 *%speed of light in free space=*pi*; %permeability of free space =*; %permittivity of free space%* % Grid parameters %*=get,Value); %location of z-directed current source=get,Value); %location of z-directed current source=floor(ke/2);

13、%Surface of observation=get,Value); %space increment of cubic lattice=*; %time step=get,Value); %total number of time steps%*% Differentiated Gaussian pulse excitationv1.0可编辑可修改%*=get,Value)*100e-12;=；=3*;=get,Value)*10e11;%*Material parameters %*=get,Value);=get,Value);%* % Updating coefficients %*

14、=*/*)/+*/*) =/+*/*);%* % Calling FDTD Algorithm %* ex=zeros(ib,jb,kb);v1.0可编辑可修改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,1);Panl = zeros,1);if ( = 1)& = 0)x = ceil(ie/=1:x-1:ie-x;=x

15、+1:x-1:ie;=1,1;=je,je;m2 = 1;for n=1:1:for m1=1:1:ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p) ;hx,hy,hzHfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endPsim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n);10v1.0可编辑可修改field_viz(param,handles,ex,ey,

16、ez,X,Y,Z,n,Psim,Panl,p);Dyn_FFTpause;endelseif ( = 0)& = 1)y = ceil(je/;=1：y-1：je-y;=y+1：y-1：je;=1,1;=ie,ie;m1 = 1;for n=1:1:for m2=1:1:ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p) ;hx,hy,hzHfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endPsim(n),Panl(

17、n) = Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause;endelseif ( = 1)& = 1)x = ceil(ie/;=1:x-1:ie-x;=x+1:x-1:ie;11v1.0可编辑可修改y = ceil(je/;=1：y-1：je-y;=y+1：y-1：je;for n=1:1:for m2=1:1:for m1=1:1:ex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,k

18、e,ib,jb,kb,n,m1,m2,p) ;hx,hy,hzHfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endendPsim(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;endelse=1;=ie;=1;=je;m1 = 1;m2=1;for n=1:1:ex,ey,ez=Efields(param,handles,ex,ey,ez,hx

19、,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p)12v1.0可编辑可修改5hx,hy,hzHfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);SimRunStop = get,value);if SimRunStop = 1h= warndlg(Simulation Run is Stopped by User !,!Warning !);waitfor(h);break;endPsim(n),Panl(n) = Cavity_Power(param,handles,ex,ey,ez,n);if n=2

20、Panl(n)=Panl(n) + Panl(n-1);endfield_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause;endend。以2 2% Cavity Field Viz. %function = field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p)%*% Cross-Section initialization13v1.0可编辑可修改%*tview = squeeze(ez(:,:,);sview = squeeze(ez(:,:);ax1 =;ax2 =;ax3 =;

21、%*% Visualize fields%*timestep=int2str(n);ezview=squeeze(ez(:,:,);exview=squeeze(ex(:,:,);eyview=squeeze(ey(:,:,);xmin = min(X(:);xmax = max(X(:);ymin = min(丫(:)；ymax = max(丫(:)；zmin = min(Z(:);daspect(2,2,1)xrange = linspace(xmin,xmax,8);yrange = linspace(ymin,ymax,8);zrange = 3:4:15;14v1.0可编辑可修改cx

22、 cy cz = meshgrid(xrange,yrange,zrange);% sview=squeeze(ez(:,:);rect = -50 -35 360 210;rectp = -50 -40 350 260;axes(ax1);axis tightset(gca,nextplot,replacechildren);E_total = sqrt(ex.A2 + ey.A2 + ez.A2);etview=squeeze(E_total(:,:,);sview = squeeze(E_total(:,:);surf(etview);shading interp;caxis();col

23、orbar;axis image;title(TotalE-Field, time step = ,timestep,fontname,TimesNewRoman,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

24、_ETotal,num2str(n),.jpeg15v1.0可编辑可修改);% imwrite(M1,filename,jpeg)axes(ax2);axis tightset(gca,nextplot,replacechildren);surf(sview);shading interp;caxis();colorbar;axis image;title(Ez(i,j=13,k),time step = ,timestep,fontname,TimesNewRoman,fontsize,12,FontWeight,BOLD);xlabel(i coordinate);ylabel(k coo

25、rdinate);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%*16v1.0可编辑可修改axes(ax3);% axis tight% set(gca,nextplo

26、t,replacechildren);t = 1:1:;Psim = 1e3*Psim./max(Psim);Panl = 1e3*Panl./max(Panl);semilogy(t,Psim,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

27、(Simulation,Analytic,location,SouthEast);set(gca,fontname,Times New Roman,fontsize,10);% F3 = getframe(gca,rectp);% M3 = frame2im(F3);%filenamefullfile(C:MATLAB7workcavityfigures,strcat(Power,num2str(n),.jpeg);% imwrite(M3,filename,jpeg)。汉3 3% Source waveform functionfunction source=waveform(param,h

28、andles,n)ax1 =;17v1.0可编辑可修改type = get,value);amp = get,value)*1e4;f = get,value)*1e9;switch typecase 2source = *exp(-(A2/A2);case 1source = *sin*n*2*pi*f);case 3source = *exp(-(A2/A2)*sin(2*pi*f*;otherwisesource = *exp(-(A2/A2);end。汉4 4function hx,hy,hz = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,

29、jb,kb,n,x,y,p)a = (x);b = (x);c = (y);d = (y);hx(a+1:b,c:d,1:ke)=.hx(a+1:b,c:d,1:ke)+.*(ey(a+1:b,c:d,2:kb)-.ey(a+1:b,c:d,1:ke)+.ez(a+1:b,c:d,1:ke)-.18v1.0可编辑可修改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)+.*(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:

30、d,1:ke);hz(a:b,c:d,2:ke)=.hz(a:b,c:d,2:ke)+.*(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);。以5 5functionex,ey,ez=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,x,y,p)a = (x);b = (x);c = (y);d = (y);if =a)&=c)&=d)source = waveform(param,handles,n);else19

31、v1.0可编辑可修改source = 0;endex(a:b,c+1:d,2:ke)=.*ex(a:b,c+1:d,2:ke)+.*(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)=.*ey(a+1:b,c:d,2:ke)+.*(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)=.*ez(a+

32、1:b,c+1:d,1:ke)+.*(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,1:ke)=ez,1:ke)+source;20v1.0可编辑可修改function FDTDonedimensionpipei(L,d,T)%终端匹配%FDTDonedimensionpipei(6,t0=3*T;c=3e8;u=4*pi*1e-7;e=;dz=T*c/10;Nz=L/dz;Nz=fix(Nz);dt=dz/2/c;Ex=zeros(1,Nz+1);B=zero

33、s(1,Nz+1);Hy=zeros(1,Nz);Nt=2*Nz;for n=0:Ntt=n*dt;F=exp(-(t-t0).A2./TA2);Ex(1)=F;for k=1:NzHy(k尸Hy(k)+dt./u.*(Ex(k)-Ex(k+1)./dz;end21v1.0可编辑可修改for k=1:Nz-1Ex(k+1)=Ex(k+1)+dt./e.*(Hy(k)-Hy(k+1)./dz;endEx(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);Vr

34、ef1=d.*Ex(Nz-300);Vref2=d.*Ex(Nz-100);plot(t,Vref1,s);hold on;plot(t,Vref2,rx);hold on;B=Ex;function FDTD_debug%constantsc_0 = 3E8;Nz = 100;Nt = 200;N = Nz;E = zeros(Nz,1);E2 = zeros(Nz,1);end% Speed of light in free space% Number of cells in z-direction% Number of time steps% Global number of cells

35、% E component22v1.0可编辑可修改Es = zeros(Nz,Nz);% Es is the output for surf function%H = zeros(Nz,1);% H componentpulse = 0;% Pulse%to be setmu_0 = *pi*;% Permeability of free spaceeps_0 = (c_0*c_0*mu_0);%c_ref = c_0;eps_ref = eps_0;mu_ref = mu_0;f_0 = 10e9;f_ref = f_0;omega_0 = *pi*f_0;omega_ref = omega

36、_0;lambda_ref = c_ref/f_ref;Dx_ref = lambda_ref/20;Dz = Dx_ref;nDz = Dz/Dx_ref;Z = Nz*Dz;r = 1;Dt_ref = r*Dx_ref/c_ref;Dt = Dt_ref;% Source positionNz_Source = *Nz;N_Source = Nz_Source;for n = 1:Nt-1Permittivty of free space% Reference velocity% Excitation frequency% Reference frequency% Excitation

37、circular frequency% Reference wavelength% Reference cells width% Normalized time step% Reference time stept = Dt_ref*r*(n-1);% Actual time23E1 = zeros(Nz,1);v1.0可编辑可修改%Source function seriesSource_type = 1;switch Source_type case 1% modified source functionncycles = 2;if t ncycles*pi/(omega_0)pulse

38、= *( - cos(omega_0*t/ncycles) ) * sin(omega_0*t); elsepulse = 0; end% sigle cos source functionif t *pi/(omega_0)pulse8*c_0A2*Dt_refA2*mu_ref*omega_0*cos(omega_0*t); elsepulse = 0; end% Gaussian pulseif t Dt_ref*r*50pulse-40*c_0*(t-t*25/(n-1)*exp(-(t-t*25/(n-1)A2/2/(50/A2/(t/(n-1)A2); elsepulse = 0;

39、 endendE2(N_Source) = E2(N_Source) - r*pulse;24v1.0可编辑可修改E1 = E2;E2 = E;m = 2:Nz-1;E(m) = rA2*(E2(m+1) - 2*E2(m) + E2(m-1) + 2*E2(m) - E1(m);%Boundary% Neumann% Mur ABCright side%Mur ABCleftE(1) = 0; E(Nz) = 0;Dirichlet%E(1) = E(2);E(Nz) = E(Nz-1);%E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz);z =

40、Nz%E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1);side z = 0%display%choice*display = 3; % choice: 1 = line plot; 2 = imagesc; 3 = surfswitch display= 1:Nz;plot(i, E(i);axis(0 Nz -4 4);A = rot90(E);imagesc(A);25v1.0可编辑可修改i = 1:Nz;for j = 1:NzEs(i,j) = E(i);endEs = rot90(Es);j = 1:Nz;surf(i,j,Es);axis(0 Nz 0

41、 Nz -4 4)otherwiseA = rot90(E);imagesc(A);endpause;end26v1.0可编辑可修改%*% 3-D FDTD code with PEC boundaries%*% Program author: Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% 608-265-5739% Date of this

42、version: February 2000% 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 27v1.0可编辑可修改% resonator

43、 is modeled. The length, width, and height of the% cavity are cm (x-direction), cm (y-direction), and% 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=

44、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

45、(-(t-t0)A2/tauA2),% 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-fi

46、le 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.%o片 *28v1.0可编辑可修改clear0片 *% Fundamental constantso片 *cc=; %speed of light in free spacemuz=*pi*; %permeability of fr

47、ee spaceepsz=(cc*cc*muz); %permittivity of free space%*% Grid parameters%*ie=50; %number of grid cells in x-directionje=24; %number of grid cells in y-directionke=10; %number of grid cells in z-directionib=ie+1;jb=je+1;kb=ke+1;is=26; %location of z-directed current sourcejs=13; %location of z-direct

48、ed current sourcekobs=5;29v1.0可编辑可修改dx=; %space increment of cubic latticedt=dx/*cc); %time stepnmax=500; %total number of time steps0%1tzmm*% Differentiated Gaussian pulse excitation%* rtau=;tau=rtau/dt;ndelay=3*tau;srcconst=-dt*+11;%* % Material parameters %* eps=;sig=;%*% Updating coefficients%*c

49、a=(dt*sig)/*epsz*eps)/+(dt*sig)/*epsz*eps);cb=(dt/epsz/eps/dx)/+(dt*sig)/*epsz*eps);30v1.0可编辑可修改da=;db=dt/muz/dx;0片 *% Field arrayso片 *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);%*% Movie initialization%*tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:);subplot(position, ),pcolor(tview);shading flat;caxis();colorbar;axis image;title(Ez(i,j,k=5), time step = 0);31v1.0可编辑可修改xlabel(i coordinate);ylabel(j coordinate);subplot(position, ),pcolor(sview)