toyfdtd1.f90

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

!! ToyFDTD1, version 1.03 (F90)
!!      The if-I-can-do-it-you-can-do-it FDTD! 
!! Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe 
!!               1999 Max Smirnoff, Matt Rundquist (F90 translation)
!! 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,
!! Max Smirnoff, Matt Rundquist
!! 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.03 from version 1.02: Updating contact & web info.
!! Changes to version 1.02 from version 1.0:  There is no 1.0 version of the 
!!      FORTRAN code.  For a listing of the changes to the C code, see 
!!      the changelog file.  
!! For some notes on the F90 translation, see the README file.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!! This is a very simple Yee algorithm 3D FDTD code in F90 implementing
!!      the free space form of Maxwell's equations on a Cartesian grid.
!! There are no internal materials or geometry.  
!! The code 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 homepage 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/


!! 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.0d9
          !! frequency of the stimulus in Hertz
#define GUIDE_WIDTH 0.0229d0
          !! meters
#define GUIDE_HEIGHT 0.0102d0
          !! meters
#define LENGTH_IN_WAVELENGTHS 5.0d0
          !! length of the waveguide in wavelengths of the stimulus wave
#define CELLS_PER_WAVELENGTH 25.0d0
          !! minimum number of grid cells per wavelength in the x, y, and
          !!     z directions

!! physical constants
#define LIGHT_SPEED     299792458.0d0
          !! speed of light in a vacuum in meters/second
#define LIGHT_SPEED_SQUARED 89875517873681764.0d0    
          !! m^2/s^2
#define MU_0 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271d-6
          !! permeability of free space in henry/meter
#define EPSILON_0 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932d-12
          !! permittivity of free space in farad/meter

!!
!! These parameters are allowed in C through include files in the compilation.
!!
#define M_PI 3.14159265358979323846d0
          !! pi as defined in /usr/include/math.h on SGI IRIX 6.2
#define FLT_MAX	3.40282347d+38
          !! max float as defined in /usr/include/float.h on SGI IRIX 6.2
#define	DBL_EPSILON	2.2204460492503131d-16
          !! DBL_EPSILON as defined in /usr/include/float.h on SGI IRIX 6.2
#define SIZEOF_DOUBLE 8
          !! bytes in a double-precision floating point number


program ToyFDTD1
implicit none !! All undeclared variables will be reported as errors;
              !! this prevents some major problems.

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! variable declarations
  integer :: i,j,k
           !! indices of the 3D array cells
  integer :: nx,ny,nz
           !! total number of cells along the x,y,z axes, respectively
  integer :: allocatedBytes = 0;          
           !! a counter to track number of bytes allocated

  integer :: iteration = 0
           !! counter to track how many timesteps have been computed
  real*8 :: stimulus = 0.0
           !! value of the stimulus at current timestep
  real*8 :: currentSimulatedTime = 0.0
           !! time in simulated seconds that the simulation has progressed
  real*8 :: totalSimulatedTime = 0.0
           !! time in seconds that will be simulated by the program
  real*8 :: omega
           !! angular frequency in radians/second
  real*8 :: lambda
           !! wavelength of the stimulus in meters
  real*8 :: dx,dy,dz;
           !! space differentials (or dimensions of a single cell) in meters
  real*8 :: dt;
           !! time differential (how much time between timesteps) in seconds
  real*8 :: dtmudx,dtepsdx
           !! physical constants used in the field update equations  
  real*8 :: dtmudy,dtepsdy
           !! physical constants used in the field update equations  
  real*8 :: dtmudz,dtepsdz
           !! physical constants used in the field update equations  

  real*8,target,allocatable,dimension(:,:,:) :: ex,ey,ez
           !! pointers to the arrays of ex, ey, and ez values
  real*8,target,allocatable,dimension(:,:,:) :: hx,hy,hz
           !! pointers to the arrays of hx, hy, and hz values

  !! bob output routine variables:
  character(LEN=1024) :: filename 
           !! filename variable for 3D bob files
  real*8 :: simulationMin = FLT_MAX
           !! tracks minimum value output by the entire simulation
  real*8 :: simulationMax = -FLT_MAX
           !! tracks maximum value output by the entire simulation
  real*8 :: min, max
           !! these track minimum and maximum values output in one timestep
  real*8 :: norm
           !! norm is set each iteration to be max or min, whichever is 
           !!     greater in magnitude
  real*8 :: scalingValue
           !! multiplier used in output scaling, calculated every timestep

  character(LEN=10) :: numbers = "0123456789" 
           !! this is a handy string with all the decimal digits

  integer :: temp
  integer :: ios, record
           !! these are just dummy variables used during file writes 

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! 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.0d0*M_PI*FREQUENCY

  !! set ny and dy:
  !! start with 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:
  do while (dy >= lambda/CELLS_PER_WAVELENGTH)
    ny = ny + 1
    dy = GUIDE_WIDTH/ny
  enddo

  !! start nz and dz:
  !! start with 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:
  do while (dz >= lambda/CELLS_PER_WAVELENGTH)
    nz = nz + 1
    dz = GUIDE_HEIGHT/nz
  enddo

  !! set dx, nx, and dt:
  !! set dx equal to dy or dz, whichever is smaller:
  if(dy < dz) then
    dx = dy