readconfig.c

a full 3D simulation of electromagnetic waves with efficient absorbing boundary a full 3D simulation

/* Copyright notice:
   radarFDTD - A free 3D-FDTD simulation for EM-waves with GPML-ABCs ;-)
   Copyright (C) 2000 Carsten Aulbert
   ( 
   I used the following program as a 'manual', so in my opinion this needs to be quoted:
   ToyFDTD, version 1.02 
   The if-I-can-do-it-you-can-do-it FDTD! 
   Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe 
   )

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

void ReadNextLine (FILE *file, char s[]) {
  /*  Read next line, which is not a comment, i.e. there is no '#' at the beginning */
   s[0] = '#';
   while ((s[0] == '#') && (!feof(file)))
      fscanf (file, "%[^\n]\n", s);
};

int currentConfigFileVersion = 3; /* this is needed to detect older config-files */

int nx, ny, nz; /* number of cells in that direction */
PRECISION dim_X, dim_Y, dim_Z; /* overall dimension of simulation space */
int pmlWidth; /* width of absobing pmls */
int maxAllowedMemory; /* Maximum number of bytes allowd for storing */
int maxIteration; /* max number of iteration steps */
int timeToPlotReceiver, timeToPlotSnapshot; /* decide when to plot */
int maxMaterials; /* maximum number of materials */
int transX, transY, transZ; /* position of transmitter */
int maxReceiver; /* max number of receiver boxes*/
int **receiver; /* receiver positions */
int **tmpReceiver; /* temp. receiver array */
int **tmpTransmitter; /* temp. transmitter positions */
int **transmitter; /* transmitter positions */
int maxTransmitter; /* max. number of transmitter boxes */
PRECISION shortestWavelength; /* shortest desired wavelength [m] */
PRECISION theoreticalReflection; 
char *materialFilename;
int numberOfBoxes, currentMaterial, numberOfReceiverBoxes, numberOfTransmitterBoxes, tmpCounter;
int snapshotMode, snapshotPosition, snapshotPlane;
PRECISION transmitterAmplitude;
PRECISION transmitterCharacteristicTime;
PRECISION transmitterFrequency;
PRECISION transmitterParameter;
PRECISION simulationLength;
PRECISION maxSigma;
PRECISION maxStretching;
PRECISION minWavelength;
int stretchSteepness;
int transmitterMode;
int transmitterStimulusComponent;
int storage = 0;
void ReadConfigfile(char *directory, char *filename)
{
  FILE *configFile;
  float a,b,c,d;
  int h,i,j,k;
  PRECISION maxEpsilon;
  PRECISION tmpdt;
  PRECISION tmpSigma, tmpEpsilon;
  char string[4000];
  PRECISION maxStep;
  PRECISION tmpReflection;
  int boxStartX, boxStartY, boxStartZ, boxEndX, boxEndY, boxEndZ;
  int warningStop = 0;
  
  /* open configfile */

  fprintf(stdout, "Opening configfile %s\n\n", filename);
 
  configFile = fopen(filename, "r");
  
  if (configFile == NULL)
    {
      fprintf(stderr, "Cannot open file %s. Please check this!\n", filename);
      exit(1);
    }

  /* ok, file is opened now read that cryptic file */

  /* check version of config-file */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &h);

  if (h != currentConfigFileVersion)
    {
      fprintf(stderr, "WARNING!!\n\nThis configfile (%s) is NOT suitable for this program.\n\n", filename);
      exit(1);
    }
    
  /* max memory [bytes] */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &maxAllowedMemory);
  /* range checking: must not be negative */
  if (maxAllowedMemory < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: max. allowed memory is less than 0!\n");
    }

  /* maximum Sigma [S/m] */
  ReadNextLine(configFile, string);
  sscanf(string, "%f\n", &a);
  maxSigma = (PRECISION) a;
  /* range checking: must not be negative */
  if (maxSigma < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: max. conductivity is less than 0!\n");
    }

  /* maximum Stretching [1] */
  ReadNextLine(configFile, string);
  sscanf(string, "%f\n", &a);
  maxStretching = (PRECISION) a;
  /* range checking: must not be negative */
  if (maxStretching < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: max. stretching is less than 0!\n");
    }
  
  /* steepness of stretch parameter */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &stretchSteepness);
  /* range checking: must not be zero */
  if (fabs(stretchSteepness) < LIMIT)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: stretching-exponent must not be 0!\n");
    }
  
  /* length of simulation */
  ReadNextLine(configFile, string);
  sscanf(string, "%g\n", &a);
  simulationLength = (PRECISION) a;
  /* range checking: must not be negative */
  if (simulationLength < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: length of simulation is less than 0!\n");
    }
  
  /* number of cells in x,y,z-direction */
  ReadNextLine(configFile, string);
  sscanf(string, "%d %d %d\n", &nx, &ny, &nz);
  /* range checking: must not be negative */
  if (nx < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in x-direction is less than 0!\n");
    }
  if (ny < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in y-direction is less than 0!\n");
    }
  if (nz < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in z-direction is less than 0!\n");
    }

  /* cell dimensions (x,y,z) */
  ReadNextLine(configFile, string);
  sscanf(string, "%g %g %g\n", &a, &b, &c);
  dx = (PRECISION) a;
  dy = (PRECISION) b;
  dz = (PRECISION) c;
  /* range checking: must not be negative */
  if (dx < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in x-direction is less than 0!\n");
    }
  if (dy < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in y-direction is less than 0!\n");
    }
  if (dz < 0.0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of cells in z-direction is less than 0!\n");
    }
  dim_X = nx * dx;
  dim_Y = ny * dy;
  dim_Z = nz * dz;

  /* time */
  ReadNextLine(configFile, string);
  sscanf(string, "%g\n", &a);
  dt = (PRECISION) a;
  /* range checking: must not be 0 */
  if (fabs(dt) < LIMIT)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: timestep is 0!\n");
    }
  
 /* when to call the output routine for 'output.dat' */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &timeToPlotReceiver);
  /* range checking: must not be negative */
  if (timeToPlotReceiver < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: receiver plot divider is less than 0!\n");
    }
  
  /* when to call the output routine for a snapshot */
  ReadNextLine(configFile, string);
  sscanf(string, "%d %d %d %d\n", &timeToPlotSnapshot, &snapshotPlane, &snapshotPosition, &snapshotMode);
  /* range checking: must not be negative */
  if (timeToPlotSnapshot < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: snapshot plot divider is less than 0!\n");
    }
  if ((snapshotPlane < 1) || (snapshotPlane > 3))
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: snapshot plane number is not valid!\n");
    }
  switch(snapshotPlane)
    {
    case 1:
      if ((snapshotPosition < 0) || (snapshotPosition >= nx))
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: snapshot plane position is out of range!\n");
	}
      break;
    case 2:
      if ((snapshotPosition < 0) || (snapshotPosition >= ny))
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: snapshot plane position is out of range!\n");
	}
      break;
    case 3:
      if ((snapshotPosition < 0) || (snapshotPosition >= nz))
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: snapshot plane position is out of range!\n");
	}
      break;
    default:;
    }
  if ((snapshotMode < 0) || (snapshotMode > 7))
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: snapshot mode is not valid!\n");
    }

  /* width of pmls */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &pmlWidth);
  /* range checking: must not be negative */
  if (pmlWidth < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number of PM-Layers is less than 0!\n");
    }
  
  /* materials, first maximum number of materials */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &maxMaterials);
  /* range checking: must not be negative */
  if (maxMaterials < 0)
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: number different materials is less than 0!\n");
    }
  
  /* allocate that much memory */

  materialConstants = (PRECISION **) malloc((maxMaterials+2) * sizeof(PRECISION *));
  if (materialConstants == NULL)
    {
      fprintf(stderr, "I cannot allocate memory for the constants.\n");
      exit(1);
    }
  for (i=0; i<maxMaterials+2; i++)
    {
      materialConstants[i] = (PRECISION *) malloc(9 * sizeof(PRECISION));
      if (materialConstants[i] == NULL)
	{
	  fprintf(stderr, "I cannot allocate memory for the constants.\n");
	  exit(1);
	}
      for(j=0; j<9; j++)
	materialConstants[i][j] = 0.0;
    }
  
 
  /* read that number of materials */

  maxEpsilon = 1;
  
  for(i=0; i<maxMaterials; i++)
    {
      ReadNextLine(configFile, string);
      sscanf(string, "%g %g\n", &a, &b);
      materialConstants[i+2][0] = (PRECISION) a;
      materialConstants[i+2][1] = (PRECISION) b;
      if (maxEpsilon < a) maxEpsilon = a;	
      /* range checking: must not be negative */
      if (a < 0.0)
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: epsilon for material %d is less than 0!\n", i);
	}
      if (b < 0.0)
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: sigma for material %d is less than 0!\n", i);
	}
    }

  /* allocate memory for materials */
  material = MyAllocShortInt(nx, ny, nz,&storage);
  
  /* number of boxes  */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &numberOfBoxes);

  for(h=0; h<numberOfBoxes; h++)
    {
      ReadNextLine(configFile, string);
      sscanf(string, "%d %d %d %d %d %d %d\n", &boxStartX, &boxStartY, &boxStartZ, &boxEndX, &boxEndY, &boxEndZ, &currentMaterial);
      /* range checking: must be within simspace */
      if (\
	  (boxStartX >= 0) && (boxStartX < nx) && \
	  (boxStartY >= 0) && (boxStartY < ny) && \
	  (boxStartZ >= 0) && (boxStartZ < nz) && \
	  (boxEndX >= 0) && (boxEndX < nx) && \
	  (boxEndY >= 0) && (boxEndY < ny) && \
	  (boxEndZ >= 0) && (boxEndZ < nz) && \