#pragma once

#ifndef NULL
	#define NULL 0
#endif

#include "FDTD_1D_EH_PML_LOSS.h"

class CFDTD_3D
{
private:
	double pi, eps_0, mu_0;
	
	long ***Ind;  //indicates the material type in each FDTD cell
	
	double **Mat, *mu_r; //coefficients of the Lorentz materials
				  // - Mat[l][0] - nr of Lorentz terms
	              // - Mat[l][1] - eps_r_inf
				  // - Mat[l][2] - a(l) ; Mat[i][3] - delta(l); Mat[i][4] - omega_0(l)
				  // ... 
	              // - Mat[l][3*Mat[l][0]+3] - mu_r
	long n_Mat;   //number of materials present in the investigated geometry

	//MPI parameters
	int myrank, myrank_i, myrank_j, myrank_k;
	int iprocs, jprocs, kprocs;
	int iprocsMIN1, jprocsMIN1, kprocsMIN1;

	//OpenMP parameters
	int nr_threads;

	//the full size of the computational volume
	long nx, ny, nz;

	double time;
	double dx, dy, dz, dt;
	double inv_dx, inv_dy, inv_dz;

	long iter;
	long switch_off_time;

	//Decompositon of the geometry
	long ista, iend;
	long jsta, jend;
	long ksta, kend;
	long nlx, nly, nlz;

	//Initialization of the FDTD algorithm parameters
	double ***Ex, ***Gx, ***Fx;
	double **Hz_recv_j, **Hy_recv_k, **Ex_send_j, **Ex_send_k;

	double ***Ey, ***Gy, ***Fy;
	double **Hx_recv_k, **Hz_recv_i, **Ey_send_i, **Ey_send_k;
	
	double ***Ez, ***Gz, ***Fz;
	double **Hy_recv_i, **Hx_recv_j, **Ez_send_i, **Ez_send_j;
	
	double ***Hx, ***Bx; 
	double **Ey_recv_k, **Ez_recv_j, **Hx_send_j, **Hx_send_k;
	
	double ***Hy, ***By;
	double **Ez_recv_i, **Ex_recv_k, **Hy_send_i, **Hy_send_k;
	
	double ***Hz, ***Bz;
	double **Ex_recv_j, **Ey_recv_i, **Hz_send_i, **Hz_send_j;

	long ista_Ex, iend_Ex, jsta_Ex, jend_Ex, ksta_Ex, kend_Ex;
	long ista_Ey, iend_Ey, jsta_Ey, jend_Ey, ksta_Ey, kend_Ey;
	long ista_Ez, iend_Ez, jsta_Ez, jend_Ez, ksta_Ez, kend_Ez;
	long ista_Hx, iend_Hx, jsta_Hx, jend_Hx, ksta_Hx, kend_Hx;
	long ista_Hy, iend_Hy, jsta_Hy, jend_Hy, ksta_Hy, kend_Hy;
	long ista_Hz, iend_Hz, jsta_Hz, jend_Hz, ksta_Hz, kend_Hz;
	long nlx_Ex, nly_Ex, nlz_Ex, nlx_Ey, nly_Ey, nlz_Ey, nlx_Ez, nly_Ez, nlz_Ez;
	long nlx_Hx, nly_Hx, nlz_Hx, nlx_Hy, nly_Hy, nlz_Hy, nlx_Hz, nly_Hz, nlz_Hz;
	long nly_HxMIN1, nlz_HxMIN1, nlx_HyMIN1, nlz_HyMIN1, nlx_HzMIN1, nly_HzMIN1;

	//Coefficients containing the PML boundary parameters
	//Electric field
	double *K_Gx_a, *K_Gx_b;
	double *K_Ex_a, *K_Ex_b, *K_Ex_c, *K_Ex_d; 
	double *K_Gy_a, *K_Gy_b;
	double *K_Ey_a, *K_Ey_b, *K_Ey_c, *K_Ey_d; 
	double *K_Gz_a, *K_Gz_b;
	double *K_Ez_a, *K_Ez_b, *K_Ez_c, *K_Ez_d; 
	
	//Magnetic field
	double *K_Bx_a, *K_Bx_b;
	double *K_Hx_a, *K_Hx_b, *K_Hx_c, *K_Hx_d; 
	double *K_By_a, *K_By_b;
	double *K_Hy_a, *K_Hy_b, *K_Hy_c, *K_Hy_d; 
	double *K_Bz_a, *K_Bz_b;
	double *K_Hz_a, *K_Hz_b, *K_Hz_c, *K_Hz_d; 

	//Material parameters
	double *K_a, *K_b;

	//PML parametrers
	long nPML_x_1, nPML_x_2, nPML_y_1, nPML_y_2, nPML_z_1, nPML_z_2;

	int jel_plane_wave; //0 - point source; 1- plane wave source

	//number of point sources present in the computational domain
	double **Param_ptSource;
	long **Coord_ptSource;
	long n_Coord_ptSource;
	int source_type;
	int pt_source_Ex, pt_source_Ey, pt_source_Ez;
	int pt_source_Hx, pt_source_Hy, pt_source_Hz;
	double constORalfa;

