// FDTD_1D_HzEy_LORENTZ.cpp: implementation of the CFDTD_1D_HzEy_LORENTZ class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "fdtd_2D_TE_PML_Lorentz_period.h"
#include "FDTD_1D_HzEy_LORENTZ.h"
#include "Math.h"
#include "Matrix.h"


#ifdef _DEBUG
#undef THIS_FILE
static char THIS_FILE[]=__FILE__;
#define new DEBUG_NEW
#endif


//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_HzEy_LORENTZ::CFDTD_1D_HzEy_LORENTZ()
{
	Hz_1D = NULL; Dy = NULL; Ey_1D = NULL;
	Sy = NULL; Sy_2 = NULL;
	K_E1_a = NULL; K_E1_b = NULL; 
	K_E2_a = NULL; K_E2_b = NULL;
	
	pi = 3.1415926535897932384626433832795;
	eps_0 = 8.854e-12; // [F/m]
	mu_0 = 4*pi*1e-7; // [H/m]
}

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

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
BOOL CFDTD_1D_HzEy_LORENTZ::Init_Main_Param_1D(int n_x, int n_pml, double *&par_mat, 
											   double d_t, double d_x)
{
	//dimension of the computational space
    nx = n_x;
	dt = d_t;
	dx = d_x;
	n_PML = n_pml;

	n_L_mat = par_mat[0];//number of Lorentz terms 
	eps_inf = par_mat[1];
	mu_r = par_mat[2+3*n_L_mat];

	Hz_1D = (double *) calloc(nx,sizeof(double));
	if(!Hz_1D)
	{
		return FALSE;
	}
	
	Dy = (double *) calloc(nx-1,sizeof(double));
	if(!Dy)
	{
		Free_Mem_1D();
		return FALSE;
	}

	Ey_1D = (double *) calloc(nx-1,sizeof(double));
	if(!Ey_1D)
	{
		Free_Mem_1D();
		return FALSE;
	}

	Sy = Init_Matrix_2D<double>(nx-1,n_L_mat);
	if(!Sy)
	{
		Free_Mem_1D();
		return FALSE;
	}

	Sy_2 = Init_Matrix_2D<double>(nx-1,n_L_mat);
	if(!Sy_2)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	//////////////////////////
	K_E1_a = (double *) calloc(nx,sizeof(double));
	if(!K_E1_a)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	K_E1_b = (double *) calloc(nx,sizeof(double));
	if(!K_E1_b)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	K_E2_a = (double *) calloc(nx-1,sizeof(double));
	if(!K_E2_a)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	K_E2_b = (double *) calloc(nx-1,sizeof(double));
	if(!K_E2_b)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	//////////////////////////
	K_a = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_a)
	{
		Free_Mem_1D();
		return FALSE;
	}
	 
	K_b = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_b)
	{
		Free_Mem_1D();
		return FALSE;
	}
	
	K_c = (double *) calloc(n_L_mat,sizeof(double));
	if(!K_c)
	{
		Free_Mem_1D();
		return FALSE;
	}

	double am, delta_0, omega_0;
	
	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 TRUE;
}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Gaussian pulse
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Init_Gauss_1D(double H_0, double t_0, double t_w)
{
	source_type = 1;
	H0 = H_0;
	t0 = t_0;
	tw = t_w;
}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Sinusoidal plane wave
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Init_Sinus_1D(double H_0, double om, double Phi)
{
	source_type = 2;
	H0 = H_0;
	omega = om;
	phi = Phi;
}


