//////////////////////////////////////////////////////////////////////
// FDTD_1D_EH_LORENTZ.cpp: implementation of the CFDTD_1D_EH_LORENTZ class.
// H_1D = Hy
// E_1D = Ez
// the field is propagating in x direction
//////////////////////////////////////////////////////////////////////
#include "FDTD_1D_EH_LORENTZ.h"
#include "run_enviroment.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_EH_LORENTZ::CFDTD_1D_EH_LORENTZ()
{
	H_1D = NULL;   D = NULL;     E_1D = NULL;
	S = NULL;      S_2 = NULL;
	K_D_a = NULL;  K_D_b = NULL; 
	K_H_a = NULL;  K_H_b = NULL;
	K_a = NULL;    K_b = NULL;    K_c = NULL;

	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]

	nr_threads = 1;
}

CFDTD_1D_EH_LORENTZ::~CFDTD_1D_EH_LORENTZ()
{
	Free_Mem_1D();
}

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

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_EH_LORENTZ::Init_Main_Param_1D(long n1d, long n_pml, double *&par_mat,
		                                    double d_t, double d_d)
{
	//dimension of the computational space
    n_1D  = n1d;
	dt    = d_t;
	d     = d_d;
	n_PML = n_pml;

	 //auxiliary parameters to reduce computational cost
	n_L_mat = (long) par_mat[0];//number of Lorentz terms 
	Inv_eps = 1.0/(eps_0*par_mat[1]);
	mu_r = par_mat[2+3*n_L_mat]; //magnetic permittivity
	n_1DMIN1 = n_1D - 1;

	//start the memory allocation
	D = (double *) calloc(n_1D,sizeof(double));
	if(!D)
	{
		Free_Mem_1D();
		return 1;
	}

	E_1D = (double *) calloc(n_1D,sizeof(double));
	if(!E_1D)
	{
		Free_Mem_1D();
		return 1;
	}

	S = Init_Matrix_2D<double>(n_1D,n_L_mat);
	if(!S)
	{
		Free_Mem_1D();
		return 1;
	}

	S_2 = Init_Matrix_2D<double>(n_1D,n_L_mat);
	if(!S_2)
	{
		Free_Mem_1D();
		return 1;
	}
	
	H_1D = (double *) calloc(n_1DMIN1,sizeof(double));
	if(!H_1D)
	{
		return 1;
	}

	//////////////////////////
	K_D_a = (double *) calloc(n_1D,sizeof(double));
	if(!K_D_a)
	{
		Free_Mem_1D();
		return 1;
	}
	
	K_D_b = (double *) calloc(n_1D,sizeof(double));
	if(!K_D_b)
	{
		Free_Mem_1D();
		return 1;
	}
	
	K_H_a = (double *) calloc(n_1DMIN1,sizeof(double));
	if(!K_H_a)
	{
		Free_Mem_1D();
		return 1;
	}
	
	K_H_b = (double *) calloc(n_1DMIN1,sizeof(double));
	if(!K_H_b)
	{
		Free_Mem_1D();
		return 1;
	}
	
	//////////////////////////
	K_a = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_a)
	{
		Free_Mem_1D();
		return 1;
	}
	 
	K_b = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_b)
	{
		Free_Mem_1D();
		return 1;
	}
	
	K_c = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_c)
	{
		Free_Mem_1D();
		return 1;
	}

	double am, delta_0, omega_0;
	
	long i;
	for (i = 0; i<n_L_mat; i++)
	{
		am = par_mat[2+3*i];
		delta_0 = par_mat[2+3*i+1];
		omega_0 = par_mat[2+3*i+2];
		
		K_a[i] = 2*(2 - omega_0*omega_0*dt*dt)/(2 + delta_0*omega_0*dt);
		K_b[i] = (-2 + delta_0*omega_0*dt)/(2 + delta_0*omega_0*dt);
		K_c[i] = 2*am*omega_0*omega_0*dt*dt/(2 + delta_0*omega_0*dt);
	}
	
	return 0;
}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Gaussian pulse
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Init_Gauss_1D(double E_0, double t_0, double t_w)
{
	source_type = 1;
	E0 = E_0;
	t0 = t_0;
	twORtw = t_w*t_w;
}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Sinusoidal plane wave
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Init_Sinus_1D(double E_0, double om, double Phi)
{
	source_type = 2;
	E0 = E_0;
	omega = om;
	phase = Phi;

	alfa = 1.0*omega/(2.0*pi);
}


