// 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 = (double ***)malloc((nx)*sizeof(double **));  for(i=0; i<(nx); i++)    {        ex[i] = (double **)malloc((ny+1)*sizeof(double *));      for(j=0; j<(ny+1); j++)	{	  ex[i][j] = (double *)malloc((nz+1)*sizeof(double));	  for(k=0; k<(nz+1); k++)	    {	      ex[i][j][k] = 0.0;	    }	}    }  allocatedBytes += ( (nx)*(ny+1)*(nz+1) * sizeof(double));    // allocate the array of ey components:  ey = (double ***)malloc((nx+1)*sizeof(double **));  for(i=0; i<(nx+1); i++)    {        ey[i] = (double **)malloc((ny)*sizeof(double *));      for(j=0; j<(ny); j++)	{	  ey[i][j] = (double *)malloc((nz+1)*sizeof(double));	  for(k=0; k<(nz+1); k++)	    {	      ey[i][j][k] = 0.0;	    }	}    }  allocatedBytes += ( (nx+1)*(ny)*(nz+1) * sizeof(double));  // allocate the array of ez components:  ez = (double ***)malloc((nx+1)*sizeof(double **));  for(i=0; i<(nx+1); i++)    {        ez[i] = (double **)malloc((ny+1)*sizeof(double *));      for(j=0; j<(ny+1); j++)	{	  ez[i][j] = (double *)malloc((nz)*sizeof(double));	  for(k=0; k<(nz); k++)	    {	      ez[i][j][k] = 0.0;	    }	}    }  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:  hx = (double ***)malloc((nx-1)*sizeof(double **));  for(i=0; i<(nx-1); i++)    {        hx[i] = (double **)malloc((ny)*sizeof(double *));      for(j=0; j<(ny); j++)	{	  hx[i][j] = (double *)malloc((nz)*sizeof(double));	  for(k=0; k<(nz); k++)	    {	      hx[i][j][k] = 0.0;	    }	}    }  allocatedBytes += ( (nx-1)*(ny)*(nz) * sizeof(double));    // allocate the array of hy components:  hy = (double ***)malloc((nx)*sizeof(double **));  for(i=0; i<(nx); i++)    {        hy[i] = (double **)malloc((ny-1)*sizeof(double *));      for(j=0; j<(ny-1); j++)	{	  hy[i][j] = (double *)malloc((nz)*sizeof(double));	  for(k=0; k<(nz); k++)	    {	      hy[i][j][k] = 0.0;	    }	}    }  allocatedBytes += ( (nx)*(ny-1)*(nz) * sizeof(double));    // allocate the array of hz components:  hz = (double ***)malloc((nx)*sizeof(double **));  for(i=0; i<(nx); i++)    {        hz[i] = (double **)malloc((ny)*sizeof(double *));      for(j=0; j<(ny); j++)	{	  hz[i][j] = (double *)malloc((nz-1)*sizeof(double));	  for(k=0; k<(nz-1); k++)	    {	      hz[i][j][k] = 0.0;	    }	}    }  allocatedBytes += ( (nx)*(ny)*(nz-1) * sizeof(double));  /////////////////////////////////////////////////////////////////////////////  // write some progress notes to standard output  // print out some identifying information   fprintf(stdout, "\n");  fprintf(stdout, "\n");  fprintf(stdout, "\n");  fprintf(stdout, "\n");  fprintf(stdout, "ToyFDTD1 version 1.03\n");  fprintf(stdout, "Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, ");  fprintf(stdout, "Matthew O'Keefe\n");  fprintf(stdout, "\n");  fprintf(stdout, "ToyFDTD1 is free software published under the terms\n");   fprintf(stdout, "of the GNU General Public License as published by the\n");   fprintf(stdout, "Free Software Foundation.\n");    fprintf(stdout, "\n");  fprintf(stdout, "\n");  // print out a bob command line,    //     including the dimensions of the output files:   fprintf(stdout, "bob -cmap chengGbry.cmap -s %dx%dx%d *.bob\n", nx+1, ny+1, nz);   fprintf(stdout, "\n");   // print out a viz command line:   fprintf(stdout, "viz ToyFDTD1c.viz\n");   fprintf(stdout, "\n");   fprintf(stdout, "\n");   // print out how much memory has been allocated:    fprintf(stdout, "Dynamically allocated %d bytes\n", allocatedBytes);   fprintf(stdout, "\n");   // print out some simulation parameters:   fprintf(stdout, "PLOT_MODULUS = %d\n", PLOT_MODULUS);   fprintf(stdout, "rectangular waveguide\n");   fprintf(stdout, "Stimulus = %lg Hertz ", FREQUENCY);   fprintf(stdout, "continuous plane wave\n");   fprintf(stdout, "\n");   fprintf(stdout, "Meshing parameters:\n");   fprintf(stdout, "%dx%dx%d cells\n", nx, ny, nz);   fprintf(stdout, "dx=%lg, dy=%lg, dz=%lg meters\n", dx, dy, dz);   fprintf(stdout, "%lg x %lg x %lg meter^3 simulation region\n",  	  GUIDE_WIDTH, GUIDE_HEIGHT, LENGTH_IN_WAVELENGTHS*lambda);   fprintf(stdout, "\n");   fprintf(stdout, "Time simulated will be %lg seconds, %d timesteps\n",  	  totalSimulatedTime, MAXIMUM_ITERATION);   fprintf(stdout, "\n");   fprintf(stdout, "3D output scaling parameters:\n");   fprintf(stdout, "Autoscaling every timestep\n");   fprintf(stdout, "\n");   fprintf(stdout, "\n");   // print out some info on each timestep:   fprintf(stdout, "Following is the iteration number and current\n");   fprintf(stdout, "simulated time for each timestep/iteration of\n");   fprintf(stdout, "the simulation.  For each timestep that 3D data is\n");   fprintf(stdout, "output to file, the maximum and minimum data\n");   fprintf(stdout, "values are printed here with the maximum and\n");   fprintf(stdout, "minimum scaled values in parentheses.\n");   fprintf(stdout, "\n");   /////////////////////////////////////////////////////////////////////////////  // open and start writing the .viz file  //    this file will be handy to feed parameters to viz if desired  while ((vizFilePointer = fopen("ToyFDTD1c.viz", "w")) == NULL)    {      fprintf(stderr, "Difficulty opening ToyFDTD1c.viz");      perror(" ");    }  fprintf(vizFilePointer, "#Viz V1.0\n");  fprintf(vizFilePointer, "time: %lg %lg\n", currentSimulatedTime, dt);  fprintf(vizFilePointer, "color: chengGbry.cmap\n");  fprintf(vizFilePointer, "\n");  /////////////////////////////////////////////////////////////////////////////  // main loop:  for(iteration = 0; iteration < MAXIMUM_ITERATION; iteration++)    {// mainloop      /////////////////////////////////////////////////////////////////////////       // Output section:       // time in simulated seconds that the simulation has progressed:      currentSimulatedTime = dt*(double)iteration;        // print to standard output the iteration number       //     and current simulated time:      fprintf(stdout, "#%d %lgsec", iteration, currentSimulatedTime);      // 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 ( (iteration % PLOT_MODULUS) == 0)	{// bob output section	  // create the filename for this iteration, 	  //     which includes the iteration number:	  sprintf(filename, "c_%06d.bob", iteration);	  // open a new data file for this iteration:	  while ((openFilePointer = fopen(filename, "wb")) == NULL)	    {// if the file fails to open, print an error message to 	     //     standard output:	      fprintf(stderr, "Difficulty opening c_%06d.bob", iteration);	      perror(" ");	    }	  // find the max and min values to be output this timestep:  	  min = FLT_MAX;	  max = -FLT_MAX;	  for(i=0; i<(nx+1); i++)	    { 	      for(j=0; j<(ny+1); j++)		{		  for(k=0; k<(nz); k++)		    {		      if (ez[i][j][k] < min) min = ez[i][j][k];		      if (ez[i][j][k] > max) max = ez[i][j][k];		    }		}	    }	  // update the tracking variables for minimum and maximum 	  //      values for the entire simulation:	  if (min < simulationMin) simulationMin = min;	  if (max > simulationMax) simulationMax = max;	  // set norm to be max or min, whichever is greater in magnitude:	  norm = (fabs(max) > fabs(min)) ? fabs(max) : fabs(min);	  if (norm == 0.0)	    {// if everything is zero, give norm a tiny value 	     //     to avoid division by zero:	      norm = DBL_EPSILON;	    }	  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:	  fprintf(stdout, "\t%lg(%d) < ez BoB < %lg(%d)",		  min, (int)(127.0 + scalingValue*min),		  max, (int)(127.0 + scalingValue*max));	  // 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:	  for(k=0; k<(nz); k++)	    { 	      for(j=0; j<(ny+1); j++)		{		  for(i=0; i<(nx+1); i++)