三维电磁场FDTD程序(3D-FDTD matlab)PEC边界

% 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 Maxwell's 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/tau^2),

% 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:\MATLAB7\work\cavity\figures')

% mkdir 'C:\MATLAB7\work\cavity\figures'

% 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-direction

je = get(handles.yslider,'Value'); %number of grid cells in y-direction

ke = get(handles.zslider,'Value'); %number of grid cells in z-direction

ib=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

%***********************************************************************

https://www.360docs.net/doc/7011902335.html,=2.99792458e8; %speed of light in free space

param.muz=4.0*pi*1.0e-7; %permeability of free space

param.epsz = 1.0/(https://www.360docs.net/doc/7011902335.html,*https://www.360docs.net/doc/7011902335.html,*param.muz); %permittivity of free space

%***********************************************************************

% Grid parameters

%***********************************************************************

param.is = get(handles.xsource,'Value'); %location of z-directed current source param.js = get(handles.ysource,'Value'); %location of z-directed current source

param.kobs = floor(ke/2); %Surface of observation

param.dx = get(handles.CellSize,'Value');; %space increment of cubic lattice param.dt = param.dx/(2.0*https://www.360docs.net/doc/7011902335.html,); %time step

param.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.si g)/(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 coordinates

Psim = 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);

end

elseif ((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 tight

set(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:\MATLAB7\work\cavity\figures',strcat('XY_ETotal',num2str(n),'.jpeg')); % imwrite(M1,filename,'jpeg')

axes(ax2);

axis tight

set(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:\MATLAB7\work\cavity\figures',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:\MATLAB7\work\cavity\figures',strcat('Power',num2str(n),'.jpeg'));

% imwrite(M3,filename,'jpeg')

%文件3

% Source waveform function

function [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.tau^2));

case 1

source = param.srcconst*sin(param.dt*n*2*pi*f);

case 3

source =

param.srcconst*exp(-((n-param.ndelay)^2/param.tau^2))*sin(2*pi*f*(n-param.ndelay)*para m.dt);

otherwise

source =

param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)^2/param.tau^2));

end

%文件4

function [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));

%文件5

function [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<=b)|(param.js>=c)&(param.js<=d)

source = waveform(param,handles,n);

else

source = 0;

end

ex(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;