{		      putc((int)(127.0 + 				 scalingValue*ez[i][j][k]), openFilePointer);		    }		}	    }	  	  // close the output file for this iteration:	  fclose(openFilePointer);	  // write the dimensions and name of the output file for this 	  //     iteration to the viz file:	  fprintf(vizFilePointer, "%dx%dx%d %s\n", nx+1, ny+1, nz, filename);	  // write identification of the corners of the mesh and the max	  //     and min values for this iteration to the viz file:	  fprintf(vizFilePointer, "bbox: 0.0 0.0 0.0 %lg %lg %lg %lg %lg\n",		  dx*(double)nx, dy*(double)ny, dz*(double)nz, min, max);	}// end bob output section      /////////////////////////////////////////////////////////////////////////      // Compute the stimulus: a plane wave emanates from the x=0 face:      //     The length of the guide lies in the x-direction, the width of the       //     guide lies in the y-direction, and the height of the guide lies      //     in the z-direction.  So the guide is sourced by all the ez       //     components on the stimulus face.        stimulus = sin(omega*currentSimulatedTime);      for (i=0; i<(1); i++)	{ 	  for(j=0; j<(ny+1); j++)	    {	      for(k=0; k<nz; k++)		{		  ez[i][j][k] = stimulus;		}	    }	}      /////////////////////////////////////////////////////////////////////////      // Update the interior of the mesh:      //    all vector components except those on the faces of the mesh.        //      // Update all the H field vector components within the mesh:      //     Since all H vectors are internal, all H values are updated here.      //     Note that the normal H vectors on the faces of the mesh are not       //     computed here, and in fact were never allocated -- the normal       //     H components on the faces of the mesh are never used to update       //     any other value, so they are left out of the memory allocation        //     entirely.               // Update the hx values:      for(i=0; i<(nx-1); i++)	{  	  for(j=0; j<(ny); j++)	    {	      for(k=0; k<(nz); k++)		{		  hx[i][j][k] += (dtmudz*(ey[i+1][j][k+1] - ey[i+1][j][k]) - 		                  dtmudy*(ez[i+1][j+1][k] - ez[i+1][j][k]));		}	    }	}      // Update the hy values:      for(i=0; i<(nx); i++)	{  	  for(j=0; j<(ny-1); j++)	    {	      for(k=0; k<(nz); k++)		{		  hy[i][j][k] +=  (dtmudx*(ez[i+1][j+1][k] - ez[i][j+1][k]) -      		                   dtmudz*(ex[i][j+1][k+1] - ex[i][j+1][k]));		}	    }	}      // Update the hz values:      for(i=0; i<(nx); i++)	{  	  for(j=0; j<(ny); j++)	    {	      for(k=0; k<(nz-1); k++)		{		  hz[i][j][k] +=  (dtmudy*(ex[i][j+1][k+1] - ex[i][j][k+1]) -      		                   dtmudx*(ey[i+1][j][k+1] - ey[i][j][k+1]));		}	    }	}      // Update the E field vector components.        // The values on the faces of the mesh are not updated here; they       //      are handled by the boundary condition computation       //      (and stimulus computation).        // Update the ex values:      for(i=0; i<(nx); i++)	{  	  for(j=1; j<(ny); j++)	    {	      for(k=1; k<(nz); k++)		{		  ex[i][j][k] += (dtepsdy*(hz[i][j][k-1] - hz[i][j-1][k-1]) -				  dtepsdz*(hy[i][j-1][k] - hy[i][j-1][k-1]));		}	    }	}            // Update the ey values:      for(i=1; i<(nx); i++)	{  	  for(j=0; j<(ny); j++)	    {	      for(k=1; k<(nz); k++)		{		  ey[i][j][k] += (dtepsdz*(hx[i-1][j][k] - hx[i-1][j][k-1]) -				  dtepsdx*(hz[i][j][k-1] - hz[i-1][j][k-1]));		}	    }	}      // Update the ez values:      for(i=1; i<(nx); i++)	{  	  for(j=1; j<(ny); j++)	    {	      for(k=0; k<(nz); k++)		{		  ez[i][j][k] += (dtepsdx*(hy[i][j-1][k] - hy[i-1][j-1][k]) -				  dtepsdy*(hx[i-1][j][k] - hx[i-1][j-1][k]));		}	    }	}            fprintf(stdout, "\n");      /////////////////////////////////////////////////////////////////////////      // Compute the boundary conditions:      // OK, so I'm yanking your chain on this one.  The PEC condition is       // enforced by setting the tangential E field components on all the      // faces of the mesh to zero every timestep (except the stimulus       // face).  But the lazy/efficient way out is to initialize those       // vectors to zero and never compute them again, which is exactly       // what happens in this code.            }// end mainloop  /////////////////////////////////////////////////////////////////////   // Output section:    // The output routine is repeated one last time to write out    // the last data computed.       // 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 for the last timestep:   // 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++)	    {	      putc((int)(127.0 + 			 scalingValue*ez[i][j][k]), openFilePointer);	    }        }    }    // close the output file for this iteration:  fclose(openFilePointer);    // write the dimensions and name of the output file for this   //     iteration to the viz file:  fprintf(vizFilePointer, "%dx%dx%d %s\n", nx+1, ny+1, nz, filename);  // write identification of the corners of the mesh and the max  //     and min values for this iteration to the viz file:  fprintf(vizFilePointer, "bbox: 0.0 0.0 0.0 %lg %lg %lg %lg %lg\n",	  dx*(double)nx, dy*(double)ny, dz*(double)nz, min, max);  /////////////////////////////////////////////////////////////////////////////  // close the viz file for this simulation:  fclose(vizFilePointer);  /////////////////////////////////////////////////////////////////////////////  // write some progress notes to standard output:  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 was %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");   // print out simulationMin and simulationMax:   fprintf(stdout, "Minimum output value was %lg \n", simulationMin);   fprintf(stdout, "Maximum output value was %lg \n", simulationMax);   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");   fprintf(stdout, "\n");   fprintf(stdout, "\n"); }// end main