fdtd_2d_te.cpp

fdtd的2d算法

void FDTD_2D_TE::Init_Gauss_Point_Source(int n_x_p, int n_y_p, double h0, double t_0, 
										 double t_w)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   jel_source_type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Sinusoidal
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Init_Sinus_Point_Source(int n_x_p, int n_y_p, double h0, double om, 
										 double phase)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   omega = om;
   phi = phase;
   jel_source_type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Sinusoidal Modulated Gauss Pulse 
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Init_PulseGS_Point_Source(int n_x_p, int n_y_p, double h0, double t_0, 
							 			   double t_w, double om, double phase)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   omega = om;
   phi = phase;
   jel_source_type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Soft Point Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Point_Source(double time)
{
	switch (jel_source_type)
	{
		case 1: //Gaussian pulse
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] + H0*exp( -pow( (time-t0)/tw ,2) ); 
			break;
		case 2: //Sinusoidal plane wave
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] = H0*cos( omega*time + phi);
			break;
		case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] = H0*cos( omega*time + phi)*
				                                     exp( -pow( (time-t0)/tw ,2) );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Line Gauss Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Init_Gauss_Line_Source(int **&Co, int n_Co, double h0, double t_0, 
										 double t_w)
{
   Coord = Co;
   n_Coord = n_Co;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   jel_source_type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Soft Point Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Line_Source(double time)
{
	if (jel_source_type == 1) //Gauss pulse
	{
		for (i = 0; i< n_Coord; i++)
			Hz[Coord[i][0]][Coord[i][1]] = Hz[Coord[i][0]][Coord[i][1]] +
			                                       H0*exp( -pow((time-t0)/tw,2) );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Total field - Scattered field formulation
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Init_TotFScatF(int n_xa, int n_xb, int n_ya, int n_yb, double H0, 
								double t0, double tw, double omega, double phi, int jel)
{
	n_x_a = n_xa;
	n_x_b = n_xb;
	n_y_a = n_ya;
	n_y_b = n_yb;

	Init_Main_Param_1D(n_x_b+5, n_PML, Mater[0][0], Mater[0][1], dt, dx);
	switch (jel)
	{
		case 1:
			Init_Gauss_1D(H0, t0, tw);
			break;
		case 2:
			Init_Sinus_1D(H0, omega, phi);
			break;
		case 3:
			Init_Wave_Packet_1D(H0, t0, tw, omega, phi);
	}
	Set_PML_Param_1D();
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Hz
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::incident_Hz()
{
	for (j = n_y_a; j <= n_y_b; j++)
	{
		Hz[n_x_a][j] = Hz[n_x_a][j] + 
						dt/(mu_0*Mater[Index[n_x_a][j]][1]*dx)*Ey_1D[0];

		Hz[n_x_b][j] = Hz[n_x_b][j] - 
						dt/(mu_0*Mater[Index[n_x_b][j]][1]*dx)*Ey_1D[n_x_b-n_x_a+1];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Ex
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::incident_Ex()
{
	for (i = n_x_a; i <= n_x_b; i++)
	{
		Ex[i][n_y_a-1] = Ex[i][n_y_a-1] -
						  dt/(eps_0*Mater[Index[i][n_y_a-1]][0]*dy)*Hz_1D[i-n_x_a+1];
		Ex[i][n_y_b] =   Ex[i][n_y_b] +
						  dt/(eps_0*Mater[Index[i][n_y_b]][0]*dy)*Hz_1D[i-n_x_a+1];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Ey
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::incident_Ey()
{
	for (j = n_y_a; j <= n_y_b; j++)
	{
		Ey[n_x_a-1][j] = Ey[n_x_a-1][j] +
						  dt/(eps_0*Mater[Index[n_x_a-1][j]][0]*dx)*Hz_1D[1];
		Ey[n_x_b][j] = Ey[n_x_b][j] -
						  dt/(eps_0*Mater[Index[n_x_b][j]][0]*dx)*Hz_1D[n_x_b-n_x_a+1];
	}

}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Free_Mem()
{
	if(Fz)
		Fz = Free_Matrix_2D<double>(Fz);
	if(Hz)
		Hz = Free_Matrix_2D<double>(Hz);
	if(Ex)
		Ex = Free_Matrix_2D<double>(Ex);
	if(Gx)
		Gx = Free_Matrix_2D<double>(Gx);
	if(Ey)
		Ey = Free_Matrix_2D<double>(Ey);
    if(Gy)
		Gy = Free_Matrix_2D<double>(Gy);
	
	if(K_E1_a)
	{
	    free(K_E1_a); 	
		K_E1_a = NULL;
	}
	if(K_E1_b)
	{
		free(K_E1_b);
		K_E1_b = NULL;
	}
	if(K_E2_a)
	{
		free(K_E2_a); 	
		K_E2_a = NULL;
	}
	if(K_E2_b)
	{
		free(K_E2_b);
		K_E2_b = NULL;
	}
	if(K_E3_a)
	{
		free(K_E3_a); 	
		K_E3_a = NULL;
	}
	if(K_E3_b)
	{
		free(K_E3_b);
		K_E3_b = NULL;
	}
	if(K_E4_a)
	{
		free(K_E4_a); 	
		K_E4_a = NULL;
	}
	if(K_E4_b)
	{
		free(K_E4_b);
		K_E4_b = NULL;
	}
	if(K_E5_a)
	{
		free(K_E5_a);
		K_E5_a = NULL;
	}
	if(K_E5_b)
	{
		free(K_E5_b);
		K_E5_b = NULL;
	}
	if(K_E6_a)
	{
		free(K_E6_a); 
		K_E6_a = NULL;
	}
	if(K_E6_b)
	{
		free(K_E6_b);
		K_E6_b = NULL;
	}
	
	if(Hz_Foll)
		Hz_Foll = Free_Matrix_2D<double>(Hz_Foll);
	if(Ex_Foll)
		Ex_Foll = Free_Matrix_2D<double>(Ex_Foll);
	if(Ey_Foll)
		Ey_Foll = Free_Matrix_2D<double>(Ey_Foll);
}

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

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

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

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

	return TRUE;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the Followed Energy components
///////////////////////////////////////////////////////////////////////////////////////
BOOL FDTD_2D_TE::Init_Followed_W(int **Ind_Followed_W, int n_W, int n_t)
{
	Ind_Foll_W = Ind_Followed_W;
	length_Ind_Foll_W = n_W;
	n_tot = n_t;

	W_Foll = Init_Matrix_2D<double>(n_W,n_t);
	if(!W_Foll)	
	{
		Free_Mem();
		return FALSE;
	}

	return TRUE;
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Set_Data_Followed(int n_t)
{
	for (i = 0; i<length_Ind_Foll; i++)
	{
		Hz_Foll[i][n_t] = Hz[Ind_Foll[i][0]][Ind_Foll[i][1]];
		Ex_Foll[i][n_t] = Ex[Ind_Foll[i][0]][Ind_Foll[i][1]];
		Ey_Foll[i][n_t] = Ey[Ind_Foll[i][0]][Ind_Foll[i][1]];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed Energy
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Set_Data_Followed_W(int n_t)
{
	for (i = 0; i<length_Ind_Foll_W; i++)
	{
		W_Foll[i][n_t] = 0.5*( mu_0*Mater[Index[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]][1]*
			                        Hz[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]*
			                        Hz[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]] +
							   eps_0*Mater[Index[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]][0]*
									Ex[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]*
			                        Ex[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]] +
							   eps_0*Mater[Index[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]][0]*
									Ey[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]]*
			                        Ey[Ind_Foll_W[i][0]][Ind_Foll_W[i][1]] );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Get_Data_Followed(double **&A, double **&B, double **&C)
{
	A = Hz_Foll;
	B = Ex_Foll;
	C = Ey_Foll;
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed Energy
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Get_Data_Followed_W(double **&A)
{
	A = W_Foll;
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Get_Data(double **&A, double **&B, double **&C)
{
	A = Hz;
	B = Ex;
	C = Ey;
}