toyfdtd1 c source code listing.htm

toyFDTD GNU

  //     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