else    dx = dz  end if  !! 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.0d0/(LIGHT_SPEED*sqrt(1.0d0/(dx*dx) + 1.0d0/(dy*dy) + 1.0d0/(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.  !! 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:  allocate(ex(0:nx-1,0:ny,0:nz))  do i=0,(nx-1),1     do j=0,ny,1        do k=0,nz,1           ex(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx)*(ny+1)*(nz+1) * SIZEOF_DOUBLE)  !! allocate the array of ey components:  allocate(ey(0:nx,0:ny-1,0:nz))  do i=0,nx,1     do j=0,(ny-1),1        do k=0,nz,1           ey(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx+1)*(ny)*(nz+1) * SIZEOF_DOUBLE)  !! allocate the array of ez components:  allocate(ez(0:nx,0:ny,0:nz-1))  do i=0,nx,1     do j=0,ny,1        do k=0,(nz-1),1           ez(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx+1)*(ny+1)*(nz) * SIZEOF_DOUBLE)  !! 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:  allocate(hx(0:nx-2,0:ny-1,0:nz-1))  do i=0,(nx-2),1     do j=0,(ny-1),1        do k=0,(nz-1),1           hx(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx-1)*(ny)*(nz) * SIZEOF_DOUBLE)  !! allocate the array of hy components:  allocate(hy(0:nx-1,0:ny-2,0:nz-1))  do i=0,(nx-1),1     do j=0,(ny-2),1        do k=0,(nz-1),1           hy(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx)*(ny-1)*(nz) * SIZEOF_DOUBLE)  !! allocate the array of hz components:  allocate(hz(0:nx-1,0:ny-1,0:nz-2))  do i=0,(nx-1),1     do j=0,(ny-1),1        do k=0,(nz-2),1           hz(i,j,k) = 0.0d0        end do     end do  end do  allocatedBytes = allocatedBytes + ( (nx)*(ny)*(nz-1) * SIZEOF_DOUBLE)  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! write some progress notes to standard output  !! print out some identifying information   print*," "  print*," "  print*," "  print*," "  print*,"ToyFDTD1 version 1.03 (F90)"  print*,"Copyright (C) 1998, 1999 Laurie E. Miller, Paul Hayes, "  print*,"Matthew O'Keefe, Max Smirnoff, Matt Rundquist"  print*," "  print*,"ToyFDTD1 is free software published under the terms"   print*,"of the GNU General Public License as published by the"   print*,"Free Software Foundation. "   print*," "  print*," "  !! print out a bob command line,   !! including the dimensions of the output files:  print'(" bob -cmap chengGbry.cmap -s ",I5.5,"x",I5.5,"x",I5.5," *.bob")',nx+1,ny+1,nz  print*," "  !! print out a viz command line:  print*,"viz ToyFDTD1f90.viz"  print*," "  print*," "  !! print out how much memory has been allocated:   print*,"Dynamically allocated ", allocatedBytes, " bytes"  print*," "  !! print out some simulation parameters:   print*,"PLOT_MODULUS = ", PLOT_MODULUS  print*,"rectangular waveguide"  print*,"Stimulus = ", FREQUENCY, " Hertz continuous plane wave"  print*," "  print*,"Meshing parameters:"  print '(I4,"x",I3,"x",I3," cells")', nx, ny, nz  print '(" dx=",F11.8," dy=",F11.8," dz=",F11.8," meters")', dx, dy, dz  print '(F11.8, " x ", F11.8, " x ", F11.8, " meter^3 simulation region")', &       & GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda  print*," "  print '(" Time simulated will be ",E12.5E2," seconds, ",I5.5," timesteps")', &       & totalSimulatedTime, MAXIMUM_ITERATION  print*," "  print*,"3D output scaling parameters:"  print*,"Autoscaling every timestep"  print*," "  print*," "  !! print out some info on each timestep:   print*,"Following is the iteration number and current "  print*,"simulated time for each timestep/iteration of "  print*,"the simulation.  For each timestep that data is "  print*,"output to file, the maximum and minimum data "  print*,"values are printed here with the maximum and "  print*,"minimum scaled values in parentheses. "  print*," "  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! open and start writing the .viz file  !!    this file will be handy to feed parameters to viz if desired  open(unit=33,file="ToyFDTD1f90.viz",iostat=ios)  !! if the file fails to open, print an error message to   !!     standard output:  if(ios /= 0) then    print*,"Difficulty opening ToyFDTD1f90.viz"    STOP  end if  write(unit=33, fmt='("#Viz V1.0")')  write(unit=33, fmt='("time: ",F5.3,ES12.5E2)') currentSimulatedTime,dt  write(unit=33, fmt='("color: chengGbry.cmap")')  write(unit=33, fmt='("")')                                       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  !! main loop:  do iteration=0,(MAXIMUM_ITERATION-1),1           !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      !! Output section:       !! time in simulated seconds that the simulation has progressed:      currentSimulatedTime = dt*iteration      !! print to standard output the iteration number       !!     and current simulated time:      write(*,372) iteration, currentSimulatedTime372   format ("#",i5," ", ES14.5E2, "sec", $)      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!      !! 3D data output every PLOT_MODULUS timesteps:      !!     The first time through the main loop all the data written to       !!     file will be zeros.  If anything is nonzero, there's a bug.  :>       if(mod(iteration,PLOT_MODULUS) == 0) then         !! Construct the filename:         !! The numerical part is obtained by successively         !!      dividing the current iteration by ten, and the         !!      remainders of these divisions are the digits.         !!      For example, 12345 would be converted as follows:         !!      12345/10 = 1234+5/10 (5 is the first digit in the number)         !!      1234/10 = 123 + 4/10 (4 is the next digit)         !!      ....and so on..this is performed six times to create         !!      a filename out of any 6-digit number.         temp = iteration         filename(1:4) = "f90_"         filename(10:10) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(9:9) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(8:8) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(7:7) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(6:6) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(5:5) = numbers(mod(temp,10)+1:mod(temp,10)+1)         temp = int(temp/10)         filename(11:14) = ".bob"         !! open a new data file for this iteration:         open(unit=69,file=filename,access="direct",iostat=ios,recl=1)         !! if the file fails to open, print an error message to 	 !!     standard output:         if(ios /= 0) then           print*,"Difficulty opening ",filename           STOP         end if         !! find the max and min values to be output this timestep:           min = FLT_MAX         max = -FLT_MAX         do i=0,(nx),1            do j=0,(ny),1               do k=0,(nz-1),1                  if(ez(i,j,k) < min) then                     min = ez(i,j,k)                  end if                  if(ez(i,j,k) > max) then                     max = ez(i,j,k)                  end if               end do            end do         end do         !!  update the tracking variables for minimum and maximum          !!  values for the entire simulation:         if(min < simulationMin) then            simulationMin = min         end if         if(max > simulationMax) then            simulationMax = max         end if                !! set norm to be max or min, whichever is greater in magnitude:         if(abs(max) > abs(min)) then            norm = abs(max)         else            norm = abs(min)         end if         !! if everything is zero, give norm a tiny value 	 !!     to avoid division by zero:         if(norm == 0.0d0) then            norm = DBL_EPSILON         end if         scalingValue = 127.0/norm         !! write to standard output the minimum and maximum values          !!     from this iteration and the minimum and maximum values         !!     that will be written to the bob file this iteration:         write(*, 373) min, int(127.0d0 + scalingValue * min), max, int(127.0d0 + scalingValue * max)373      format(" ", f9.5, " (", i3, ") < ez BoB < ", f9.5, " (", i3, ")", $)         !! scale each ez value in the mesh to the range of integers          !!     from zero through 254 and write them to the output          !!     file for this iteration:         record = 1         do k=0,(nz-1),1            do j=0,(ny),1               do i=0,(nx),1