fdtd_2d_tm.cpp

2维时域有限差分模拟

#include "StdAfx.h"
#include ".\fdtd_2d_tm.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_2D_TM::CFDTD_2D_TM(void)
{
	pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679;
	//permittivity of free space 
	eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]

	Fz = NULL; Ez = NULL; Hx = NULL; Gx = NULL; Hy = NULL; Gy = NULL;  

	K_Fz_a = NULL; K_Fz_b = NULL; K_Ez_a = NULL; K_Ez_b = NULL;
	K_Gx_a = NULL; K_Gx_b = NULL; K_Hx_a = NULL; K_Hx_b = NULL;
	K_Gy_a = NULL; K_Gy_b = NULL; K_Hy_a = NULL; K_Hy_b = NULL;

	//the interface between the total field-scattered field zone
	ll_1D_E = NULL; ll_1D_H = NULL;
	d_f_E = NULL; d_b_E = NULL; d_l_E = NULL; d_r_E = NULL;
	d_f_H = NULL; d_b_H = NULL; d_l_H = NULL; d_r_H = NULL;
	Ez_f_1D = NULL; Ez_b_1D = NULL; Ez_l_1D = NULL; Ez_r_1D = NULL;
	Hx_f_1D = NULL; Hx_b_1D = NULL; Hy_l_1D = NULL; Hy_r_1D = NULL;

	//Contains the followed field components
	Ez_Foll = NULL;	Hx_Foll = NULL;	 Hy_Foll = NULL;
	
	//to save the field components
	path_name_Ez = NULL; 
	path_name_Hx = NULL; 
	path_name_Hy = NULL; 

	Ez_1D = NULL; Hy_1D = NULL;

	#ifdef senddata
		ep = NULL;
	#endif
}

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

#ifdef senddata
	///////////////////////////////////////////////////////////////////////////////////////
	//Gets the identifier of the Matlab engine
	///////////////////////////////////////////////////////////////////////////////////////
	void CFDTD_2D_TM::set_MATLABengine(Engine *Ep)
	{
		ep = Ep;
	}
