% 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
%***********************************************************************
param.cc=2.99792458e8; %speed of light in free space
param.muz=4.0*pi*1.0e-7; %permeability of free space
param.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 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*param.cc); %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.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 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 SimRunSt