电磁场FDTD算法以及仿真图.doc

function output_args = Untitled2( input_args )%UNTITLED2 Summary of this function goes here% Detailed explanation goes here%*% 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 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% resonator is modeled. The length, width, and height of the % cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and % 2.0 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.%这个MATLAB的m文件实现了麦克斯韦旋度方程的有限差分方法在三维笛卡尔空间点阵组成的统一的立方网格细胞。为了说明算法,建模一个充气矩形空腔谐振器。腔的长度、宽度和高度是10.0厘米(x方向)(y方向)4.8厘米,2.0厘米(z三个方向),分别为。计算域边界截断使用压电陶瓷这些压电陶瓷外无损耗腔壁边界形式。腔是由面向添加剂电流源的兴奋沿着z三个方向。源波形是一个差异化的高斯脉冲的半最大值光谱带宽zero-dc -内容脉冲大约是7 GHz。网格分辨率(dx = 2毫米)被选为至少10样品每波长通过15 GHz。执行这个m文件,输入“fdtd3D”MATLAB提示。这个m文件显示在每个其他FDTD-computed Ez字段时间步长,并记录这些框架在电影矩阵,M播放结束时模拟使用“电影”命令。clear%*% Fundamental constants%*cc=2.99792458e8; %speed of light in free space muz=4.0*pi*1.0e-7; %permeability of free space epsz=1.0/(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-directed current source kobs=5; dx=0.002; %space increment of cubic latticedt=dx/(2.0*cc); %time step nmax=500; %total number of time steps%* % Differentiated Gaussian pulse excitation%*rtau=50.0e-12;tau=rtau/dt;ndelay=3*tau;srcconst=-dt*3.0e+11;% Material parameters%*eps=1.0;sig=0.0;%* % Updating coefficients%*ca=(1.0-(dt*sig)/(2.0*epsz*eps)/(1.0+(dt*sig)/(2.0*epsz*eps); cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps);da=1.0;db=dt/muz/dx;%* % Field arrays%*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,0.15 0.45 0.7 0.45),pcolor(tview);shading flat;caxis(-1.0 1.0);colorbar;axis image;title(Ez(i,j,k=5), time step = 0);xlabel(i coordinate);ylabel(j coordinate);subplot(position,0.15 0.10 0.7 0.25),pcolor(sview);shading flat;caxis(-1.0 1.0);colorbar;axis image;title(Ez(i,j=13,k), time step = 0);xlabel(i coordinate);ylabel(k coordinate);rect=get(gcf,Position);rect(1:2)=0 0;M=moviein(nmax/2,gcf,rect);%* % BEGIN TIME-STEPPING LOOP%*for n=1:nmax%*% Update electric fields%*ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+. cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+. hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke);ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+. cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+. hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke);ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+. cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+. hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke);ez(is,js,1:ke)=ez(is,js,1:ke)+. srcconst*(n-ndelay)*exp(-(n-ndelay)2/tau2);%*% Update magnetic fields%*hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+. db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+. ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke);hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+. db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+. ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke);hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+. db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+. ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke);%*% Visualize fields%*if mod(n,2)=0; timestep=int2str(n);tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:); subplot(position,0.15 0.45 0.7 0.45),pcolor(tview);shading flat;caxis(-1.0 1.0);colorbar;axis image;title(Ez(i,j,k=5), time step = ,timestep);xlabel(i coordinate);ylabel(j coordinate);subplot(position,0.15 0.10