#endif

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_2D_TM::Init(int **&index, int n_x, int n_y, double **&mater, int n_mat, 
	  				   int n_pml, double d_t, double d_x, double d_y)
{
	Ind = index; //contains the indices of the material type. Dimension [nx][ny]

	//dimensions of the computational space
    nx = n_x;
	ny = n_y; 
	
	//contains the material parameters [eps_r mu_r]
	Mat = mater;
	n_Mat = n_mat;

	//dimension of the PML region
	n_PML = n_pml;
	
	dx = d_x;
	dy = d_y;
	dt = d_t; //the time step

	invdx = 1.0/dx;
	invdy = 1.0/dy;
	
	dtDIVdx = dt/dx;
	dtDIVdy = dt/dy;

	K_Ez_x = dt/(eps_0*dy);
	K_Ez_y = dt/(eps_0*dx);
	
	K_Hx = dt/(mu_0*dy);
	K_Hy = dt/(mu_0*dx);
	
	Fz = Init_Matrix_2D<double>(nx,ny);
	if(!Fz)
	{
		Free_Mem();
		return false;
	}

	Ez = Init_Matrix_2D<double>(nx,ny);
	if(!Ez)
	{
		Free_Mem();
		return false;
	}

	Gx = Init_Matrix_2D<double>(nx,ny-1);
	if(!Gx)
	{
		Free_Mem();
		return false;
	}

	Hx = Init_Matrix_2D<double>(nx,ny-1);
	if(!Hx)
	{
		Free_Mem();
		return false;
	}

	Gy = Init_Matrix_2D<double>(nx-1,ny);
	if(!Gy)
	{
		Free_Mem();
		return false;
	}

	Hy = Init_Matrix_2D<double>(nx-1,ny);
	if(!Hy)
	{
		Free_Mem();
		return false;
	}

	K_Fz_a = (double *) calloc(ny,sizeof(double));
	if(!K_Fz_a)
	{
		Free_Mem();
		return false;
	}

	K_Fz_b = (double *) calloc(ny,sizeof(double));
	if(!K_Fz_b)
	{
		Free_Mem();
		return false;
	}

	K_Ez_a = (double *) calloc(nx,sizeof(double));
	if(!K_Ez_a)
	{
		Free_Mem();
		return false;
	}

	K_Ez_b = (double *) calloc(nx,sizeof(double));
	if(!K_Ez_b)
	{
		Free_Mem();
		return false;
	}

	K_Gx_a = (double *) calloc(ny-1,sizeof(double));
	if(!K_Gx_a)
	{
		Free_Mem();
		return false;
	}

	K_Gx_b = (double *) calloc(ny-1,sizeof(double));
	if(!K_Gx_b)
	{
		Free_Mem();
		return false;
	}

	K_Hx_a = (double *) calloc(nx,sizeof(double));
	if(!K_Hx_a)
	{
		Free_Mem();
		return false;
	}

	K_Hx_b = (double *) calloc(nx,sizeof(double));
	if(!K_Hx_b)
	{
		Free_Mem();
		return false;
	}

	K_Gy_a = (double *) calloc(nx-1,sizeof(double));
	if(!K_Gy_a)
	{
		Free_Mem();
		return false;
	}

	K_Gy_b = (double *) calloc(nx-1,sizeof(double));
	if(!K_Gy_b)
	{
		Free_Mem();
		return false;
	}

	K_Hy_a = (double *) calloc(ny,sizeof(double));
	if(!K_Hy_a)
	{
		Free_Mem();
		return false;
	}

	K_Hy_b = (double *) calloc(ny,sizeof(double));
	if(!K_Hy_b)
	{
		Free_Mem();
		return false;
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init and set the total field scattred field matrices
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::Init_TS_Param()
{
	#ifdef senddata
		// Array to store data ready to transfer to MATLAB
		mxArray *X_ll_1D = NULL;
		mxArray *X_d_fb = NULL, *X_d_lr = NULL;
		
		X_ll_1D = mxCreateDoubleMatrix(1, l_1D, mxREAL);
		X_d_fb  = mxCreateDoubleMatrix(1, l_x, mxREAL);
		X_d_lr  = mxCreateDoubleMatrix(1, l_y, mxREAL);
	#endif

	int i, j;

	ll_1D_E = (double *) calloc(l_1D,sizeof(double));
	if (!ll_1D_E)
	{
		Free_Mem();
		return 1;
	}

	ll_1D_H = (double *) calloc(l_1D-1,sizeof(double));
	if (!ll_1D_H)
	{
		Free_Mem();
		return 1;
	}

	d_f_E = (double *) calloc(l_x,sizeof(double));
	if (!d_f_E)
	{
		Free_Mem();
		return 1;
	}
	d_b_E = (double *) calloc(l_x,sizeof(double));
	if (!d_b_E)
	{
		Free_Mem();
		return 1;
	}
	d_l_E = (double *) calloc(l_y,sizeof(double));
	if (!d_l_E)
	{
		Free_Mem();
		return 1;
	}
	d_r_E = (double *) calloc(l_y,sizeof(double));
	if (!d_r_E)
	{
		Free_Mem();
		return 1;
	}

	Ez_f_1D = (double *) calloc(l_x,sizeof(double));
	if (!Ez_f_1D)
	{
		Free_Mem();
		return 1;
	}
	Ez_b_1D = (double *) calloc(l_x,sizeof(double));
	if (!Ez_b_1D)
	{
		Free_Mem();
		return 1;
	}
	Ez_l_1D = (double *) calloc(l_y,sizeof(double));
	if (!Ez_l_1D)
	{
		Free_Mem();
		return 1;
	}
	Ez_r_1D = (double *) calloc(l_y,sizeof(double));
	if (!Ez_r_1D)
	{
		Free_Mem();
		return 1;
	}

	d_f_H = (double *) calloc(l_x,sizeof(double));
	if (!d_f_H)
	{
		Free_Mem();
		return 1;
	}
	d_b_H = (double *) calloc(l_x,sizeof(double));
	if (!d_b_H)
	{
		Free_Mem();
		return 1;
	}
	d_l_H = (double *) calloc(l_y,sizeof(double));
	if (!d_l_H)
	{
		Free_Mem();
		return 1;
	}
	d_r_H = (double *) calloc(l_y,sizeof(double));
	if (!d_r_H)
	{
		Free_Mem();
		return 1;
	}

	Hx_f_1D = (double *) calloc(l_x,sizeof(double));
	if (!Hx_f_1D)
	{
		Free_Mem();
		return 1;
	}
	Hx_b_1D = (double *) calloc(l_x,sizeof(double));
	if (!Hx_b_1D)
	{
		Free_Mem();
		return 1;
	}
	Hy_l_1D = (double *) calloc(l_y,sizeof(double));
	if (!Hy_l_1D)
	{
		Free_Mem();
		return 1;
	}
	Hy_r_1D = (double *) calloc(l_y,sizeof(double));
	if (!Hy_r_1D)
	{
		Free_Mem();
		return 1;
	}
	
	for (i = 0; i < l_1D-1; i++)
	{
		ll_1D_E[i] = (i - n_PML - 5)*d;
		ll_1D_H[i] = (i + 0.5 - n_PML - 5)*d;
	}
	ll_1D_E[l_1D-1] = (l_1D-1 - n_PML - 5)*d;

	// 0 <= phi <= 90
	if ( sin_phi >= 0 && cos_phi >= 0 )
	{
		for (i = 0; i < l_x; i++)
		{
			d_f_E[i] = i*dx*cos_phi;
			d_b_E[i] = i*dx*cos_phi + (l_y-1)*dy*sin_phi;
			d_f_H[i] = i*dx*cos_phi - 0.5*dy*sin_phi;
			d_b_H[i] = i*dx*cos_phi + (l_y-1+0.5)*dy*sin_phi;
		}

		for (j = 0; j < l_y; j++)
		{
			d_l_E[j] = j*dy*sin_phi;
			d_r_E[j] = (l_x-1)*dx*cos_phi  + j*dy*sin_phi;
			d_l_H[j] = -0.5*dx*cos_phi + j*dy*sin_phi;
			d_r_H[j] = (l_x-1+0.5)*dx*cos_phi + j*dy*sin_phi;
		}
	}

	// 90 < phi <= 180
	if ( sin_phi >= 0 && cos_phi < 0 )
	{
		for (i = 0; i < l_x; i++)
		{
			d_f_E[i] = (-l_x+1+i)*dx*cos_phi;
			d_b_E[i] = (-l_x+1+i)*dx*cos_phi + (l_y-1)*dy*sin_phi;
			d_f_H[i] = (-l_x+1+i)*dx*cos_phi - 0.5*dy*sin_phi;
			d_b_H[i] = (-l_x+1+i)*dx*cos_phi + (l_y-1+0.5)*dy*sin_phi;
		}

		for (j = 0; j < l_y; j++)
		{
			d_l_E[j] = (-l_x+1)*dx*cos_phi + j*dy*sin_phi;
			d_r_E[j] = j*dy*sin_phi;
			d_l_H[j] = (-l_x+1-0.5)*dx*cos_phi + j*dy*sin_phi;
			d_r_H[j] = 0.5*dx*cos_phi + j*dy*sin_phi;
		}
	}

	// 180 <= phi <= 270
	if ( sin_phi < 0 && cos_phi <= 0 )
	{
		for (i = 0; i < l_x; i++)
		{
			d_f_E[i] = (-l_x+1+i)*dx*cos_phi + (-l_y+1)*dy*sin_phi;
			d_b_E[i] = (-l_x+1+i)*dx*cos_phi;
			d_f_H[i] = (-l_x+1+i)*dx*cos_phi + (-l_y+1-0.5 )*dy*sin_phi;
			d_b_H[i] = (-l_x+1+i)*dx*cos_phi + 0.5*dy*sin_phi;
		}

		for (j = 0; j < l_y; j++)
		{
			d_l_E[j] = (-l_x+1)*dx*cos_phi + (-l_y+1+j)*dy*sin_phi;
			d_r_E[j] = (-l_y+1+j)*dy*sin_phi;
			d_l_H[j] = (-l_x+1-0.5)*dx*cos_phi + (-l_y+1+j)*dy*sin_phi;
			d_r_H[j] = 0.5*dx*cos_phi + (-l_y+1+j)*dy*sin_phi;
		}
	}

	// 270 < phi < 360
	if ( sin_phi < 0 && cos_phi > 0 )
	{
		for (i = 0; i < l_x; i++)
		{
			d_f_E[i] = i*dx*cos_phi + (-l_y+1)*dy*sin_phi;
			d_b_E[i] = i*dx*cos_phi;
			d_f_H[i] = i*dx*cos_phi + (-l_y+1-0.5 )*dy*sin_phi;
			d_b_H[i] = i*dx*cos_phi + 0.5*dy*sin_phi;
		}

		for (j = 0; j < l_y; j++)
		{
			d_l_E[j] = (-l_y+1+j)*dy*sin_phi;
			d_r_E[j] = (l_x-1)*dx*cos_phi + (-l_y+1+j)*dy*sin_phi;
			d_l_H[j] = -dx*0.5*cos_phi + (-l_y+1+j)*dy*sin_phi;
			d_r_H[j] = (l_x-1+0.5)*dx*cos_phi + (-l_y+1+j)*dy*sin_phi;
		}
	}

	#ifdef senddata // Transfer to MATLAB
			memcpy((void *)mxGetPr(X_ll_1D),ll_1D_E, l_1D*sizeof(double));
			engPutVariable(ep,"ll_1D_E",X_ll_1D); 

			memcpy((void *)mxGetPr(X_d_fb),d_f_E, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_f_E",X_d_fb); 

			memcpy((void *)mxGetPr(X_d_fb),d_b_E, l_y*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_b_E",X_d_fb); 

			memcpy((void *)mxGetPr(X_d_fb),d_f_H, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_f_H",X_d_fb); 

			memcpy((void *)mxGetPr(X_d_fb),d_b_H, l_y*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_b_H",X_d_fb); 

			memcpy((void *)mxGetPr(X_d_lr),d_l_E, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_l_E",X_d_lr); 

			memcpy((void *)mxGetPr(X_d_lr),d_r_E, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_r_E",X_d_lr); 

			memcpy((void *)mxGetPr(X_d_lr),d_l_H, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_l_H",X_d_lr); 

			memcpy((void *)mxGetPr(X_d_lr),d_r_H, l_x*sizeof(double));
			// Transfer to MATLAB
			engPutVariable(ep,"d_r_H",X_d_lr); 
	#endif


	#ifdef senddata
		mxDestroyArray(X_ll_1D);
		mxDestroyArray(X_d_fb);
		mxDestroyArray(X_d_lr);
	#endif

	return 0;
}
///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::Init_PML_Param()
{
	int i, j;

	for (i = 0; i<nx-1; i++)
	{
		K_Ez_a[i] = 1.0;
		K_Ez_b[i] = 1.0/eps_0;

		K_Hx_a[i] = 1.0/mu_0;
		K_Hx_b[i] = 1.0/mu_0;

		K_Gy_a[i] = 1.0;
		K_Gy_b[i] = dt/dx;
	}
	K_Ez_a[nx-1] = 1.0;
	K_Ez_b[nx-1] = 1.0/eps_0;

	K_Hx_a[nx-1] = 1.0/mu_0;
	K_Hx_b[nx-1] = 1.0/mu_0;
	

	for (j = 0; j<ny-1; j++)
	{
		K_Fz_a[j] = 1.0;
		K_Fz_b[j] = dt;

		K_Hy_a[j] = 1.0/mu_0;
		K_Hy_b[j] = 1.0/mu_0;

		K_Gx_a[j] = 1.0;
		K_Gx_b[j] = dt/dy;
	}
	K_Fz_a[ny-1] = 1.0;
	K_Fz_b[ny-1] = dt;

	K_Hy_a[ny-1] = 1.0/mu_0;
	K_Hy_b[ny-1] = 1.0/mu_0;
	

	//PML parameters
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];

	double ka_max = 2.0;
	int exponent = 6;
	double R_err = 1.0e-16;
	double eta = sqrt( mu_0*mu_r/(eps_0*eps_r) );

	double sigma_x, sigma_y, sigma_max, ka_x, ka_y;
		
	sigma_max= -(exponent+1.0)*log(R_err)/(2.0*eta*n_PML*dx);

	for (i = 0; i<n_PML; i++)
	{
		sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_x            = 1.0 + (ka_max - 1.0)*pow( (n_PML-i)/((double) n_PML) ,exponent);

		K_Ez_a[i]       = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
		K_Ez_a[nx-i-1]  = K_Ez_a[i];

		K_Ez_b[i]       = 2.0/(2.0*eps_0*ka_x + sigma_x*dt);
		K_Ez_b[nx-i-1]  = K_Ez_b[i];

		K_Hx_a[i]       = (2.0*eps_0*ka_x + sigma_x*dt)/(2.0*eps_0*mu_0);
		K_Hx_a[nx-i-1]  = K_Hx_a[i];

		K_Hx_b[i]       = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*mu_0);
		K_Hx_b[nx-i-1]  = K_Hx_b[i];

		sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_x            = 1.0 + (ka_max - 1.0)*pow( (n_PML - i - 0.5)/n_PML ,exponent);

        K_Gy_a[i]       = (2.0*eps_0*ka_x - sigma_x*dt)/(2.0*eps_0*ka_x + sigma_x*dt);
		K_Gy_a[nx-i-2]  = K_Gy_a[i];

		K_Gy_b[i]       = 2.0*eps_0*dt/(2.0*eps_0*ka_x + sigma_x*dt)/dx;
		K_Gy_b[nx-i-2]  = K_Gy_b[i];
	}