	//plane wave excitation with total field - scatterd field formulation
	CFDTD_1D_EH_PML_LOSS fdtd_1D;
	long n_1D, n_1D_MIN_1, nPML_1D;
	double d_1D;

	//the interface between the total field-scattered field zone
	double ct_Ex_1, ct_Ex_2, ct_Ey_1, ct_Ey_2, ct_Ez_1, ct_Ez_2;
	double ct_Hx_1, ct_Hx_2, ct_Hy_1, ct_Hy_2, ct_Hz_1, ct_Hz_2;

	//the incident field for total field scattered field formulation
	double *E_1D, *H_1D; 
	double *ll_1D_E, *ll_1D_H;

	//the direction of propagation and the polarization of the plane wave
	double teta;  //angle relative to +z axis 0 < teta < 180
	double phi;   //angle relative to +x axis 0 <= phi < 360
	double gamma; //polarization
	
	//the recuired incident magnetic fields to compute the electric field components
	double **Hz_i0, **Hy_i0, **Hz_i1, **Hy_i1;
	double **Hz_j0, **Hx_j0, **Hz_j1, **Hx_j1;
	double **Hy_k0, **Hx_k0, **Hy_k1, **Hx_k1;
	double **face_Hz_i0, **face_Hy_i0, **face_Hz_i1, **face_Hy_i1;
	double **face_Hz_j0, **face_Hx_j0, **face_Hz_j1, **face_Hx_j1;
	double **face_Hy_k0, **face_Hx_k0, **face_Hy_k1, **face_Hx_k1;

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

	long n_TS; //the size of the Scattered Field zone
	int *jel_TS_planes, jel_TS;
	//the position of the Total Field - Scattered Field region 
	long n_TS_xa, n_TS_xb, n_TS_ya, n_TS_yb,n_TS_za, n_TS_zb;
	//the local position of the Total Field - Scattered Field region 
	long n_TS_xa_Loc, n_TS_xb_Loc, n_TS_ya_Loc, n_TS_yb_Loc, n_TS_za_Loc, n_TS_zb_Loc;
	long n_TS_xa_Glob, n_TS_ya_Glob, n_TS_za_Glob;
	long ista_TS,  jsta_TS,  ksta_TS; //, iend_TS, jend_TS, kend_TS;
	long l_TS_x, l_TS_y, l_TS_z;
	long l_y_Loc_HzEy_i, l_z_Loc_HzEy_i, l_y_Loc_HyEz_i, l_z_Loc_HyEz_i;
    long l_x_Loc_HzEx_j, l_z_Loc_HzEx_j, l_x_Loc_HxEz_j, l_z_Loc_HxEz_j;
	long l_x_Loc_HyEx_k, l_y_Loc_HyEx_k, l_x_Loc_HxEy_k, l_y_Loc_HxEy_k;

	//Auxiliary parameters for the Total Field - Scatterd Field formulation
	long n_TS_xa_Loc_MIN_1, n_TS_ya_Loc_MIN_1, n_TS_za_Loc_MIN_1;
	double cos_teta, sin_teta, cos_phi, sin_phi, sin_gamma, cos_gamma;
	long lyLocORlzLoc_HzEy_i, lyLocORlzLoc_HyEz_i;
	long lxLocORlzLoc_HzEx_j, lxLocORlzLoc_HxEz_j;
	long lxLocORlyLoc_HyEx_k, lxLocORlyLoc_HxEy_k;
	long l_x_Loc, l_y_Loc, l_z_Loc;
	long n_TS_xb_Loc_HzEx_j, n_TS_zb_Loc_HzEx_j, n_TS_xb_Loc_HyEx_k, n_TS_yb_Loc_HyEx_k;
	long n_TS_yb_Loc_HzEy_i, n_TS_zb_Loc_HzEy_i, n_TS_xb_Loc_HxEy_k, n_TS_yb_Loc_HxEy_k;
	long n_TS_yb_Loc_HyEz_i, n_TS_zb_Loc_HyEz_i, n_TS_xb_Loc_HxEz_j, n_TS_zb_Loc_HxEz_j;
	
	//the coordinates of the followed field components
	long **Ind_Foll; 
	 //number of the followed points and the number of iterations
	long length_Ind_Foll, n_tot;
	//stores the followed field components
	double **Hx_Foll, **Hy_Foll, **Hz_Foll, **Ex_Foll, **Ey_Foll, **Ez_Foll;

	//to save the field components
	long n_x_a, n_x_b, n_y_a, n_y_b, n_z_a, n_z_b; //identify the saved volume
	long n_x_b1, n_y_b1, n_z_b1;
	long nx_yz, ny_xz, nz_xy;
	char *path_name_Ex, *path_name_Ey, *path_name_Ez;
	char *path_name_Hx, *path_name_Hy, *path_name_Hz;
	int jel_Save_Slice;
	long nr_save;
	long saveFROMinst, saveTOinst;

public:
	CFDTD_3D(void);
	~CFDTD_3D(void);

