fdtd_2d_tm.cpp

2维时域有限差分模拟

	{
		mu_r = Mat[Ind[i][n_TS_ya-1]][1];
		Hx[i][n_TS_ya-1] += K_Hx*Ez_f_1D[i-n_TS_xa]/mu_r;

		mu_r = Mat[Ind[i][n_TS_yb]][1];
		Hx[i][n_TS_yb] -=   K_Hx*Ez_b_1D[i-n_TS_xa]/mu_r;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hy on total-field scattered field boundary - x direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::TS_Hy()
{
	int j;
	double mu_r;
	
	for (j = n_TS_ya; j <= n_TS_yb; j++)
	{
		mu_r = Mat[Ind[n_TS_xa-1][j]][1];
		Hy[n_TS_xa-1][j] -= K_Hy*Ez_l_1D[j-n_TS_ya]/mu_r;

		mu_r = Mat[Ind[n_TS_xb][j]][1];
		Hy[n_TS_xb][j] += K_Hy*Ez_r_1D[j-n_TS_ya]/mu_r;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Collect the function calls for the FDTD algorithm
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::fdtd_marching(int n_max_steep)
{
	double time = 0;
	int iter;

	#ifdef senddata //Debugging with Matlab
		// Array to store data ready to transfer to MATLAB
		mxArray *X_Ez_1D = NULL, *X_Hy_1D = NULL;
		mxArray *X_fb_1D = NULL, *X_lr_1D = NULL;
		
		X_Ez_1D = mxCreateDoubleMatrix(1, l_1D, mxREAL);
		X_Hy_1D = mxCreateDoubleMatrix(1, l_1D-1, mxREAL);

		//i
		X_fb_1D = mxCreateDoubleMatrix(1, l_x, mxREAL);
		//j
		X_lr_1D = mxCreateDoubleMatrix(1, l_y, mxREAL);
	#endif

	for (iter = 0; iter < n_max_steep; iter++)
	{
		time += dt;
		
		//1D wave propagation
		if (jel_plane_wave == 1) 
		{
			fdtd_1D.calc_Hy();
			//front face i=i0...i1, j=j0 
			interp_1D(ll_1D_H,l_1D-1,Hy_1D,d_f_H,l_x,Hx_f_1D);
			//back face i=i0...i1, j=j1 
			interp_1D(ll_1D_H,l_1D-1,Hy_1D,d_b_H,l_x,Hx_b_1D);
			//left face i=i0, j=j0:j1 
			interp_1D(ll_1D_H,l_1D-1,Hy_1D,d_l_H,l_y,Hy_l_1D);
			//rigth face i=i1, j=j0:j1 
			interp_1D(ll_1D_H,l_1D-1,Hy_1D,d_r_H,l_y,Hy_r_1D);

			fdtd_1D.calc_Ez(time);
			//front face i=i0...i1, j=j0 
			interp_1D(ll_1D_E,l_1D,Ez_1D,d_f_E,l_x,Ez_f_1D);
			//back face i=i0...i1, j=j1 
			interp_1D(ll_1D_E,l_1D,Ez_1D,d_b_E,l_x,Ez_b_1D);
			//left face i=i0, j=j0:j1 
			interp_1D(ll_1D_E,l_1D,Ez_1D,d_l_E,l_y,Ez_l_1D);
			//rigth face i=i1, j=j0:j1 
			interp_1D(ll_1D_E,l_1D,Ez_1D,d_r_E,l_y,Ez_r_1D);
			
			#ifdef senddata // Transfer to MATLAB for debugging
				memcpy((void *)mxGetPr(X_Ez_1D),Ez_1D, l_1D*sizeof(double));
				engPutVariable(ep,"Ez_1D",X_Ez_1D); 

				memcpy((void *)mxGetPr(X_Hy_1D),Hy_1D, (l_1D-1)*sizeof(double));
				engPutVariable(ep,"Hy_1D",X_Hy_1D); 

				memcpy((void *)mxGetPr(X_fb_1D),Ez_f_1D, l_x*sizeof(double));
				engPutVariable(ep,"Ez_f_1D",X_fb_1D); 

				memcpy((void *)mxGetPr(X_fb_1D),Ez_b_1D, l_x*sizeof(double));
				engPutVariable(ep,"Ez_b_1D",X_fb_1D); 

				memcpy((void *)mxGetPr(X_lr_1D),Ez_l_1D, l_y*sizeof(double));
				engPutVariable(ep,"Ez_l_1D",X_lr_1D); 

				memcpy((void *)mxGetPr(X_lr_1D),Ez_r_1D, l_y*sizeof(double));
				engPutVariable(ep,"Ez_r_1D",X_lr_1D); 

				memcpy((void *)mxGetPr(X_fb_1D),Hx_f_1D, l_x*sizeof(double));
				engPutVariable(ep,"Hx_f_1D",X_fb_1D); 

				memcpy((void *)mxGetPr(X_fb_1D),Hx_b_1D, l_x*sizeof(double));
				engPutVariable(ep,"Hx_b_1D",X_fb_1D); 

				memcpy((void *)mxGetPr(X_lr_1D),Hy_l_1D, l_y*sizeof(double));
				engPutVariable(ep,"Hy_l_1D",X_lr_1D); 

				memcpy((void *)mxGetPr(X_lr_1D),Hy_r_1D, l_y*sizeof(double));
				engPutVariable(ep,"Hy_r_1D",X_lr_1D); 
			#endif

		}

		calc_Ez();
		if (jel_plane_wave == 1) 
		{
		  TS_Ez();
		}

		//point source
		if (jel_plane_wave == 0) 
		{
			ptSource_J(Ez, Par_ptSource, time);
		}
		
		calc_Hx();
		calc_Hy();
		if (jel_plane_wave == 1) 
		{
			TS_Hx();
			TS_Hy();
		}

		//the followed field components
	    Set_Data_Followed(iter);
		if (jel_Save == 1 && iter%nr_save == 1 && iter>10000 && iter<50000)
		{
			Save_Field(iter);
		}
		
		if (iter%1000 == 1)
			cout<< iter << endl;
			
	}
	
	#ifdef senddata
		mxDestroyArray(X_Ez_1D);
		mxDestroyArray(X_Hy_1D);
		mxDestroyArray(X_fb_1D);
		mxDestroyArray(X_lr_1D);
	#endif

	return iter; //the number of executed steeps
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_2D_TM::Init_Followed(int **Ind_Followed, int n, int n_t)
{
	Ind_Foll = Ind_Followed;
	length_Ind_Foll = n;
	n_tot = n_t;

	Hx_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hx_Foll)	
	{
		Free_Mem();
		return false;
	}

	Hy_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hy_Foll)	
	{
		Free_Mem();
		return false;
	}

	Ez_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ez_Foll)	
	{
		Free_Mem();
		return false;
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::Set_Data_Followed(int n_t)
{
	int i, i1, j1;
	for (i = 0; i<length_Ind_Foll; i++)
	{
		i1 = Ind_Foll[i][0];
		j1 = Ind_Foll[i][1];

		Ez_Foll[i][n_t] = Ez[i1][j1];
		Hx_Foll[i][n_t] = Hx[i1][j1];
		Hy_Foll[i][n_t] = Hy[i1][j1];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::Get_Data_Followed(double **&hx, double **&hy, double **&ez)
{
	hx = Hx_Foll;
	hy = Hy_Foll;
	ez = Ez_Foll;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init to save the field components
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_2D_TM::Init_Save_Field(int nxa, int nxb, int nya, int nyb, 
				  			      int nr_Save, char *path)
{
	n_x_a = nxa; 
	n_x_b = nxb;
	n_y_a = nya; 
	n_y_b = nyb;
	
	
	nr_save = nr_Save;
	jel_Save = 1;

	path_name_Ez_1D= (char *) calloc(512,sizeof(char)); 
	if (!path_name_Ez_1D)
	{
		Free_Mem();
		return false;
	}
	
	path_name_Ez = (char *) calloc(512,sizeof(char));
	if (!path_name_Ez)
	{
		Free_Mem();
		return false;
	}
	
	path_name_Hx = (char *) calloc(512,sizeof(char)); 
	if (!path_name_Hx)
	{
		Free_Mem();
		return false;
	}
	
	path_name_Hy = (char *) calloc(512,sizeof(char));
	if (!path_name_Hy)
	{
		Free_Mem();
		return false;
	}
	
	strcpy(path_name_Ez,path);
    strcat(path_name_Ez,"/Ez");
	
	strcpy(path_name_Hx,path);
    strcat(path_name_Hx,"/Hx");
	
	strcpy(path_name_Hy,path);
    strcat(path_name_Hy,"/Hy");

	strcpy(path_name_Ez_1D,path);
    strcat(path_name_Ez_1D,"/Ez_1D");
	
	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save the field components
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::Save_Field(int iter)
{
	/*if (jel_plane_wave ==1) 
	{
		save_1D(Ez_1D,0,l_1D,iter,path_name_Ez_1D);
	}*/

	//if (!save_2D_binary(Ez,nx,ny,iter,path_name_Ez))
	if (!save_2D(Ez,n_x_a,n_x_b,n_y_a,n_y_b,iter,path_name_Ez))
		return -1;
	
	/*if (!save_2D(Hx,n_x_a,n_x_b,n_y_a,n_y_b-1,iter,path_name_Hx))
		return -1;

	if (!save_2D(Hy,n_x_a,n_x_b-1,n_y_a,n_y_b,iter,path_name_Hy))
		return -1;*/

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Interpolate according to the ordered values of x_i
///////////////////////////////////////////////////////////////////////////////////////
/*int CFDTD_2D_TM::interp_1D(double *x, int length_x, double *f, double *x_i, 
						   int length_x_i, double *g)
{ 
  //x - the coordinates
  //f - the corresponding function values
  //x_i - the points where I will determine the values of f with linear interpolation 
	
	int i, j, jel;

	if ( (x_i[1] - x_i[0]) > 0)
	{ 
		jel = 1;
	}
	else
	{
		jel = 0;
	}

	//if (x_i[0]<x[0] || x_i[length_x_i-1]>x[length_x-1])
	//	return 0;
	//else
	//{ 	
		if (jel == 1)
		{
			int x_r_i = 0;
			for (i = 0; i<length_x_i; i++)
			{
				for(j = x_r_i; j<length_x-1; j++)
				{
					if ( x[j] <= x_i[i] && x[j+1] > x_i[i])
					{
						x_r_i = j;
						g[i] = ( f[j]*(x[j+1]-x_i[i]) + f[j+1]*(x_i[i]-x[j]) )
							/(x[j+1]-x[j]);
						break;
					}
				}
			}
		}
		else
		{
			for (i = 0; i<length_x_i; i++)
			{
				for(j = 0; j< length_x-1; j++)
				{
					if ( x[j] <= x_i[i] && x[j+1] > x_i[i])
					{
						g[i] = ( f[j]*(x[j+1]-x_i[i]) + f[j+1]*(x_i[i]-x[j]) )
							/(x[j+1]-x[j]);
						break;
					}
				}
			}
		//}
		return 1;
	//}
}
*/
int CFDTD_2D_TM::interp_1D(double *x, int length_x, double *f, double *x_i, int length_x_i, double *g)
{
    int i, j; 
	int n_x_i, n_x;
	n_x_i = length_x_i;
	n_x = length_x;

	for (i = 0; i < n_x_i; i++)
	{
		if( (x_i[i] >= x[0]) && (x_i[i] <= x[n_x-1]))
		{
			for(j = 1; j < n_x; j++)
			{
				if(x_i[i]<=x[j])
				{
					g[i] = (x_i[i]-x[j-1]) / (x[j]-x[j-1]) * (f[j]-f[j-1]) + f[j-1];
					break;
				}
			}
		}
		else
		{
			//cannot interpolate
			return 1;
		}
	}

    return 0;
} 

///////////////////////////////////////////////////////////////////////////////////////
//calculates the phase velocity with Newton_Rhapson scheme
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::phase_velocity(double dx, double dy, double dt, double eps_r,
								double mu_r, double om, double phi, double *phase_vel)
{
	//permittivity of free space 
	double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]


	double vel_light = 1.0/sqrt(eps_0*eps_r*mu_0*mu_r);
	
	double dx2 = dx*dx;
	double dy2 = dy*dy;

	double A = dx*cos(phi)/2.0;
	double B = dy*sin(phi)/2.0;
	
	double C = sin(om*dt/2.0)/(vel_light*dt);
	C = C*C;
	double c_a, c_b;

	double k_r = 0.0;
	double k = om/vel_light;
	while ( abs(k-k_r) > 1e-12 )
	{
		k_r = k;
		c_a = sin(A*k)/dx;
		c_b = sin(B*k)/dy;

		k -= ( c_a*c_a + c_b*c_b - C)/( A*sin(2.0*A*k)/dx2 + B*sin(2.0*B*k)/dy2 );
	}

	*phase_vel = om/k;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::Free_Mem()
{
	if(Fz)
		Fz = Free_Matrix_2D<double>(Fz);
	if(Ez)
		Ez = Free_Matrix_2D<double>(Ez);
	if(Hx)
		Hx = Free_Matrix_2D<double>(Hx);
	if(Gx)
		Gx = Free_Matrix_2D<double>(Gx);
	if(Hy)
		Hy = Free_Matrix_2D<double>(Hy);
	if(Gy)
		Gy = Free_Matrix_2D<double>(Gy);
	
	
	if (path_name_Ez)
	{
		free(path_name_Ez);
		path_name_Ez = NULL;
	}

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

	if (path_name_Hy)
	{
		free(path_name_Hy);
		path_name_Hy = NULL;
	}
	
	//The followed field components
	if(Hx_Foll)
		Hx_Foll = Free_Matrix_2D<double>(Hx_Foll);
	if(Hy_Foll)
		Hy_Foll = Free_Matrix_2D<double>(Hy_Foll);
	if(Ez_Foll)
		Ez_Foll = Free_Matrix_2D<double>(Ez_Foll);

	//total field - scattered field formulation with arbitrary incident angle
	if (ll_1D_E)
	{
		free(ll_1D_E);
		ll_1D_E = NULL;
	}

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

	if (d_f_E)
	{
		free(d_f_E);
		d_f_E = NULL;
	}
	if (d_b_E)
	{
		free(d_b_E);
		d_b_E = NULL;
	}
	if (d_l_E)
	{
		free(d_l_E);
		d_l_E = NULL;
	}
	if (d_r_E)
	{
		free(d_r_E);
		d_r_E = NULL;
	}

	if (Ez_f_1D)
	{
		free(Ez_f_1D);
		Ez_f_1D = NULL;
	}
	if (Ez_b_1D)
	{
		free(Ez_b_1D);
		Ez_b_1D = NULL;
	}
	if (Ez_l_1D)
	{
		free(Ez_l_1D);
		Ez_l_1D = NULL;
	}
	if (Ez_r_1D)
	{
		free(Ez_r_1D);
		Ez_r_1D = NULL;
	}

	if (d_f_H)
	{
		free(d_f_H);
		d_f_H = NULL;
	}
	if (d_b_H)
	{
		free(d_b_H);
		d_b_H = NULL;
	}
	if (d_l_H)
	{
		free(d_l_H);
		d_l_H = NULL;
	}
	if (d_r_H)
	{
		free(d_r_H);
		d_r_H = NULL;
	}

	if (Hx_f_1D)
	{
		free(Hx_f_1D);
		Hx_f_1D = NULL;
	}
	if (Hx_b_1D)
	{
		free(Hx_b_1D);
		Hx_b_1D = NULL;
	}
	if (Hy_l_1D)
	{
		free(Hy_l_1D);
		Hy_l_1D = NULL;
	}
	if (Hy_r_1D)
	{
		free(Hy_r_1D);
		Hy_r_1D = NULL;
	}
}