fdtd3d_main[1].m

fdtd一维到三维算法

% 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