#include "stdafx.h"
#include "FDTD_1D_ExHy.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFDTD_1D_ExHy::CFDTD_1D_ExHy(void)
{
	Ex = NULL; 
	Hy = NULL;
	
	Hy_n_1[0] = 0.0;   Hy_n_1[1] = 0.0;
	Hy_n[0]   = 0.0;   Hy_n[1]   = 0.0;

	nr_threads = 1;
}

CFDTD_1D_ExHy::~CFDTD_1D_ExHy(void)
{
	if (Ex)
	{
		free(Ex);
		Ex = NULL;
	}

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

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

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_ExHy::Init(long n_z, double dt, double dz, double eps_r, double mu_r)
{
	nz = n_z;
	nzMIN1 = nz -1;

	Ex = (double *)calloc(nz,sizeof(double));
	if (!Ex)
	{
		return 1;
	}

	Hy = (double *)calloc(nz,sizeof(double));
	if (!Hy)
	{
		free(Ex);
		Ex = NULL;
		return 1;
	}

	double pi = 180.0*atan(1.0)/45.0;
	//permittivity of free space 
	double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]
	
	K_E = dt/(eps_0*eps_r*dz);
	K_H = dt/(mu_0*mu_r*dz);

	double speed = 1.0/sqrt(eps_0*eps_r*mu_0*mu_r); //the speed of the light in the vacuum

	K_MUR_1 = ( speed*dt - dz )/( speed*dt + dz );
	K_MUR_2 =  2.0*dz/( speed*dt + dz );

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::Init_Sin_Source(double E_0, double om, double Phi)
{
	double pi = 180.0*atan(1.0)/45.0;

	E0 = E_0;
	omega = om;
	phi = Phi;

	alfa = 3.0*omega/(2.0*pi);

	Source_Type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize  a Gaussian source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::Init_Gauss_Source(double E_0, double t_0, double t_w)
{
	E0 = E_0;
	t0 = t_0;
	tw = t_w;
	
	Source_Type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal modulated Gaussian pulse
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::Init_GaussSin(double E_0, double om, double Phi, double t_0, double t_w)
{
	E0 = E_0;
	omega = om;
	phi = Phi;
	t0 = t_0;
	tw = t_w;

	Source_Type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ex field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::calc_Ex(double time)
{
	long k;

	#pragma omp parallel for default(shared) private(k) schedule(dynamic,nr_threads)
	for (k = 1; k < nz; k++)
	{
		Ex[k] += K_E*( Hy[k-1] - Hy[k] );
	}

	switch (Source_Type)
	{
		case 1:
			Ex[0] = (1.0-exp(-alfa*time))*E0*cos(omega*time + phi);
			break;
		case 2:
			Ex[0] = E0*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
		case 3:
			Ex[0] = E0*cos(omega*(time-t0) + phi)*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
	}
	
	
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::calc_Hy()
{
	long k;

	#pragma omp parallel for default(shared) private(k) schedule(dynamic,nr_threads)
	for (k = 0; k < nzMIN1; k++)
	{
		Hy[k] += K_H*( Ex[k] - Ex[k+1] );
	}

	Hy[nz-1] = -Hy_n_1[1] + K_MUR_1*(Hy_n_1[0] + Hy[nz-2]) + K_MUR_2*(Hy_n[0] + Hy_n[1]);
	
	Hy_n_1[0] = Hy_n[0]; 
	Hy_n_1[1] = Hy_n[1];
	Hy_n[0]   = Hy[nz-1]; 
	Hy_n[1]   = Hy[nz-2];
}

///////////////////////////////////////////////////////////////////////////////////////
//Reinitialize with zero all the field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::Reset_Field_Comp(void)
{
	long k;
	for (k = 0; k < nz; k++)
	{
		Ex[k] = 0.0;
		Hy[k] = 0.0;
	}

	Hy_n_1[0] = 0.0; 
	Hy_n_1[1] = 0.0;
	Hy_n[0]   = 0.0;
	Hy_n[1]   = 0.0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the address of Ex and Hy
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_ExHy::Get_Data(double *&X, double *&Y)
{
	X = Ex;
	Y = Hy;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save out the workspace data - the possibility to continue computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_ExHy::Save_FDTD_1D_Workspace(char *path_name)
{	
	char *file_name = NULL;
	file_name =(char *) calloc(256,sizeof(char));
	if (!file_name)
	{
		if (file_name)
			free(file_name);
		return 1;
	}
	//save Ex_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/Ex_1D");
	if (save_1D_binary(Ex, nz, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save Hy_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/Hy_1D");
	if(save_1D_binary(Hy, nz, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save Hy_n
	strcpy(file_name,path_name);
	strcat(file_name,"/Hy_n_1D");
	if(save_1D_binary(Hy_n, 2, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save Hy_n_1
	strcpy(file_name,path_name);
	strcat(file_name,"/Hy_n_1_1D");
	if(save_1D_binary(Hy_n_1, 2, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	free(path_name);
	free(file_name);

	return 0;
}

int CFDTD_1D_ExHy::Load_FDTD_1D_Workspace(char *path)
{
	char *file_name = NULL;
	file_name =(char *) calloc(256,sizeof(char));
	if (!file_name)
	{
		return 1;
	}
	
	//load Ex_1D
	strcpy(file_name,path);
	strcat(file_name,"/Ex_1D_0.dat");
	if (load_1D_binary(Ex, nz, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hy_1D
	strcpy(file_name,path);
	strcat(file_name,"/Hy_1D_0.dat");
	if (load_1D_binary(Hy, nz, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hy_n
	strcpy(file_name,path);
	strcat(file_name,"/Hy_n_1D_0.dat");
	if (load_1D_binary(Hy_n, 2, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hy_n_1
	strcpy(file_name,path);
	strcat(file_name,"/Hy_n_1_1D_0.dat");
	if (load_1D_binary(Hy_n_1, 2, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	free(file_name);

	return 0;
}