fdtd3d_main.m

３Ｄ的ＦＤＴＤ模拟程序

% 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 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)*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(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;