//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Wave Packet (sinusoidal modulated Gauss pulse)
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Init_Wave_Packet_1D(double E_0, double t_0, double t_w, 
									          double om, double Phi)
{
	source_type = 3;
	E0 = E_0;
	t0 = t_0;
	omega = om;
	phase = Phi;

	twORtw = t_w*t_w;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Set_PML_Param_1D(double eps_r, double mu_r)
{
	long i;
	for (i = 0; i<n_1D; i++)
	{
		K_D_a[i] = 1.0;
		K_D_b[i] = dt/d;

		if (i < n_1DMIN1)
		{
			K_H_a[i] = 1.0;
			K_H_b[i] = dt/(mu_0*mu_r*d);
		}			
	}

	//PML parameters
	double ka_max = 1;
	int exponent = 4;
	double R_err = 1e-16;
	
	double eta = sqrt(mu_0*mu_r/eps_0/eps_r);

	double sigma_PML, sigma_max, ka_PML;
		
	sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*d);
	//sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*d*sqrt(mu_0*mu_r*eps_0*eps_r));

	long cik = 0, ii;
	for (i = 0; i<n_PML; i++)
	{
		//i
		ii = n_1D - 1 - cik;
		sigma_PML  = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_PML     = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent);
		
		K_D_a[ii]  = (2*eps_0*ka_PML - sigma_PML*dt)/(2*eps_0*ka_PML + sigma_PML*dt);
		K_D_b[ii]  = 2*eps_0*dt/( (2*eps_0*ka_PML + sigma_PML*dt)*d );
		
		//i+0.5
		ii = n_1D - 2 - cik;
		sigma_PML  = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_PML     = 1.0 + (ka_max - 1.0)*pow( (n_PML - i - 0.5)/n_PML ,exponent);
        
		K_H_a[ii]  = (2.0*eps_0*ka_PML - sigma_PML*dt)/(2.0*eps_0*ka_PML + sigma_PML*dt);
		K_H_b[ii]  = 2.0*eps_0*dt/( (2.0*eps_0*ka_PML + sigma_PML*dt)*d*mu_0*mu_r);

		cik++;
	}
}

//////////////////////////////////////////////////////////////////////
//Compute the Hz component
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Calc_H_1D()
{
	long i;

	#pragma omp parallel for default(shared) private(i) num_threads(nr_threads)
	for (i = 0; i < n_1DMIN1; i++)
	{
		H_1D[i] = K_H_a[i]*H_1D[i] + K_H_b[i]*(E_1D[i] - E_1D[i+1]);
	}
}

