toyfdtd1.m

This a translation of the ToyFDTD c code available from the web site http://www.borg.umn.edu/toyfdt

function toyFDTD1(interactive)
%toyFDTD1(interactive)
% nargin==0 : run without using sliceomatic.
% interactive=0 : stop at every saved time step and run sliceomatic.m for viewing the ez field.
% interactive=1 :stop only at the end and run sliceomatic.m for viewing of the ez field
%
% This is a translated version of the work by the authors below:
% It has been modified to work with MATLAB R13
% You will need to download sliceomatic.m from MATLAB EXCHANGE in order to view the results
% interactively. I dont know how to use VIZ or BOB which is the format of the files used in 
% the original C code. The original code is included in the ZIP file. The authors have another
% version which has additional features. If someone wants to submit that version - please do!
%
% See below for credits...
% Have fun! 
% Patrick Moran, Airsprite Technologies Inc., 2003


%  ToyFDTD1, version 1.03 
%       The if-I-can-do-it-you-can-do-it FDTD! 
%  Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe 

%  This program is free software; you can redistribute it and/or 
%      modify it under the terms of the GNU General Public License 
%      as published by the Free Software Foundation; either version 2
%      of the License, or any later version, with the following conditions
%      attached in addition to any and all conditions of the GNU
%      General Public License:
%      When reporting or displaying any results or animations created
%      using this code or modification of this code, make the appropriate
%      citation referencing ToyFDTD1 by name and including the version
%      number.  
% 
%  This program is distributed in the hope that it will be useful,
%      but WITHOUT ANY WARRANTY; without even the implied warranty 
%      of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
%      See the GNU General Public License for more details.
% 
%  You should have received a copy of the GNU General Public License
%      along with this program; if not, write to the Free Software
%      Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  
%      02111-1307  USA

%  Contacting the authors:
% 
%  Laurie E. Miller, Paul Hayes, Matthew O'Keefe
%  Department of Electrical and Computer Engineering
%       200 Union Street S. E.
%       Minneapolis, MN 55455
% 
%  [e-mail]
%  
%  http:% www.borg.umn.edu/toyfdtd/
%  http:% www.borg.umn.edu/toyfdtd/ToyFDTD1.html
% http:%www.toyfdtd.org/

% This code is here for everyone, but not everyone will need something 
%      so simple, and not everyone will need to read all the comments.  
% This file is over 700 lines long, but less than 400 of that is actually
%      code. 

% This ToyFDTD1 is a stripped-down, minimalist 3D FDTD code.  It 
%      illustrates the minimum factors that must be considered to 
%      create a simple FDTD simulation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Changes to version 1.02 from version 1.03: updating contact & web info.
%      For further notes on the revision history, see the changelog file.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% This is a very simple Yee algorithm 3D FDTD code in C implementing
%      the free space form of Maxwell's equations on a Cartesian grid.
% There are no internal materials or geometry.  
% The code as delivered simulates an idealized rectangular waveguide 
%      by treating the interior of the mesh as free space/air and enforcing
%      PEC (Perfect Electric Conductor) conditions on the faces of the mesh.
% The problem is taken from Field and Wave Electromagnetics, 2nd ed., by 
%      David K. Cheng, pages 554-555.  It is a WG-16 waveguide 
%      useful for X-band applications, interior width = 2.29cm, 
%      interior height = 1.02cm.  The frequency (10 GHz) is chosen to be 
%      in the middle of the frequency range for TE10 operation.  
% Boundaries: PEC (Perfect Electric Conductor).
% Stimulus: A simplified sinusoidal plane wave emanates from x = 0 face.
% 3D output: The electric field intensity vector components in the direction 
%      of the height of the guide (ez) are output to file every 
%      PLOT_MODULUS timesteps (and for the last timestep), scaled to 
%      the range of integers from zero through 254.  (The colormap 
%      included in the tar file assigns rgb values to the range zero through
%      255.)   Scaling is performed on each timestep individually, since 
%      it is not known in advance what the maximum and minimum 
%      values will be for the entire simulation.  The integer value 127 is 
%      held to be equal to the data zero for every timestep.  This method 
%      of autoscaling every timestep can be very helpful in a simulation 
%      where the intensities are sometimes strong and sometimes faint, 
%      since it will highlight the presence and structure of faint signals 
%      when stronger signals have left the mesh.  
%      Each timestep has it's own output file.  This data output file format
%      can be used in several visualization tools, such as animabob and viz. 
% Other output: Notes on the progress of the simulation are written to standard
%      output as the program runs.  
%      A .viz file is output to feed parameters to viz, should viz later be
%      used to view the data files.