//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Wave Packet (sinusoidal modulated Gauss pulse)
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Init_Wave_Packet_1D(double H_0, double t_0, double t_w, 
									   double om, double Phi)
{
	source_type = 3;
	H0 = H_0;
	t0 = t_0;
	tw = t_w;
	omega = om;
	phi = Phi;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Set_PML_Param_1D(double eps_r, double mu_r)
{
	for (i = 0; i<nx; i++)
	{
		K_E1_a[i] = 1.0;
		K_E1_b[i] = dt/(mu_0*mu_r*dx);

		if (i < nx-1)
		{
			K_E2_a[i] = 1.0;
			K_E2_b[i] = dt/dx;
		}			
	}

	//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_x, sigma_max, ka_x;
		
	sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx);
	//sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*dx*sqrt(mu_0*mu_r*eps_0*eps_r));

	for (i = 0; i<n_PML; i++)
	{
		sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_x            = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent);
		K_E1_a[nx-i-1]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		//K_E1_a[i]       = K_E1_a[nx-i-1];
		K_E1_b[nx-i-1]  = 2*eps_0*dt/(2*eps_0*ka_x+sigma_x*dt)/(mu_0*mu_r*dx);
		//K_E1_b[i]       = K_E1_b[nx-i-1];
		
		sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_x            = 1 + (ka_max - 1)*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		K_E2_a[nx-i-2]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		//K_E2_a[i]       = K_E2_a[nx-i-2];
		K_E2_b[nx-i-2]  = 2*eps_0*dt/(2*eps_0*ka_x + sigma_x*dt)/dx;
		//K_E2_b[i]       = K_E2_b[nx-i-2];
	}

}

//////////////////////////////////////////////////////////////////////
//Compute the Hz component
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::calc_Hz_1D(double t)
{
	switch (source_type)
	{
		case 1: //Gaussian pulse
			Hz_1D[0] = H0*exp( -pow( (t-t0)/tw ,2) ); 
			break;
		case 2: //Sinusoidal plane wave
			Hz_1D[0] = H0*cos( omega*t + phi);
			break;
		case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
			Hz_1D[0] = H0*cos( omega*t + phi)*exp( -pow( (t-t0)/tw ,2) );
	}
	 
	for (i = 1; i < nx - 1; i++)
	{
		Hz_1D[i] = K_E1_a[i]*Hz_1D[i] - K_E1_b[i]*(Ey_1D[i] - Ey_1D[i-1]);
	}

}

//////////////////////////////////////////////////////////////////////
//Compute the Ey component
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::calc_Ey_1D()
{
	int kk;
	double Sum_S;

	for (i = 0; i < nx - 1; i++)
	{
		double Dy_r = Dy[i];

		Dy[i] = K_E2_a[i]*Dy[i] - K_E2_b[i]*(Hz_1D[i+1] - Hz_1D[i]);

		Sum_S = 0;
		for (kk = 0; kk<n_L_mat; kk++)
		{
			 Sum_S = Sum_S + Sy[i][kk];
		}
		Ey_1D[i] = (Dy[i] - eps_0*Sum_S)/(eps_0*eps_inf);

		for (kk = 0; kk<n_L_mat; kk++)
		{
			Sy_r = Sy[i][kk]; 
			Sy[i][kk]   = K_a[kk] * Sy_r +
				 		  K_b[kk] * Sy_2[i][kk] + 
						  K_c[kk] * Ey_1D[i];
			Sy_2[i][kk] = Sy_r;
		}
	}

}

//////////////////////////////////////////////////////////////////////
//Free the allocated memory
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Free_Mem_1D()
{
	if (Hz_1D)
		free(Hz_1D);
	if (Dy)
		free(Dy);
	if (Ey_1D)
		free(Ey_1D);
	if (Sy)
		Sy = Free_Matrix_2D<double>(Sy);
	if (Sy_2)
		Sy_2 = Free_Matrix_2D<double>(Sy_2);
	if (K_E1_a)
		free(K_E1_a);
	if (K_E1_b)
		free(K_E1_b);
	if (K_E2_a)
		free(K_E2_a);
	if (K_E2_b)
		free(K_E2_b);
	if (K_a)
		free(K_a);
	if (K_b)
		free(K_b);
	if (K_c)
		free(K_c);
}

//////////////////////////////////////////////////////////////////////
//Access Hz_1D and Ey_1D
//////////////////////////////////////////////////////////////////////
void CFDTD_1D_HzEy_LORENTZ::Get_Data_1D(double *&X, double *&Y)
{
	X = Hz_1D;
	Y = Ey_1D;
}