/* this is the program for 2D ADI-FDTD method (TE wave) - a point source in a cavity
(dimensions of the cavity are defined with a and b) that is closed with PECs at all 
the four outer boundaries. The 1st and 2nd resonance frequencies of the cavity 
are 150GHz and 300 GHz*/
/*the final output of this program is the resonance frequency of the 2-D cavity*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <iostream>
typedef double *PTR;
typedef double **DPTR;
typedef double **TPTR;

int main ()
{ 
  unsigned int myt = time(NULL); 
// #define CFLN 5   /*you can change this number to 1,2,3,4,5,6,...*/
// #define imax 100 /*number of cells in x direction*/
// #define jmax 50 /*number of cells in y direction*/
// #define itmax 200000 /*total iterations*/
// #define snapStep 1000
  double eps0=8.854e-12;
  double xmu0,xx,yy;
//  double a=1.0e-3,b=0.5e-3; /*a is the the length of the cavity and b is the width*/
  double e,m,ttt,tt,t,t0,c0,pi,t2e,t2h;
  int i,j,i0,j0,it,k, nobsX, nobsY;
  int nstop;
  float a, b, fmax, sourceX, sourceY, obsX, obsY;

  float leftX, leftY, rightX, rightY, thetaA;
  int nleftX, nleftY, nrightX, nrightY;
  
  int nBool;

  double xre1,xim1; /*used for F-transform*/
  double w,rr1[1500],ff[1500];/*used for F-transform*/

  float CFLN;
  int imax, jmax, itmax, snapStep;

  printf("==================Object Grid Configuration==============\n");
  printf("The X size of object(mm): ");
  scanf("%f", &a);
  
  printf("The Y size of object(mm): ");
  scanf("%f", &b);

  
  printf("The  number of cells in X: ");
  scanf("%d", &imax);

  printf("The  number of cells in Y:");
  scanf("%d", &jmax);

  xx = a*(1e-3)/imax;
  yy = b*(1e-3)/jmax;
  
  /////////////////////////////////
  printf("If the domain contain the inner conductor (1 or 0): ");
  scanf("%d", &nBool);

  if(nBool == 1)
  {
	  /// conductor object
	  
	  printf("The conductivity of the inner conductor: ");
	  scanf("%f", &thetaA);

	  printf("The inner conductor coordiation, bottom left X (mm): ");
	  scanf("%f", &leftX);

	  printf("The inner conductor coordiation, bottom left Y (mm): : ");
	  scanf("%f", &leftY);

	  printf("The inner conductor coordiation, top right X (mm): ");
	  scanf("%f", &rightX);

	  printf("The inner conductor coordiation, top right Y (mm): : ");
	  scanf("%f", &rightY);

	  nleftX = int(leftX*(1e-3)/xx);
	  nleftY = int(leftY*(1e-3)/yy);
	  nrightX = int(rightX*(1e-3)/xx);
	  nrightY = int(rightY*(1e-3)/yy);
  }

  printf("==============Excitation and Time Configuration==========\n");
  printf("The source location (mm): X = ");
  scanf("%f", &sourceX);
  printf("The source location (mm): Y = ");
  scanf("%f", &sourceY);
  i0 = int(sourceX*(1e-3)/xx);
  j0 = int(sourceY*(1e-3)/yy);

  printf("The CFLN: ");
  scanf("%f", &CFLN);

  printf("The total iterations in time: ");
  scanf("%d", &itmax);

  printf("The max frequency(GHz): ");
  scanf("%f", &fmax);
  fmax = fmax*1e9;

  printf("==============Observation Point==========\n");
  printf("The number of snapstep: ");
  scanf("%d", &snapStep);
  printf("The observation point (mm): X = ");
  scanf("%f", &obsX);
  printf("The observation point (mm): Y = ");
  scanf("%f", &obsY);
  nobsX = int(obsX*(1e-3)/xx);
  nobsY = int(obsY*(1e-3)/yy);

  FILE *fp1,*fp2, *fp3;
  char stimulusFilename1[40], stimulusFilename2[40], stimulusFilename3[40]; 
  char addString1[40], addString2[40], addString3[40];

  if(nBool == 1)
	  sprintf(addString1, "DoolittleField_Cond_CFLN=%f.csv", CFLN);
  else
	  sprintf(addString1, "DoolittleField_CFLN=%f.csv", CFLN);
  strcpy(stimulusFilename1, addString1);
  fp1=fopen(stimulusFilename1,"w");

  if(nBool == 1)
	  sprintf(addString2, "DoolittleFreq_Cond_CFLN=%f.csv", CFLN);
  else
	  sprintf(addString2, "DoolittleFreq_CFLN=%f.csv", CFLN);
  strcpy(stimulusFilename2, addString2);
  fp2=fopen(stimulusFilename2,"w");

  if(nBool == 1)
	  sprintf(addString3, "DoolittleConfig_Cond_CFLN=%f.csv", CFLN);
  else
	  sprintf(addString3, "DoolittleConfig_CFLN=%f.csv", CFLN);
  strcpy(stimulusFilename3, addString3);
  fp3=fopen(stimulusFilename3,"w");
 
  fprintf(fp3, "The size of object: X = %f mm, Y = %f mm\n",a, b);
  fprintf(fp3, "imax = %d, jmax = %d\n", imax,jmax);
  fprintf(fp3, "xx = %f(mm),yy = %(mm), ITMAX = %d\n",xx/1.0e-3,yy/1.0e-3,itmax);
  fprintf(fp3, "The max frequency(GHz): %10.7f GHz\n", fmax/1e9);
  fprintf(fp3, "The CFLN = %f\n", CFLN);
  fprintf(fp3, "The source location (mm): X = %f mm , Y = %f mm\n", sourceX, sourceY);
  fprintf(fp3, "The observation point (mm): X = %f mm , Y = %f mm\n", obsX, obsY);
  fclose(fp3);

  printf("imax=%5d,jmax=%5d\n",imax,jmax);
  printf("xx=%12.6f(mm),yy=%12.6f(mm), ITMAX=%d\n",xx/1.0e-3,yy/1.0e-3,itmax);
  c0=3.0e8;
  xmu0=1.0/(eps0*c0*c0);
  e=eps0;
  m=xmu0;
  
