#include "FDTD_3D.h"
#include "run_enviroment.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_3D::CFDTD_3D(void)
{
	pi = 180.0*atan(1.0)/45.0;
	//permittivity of free space 
	eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]

	Ind = NULL;

	Mat = NULL;
	mu_r = NULL;

	Ex = NULL; Gx = NULL; Fx = NULL;
	Hz_recv_j = NULL; Hy_recv_k = NULL;
	Ex_send_j = NULL; Ex_send_k = NULL;

	Ey = NULL; Gy = NULL; Fy = NULL;
	Hx_recv_k = NULL; Hz_recv_i = NULL;
	Ey_send_i = NULL; Ey_send_k = NULL;
	
	Ez = NULL; Gz = NULL; Fz = NULL;
	Hy_recv_i = NULL; Hx_recv_j = NULL;
	Ez_send_i = NULL; Ez_send_j = NULL;
	
	Hx = NULL; Bx = NULL; 
	Ey_recv_k = NULL; Ez_recv_j = NULL;
	Hx_send_j = NULL; Hx_send_k = NULL;
	
	Hy = NULL; By = NULL;
	Ez_recv_i = NULL; Ex_recv_k = NULL;
	Hy_send_i = NULL; Hy_send_k = NULL;
	
	Hz = NULL; Bz = NULL;
	Ex_recv_j = NULL; Ey_recv_i = NULL;
	Hz_send_i = NULL; Hz_send_j = NULL;
	
	//Coefficients containing the PML boundary parameters
	//Electric field
	//Electric field
	K_Gx_a = NULL; K_Gx_b = NULL;
	K_Ex_a = NULL; K_Ex_b = NULL; K_Ex_c = NULL; K_Ex_d = NULL; 
	K_Gy_a = NULL; K_Gy_b = NULL;
	K_Ey_a = NULL; K_Ey_b = NULL; K_Ey_c = NULL; K_Ey_d = NULL; 
	K_Gz_a = NULL; K_Gz_b = NULL;
	K_Ez_a = NULL; K_Ez_b = NULL; K_Ez_c = NULL; K_Ez_d = NULL; 

	
	//Magnetic field
	K_Bx_a = NULL; K_Bx_b = NULL;
	K_Hx_a = NULL; K_Hx_b = NULL; K_Hx_c = NULL; K_Hx_d = NULL; 
	K_By_a = NULL; K_By_b = NULL;
	K_Hy_a = NULL; K_Hy_b = NULL; K_Hy_c = NULL; K_Hy_d = NULL; 
	K_Bz_a = NULL; K_Bz_b = NULL;
	K_Hz_a = NULL; K_Hz_b = NULL; K_Hz_c = NULL; K_Hz_d = NULL; 

	K_a = NULL; K_b  = NULL;

	jel_TS = 0;

	//the incident field for total field scattered field formulation
	jel_TS_planes = NULL;

	E_1D = NULL; H_1D = NULL; 
	ll_1D_E = NULL; ll_1D_H = NULL;

	//the recuired incident magnetic fields to compute the electric field components
	Hz_i0 = NULL; Hy_i0 = NULL; Hz_i1 = NULL; Hy_i1 = NULL;
	Hz_j0 = NULL; Hx_j0 = NULL; Hz_j1 = NULL; Hx_j1 = NULL;
	Hy_k0 = NULL; Hx_k0 = NULL; Hy_k1 = NULL; Hx_k1 = NULL;
	face_Hz_i0 = NULL; face_Hy_i0 = NULL; face_Hz_i1 = NULL; face_Hy_i1 = NULL;
	face_Hz_j0 = NULL; face_Hx_j0 = NULL; face_Hz_j1 = NULL; face_Hx_j1 = NULL;
	face_Hy_k0 = NULL; face_Hx_k0 = NULL; face_Hy_k1 = NULL; face_Hx_k1 = NULL;

	//the recuired incident electric fields to compute the magnetic field components
	Ey_i0 = NULL; Ez_i0 = NULL; Ey_i1 = NULL; Ez_i1 = NULL;
	Ex_j0 = NULL; Ez_j0 = NULL; Ex_j1 = NULL; Ez_j1 = NULL;
	Ex_k0 = NULL; Ey_k0 = NULL; Ex_k1 = NULL; Ey_k1 = NULL;
	face_Ey_i0 = NULL; face_Ez_i0 = NULL; face_Ey_i1 = NULL; face_Ez_i1 = NULL;
	face_Ex_j0 = NULL; face_Ez_j0 = NULL; face_Ex_j1 = NULL; face_Ez_j1 = NULL;
	face_Ex_k0 = NULL; face_Ey_k0 = NULL; face_Ex_k1 = NULL; face_Ey_k1 = NULL;

	jel_plane_wave = 0;

	pt_source_Ex = 0;
	pt_source_Ey = 0;
	pt_source_Ez = 0;
	pt_source_Hx = 0;
	pt_source_Hy = 0;
	pt_source_Hz = 0;

	Param_ptSource = NULL;
	Coord_ptSource = NULL;

	path_name_Ex = NULL; path_name_Ey = NULL; path_name_Ez = NULL;
	path_name_Hx = NULL; path_name_Hy = NULL; path_name_Hz = NULL;

	nx_yz = -1; ny_xz = -1; nz_xy = -1;

	Ex_Foll = NULL; Ey_Foll = NULL; Ez_Foll = NULL;
	Hx_Foll = NULL; Hy_Foll = NULL; Hz_Foll = NULL; 
}