	void Init_nr_THR(int nr_Threads);
	
	int 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);

	int Init_Followed(long **Ind_Followed, long nFoll, long n_t);

	void Set_Data_Followed(long n_t);

	void Get_Data_Followed(double **&e_x, double **&e_y, double **&e_z, 
						   double **&h_x, double **&h_y, double **&h_z);

	int Init_Save_FieldSlices(long nX_yz, long nY_xz, long nZ_xy,
							  char *path_name_ex, char *path_name_ey, char *path_name_ez,
							  char *path_name_hx, char *path_name_hy, char *path_name_hz);

	void Init_ptSource(long **Coord, long nCoord, double **Par_Excit, int s_type,
					   int Pt_Source_Ex, int Pt_Source_Ey, int Pt_Source_Ez,
	                   int Pt_Source_Hx, int Pt_Source_Hy, int Pt_Source_Hz,
					   long switch_OFF_time);


	int Init_TS(double Teta, double Phi, double Gamma, long N_TS,
				long nPML_1D, int source_type, double X0, double t0,
				double tw, double omega, double phase, double PML_eps_r_1D,
				double PML_mu_r_1D);

	void Get_SendRecv(double **&hz_recv_j, double **&hy_recv_k, 
					  double **&ex_send_j, double **&ex_send_k,
					  double **&hx_recv_k, double **&hz_recv_i, 
					  double **&ey_send_i, double **&ey_send_k,
					  double **&hy_recv_i, double **&hx_recv_j,
					  double **&ez_send_i, double **&ez_send_j,
					  double **&ey_recv_k, double **&ez_recv_j,
					  double **&hx_send_j, double **&hx_send_k,
					  double **&ez_recv_i, double **&ex_recv_k, 
					  double **&hy_send_i, double **&hy_send_k,
					  double **&ex_recv_j, double **&ey_recv_i, 
					  double **&hz_send_i, double **&hz_send_j);

	int Interp_H_TS();
	int Interp_E_TS();

	int Init_PML_Par_x(double eps_r_x_1, double mu_r_x_1, double eps_r_x_2, 
					   double mu_r_x_2);
	int Init_PML_Par_y(double eps_r_y_1, double mu_r_y_1, double eps_r_y_2, 
					   double mu_r_y_2);
	int Init_PML_Par_z(double eps_r_z_1, double mu_r_z_1, double eps_r_z_2, 
					   double mu_r_z_2);

	void Get_jel_TS(int *jelTS); 

	void Set_Time(double Time, long Iter); 

	void Calc_Ex(long  nlx, long  nly, long  nlz);
	void Calc_Ey(long  nlx, long  nly, long  nlz);
	void Calc_Ez(long  nlx, long  nly, long  nlz);
	void Calc_Hx(long  nlx, long  nly, long  nlz);
	void Calc_Hy(long  nlx, long  nly, long  nlz);
	void Calc_Hz(long  nlx, long  nly, long  nlz);

	void Update_Ex_send();
	void Update_Ey_send();
	void Update_Ez_send();
	void Update_Hx_send();
	void Update_Hy_send();
	void Update_Hz_send();

	int Save_FieldSlices(long iter);

	int Load_FGE_BH(char *Path, long *num_iter_0, double *time, long *iter_nr_Save);
	int Save_FGE_BH(char *Path, long iter, double time, long iter_nr_Save);

private:
	//Point Source
	void PtSource_J(double ***&xx, double time);

	int Interp_1D(double *x, int length_x, double *f, double *x_i, int length_x_i, double *g);
	
	//TF_SF Ex
	void TS_Ex_j0();
	void TS_Ex_j1();
	void TS_Ex_k0();
	void TS_Ex_k1();

	//TF_SF Ey
	void TS_Ey_i0(); 
	void TS_Ey_i1(); 
	void TS_Ey_k0(); 
	void TS_Ey_k1();

	//TF_SF Ez
	void TS_Ez_i0(); 
	void TS_Ez_i1(); 
	void TS_Ez_j0(); 
	void TS_Ez_j1();

	//TF_SF Hx
	void TS_Hx_j0(); 
	void TS_Hx_j1(); 
	void TS_Hx_k0(); 
	void TS_Hx_k1();

	//TF_SF Hy
	void TS_Hy_i0();
	void TS_Hy_i1(); 
	void TS_Hy_k0(); 
	void TS_Hy_k1();

	//TF_SF Hz
	void TS_Hz_i0();
	void TS_Hz_i1(); 
	void TS_Hz_j0(); 
	void TS_Hz_j1();

	//Compute the phase velocity
	int Phase_velocity(double dx, double dy, double dz, double dt, 
				       double eps_r, double mu_r, double om, 
				       double teta, double phi, double *phase_vel);

	void Free_Mem();
	
	void ErrorMessage(int er, int myrank, char* text);
};