toyfdtd1.c

toyFDTD一个简单的三维FDTD的程序用Fortran编的

// 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 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_MAX

main()	
{// 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)));