CFDTD_3D::~CFDTD_3D(void)
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the number of threads
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Init_nr_THR(int nr_Threads)
{
	nr_threads = nr_Threads;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the time
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_3D::Set_Time(double Time, long Iter)
{
	time = Time;
	iter = Iter;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the FDTD algorithm
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_3D::Init(long ***&Index, long iSta, long iEnd, long jSta, long jEnd, long kSta, 
				   long kEnd, long nLx, long nLy, long nLz, long Nx, long Ny, long Nz,
				   int myRank, int myRanki, int myRankj, int myRankk, int iProcs, 
				   int jProcs, int kProcs, double **&Mater, long nMat, long npml_x1, 
				   long npml_x2, long npml_y1, long npml_y2, long npml_z1, 
				   long npml_z2, double d_x, double d_y, double d_z, double d_t)
{
	Ind = Index;  //Dimensions [nlx][nly][nlz]

	Mat = Mater;  //coefficients of the Lorentz materials
	n_Mat = nMat; //number of materials presents in the investigated geometry
	
	//the global indices of the computational subvolume
	ista = iSta;  
	iend = iEnd; 
	jsta = jSta;  
	jend = jEnd; 
	ksta = kSta;  
	kend = kEnd; 
	
	//dimensions of the full computational space
    nx = Nx;
	ny = Ny; 
	nz = Nz;

	//the size of the subvolume
	nlx = nLx; 
	nly = nLy; 
	nlz = nLz; 

	myrank  = myRank;
	myrank_i = myRanki;
	myrank_j = myRankj;
	myrank_k = myRankk;
						    
	iprocs = iProcs;
	jprocs = jProcs;
	kprocs = kProcs;

	iprocsMIN1 = iprocs - 1;
	jprocsMIN1 = jprocs - 1;
	kprocsMIN1 = kprocs - 1;

	//time increment
	dt = d_t;

	//elementary cell sizes
	dx = d_x;
	dy = d_y;
	dz = d_z;

	nPML_x_1 = npml_x1;
	nPML_x_2 = npml_x2;
	nPML_y_1 = npml_y1;
	nPML_y_2 = npml_y2;
	nPML_z_1 = npml_z1;
	nPML_z_2 = npml_z2;

	inv_dx = 1.0/dx;
	inv_dy = 1.0/dy;
	inv_dz = 1.0/dz;

	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Ex 
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Ex = ista;	   iend_Ex = iend;    jsta_Ex = jsta;
	jend_Ex = jend;	   ksta_Ex = ksta;	  kend_Ex = kend;
	nlx_Ex = nlx;	   nly_Ex = nly;	  nlz_Ex = nlz;
	if (myrank_i == iprocs-1)
	{
		iend_Ex--;
		nlx_Ex--;
	}
	
	Ex = Init_Matrix_3D<double>(nlx_Ex, nly_Ex, nlz_Ex);
	if (!Ex)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex - ");
		return 1;
	}
	
	Gx = Init_Matrix_3D<double>(nlx_Ex, nly_Ex, nlz_Ex);
	if (!Gx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Gx - ");
		return 1;
	}

	Fx = Init_Matrix_3D<double>(nlx_Ex, nly_Ex, nlz_Ex);
	if (!Fx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Fx - ");
		return 1;
	}

	Hz_recv_j = Init_Matrix_2D<double>(nlx_Ex, nlz_Ex);
	if (!Hz_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_recv_j - ");
		return 1;
	}

	Hy_recv_k = Init_Matrix_2D<double>(nlx_Ex, nly_Ex);
	if (!Hy_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_recv_k - ");
		return 1;
	}

	Ex_send_j = Init_Matrix_2D<double>(nlx_Ex, nlz_Ex);
	if (!Ex_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_send_j - ");
		return 1;
	}

	Ex_send_k = Init_Matrix_2D<double>(nlx_Ex, nly_Ex);
	if (!Ex_send_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_send_k - ");
		return 1;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Ey
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Ey = ista;	   iend_Ey = iend;    jsta_Ey = jsta;
	jend_Ey = jend;	   ksta_Ey = ksta;	  kend_Ey = kend;
	nlx_Ey = nlx;	   nly_Ey = nly;	  nlz_Ey = nlz;
	if (myrank_j == jprocs-1)
	{
		jend_Ey--;
		nly_Ey--;
	}

	Ey = Init_Matrix_3D<double>(nlx_Ey, nly_Ey, nlz_Ey);
	if (!Ey)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey - ");
		return 1;
	}

	Gy = Init_Matrix_3D<double>(nlx_Ey, nly_Ey, nlz_Ey);
	if (!Gy)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Gy - ");
		return 1;
	}

	Fy = Init_Matrix_3D<double>(nlx_Ey, nly_Ey, nlz_Ey);
	if (!Fy)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Fy - ");
		return 1;
	}
	
	Hx_recv_k = Init_Matrix_2D<double>(nlx_Ey, nly_Ey);
	if (!Hx_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_recv_k - ");
		return 1;
	}

	Hz_recv_i = Init_Matrix_2D<double>(nly_Ey, nlz_Ey);
	if (!Hz_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hz_recv_i - ");
		return 1;
	}

	Ey_send_i = Init_Matrix_2D<double>(nly_Ey, nlz_Ey);
	if (!Ey_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_send_i - ");
		return 1;
	}

	Ey_send_k = Init_Matrix_2D<double>(nlx_Ey, nly_Ey);
	if (!Ey_send_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_send_k - ");
		return 1;
	}

	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Ez
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Ez = ista;	   iend_Ez = iend;    jsta_Ez = jsta;
	jend_Ez = jend;	   ksta_Ez = ksta;	  kend_Ez = kend;
	nlx_Ez = nlx;	   nly_Ez = nly;	  nlz_Ez = nlz;
	if (myrank_k == kprocs-1)
	{
		kend_Ez--;
		nlz_Ez--;
	}

	Ez = Init_Matrix_3D<double>(nlx_Ez, nly_Ez, nlz_Ez);
	if (!Ez)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez - ");
		return 1;
	}

	Gz = Init_Matrix_3D<double>(nlx_Ez, nly_Ez, nlz_Ez);
	if (!Gz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Gz - ");
		return 1;
	}

	Fz = Init_Matrix_3D<double>(nlx_Ez, nly_Ez, nlz_Ez);
	if (!Fz)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Fz - ");
		return 1;
	}

	Hy_recv_i = Init_Matrix_2D<double>(nly_Ez, nlz_Ez);
	if (!Hy_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_recv_i - ");
		return 1;
	}

	Hx_recv_j = Init_Matrix_2D<double>(nlx_Ez, nlz_Ez);
	if (!Hx_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_recv_j - ");
		return 1;
	}

	Ez_send_i = Init_Matrix_2D<double>(nly_Ez, nlz_Ez);
	if (!Ez_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_send_i - ");
		return 1;
	}

	Ez_send_j = Init_Matrix_2D<double>(nlx_Ez, nlz_Ez);
	if (!Ez_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_send_j - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Hx
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hx = ista;	   iend_Hx = iend;    jsta_Hx = jsta;
	jend_Hx = jend;	   ksta_Hx = ksta;	  kend_Hx = kend;
	nlx_Hx = nlx;	   nly_Hx = nly;	  nlz_Hx = nlz;
	if (myrank_j == jprocs-1)
	{
		jend_Hx--;
		nly_Hx--;
	}
	if (myrank_k == kprocs-1)
	{
		kend_Hx--;
		nlz_Hx--;
	}

	nly_HxMIN1 = nly_Hx - 1;
	nlz_HxMIN1 = nlz_Hx - 1;

	Hx = Init_Matrix_3D<double>(nlx_Hx, nly_Hx, nlz_Hx);
	if (!Hx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx - ");
		return 1;
	}

	Bx = Init_Matrix_3D<double>(nlx_Hx, nly_Hx, nlz_Hx);
	if (!Bx)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Bx - ");
		return 1;
	}

	Ey_recv_k = Init_Matrix_2D<double>(nlx_Hx, nly_Hx);
	if (!Ey_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ey_recv_k - ");
		return 1;
	}

	Ez_recv_j = Init_Matrix_2D<double>(nlx_Hx, nlz_Hx);
	if (!Ez_recv_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_recv_j - ");
		return 1;
	}

	Hx_send_j = Init_Matrix_2D<double>(nlx_Hx, nlz_Hx);
	if (!Hx_send_j)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_send_j - ");
		return 1;
	}

	Hx_send_k = Init_Matrix_2D<double>(nlx_Hx, nly_Hx);
	if (!Hx_send_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hx_send_k - ");
		return 1;
	}
	
	///////////////////////////////////////////////////////////////////////////////////////
	//memory allocations to compute Hy
	///////////////////////////////////////////////////////////////////////////////////////
	ista_Hy = ista;	   iend_Hy = iend;    jsta_Hy = jsta;
	jend_Hy = jend;	   ksta_Hy = ksta;	  kend_Hy = kend;
	nlx_Hy = nlx;	   nly_Hy = nly;	  nlz_Hy = nlz;
	if (myrank_i == iprocs-1)
	{
		iend_Hy--;
		nlx_Hy--;
	}
	if (myrank_k == kprocs-1)
	{
		kend_Hy--;
		nlz_Hy--;
	}

	nlx_HyMIN1 = nlx_Hy - 1;
	nlz_HyMIN1 = nlz_Hy - 1;

	Hy = Init_Matrix_3D<double>(nlx_Hy, nly_Hy, nlz_Hy);
	if (!Hy)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy - ");
		return 1;
	}

	By = Init_Matrix_3D<double>(nlx_Hy, nly_Hy, nlz_Hy);
	if (!By)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - By - ");
		return 1;
	}

	Ez_recv_i = Init_Matrix_2D<double>(nly_Hy, nlz_Hy);
	if (!Ez_recv_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ez_recv_i - ");
		return 1;
	}

	Ex_recv_k = Init_Matrix_2D<double>(nlx_Hy, nly_Hy);
	if (!Ex_recv_k)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Ex_recv_k - ");
		return 1;
	}

	Hy_send_i = Init_Matrix_2D<double>(nly_Hy, nlz_Hy);
	if (!Hy_send_i)
	{
		ErrorMessage(1, myrank, " -- Memory allocation problem - Hy_send_i - ");
		return 1;
	}

	Hy_send_k = Init_Matrix_2D<double>(nlx_Hy, nly_Hy);
	if (!Hy_send_k)
	{