電磁場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文件實現(xiàn)了麥克斯韋旋度方程的有限差分方法在三維笛卡爾空間點陣組成的統(tǒng)一的立方網(wǎng)格細(xì)胞。為了說明算法,建模一個充氣矩形空腔諧振器。腔的長度、寬度和高度是10.0厘米(x方向)(y方向)4.8厘米,2.0厘米(z三個方向),分別為。計算域邊界截斷使用壓電陶瓷這些壓電陶瓷外無損耗腔壁邊界形式。腔是由面向添加劑電流源的興奮沿著z三個方向。源波形是一個差異化的高斯脈沖的半最大值光譜帶寬zero-dc -內(nèi)容脈沖大約是7 GHz。網(wǎng)格分辨率(dx = 2毫米)被選為至少10樣品每波長通過15 GHz。執(zhí)行這個m文件,輸入“fdtd3D”MATLAB提示。這個m文件顯示在每個其他FDTD-computed Ez字段時間步長,并記錄這些框架在電影矩陣,M播放結(jié)束時模擬使用“電影”命令。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