readconfig.c

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

	  (boxStartX <= boxEndX) && (boxStartY <= boxEndY) && (boxStartZ <= boxEndZ))
	{
	  for(i=boxStartX; i<=boxEndX; i++)
	    for(j=boxStartY; j<=boxEndY; j++)
	      for(k=boxStartZ; k<=boxEndZ; k++)
		material[i][j][k] = currentMaterial;
	}
      else
	{ 
	  warningStop = 1;
	  fprintf(stderr, "WARNING: error in material box #%d coordinated\n", h);
	}
    }
  
 /* how many transmitter boxes are used? */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &numberOfTransmitterBoxes);
  /* range checking: must not be negative */
      if (numberOfTransmitterBoxes < 0)
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: number of receiver boxes is less than 0!\n");
	}

  /* get some temporary memory */

  tmpTransmitter = (int **) malloc(numberOfTransmitterBoxes * sizeof(int *));
  if (tmpTransmitter == NULL)
    {
      fprintf(stderr, "Not enough memory available!\n");
      exit(1);
    }
  for(h=0; h<numberOfTransmitterBoxes; h++)
    {
      tmpTransmitter[h] = (int *) malloc(8 * sizeof(int));
      if (tmpTransmitter[h] == NULL)
	{
	  fprintf(stderr, "Not enough memory available!\n");
	  exit(1);
	}
    }
 
  maxTransmitter = 0;

  for(h=0; h<numberOfTransmitterBoxes; h++)
    {
      ReadNextLine(configFile, string);
      sscanf(string, "%d %d %d %d %d %d %d\n", &tmpTransmitter[h][0], &tmpTransmitter[h][1], &tmpTransmitter[h][2], &tmpTransmitter[h][3], &tmpTransmitter[h][4], &tmpTransmitter[h][5], &tmpTransmitter[h][6]);
    }
   
  for(h=0; h<numberOfTransmitterBoxes; h++)
    {
      if (\
	  (tmpTransmitter[h][0] >= 0) && (tmpTransmitter[h][0] < nx) && \
	  (tmpTransmitter[h][1] >= 0) && (tmpTransmitter[h][1] < ny) && \
	  (tmpTransmitter[h][2] >= 0) && (tmpTransmitter[h][2] < nz) && \
	  (tmpTransmitter[h][3] >= 0) && (tmpTransmitter[h][3] < nx) && \
	  (tmpTransmitter[h][4] >= 0) && (tmpTransmitter[h][4] < ny) && \
	  (tmpTransmitter[h][5] >= 0) && (tmpTransmitter[h][5] < nz) && \
	  (tmpTransmitter[h][0] <= tmpTransmitter[h][3]) && (tmpTransmitter[h][1] <= tmpTransmitter[h][4]) && (tmpTransmitter[h][2] <= tmpTransmitter[h][5]) && \
	  (tmpTransmitter[h][6] > 0) \
	  )
	{
	  for(i=tmpTransmitter[h][0]; i<=tmpTransmitter[h][3]; i+=tmpTransmitter[h][6])
	    for(j=tmpTransmitter[h][1]; j<=tmpTransmitter[h][4]; j+=tmpTransmitter[h][6])
	      for(k=tmpTransmitter[h][2]; k<=tmpTransmitter[h][5]; k+=tmpTransmitter[h][6])
		maxTransmitter++;
	  
	  transmitter = (int **) malloc(maxTransmitter * sizeof(int *));
	  if (transmitter == NULL)
	    {
	      fprintf(stderr, "Not enough memory available!\n");
	      exit(1);
	    }
	  for(i=0; i<maxTransmitter; i++)
	    {
	      transmitter[i] = (int *) malloc(3 * sizeof(int));
	      if (transmitter[i] == NULL)
		{
		  fprintf(stderr, "Not enough memory available!\n");
		  exit(1);
		}
	    }
	  
	  tmpCounter = 0;
	  
	  for(i=tmpTransmitter[h][0]; i<=tmpTransmitter[h][3]; i+=tmpTransmitter[h][6])
	    for(j=tmpTransmitter[h][1]; j<=tmpTransmitter[h][4]; j+=tmpTransmitter[h][6])
	      for(k=tmpTransmitter[h][2]; k<=tmpTransmitter[h][5]; k+=tmpTransmitter[h][6])
		{
		  transmitter[tmpCounter][0] = i;
		  transmitter[tmpCounter][1] = j;
		  transmitter[tmpCounter][2] = k;
		  tmpCounter++;
		}
	}
      else
	{ 
	  warningStop = 1;
	  fprintf(stderr, "WARNING: error in transmitter box #%d coordinates \n", h);
	}
    }
  
  for(h=0; h<numberOfTransmitterBoxes; h++)
    free(tmpTransmitter[h]);
  free(tmpTransmitter);
  
  /* which transmitter mode is used? */
  ReadNextLine(configFile, string);
  sscanf(string, "%d %f %f %f %f %d\n", &transmitterMode, &a, &b, &c, &d, &transmitterStimulusComponent);
  transmitterAmplitude          = (PRECISION) a;
  transmitterCharacteristicTime = (PRECISION) b;
  transmitterFrequency          = (PRECISION) c;
  transmitterParameter          = (PRECISION) d;

  /* range checking: ONLY THE LAST COMPONENT IS CURRENTLY CHECKED!!! */
  if ((transmitterStimulusComponent < 0) || (transmitterStimulusComponent > 2))
    {
      warningStop = 1;
      fprintf(stderr, "WARNING: stimulus component of the transmitter is out of range\n");
    }
  
  /* how many receiving boxes are used? */
  ReadNextLine(configFile, string);
  sscanf(string, "%d\n", &numberOfReceiverBoxes);
  /* range checking: must not be negative */
      if (numberOfReceiverBoxes < 0)
	{
	  warningStop = 1;
	  fprintf(stderr, "WARNING: number of receiver boxes is less than 0!\n");
	}

  /* get some temporary memory */

  tmpReceiver = (int **) malloc(numberOfReceiverBoxes * sizeof(int *));
  if (tmpReceiver == NULL)
    {
      fprintf(stderr, "Not enough memory available!\n");
      exit(1);
    }
  for(h=0; h<numberOfReceiverBoxes; h++)
    {
      tmpReceiver[h] = (int *) malloc(8 * sizeof(int));
      if (tmpReceiver[h] == NULL)
	{
	  fprintf(stderr, "Not enough memory available!\n");
	  exit(1);
	}
    }
 
  maxReceiver = 0;

  for(h=0; h<numberOfReceiverBoxes; h++)
    {
      ReadNextLine(configFile, string);
      sscanf(string, "%d %d %d %d %d %d %d %d\n", &tmpReceiver[h][0], &tmpReceiver[h][1], &tmpReceiver[h][2], &tmpReceiver[h][3], &tmpReceiver[h][4], &tmpReceiver[h][5], &tmpReceiver[h][6], &tmpReceiver[h][7]);
    }
   
  for(h=0; h<numberOfReceiverBoxes; h++)
    {
      if (\
	  (tmpReceiver[h][0] >= 0) && (tmpReceiver[h][0] < nx) && \
	  (tmpReceiver[h][1] >= 0) && (tmpReceiver[h][1] < ny) && \
	  (tmpReceiver[h][2] >= 0) && (tmpReceiver[h][2] < nz) && \
	  (tmpReceiver[h][3] >= 0) && (tmpReceiver[h][3] < nx) && \
	  (tmpReceiver[h][4] >= 0) && (tmpReceiver[h][4] < ny) && \
	  (tmpReceiver[h][5] >= 0) && (tmpReceiver[h][5] < nz) && \
	  (tmpReceiver[h][0] <= tmpReceiver[h][3]) && (tmpReceiver[h][1] <= tmpReceiver[h][4]) && (tmpReceiver[h][2] <= tmpReceiver[h][5]) && \
	  (tmpReceiver[h][6] >= 0) && (tmpReceiver[h][6] <= 7) && \
	  (tmpReceiver[h][7] > 0) \
	  )
	{
	  for(i=tmpReceiver[h][0]; i<=tmpReceiver[h][3]; i+=tmpReceiver[h][7])
	    for(j=tmpReceiver[h][1]; j<=tmpReceiver[h][4]; j+=tmpReceiver[h][7])
	      for(k=tmpReceiver[h][2]; k<=tmpReceiver[h][5]; k+=tmpReceiver[h][7])
		maxReceiver++;
	  
	  receiver = (int **) malloc(maxReceiver * sizeof(int *));
	  if (receiver == NULL)
	    {
	      fprintf(stderr, "Not enough memory available!\n");
	      exit(1);
	    }
	  for(i=0; i<maxReceiver; i++)
	    {
	      receiver[i] = (int *) malloc(4 * sizeof(int));
	      if (receiver[i] == NULL)
		{
		  fprintf(stderr, "Not enough memory available!\n");
		  exit(1);
		}
	    }
	  
	  tmpCounter = 0;
	  
	  for(i=tmpReceiver[h][0]; i<=tmpReceiver[h][3]; i+=tmpReceiver[h][7])
	    for(j=tmpReceiver[h][1]; j<=tmpReceiver[h][4]; j+=tmpReceiver[h][7])
	      for(k=tmpReceiver[h][2]; k<=tmpReceiver[h][5]; k+=tmpReceiver[h][7])
		{
		  receiver[tmpCounter][0] = i;
		  receiver[tmpCounter][1] = j;
		  receiver[tmpCounter][2] = k;
		  receiver[tmpCounter][3] = tmpReceiver[h][6]; /* mode */
		  tmpCounter++;
		}
	}
      else
	{ 
	  warningStop = 1;
	  fprintf(stderr, "WARNING: error in receiver box #%d coordinates \n", h);
	}
    }
  
  for(h=0; h<numberOfReceiverBoxes; h++)
    free(tmpReceiver[h]);
  free(tmpReceiver);
  
  /* NOW CHECK IF THE PROGRAM CAN CONTINUE BEYOND THIS POINT */
  
  if (warningStop == 0)
    {
      fprintf(stdout, "Settings seem to be alright. Starting the computation.\n\n");
    }
  else
    {
      fprintf(stderr, "Since there is something wrong, please correct your configfile (%s/%s) and restart the program.\n", directory, filename);
    }
	  
 
  /*********** Now the computations ***********/

  /* Courant-Criterion */

  tmpdt =  1 / (LIGHT_SPEED * sqrt(1/(dx*dx) + 1/(dy*dy) + 1/(dz*dz)));
  
  if (dt < 0) dt = -tmpdt * (PRECISION) dt;
  
  /* material constants */

  /*  set values for vacuum */

  materialConstants[0][DTMUDX]   = dt/(MU_0*dx);
  materialConstants[0][DTEPSDX]  = dt/(EPSILON_0*dx);
  materialConstants[0][DTMUDY]   = dt/(MU_0*dy);
  materialConstants[0][DTEPSDY]  = dt/(EPSILON_0*dy);
  materialConstants[0][DTMUDZ]   = dt/(MU_0*dz);
  materialConstants[0][DTEPSDZ]  = dt/(EPSILON_0*dz);
  materialConstants[0][DTSIGEPS] = 1.0;
  materialConstants[0][EPSILON]  = 1.0;
  materialConstants[0][SIGMA]    = 0.0;

  /* for ideal conductors (this will _not_ be used, but just in case) */

  materialConstants[1][DTMUDX]   = dt/(MU_0*dx);
  materialConstants[1][DTEPSDX]  = 0.0;
  materialConstants[1][DTMUDY]   = dt/(MU_0*dy);
  materialConstants[1][DTEPSDY]  = 0.0;
  materialConstants[1][DTMUDZ]   = dt/(MU_0*dz);
  materialConstants[1][DTEPSDZ]  = 0.0;
  materialConstants[1][DTSIGEPS] = 0.0;
  materialConstants[1][EPSILON]  = 0.0;
  materialConstants[1][SIGMA]    = 0.0;


  /*and for the various other materials */

  for(i=0; i<maxMaterials; i++)
    {
      tmpEpsilon = materialConstants[i+2][0];
      tmpSigma   = materialConstants[i+2][1];
      materialConstants[i+2][DTMUDX]   = dt/(MU_0 * dx);
      materialConstants[i+2][DTEPSDX]  = dt / (dx * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
      materialConstants[i+2][DTMUDY]   = dt/(MU_0 * dy);
      materialConstants[i+2][DTEPSDY]  = dt / (dy * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
      materialConstants[i+2][DTMUDZ]   = dt/(MU_0 * dz);
      materialConstants[i+2][DTEPSDZ]  = dt / (dz * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
      materialConstants[i+2][DTSIGEPS] = (1 - tmpSigma * dt / (2.0 * EPSILON_0 * tmpEpsilon)) / (1 + tmpSigma * dt / (2.0 * EPSILON_0 * tmpEpsilon));
      materialConstants[i+2][EPSILON]  = tmpEpsilon;
      materialConstants[i+2][SIGMA]    = tmpSigma;
    }

  /* Calculate the maximum space - step size */
  
  maxStep = dx;
  if (maxStep < dy) maxStep = dy;
  if (maxStep < dz) maxStep = dz;
  
   /* Calculate the number of iteration steps needed */
  
  maxIteration = (int) (1.0 + simulationLength / dt);

  /* Calculate theoretical reflection from the absorber */
 
  tmpReflection = 0.0;
  for(i=0; i<pmlWidth; i++)
    tmpReflection += ReturnStretching((i+1) * maxStep, pmlWidth * maxStep, maxStretching, stretchSteepness) * ReturnConductivity((i+1) * maxStep, pmlWidth * maxStep, maxSigma);

  theoreticalReflection = exp(-2.0 / (maxEpsilon * EPSILON_0 * LIGHT_SPEED) * tmpReflection * maxStep);

  minWavelength = (1.0 + maxStretching) * 3.0 * maxStep;

  /***************** Output section ******************/

  fprintf(stdout, "\n******************************************************************\n");
  fprintf(stdout, VERSION_INFO);
  fprintf(stdout, "\n******************************************************************\n");
  fprintf(stdout, "\nConfiguration file: %s\n\n", filename);
  fprintf(stdout, "Simulation steps: %d (%.4gns), dt = %.4gns, output every %d steps\n\n", maxIteration, 1e9 * (float) maxIteration * dt, 1e9*dt,timeToPlotReceiver);
  fprintf(stdout, "Simulation space: %.4f x %.4f x %.4f m砛n", dx*nx, dy*ny, dz*nz);
  fprintf(stdout, "Number of cells:\t%d\t%d\t%d\n", nx,ny,nz);
  fprintf(stdout, "Cellsizes:\t\t%.4gm\t%.4gm\t%.4gm\n\n", dx,dy,dz);
  fprintf(stdout, "Number of absorbing boundary layers: %d\nMaximum conductivity at edge: %.4g, maximum Stretching: %.4g\n\n", pmlWidth, (float) maxSigma, (float) maxStretching);
  
  fprintf(stdout, "\nThe derived theoretical reflection  (after Fang and Wu) is: %e,\nand the minimum absorber wavelength is", theoreticalReflection);
  PrintCorrectOrder(minWavelength, stdout);
  fprintf(stdout, "m (freq: ");
  PrintCorrectOrder(LIGHT_SPEED / minWavelength, stdout);
  fprintf(stdout, "Hz)\n\n");

  fprintf(stdout, "Number of materials: %d\n", (maxMaterials+2));
  
  for(i=0; i<maxMaterials+2; i++)
    if (i != 1)
      fprintf(stdout, "%d:\tEpsilon: %.6g\tSigma: %.6g\n", i,  materialConstants[i][7], materialConstants[i][8]);
  else
    fprintf(stdout, "1:\tperfect conductor\n");
  
  fprintf(stdout, "\nNumber of tranmsitters: %d\n", maxTransmitter);
  fprintf(stdout, "\nNumber of receivers: %d\n\n", maxReceiver);
  fprintf(stdout, "Remarks:\nCourant-criterion ");
  
  if (dt < tmpdt) 
    fprintf(stdout, "OK.\n");
  else 
    fprintf(stdout, "NOT OK! (dt / good) < %6.2g\n", (tmpdt/dt));

  fprintf(stdout, "Cellsize: ");
  
  if (maxStep < LIGHT_SPEED/(10.0 * transmitterFrequency * sqrt(maxEpsilon)))
    fprintf(stdout, "OK\n");
  else
    {
      fprintf(stdout, "NOT OK! Recommended value is at least:");
      PrintCorrectOrder(LIGHT_SPEED/(10 * transmitterFrequency * sqrt(maxEpsilon)), stdout);
      fprintf(stdout, "m\n\n");
    }
  
  /* close file */

  fclose(configFile);
}