% Some terminology used here:
%
% This code implements a Cartesian mesh with space differentials 
%     of dx, dy, dz.
% This means that a point in the mesh has neighboring points dx meters 
%     away in the direction of the positive and negative x-axis,
%     and neighboring points dy meters away in the directions 
%     of the +- y-axis, and neighboring points dz meters away 
%     in the directions of the +- z-axis,
% The mesh has nx cells in the x direction, ny in the y direction, 
%     and nz in the z direction.
% ex, ey, and ez refer to the arrays of electric field intensity vectors 
%     -- for example, ex is a 3-dimensional array consisting of the 
%     x component of the E field intensity vector for every point in the 
%     mesh.  ex[i][j][k] refers to the x component of the E field intensity 
%     vector at point [i][j][k].  
% hx, hy, and hz refer to the arrays of magnetic field intensity vectors.
%
% dt is the time differential -- the length of each timestep in seconds.
%
% bob is a file format that stands for 'brick of bytes', meaning a string 
%     of bytes that can be interpreted as a 3-dimensional array of byte 
%     values (integers from zero through 255).  animabob is a free 
%     visualization tool that displays and animates a sequence of bricks 
%     of bytes.  For more information on animabob or to download a copy, 
%     see the ToyFDTD website at http:%www.borg.umn.edu/toyfdtd/
%
% viz is another free visualization tool that displays and animates 
%     brick-of-byte files. For more information on viz or to download a copy, 
%     see the ToyFDTD website at http:%www.borg.umn.edu/toyfdtd/

% MATLAB fid constants
stdout=1;
stderr=2;

% program control constants
MAXIMUM_ITERATION=1000;% total number of timesteps to be computed
PLOT_MODULUS =5;% The program will output 3D data every PLOT_MODULUS timesteps,
%     except for the last iteration computed, which is always
%     output.  So if MAXIMUM_ITERATION is not an integer
%     multiple of PLOT_MODULUS, the last timestep output will
%     come after a shorter interval than that separating
%     previous outputs.  
FREQUENCY =10.0e9  ; % frequency of the stimulus in Hertz
GUIDE_WIDTH =0.0229; % meters
GUIDE_HEIGHT= 0.0102;% meters
LENGTH_IN_WAVELENGTHS =5.0; % length of the waveguide in wavelengths of the stimulus wave
CELLS_PER_WAVELENGTH =25.0; % minimum number of grid cells per wavelength in the x, y, and  z directions
% physical constants
LIGHT_SPEED =299792458.0 ;      
% speed of light in a vacuum in meters/second
LIGHT_SPEED_SQUARED =89875517873681764.0 ;       
% m^2/s^2
MU_0 =1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;
% permeability of free space in henry/meter
EPSILON_0=8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12;
% permittivity of free space in farad/meter
% The value used for pi is M_PI as found in /usr/include/math.h on SGI 
%     IRIX 6.2.  Other such constants used here are DBL_EPSILON and REALMAX
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% variable declarations
% ii,jj,kk;                  % indices of the 3D array of cells
% nx, ny, nz;                % total number of cells along the x, y, and z axes, respectively
allocatedBytes = 0;          % a counter to track number of bytes allocated
iteration = 0;               % counter to track how many timesteps have been computed
stimulus = 0.0;              % value of the stimulus at a given timestep
currentSimulatedTime = 0.0;  % time in simulated seconds that the simulation has progressed
totalSimulatedTime = 0.0;    % time in seconds that will be simulated by the program
%  omega;              % angular frequency in radians/second
%  lambda;             % wavelength of the stimulus in meters
%  dx, dy, dz;         % space differentials (or dimensions of a single cell) in meters
%  dt;                 % time differential (how much time between timesteps) in seconds
%  dtmudx, dtepsdx;    % physical constants used in the field update equations 
%  dtmudy, dtepsdy;    % physical constants used in the field update equations 
%  dtmudz, dtepsdz;    % physical constants used in the field update equations 
%  ex, ey, ez;         % pointers to the arrays of ex, ey, and ez values
%  hx, hy, hz;         % pointers to the arrays of hx, hy, and hz values

