radarfdtd.c

The purpose of this program is to enable building a config file to the radarFDTD package using a mo

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


/* Very first implementation of a FDTD-simulation code 
   written by Carsten Aulbert ([e-mail])
   (using the sources from ToyFDTD.c as a 'manual')
   HISTORY
   V1.0 - 1999/08/30 - coding started with template ToyFDTD
   V2.0 - 1999/09/04 - different materials are possible
   V3.0 - 1999/11/19 - mainly code cleaning and PML boundaries
   V4.0 - 2000/01/03 - implementing configfile
   V4.1 - 2000/01/04 - mainly bug tracking and hunting them down
   V4.2 - 2000/01/28 - implementing GPML, bugfixing
   V5.0pre - 2000/02/23 - tidy up code, quicken code
   V5.0 - 2000/03/10 - changing input scheme, correct GPMLs 
   V6.0 - 2000/04/09 - changed last bits, changing input scheme again 
        - 2000/05/04 - made final changes, release time :-)
   V7.0 - 2004/09/24 - Ron Daisy <[e-mail]> sent a patch for ReadConfig.c, new
                       major release since he also addded a Qt-GUI 
*/


#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <float.h>
#include <time.h> 

#define VERSION_INFO "RadarFDTD Version 7.0 by Carsten Aulbert ([e-mail])"

/* Define the precision */
#define PRECISION float

/* Define some common acronyms for the materialConstants array */
#define DTMUDX 0
#define DTEPSDX 1
#define DTMUDY 2
#define DTEPSDY 3
#define DTMUDZ 4
#define DTEPSDZ 5
#define DTSIGEPS 6
#define EPSILON 7
#define SIGMA 8

/* define pi if not defined */
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* under which number a value is considered zero */
#define LIMIT 1e-50

/* Now there are some physical constants, don't fiddle around with them too mucyh ;) */
#define LIGHT_SPEED 299792458.0
#define LIGHT_SPEED_SQUARED 89875517873681764.0
#define MU_0 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6
#define EPSILON_0 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12

/* Constant for plotting, every TIME_FOR_PLOT iteration there will be an output */

#ifdef DEBUG
     FILE *loggingFile;
#endif

/* Usually I do not like global variables, but if I had not chosen to use these, the function prototypes would have been quite long! */

/* global variables go here */

int i,j,k;                          /* indices of 3D array of cells */
int nx, ny, nz;                     /* number of cells along axes */
int posX1, posY1, posX2, posY2;
PRECISION allocatedBytes = 0.0;        /* tracker of how much memory is used */
int iteration = 0;                  /* iteration counter */

PRECISION mytmp;
PRECISION stimulus = 0.0;              /* stimulus value */
PRECISION currentSimulationTime = 0.0; /* current simulation time [s] */
PRECISION totalSimulationTime = 0.0;   /* total simulation time [s] */
PRECISION omega;                       /* angular frequency for stimulus [radians/second] */
PRECISION lambda;                      /* wavelength of stimulus [m] */
PRECISION dx, dy, dz;                  /* space differentials, i.e. size of a cell [m]*/
PRECISION dt,tmp;                          /* time differential [s] */
PRECISION ***ex, ***ey, ***ez;         /* pointers to 3D-arrays */
PRECISION ***hx, ***hy, ***hz;         /* pointers to 3D-arrays */

unsigned short int ***material;        /* layout of material */
PRECISION **materialConstants;         /* constants for various materials */
PRECISION tmpout;                      /* tmp output */

long int programStartTime;             /* number of seconds since 01/01/1970*/
long int tmpTime, tmpTimeRemaining;
float tmpTimePerStep;
int tmpRemainingDays,tmpRemainingHours,tmpRemainingMinutes,tmpRemainingSeconds;
char timeString[1024];
char addString[1024];
 
PRECISION ***pmlEXYBottomYZ, ***pmlEXZBottomYZ; /* pml values */
PRECISION ***pmlEYXBottomYZ, ***pmlEYZBottomYZ; /* pml values */
PRECISION ***pmlEZXBottomYZ, ***pmlEZYBottomYZ; /* pml values */
PRECISION ***pmlHXYBottomYZ, ***pmlHXZBottomYZ; /* pml values */
PRECISION ***pmlHYXBottomYZ, ***pmlHYZBottomYZ; /* pml values */
PRECISION ***pmlHZXBottomYZ, ***pmlHZYBottomYZ; /* pml values */

