// 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//// lemiller@borg.umn.edu// // 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 other notes on the revision history, see the changelog file//      which should have come with the code and can also be found at//      the ToyFDTD website.///////////////////////////////////////////////////////////////////////////////// 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/#include <math.h> #include <stdio.h>         #include <float.h>// program control constants#define MAXIMUM_ITERATION 1000          // total number of timesteps to be computed#define 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.  #define FREQUENCY 10.0e9               // frequency of the stimulus in Hertz#define GUIDE_WIDTH 0.0229          // meters#define GUIDE_HEIGHT 0.0102          // meters#define LENGTH_IN_WAVELENGTHS 5.0          // length of the waveguide in wavelengths of the stimulus wave#define CELLS_PER_WAVELENGTH 25.0          // minimum number of grid cells per wavelength in the x, y, and          //     z directions// physical constants#define LIGHT_SPEED 299792458.0                 // speed of light in a vacuum in meters/second#define LIGHT_SPEED_SQUARED 89875517873681764.0                  // m^2/s^2#define MU_0 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6          // permeability of free space in henry/meter#define 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 FLT_MAXmain()	{// main  /////////////////////////////////////////////////////////////////////////////  // variable declarations  int i,j,k;                // indices of the 3D array of cells  int nx, ny, nz;                    // total number of cells along the x, y, and z axes, respectively  int allocatedBytes = 0;                     // a counter to track number of bytes allocated  int iteration = 0;                     // counter to track how many timesteps have been computed  double stimulus = 0.0;                  // value of the stimulus at a given timestep  double currentSimulatedTime = 0.0;           // time in simulated seconds that the simulation has progressed  double totalSimulatedTime = 0.0;                  // time in seconds that will be simulated by the program  double omega;                  // angular frequency in radians/second  double lambda;                  // wavelength of the stimulus in meters  double dx, dy, dz;            // space differentials (or dimensions of a single cell) in meters  double dt;           // time differential (how much time between timesteps) in seconds  double dtmudx, dtepsdx;                  // physical constants used in the field update equations   double dtmudy, dtepsdy;                  // physical constants used in the field update equations   double dtmudz, dtepsdz;                  // physical constants used in the field update equations   double ***ex, ***ey, ***ez;             // pointers to the arrays of ex, ey, and ez values  double ***hx, ***hy, ***hz;           // pointers to the arrays of hx, hy, and hz values  // bob output routine variables:  char filename[1024];           // filename variable for 3D bob files  double simulationMin = FLT_MAX;           // tracks minimum value output by the entire simulation  double simulationMax = -FLT_MAX;           // tracks maximum value output by the entire simulation  double min, max;           // these track minimum and maximum values output in one timestep  double norm;           // norm 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*M_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++;      dy = GUIDE_WIDTH/ny;    }  // 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++;      dz = GUIDE_HEIGHT/nz;    }  // set dx, nx, and dt:  // set dx equal to dy or dz, whichever is smaller:  dx = (dy < dz) ? dy : dz;  // choose nx to make the guide LENGTH_IN_WAVELENGTHS   //     wavelengths long:    nx = (int)(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)));