/*position of the excitation, you can use almost arbitrary number for i0 and j0*/
//  i0 = int(imax/10);
 // j0 = int(jmax/2);

  ttt=1.0/(c0*sqrt(1.0/(xx*xx)+1.0/(yy*yy)));
  tt=CFLN*ttt;/*CFLN number applied here*/
  printf("ttt=%e,tt=%e,tt/ttt=%10.5f\n",ttt,tt,tt/ttt);
 
  //fmax=4e9;/*maximum frequcncy used for the gaussian pulse*/
 
  t=1.0/(2.0*fmax);/*used for gaussian pulse*/
  t0=4.0*t;/*used for gaussian pulse*/
  nstop=(int)(2*t0/tt);
  pi=4.0*atan(1.0);
  printf("c0=%e,nstop=%d\n",c0,nstop);

  /*some coefficients*/
   t2e=tt/(2.0*e);
   t2h=tt/(2.0*m);

    /*stuff related to C++*/
  // all required variables
  DPTR EX,EY,HZ,HZINC1,HZINC2,IEX,IEY,IHZ;
  /*note- EX, EY and HZ are used to store the field components at time step n and n+1; 
  IEX,IEY,and IHZ are used to store the field components at time step n+1/2;
  HZINC1 and HZINC2 are used for the excitation*/
  
  double bet1, bet2;
  DPTR muArray, thetaArray, epsArray;  /*used for materials array*/

  DPTR aa, bb, cc, r, g1, g2;  /*used for solving the tri-diagonal matrix*/
  PTR H1;/*used to store the recording field components at a certain point to calculate
  the resonance frequency of the cavity*/
  
  // allocate memories for each variable
  H1=new double [itmax+1];
  for (i=0;i<itmax+1;i++){
  H1[i]=0.0;
  }

  //// EX, EY, HZ ////////
  EX = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    EX[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      EX[i][j] = 0.0;
	}
  }

  EY = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    EY[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      EY[i][j] = 0.0;
	}
  }

  HZ = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    HZ[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      HZ[i][j] = 0.0;
	}
  }

  IEX = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    IEX[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      IEX[i][j] = 0.0;
	}
  }

  IEY = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    IEY[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      IEY[i][j] = 0.0;
	}
  }

  IHZ = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    IHZ[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      IHZ[i][j] = 0.0;
	}
  }


  HZINC1 = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    HZINC1[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      HZINC1[i][j] = 0.0;
	}
  }

 HZINC2 = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    HZINC2[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      HZINC2[i][j] = 0.0;
	}
  }

  //////materials //////////
  muArray = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    muArray[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      muArray[i][j] = xmu0;
	}
  }

  thetaArray = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    thetaArray[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      thetaArray[i][j] = 0.0;
	}
  }


  epsArray = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    epsArray[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      epsArray[i][j] = eps0;
	}
  }

  if(nBool == 1)
  {
	  for(i=nleftX; i<=nrightX; i ++)
		  for(j=nleftY; j<=nrightY;j++)
		  {
			  thetaArray[i][j] = thetaA;
			  epsArray[i][j] = eps0;
		  }
  }
  ////triagonal elements /////
  aa = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    aa[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      aa[i][j] = 0.0;
	}
  }

  bb = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    bb[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      bb[i][j] = 0.0;
	}
  }

  cc = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    cc[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      cc[i][j] = 0.0;
	}
  }

  r = new PTR[(imax+1)];
  for(i=0;i<=imax; i++ )
  {
    r[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++){
      r[i][j] = 0.0;
	}
  }

  g1 = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    g1[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      g1[i][j] = 0.0;
	}
  }

  g2 = new PTR[(imax+1)];
  for(i=0;i<=imax;i++)
  {
    g2[i] = new double[jmax+1];
    for(j=0;j<=jmax;j++)
	{
      g2[i][j] = 0.0;
	}
  }