PRECISION ***pmlEXYBottomXY, ***pmlEXZBottomXY; /* pml values */
PRECISION ***pmlEYXBottomXY, ***pmlEYZBottomXY; /* pml values */
PRECISION ***pmlEZXBottomXY, ***pmlEZYBottomXY; /* pml values */
PRECISION ***pmlHXYBottomXY, ***pmlHXZBottomXY; /* pml values */
PRECISION ***pmlHYXBottomXY, ***pmlHYZBottomXY; /* pml values */
PRECISION ***pmlHZXBottomXY, ***pmlHZYBottomXY; /* pml values */

PRECISION ***pmlEXYBottomXZ, ***pmlEXZBottomXZ; /* pml values */
PRECISION ***pmlEYXBottomXZ, ***pmlEYZBottomXZ; /* pml values */
PRECISION ***pmlEZXBottomXZ, ***pmlEZYBottomXZ; /* pml values */
PRECISION ***pmlHXYBottomXZ, ***pmlHXZBottomXZ; /* pml values */
PRECISION ***pmlHYXBottomXZ, ***pmlHYZBottomXZ; /* pml values */
PRECISION ***pmlHZXBottomXZ, ***pmlHZYBottomXZ; /* pml values */

PRECISION ***pmlEXYTopYZ, ***pmlEXZTopYZ; /* pml values */
PRECISION ***pmlEYXTopYZ, ***pmlEYZTopYZ; /* pml values */
PRECISION ***pmlEZXTopYZ, ***pmlEZYTopYZ; /* pml values */
PRECISION ***pmlHXYTopYZ, ***pmlHXZTopYZ; /* pml values */
PRECISION ***pmlHYXTopYZ, ***pmlHYZTopYZ; /* pml values */
PRECISION ***pmlHZXTopYZ, ***pmlHZYTopYZ; /* pml values */

PRECISION ***pmlEXYTopXY, ***pmlEXZTopXY; /* pml values */
PRECISION ***pmlEYXTopXY, ***pmlEYZTopXY; /* pml values */
PRECISION ***pmlEZXTopXY, ***pmlEZYTopXY; /* pml values */
PRECISION ***pmlHXYTopXY, ***pmlHXZTopXY; /* pml values */
PRECISION ***pmlHYXTopXY, ***pmlHYZTopXY; /* pml values */
PRECISION ***pmlHZXTopXY, ***pmlHZYTopXY; /* pml values */

PRECISION ***pmlEXYTopXZ, ***pmlEXZTopXZ; /* pml values */
PRECISION ***pmlEYXTopXZ, ***pmlEYZTopXZ; /* pml values */
PRECISION ***pmlEZXTopXZ, ***pmlEZYTopXZ; /* pml values */
PRECISION ***pmlHXYTopXZ, ***pmlHXZTopXZ; /* pml values */
PRECISION ***pmlHYXTopXZ, ***pmlHYZTopXZ; /* pml values */
PRECISION ***pmlHZXTopXZ, ***pmlHZYTopXZ; /* pml values */

PRECISION ***pmlComplexEXYBottomYZ, ***pmlComplexEXZBottomYZ; /* pml values */
PRECISION ***pmlComplexEYXBottomYZ, ***pmlComplexEYZBottomYZ; /* pml values */
PRECISION ***pmlComplexEZXBottomYZ, ***pmlComplexEZYBottomYZ; /* pml values */

PRECISION ***pmlComplexEXYBottomXY, ***pmlComplexEXZBottomXY; /* pml values */
PRECISION ***pmlComplexEYXBottomXY, ***pmlComplexEYZBottomXY; /* pml values */
PRECISION ***pmlComplexEZXBottomXY, ***pmlComplexEZYBottomXY; /* pml values */

PRECISION ***pmlComplexEXYBottomXZ, ***pmlComplexEXZBottomXZ; /* pml values */
PRECISION ***pmlComplexEYXBottomXZ, ***pmlComplexEYZBottomXZ; /* pml values */
PRECISION ***pmlComplexEZXBottomXZ, ***pmlComplexEZYBottomXZ; /* pml values */

PRECISION ***pmlComplexEXYTopYZ, ***pmlComplexEXZTopYZ; /* pml values */
PRECISION ***pmlComplexEYXTopYZ, ***pmlComplexEYZTopYZ; /* pml values */
PRECISION ***pmlComplexEZXTopYZ, ***pmlComplexEZYTopYZ; /* pml values */

PRECISION ***pmlComplexEXYTopXY, ***pmlComplexEXZTopXY; /* pml values */
PRECISION ***pmlComplexEYXTopXY, ***pmlComplexEYZTopXY; /* pml values */
PRECISION ***pmlComplexEZXTopXY, ***pmlComplexEZYTopXY; /* pml values */

