toyfdtd1 c source code listing.htm

toyFDTD GNU

      ///////////////////////////////////////////////////////////////////////// 
      // 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++)
                    {
                      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]) -