//////////////////////////////////////////////////////////////////////
//Compute the Ey component
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Calc_E_1D(double time)
{
	int kk;
	double Sum_S, Dy_r, S_r, timeMINt0;

	switch (source_type)
	{
		case 1:
			timeMINt0 = time - t0;
			E_1D[0] = E0*exp( - timeMINt0*timeMINt0/twORtw );
			break;
		case 2:
			E_1D[0] = (1.0-exp(-alfa*time))*E0*cos(omega*time + phase);
			break;
		case 3:
			timeMINt0 = time - t0;
			E_1D[0] = E0*cos(omega*timeMINt0 + phase)*exp( - timeMINt0*timeMINt0/twORtw );
			break;
	}
	
	long i;

	#pragma omp parallel for default(shared) private(i,Sum_S,Dy_r,S_r,kk) num_threads(nr_threads)
	for (i = 1; i < n_1D; i++)
	{
		Dy_r = D[i];

		D[i] = K_D_a[i]*D[i] + K_D_b[i]*(H_1D[i-1] - H_1D[i]);

		Sum_S = 0;
		for (kk = 0; kk<n_L_mat; kk++)
		{
			 Sum_S = Sum_S + S[i][kk];
		}
		E_1D[i] = (D[i] - eps_0*Sum_S)*Inv_eps;

		for (kk = 0; kk<n_L_mat; kk++)
		{
			S_r = S[i][kk]; 
			S[i][kk]   = K_a[kk] * S_r +
				 		 K_b[kk] * S_2[i][kk] + 
						 K_c[kk] * E_1D[i];
			S_2[i][kk] = S_r;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Save out the workspace data - the possibility to continue computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_EH_LORENTZ::Save_FDTD_1D_Workspace(char *path, int myrank)
{	
	char *path_name = NULL, *file_name = NULL;
	path_name =(char *) calloc(256,sizeof(char));
	file_name =(char *) calloc(256,sizeof(char));
	if (!path_name && !file_name)
	{
		if (path_name)
			free(path_name);
		if (file_name)
			free(file_name);
		return 1;
	}
	strcpy(path_name,path);
	strcat(path_name,"/data_fdtd_1D");

	//make directory to save F G E
	#ifdef run_enviroment 
		CreateDirectory(path_name,0);
	#else
		mode_t Mode = S_IRWXU;
		mkdir(path_name,Mode);//for UNIX-AIX
	#endif

	//save D
	strcpy(file_name,path_name);
	strcat(file_name,"/D_1D");
	if (save_1D_binary(D, n_1D, myrank, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save S
	strcpy(file_name,path_name);
	strcat(file_name,"/S_1D");
	if (save_2D_binary(S, n_1D, n_L_mat, myrank, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save S_2
	strcpy(file_name,path_name);
	strcat(file_name,"/S_2_1D");
	if (save_2D_binary(S_2, n_1D, n_L_mat, myrank, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save E_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/E_1D");
	if (save_1D_binary(E_1D, n_1D, myrank, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save H_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/H_1D");
	if(save_1D_binary(H_1D, n_1DMIN1, myrank, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}


	if (path_name)
	{
		free(path_name);
		path_name = NULL;
	}

	if (file_name)
	{
		free(file_name);
		file_name = NULL;
	}

	return 0;
}

int CFDTD_1D_EH_LORENTZ::Load_FDTD_1D_Workspace(char *path, int myrank)
{
	char *file_name = NULL;
	file_name =(char *) calloc(256,sizeof(char));
	if (!file_name)
	{
		return 1;
	}
	
	//load D
	sprintf(file_name,"%s/D_1D_%d.dat",path,myrank);
	if (load_1D_binary(D, n_1D, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//save S
	sprintf(file_name,"%s/S_1D_%dsz%d_%d.dat",path,n_1D,n_L_mat,myrank);
	if (load_2D_binary(S, n_1D, n_L_mat, file_name))
	{
		free(file_name);
		return 2; //faild to save data
	}

	//save S_2
	sprintf(file_name,"%s/S_2_1D_%dsz%d_%d.dat",path,n_1D,n_L_mat,myrank);
	if (load_2D_binary(S_2, n_1D, n_L_mat, file_name))
	{
		free(file_name);
		return 2; //faild to save data
	}

	//load E_1D
	sprintf(file_name,"%s/E_1D_%d.dat",path,myrank);
	if (load_1D_binary(E_1D, n_1D, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hx_1D
	sprintf(file_name,"%s/H_1D_%d.dat",path,myrank);
	if (load_1D_binary(H_1D, n_1DMIN1, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	if (file_name)
	{
		free(file_name);
		file_name = NULL;
	}

	return 0;
}

//////////////////////////////////////////////////////////////////////
//Access H_1D and E_1D
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Get_Data_1D(double *&e1D, double *&h1D)
{
	e1D = E_1D;
	h1D = H_1D;
}

//////////////////////////////////////////////////////////////////////
//Free the allocated memory
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_EH_LORENTZ::Free_Mem_1D()
{
	if (H_1D)
	{
		free(H_1D);
		H_1D = NULL;
	}

	if (D)
	{
		free(D);
		D = NULL;
	}

	if (E_1D)
	{
		free(E_1D);
		E_1D = NULL;
	}

	if (S)
		S = Free_Matrix_2D<double>(S);

	if (S_2)
		S_2 = Free_Matrix_2D<double>(S_2);

	if (K_D_a)
	{
		free(K_D_a);
		K_D_a = NULL;
	}

	if (K_D_b)
	{
		free(K_D_b);
		K_D_b = NULL;
	}

	if (K_H_a)
	{
		free(K_H_a);
		K_H_a = NULL;
	}

	if (K_H_b)
	{
		free(K_H_b);
		K_H_b = NULL;
	}

	if (K_a)
	{
		free(K_a);
		K_a = NULL;
	}

	if (K_b)
	{
		free(K_b);
		K_b = NULL;
	}

	if (K_c)
	{
		free(K_c);
		K_c = NULL;
	}
}