/* Define the precision */
#include "stdafx.h"
#include <iostream>
#define PRECISION double
#include "radarFDTD.h"

CradarFDTD::CradarFDTD()
{
	allocatedBytes = 0.0;        
	iteration = 0;     
	stimulus = 0.0;              
	currentSimulationTime = 0.0; 
	totalSimulationTime = 0.0; 
	currentConfigFileVersion = 3; 
	storage = 0;
}

CradarFDTD::~CradarFDTD()
{
	
}

PRECISION*** CradarFDTD::MyAllocPrecision(int dimensionX, int dimensionY, int dimensionZ, PRECISION value, int *storage)
{
	PRECISION ***returnArray;
	int i,j,k;
  
	*storage += sizeof(PRECISION) * dimensionX * dimensionY * dimensionZ;
	*storage += sizeof(PRECISION *) * dimensionX * dimensionY;
	*storage += sizeof(PRECISION **) * dimensionX;
  

	returnArray = (PRECISION ***) malloc(dimensionX * sizeof(PRECISION **));
	if (returnArray == NULL)
	{
		fprintf(stdout, "Sorry, but you don't have enough memory installed!");
		exit(1);
	}
  
	 for(i=0; i<dimensionX; i++)
	{
		returnArray[i] = (PRECISION **) malloc(dimensionY * sizeof(PRECISION *));
		if (returnArray[i] == NULL)
		{
		    fprintf(stdout, "Sorry, but you don't have enough memory installed!");
		    exit(1);
		}
		    for(j=0; j<dimensionY; j++)
		{
		    returnArray[i][j] = (PRECISION *) malloc(dimensionZ * sizeof(PRECISION));
		    if (returnArray[i][j] == NULL)
			{
			    fprintf(stdout, "Sorry, but you don't have enough memory installed!");
			    exit(1);
			}
			for(k=0; k<dimensionZ; k++)
			{
				returnArray[i][j][k] = value;
			}
		}
	}

	return returnArray;

}

unsigned short int*** CradarFDTD::MyAllocShortInt(int dimensionX, int dimensionY, int dimensionZ, int *storage)
{
	unsigned short int ***returnArray;
    int i,j,k;

    *storage += sizeof(unsigned short int) * dimensionX * dimensionY * dimensionZ;
    *storage += sizeof(unsigned short int *) * dimensionX * dimensionY;
    *storage += sizeof(unsigned short int **) * dimensionX;

    returnArray = (unsigned short int ***) malloc(dimensionX * sizeof(unsigned short int **));
    if (returnArray == NULL)
    {
        fprintf(stdout, "Sorry, but you don't have enough memory installed!");
        exit(1);
    }
  
    for(i=0; i<dimensionX; i++)
    {
        returnArray[i] = (unsigned short int **) malloc(dimensionY * sizeof(unsigned short int *));
        if (returnArray[i] == NULL)
		{
			fprintf(stdout, "Sorry, but you don't have enough memory installed!");
			exit(1);
		}

		for(j=0; j<dimensionY; j++)
		{
			returnArray[i][j] = (unsigned short int *) malloc(dimensionZ * sizeof(unsigned short int));
			if (returnArray[i][j] == NULL)
			{
				fprintf(stdout, "Sorry, but you don't have enough memory installed!");
				exit(1);
			}

			for(k=0; k<dimensionZ; k++)
			{
			   returnArray[i][j][k] = 0.0;
			}
		}
    }

  return returnArray;
  
}

int CradarFDTD::MyFreePrecision(PRECISION ***pointer, int x, int y, int z)
{
	int i,j;
  
	for(i=0; i<x; i++)
		for(j=0; j<y; j++)
			free(pointer[i][j]);
  
	 for(i=0; i<x; i++)
		free(pointer[i]);
  
	free(pointer);
	return 0;
}

int CradarFDTD::MyFreeShortInt(unsigned short int ***pointer, int x, int y, int z)
{
    int i,j;
  
    for(i=0; i<x; i++)
        for(j=0; j<y; j++)
			free(pointer[i][j]);
  
	for(i=0; i<x; i++)
		free(pointer[i]);
  
	free(pointer);
	return 0;
	
}

void CradarFDTD::PrintCorrectOrder(PRECISION number, FILE *output)
{
    int order = (int) log(number) / log(10) - 0.5;
  
	switch(order)
	{
		case -15: fprintf(output, "%5.2f f", number*1.0e+15);break;
		case -14: fprintf(output, "%5.2f f", number*1.0e+15);break;
		case -13: fprintf(output, "%5.2f f", number*1.0e+15);break;
		case -12: fprintf(output, "%5.2f p", number*1.0e+12);break;
		case -11: fprintf(output, "%5.2f p", number*1.0e+12);break;
		case -10: fprintf(output, "%5.2f p", number*1.0e+12);break;
		case -9: fprintf(output, "%5.2f n", number*1.0e+9);break;
		case -8: fprintf(output, "%5.2f n", number*1.0e+9);break;
		case -7: fprintf(output, "%5.2f n", number*1.0e+9);break;
		case -6: fprintf(output, "%5.2f t", number*1.0e+6);break;
		case -5: fprintf(output, "%5.2f t", number*1.0e+6);break;
		case -4: fprintf(output, "%5.2f t", number*1.0e+6);break;
		case -3: fprintf(output, "%5.2f m", number*1.0e+3);break;
		case -2: fprintf(output, "%5.2f m", number*1.0e+3);break;
		case -1: fprintf(output, "%5.2f m", number*1.0e+3);break;
		case 0: fprintf(output, "%5.2f", number);break;
		case 1: fprintf(output, "%5.2f", number);break;
		case 2: fprintf(output, "%5.2f", number);break;
		case 3: fprintf(output, "%5.2f k", number*1.0e-3);break;
		case 4: fprintf(output, "%5.2f k", number*1.0e-3);break;
		case 5: fprintf(output, "%5.2f k", number*1.0e-3);break;
		case 6: fprintf(output, "%5.2f M", number*1.0e-6);break;
		case 7: fprintf(output, "%5.2f M", number*1.0e-6);break;
		case 8: fprintf(output, "%5.2f M", number*1.0e-6);break;
		case 9: fprintf(output, "%5.2f G", number*1.0e-9);break;
		case 10: fprintf(output, "%5.2f G", number*1.0e-9);break;
		case 11: fprintf(output, "%5.2f G", number*1.0e-9);break;
		case 12: fprintf(output, "%5.2f T", number*1.0e-12);break;
		case 13: fprintf(output, "%5.2f T", number*1.0e-12);break;
		case 14: fprintf(output, "%5.2f T", number*1.0e-12);break;
		default: fprintf(output, "%5.2e", number);
  }
	
}


void CradarFDTD::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);
}

void CradarFDTD::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;