PRECISION ***pmlComplexEXYTopXZ, ***pmlComplexEXZTopXZ; /* pml values */
PRECISION ***pmlComplexEYXTopXZ, ***pmlComplexEYZTopXZ; /* pml values */
PRECISION ***pmlComplexEZXTopXZ, ***pmlComplexEZYTopXZ; /* pml values */

PRECISION ***pmlSigmaXBottomYZ, ***pmlSigmaXStarBottomYZ;
PRECISION ***pmlSigmaXBottomXY, ***pmlSigmaXStarBottomXY;
PRECISION ***pmlSigmaXBottomXZ, ***pmlSigmaXStarBottomXZ;
PRECISION ***pmlSigmaXTopYZ, ***pmlSigmaXStarTopYZ;
PRECISION ***pmlSigmaXTopXY, ***pmlSigmaXStarTopXY;
PRECISION ***pmlSigmaXTopXZ, ***pmlSigmaXStarTopXZ;

PRECISION ***pmlSigmaYBottomYZ, ***pmlSigmaYStarBottomYZ;
PRECISION ***pmlSigmaYBottomXY, ***pmlSigmaYStarBottomXY;
PRECISION ***pmlSigmaYBottomXZ, ***pmlSigmaYStarBottomXZ;
PRECISION ***pmlSigmaYTopYZ, ***pmlSigmaYStarTopYZ;
PRECISION ***pmlSigmaYTopXY, ***pmlSigmaYStarTopXY;
PRECISION ***pmlSigmaYTopXZ, ***pmlSigmaYStarTopXZ;

PRECISION ***pmlSigmaZBottomYZ, ***pmlSigmaZStarBottomYZ;
PRECISION ***pmlSigmaZBottomXY, ***pmlSigmaZStarBottomXY;
PRECISION ***pmlSigmaZBottomXZ, ***pmlSigmaZStarBottomXZ;
PRECISION ***pmlSigmaZTopYZ, ***pmlSigmaZStarTopYZ;
PRECISION ***pmlSigmaZTopXY, ***pmlSigmaZStarTopXY;
PRECISION ***pmlSigmaZTopXZ, ***pmlSigmaZStarTopXZ;


/* stretching parameters */

PRECISION ***pmlSXBottomYZ, ***pmlSYBottomYZ, ***pmlSZBottomYZ;
PRECISION ***pmlSXTopYZ, ***pmlSYTopYZ, ***pmlSZTopYZ;
PRECISION ***pmlSXBottomXY, ***pmlSYBottomXY, ***pmlSZBottomXY;
PRECISION ***pmlSXTopXY, ***pmlSYTopXY, ***pmlSZTopXY;
PRECISION ***pmlSXBottomXZ, ***pmlSYBottomXZ, ***pmlSZBottomXZ;
PRECISION ***pmlSXTopXZ, ***pmlSYTopXZ, ***pmlSZTopXZ;

/* materials within pmls */

unsigned short int ***pmlMaterialBottomYZ, ***pmlMaterialTopYZ;
unsigned short int ***pmlMaterialBottomXY, ***pmlMaterialTopXY;
unsigned short int ***pmlMaterialBottomXZ, ***pmlMaterialTopXZ;


/* include module for updating components */
#include "MyAlloc.c"
#include "PMLFunctions.c"
#include "Order.c"
#include "ReadConfig.c"
#include "ComponentUpdate.c"
#include "ComputeBoundarySigma.c"
#include "Output.c"
#include "Whirl.c"

 /* Start main program here */

int main (int argc, char **argv) {

  /* variables for plotting */

  FILE *foutput;                         /* pointer to output file */
  char filename[1024];                   /* dummy filename */
  char directory[1024];                  /* work directory */
  char outputFilename[1024];             /* output filename for receiver data*/
  char stimulusFilename[1024];           /* output filename for wave data*/
  FILE *fstimulus;                       /* pointer to stimulus output file */
  long int allocatedBytesSimSpace, allocatedBytesPMLSpace;
  PRECISION maxAmplitude;                /* Maximum amplitude of the stimulus */
  
  
  if (argc < 2)
    {
      fprintf(stderr, "WARNING!!\n\nPlease ensure you provide a valid directory-file name as an argument.\n\n");
      exit(-1);
    }
  
 
  strcpy(directory, argv[1]);
  
  if (argc > 2)
    {
      strcpy(filename, argv[1]);
      strcat(filename, argv[2]);
    }
  else
    {
      fprintf(stdout, "There was no configname given, using default 磗im.cfg