Help needed with fdtd matlab script

[kabukicho2001]

Here it is:
% 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:1aram.nmax
for m1=1:1.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,Panl] = 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:1aram.nmax
for m2=1:1.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,Panl] = 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:1aram.nmax
for m2=1:1.PN
for m1=1:1.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,Panl] = 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:1aram.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,Panl] = Cavity_Power(param,handles,ex,ey,ez,n);
if n>=2
Panl=Panl + 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;
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,'.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,'.jpeg'));
% imwrite(M2,filename,'jpeg')

%***********************************************************************
% Cavity Power - Analitic expression
%***********************************************************************

axes(ax3);
% axis tight
% set(gca,'nextplot','replacechildren');

t = [1:1aram.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,'.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)*param.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;
d = param.d;

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;
d = param.d;

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;