% bob output routine variables:
% filename;                    % filename variable for 3D bob files
simulationMin = REALMAX;       % tracks minimum value output by the entire simulation
simulationMax = -REALMAX;      % tracks maximum value output by the entire simulation
% double min, max;             % these track minimum and maximum values output in one timestep
% double nnorm;                % nnorm is set each iteration to be max or min, whichever is greater in magnitude
% double scalingValue;         % multiplier used in output scaling, calculated every timestep
% FILE *openFilePointer;       % pointer to 3D bob output file
% FILE *vizFilePointer;        % pointer to viz file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% setting up the problem to be modeled
%
% David K. Cheng, Field and Wave Electromagnetics, 2nd ed., 
%     pages 554-555. 
% Rectangular waveguide, interior width = 2.29cm, interior height = 1.02cm.
% This is a WG-16 waveguide useful for X-band applications.
%
% Choosing nx, ny, and nz:
% There should be at least 20 cells per wavelength in each direction, 
%     but we'll go with 25 so the animation will look prettier.   
%     (CELLS_PER_WAVELENGTH was set to 25.0 in the global 
%     constants at the beginning of the code.)
% The number of cells along the width of the guide and the width of 
%     those cells should fit the guide width exactly, so that ny*dy 
%     = GUIDE_WIDTH meters.  
%     The same should be true for nz*dz = GUIDE_HEIGHT meters.  
% dx is chosen to be dy or dz -- whichever is smaller
% nx is chosen to make the guide LENGTH_IN_WAVELENGTHS 
%     wavelengths long.  
% 
% dt is chosen for Courant stability; the timestep must be kept small 
%     enough so that the plane wave only travels one cell length 
%     (one dx) in a single timestep.  Otherwise FDTD cannot keep up 
%     with the signal propagation, since FDTD computes a cell only from 
%     it's immediate neighbors.  

% wavelength in meters:
lambda = LIGHT_SPEED/FREQUENCY; 
% angular frequency in radians/second:
omega = 2.0*pi*FREQUENCY; 
% set ny and dy:
% start with a small ny:
ny = 3;  
% calculate dy from the guide width and ny:
dy = GUIDE_WIDTH/ny;
% until dy is less than a twenty-fifth of a wavelength,
%     increment ny and recalculate dy:
while (dy >= lambda/CELLS_PER_WAVELENGTH)     
    ny=ny+1;
    dy = GUIDE_WIDTH/ny;
end
% set nz and dz:
% start with a small nz:
nz = 3;  
% calculate dz from the guide height and nz:
dz = GUIDE_HEIGHT/nz;
% until dz is less than a twenty-fifth of a wavelength,
%     increment nz and recalculate dz:
while(dz >= lambda/CELLS_PER_WAVELENGTH)
    nz=nz+1;
    dz = GUIDE_HEIGHT/nz;
end
% set dx, nx, and dt:

% set dx equal to dy or dz, whichever is smaller:
dx=min(dy,dz);
% choose nx to make the guide LENGTH_IN_WAVELENGTHS wavelengths long:  
nx =  ceil(LENGTH_IN_WAVELENGTHS*lambda/dx);
% chose dt for Courant stability:
dt = 1.0/(LIGHT_SPEED*sqrt(1.0/(dx*dx) + 1.0/(dy*dy) + 1.0/(dz*dz)));
% time in seconds that will be simulated by the program:
totalSimulatedTime = MAXIMUM_ITERATION*dt;
% constants used in the field update equations:
dtmudx = dt/(MU_0*dx);
dtepsdx = dt/(EPSILON_0*dx);
dtmudy = dt/(MU_0*dy);
dtepsdy = dt/(EPSILON_0*dy);
dtmudz = dt/(MU_0*dz);
dtepsdz = dt/(EPSILON_0*dz);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%/
% memory allocation for the FDTD mesh:
% There is a separate array for each of the six vector components, 
%      ex, ey, ez, hx, hy, and hz.
% The mesh is set up so that tangential E vectors form the outer faces of 
%     the simulation volume.  There are nx*ny*nz cells in the mesh, but 
%     there are nx*(ny+1)*(nz+1) ex component vectors in the mesh.
%     There are (nx+1)*ny*(nz+1) ey component vectors in the mesh.
%     There are (nx+1)*(ny+1)*nz ez component vectors in the mesh.
% If you draw out a 2-dimensional slice of the mesh, you'll see why
%     this is.  For example if you have a 3x3x3 cell mesh, and you 
%     draw the E field components on the z=0 face, you'll see that
%     the face has 12 ex component vectors, 3 in the x-direction
%     and 4 in the y-direction.  That face also has 12 ey components,
%     4 in the x-direction and 3 in the y-direction.  

% Allocate memory for the E field arrays:

% allocate the array of ex components:
ex=zeros([nx,ny+1,nz+1]);
allocatedBytes =prod(size(ex)); % sizeof(double)

% allocate the array of ey components:
ey=zeros([nx+1,ny,nz+1]);
allocatedBytes =allocatedBytes +prod(size(ey));

% allocate the array of ez components:
ez=zeros([nx+1,ny+1,nz]);
allocatedBytes =allocatedBytes +prod(size(ez));

% Allocate the H field arrays:
% Since the H arrays are staggered half a step off 
%     from the E arrays in every direction, the H 
%     arrays are one cell smaller in the x, y, and z 
%     directions than the corresponding E arrays. 
% By this arrangement, the outer faces of the mesh
%     consist of E components only, and the H 
%     components lie only in the interior of the mesh.  

% allocate the array of hx components:
hx=zeros([nx-1,ny,nz]);
allocatedBytes = allocatedBytes +prod(size(hx));

% allocate the array of hy components: