fdtd_2d_tm.cpp

2维时域有限差分模拟

	sigma_max = -(exponent+1.0)*log(R_err)/(2.0*eta*n_PML*dy);
	for (j = 0; j<n_PML; j++)
	{
		sigma_y         = sigma_max*pow( (n_PML - j)/((double) n_PML) ,exponent);
		ka_y            = 1.0 + (ka_max - 1.0)*pow( (n_PML-j)/((double) n_PML) ,exponent);

		K_Fz_a[j]       = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		K_Fz_a[ny-j-1]  = K_Fz_a[j];

		K_Fz_b[j]       = 2.0*eps_0*dt/(2.0*eps_0*ka_y + sigma_y*dt);
		K_Fz_b[ny-j-1]  = K_Fz_b[j];

		K_Hy_a[j]       = (2.0*eps_0*ka_y + sigma_y*dt)/(2.0*eps_0*mu_0);      
		K_Hy_a[ny-j-1]  = K_Hy_a[j];

		K_Hy_b[j]       = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*mu_0);      
		K_Hy_b[ny-j-1]  = K_Hy_b[j];

		sigma_y         = sigma_max*pow( (n_PML - j - 0.5)/n_PML, exponent);
		ka_y            = 1.0 + (ka_max - 1.0)*pow( (n_PML - j - 0.5)/n_PML, exponent);

		K_Gx_a[j]       = (2.0*eps_0*ka_y - sigma_y*dt)/(2.0*eps_0*ka_y + sigma_y*dt);
		K_Gx_a[ny-j-2]  = K_Gx_a[j];

		K_Gx_b[j]       = 2.0*eps_0*dt/( (2.0*eps_0*ka_y + sigma_y*dt)*dy );
		K_Gx_b[ny-j-2]  = K_Gx_b[j];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::calc_Ez()
{
	int i, j;
	double Fz_r, eps_r;

	/*for (i = 1; i<nx-1; i++)
	{
		for (j = 1; j<ny-1; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Fz_r = Fz[i][j];

			Fz[i][j] = K_Fz_a[j]*Fz[i][j] + K_Fz_b[j]*( (Hy[i][j] - Hy[i-1][j])/dx -
														(Hx[i][j] - Hx[i][j-1])/dy );

			Ez[i][j] = K_Ez_a[i]*Ez[i][j] + K_Ez_b[i]*(Fz[i][j] - Fz_r)/eps_r;
		}
	}*/

	
	for (i = 1; i<n_PML; i++)
	{
		for (j = 1; j<ny-1; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Fz_r = Fz[i][j];

			Fz[i][j] = K_Fz_a[j]*Fz[i][j] + K_Fz_b[j]*( (Hy[i][j] - Hy[i-1][j])*invdx -
														(Hx[i][j] - Hx[i][j-1])*invdy );

			Ez[i][j] = K_Ez_a[i]*Ez[i][j] + K_Ez_b[i]*(Fz[i][j] - Fz_r)/eps_r;
		}
	}

	for (i = n_PML; i<nx-n_PML; i++)
	{
		for (j = 1; j<n_PML; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Fz_r = Fz[i][j];

			Fz[i][j] = K_Fz_a[j]*Fz[i][j] + K_Fz_b[j]*( (Hy[i][j] - Hy[i-1][j])*invdx -
														(Hx[i][j] - Hx[i][j-1])*invdy );

			Ez[i][j] = K_Ez_a[i]*Ez[i][j] + K_Ez_b[i]*(Fz[i][j] - Fz_r)/eps_r;
		}
	}

	//outside PML
	for (i = n_PML; i<nx-n_PML; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Ez[i][j] = Ez[i][j] + ( K_Ez_y*(Hy[i][j] - Hy[i-1][j]) -
									K_Ez_x*(Hx[i][j] - Hx[i][j-1]) )/eps_r;
		}
	}

	for (i = n_PML; i<nx-n_PML; i++)
	{
		for (j = ny-n_PML; j<ny-1; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Fz_r = Fz[i][j];

			Fz[i][j] = K_Fz_a[j]*Fz[i][j] + K_Fz_b[j]*( (Hy[i][j] - Hy[i-1][j])*invdx -
														(Hx[i][j] - Hx[i][j-1])*invdy );

			Ez[i][j] = K_Ez_a[i]*Ez[i][j] + K_Ez_b[i]*(Fz[i][j] - Fz_r)/eps_r;
		}
	}

	for (i = nx-n_PML; i<nx-1; i++)
	{
		for (j = 1; j<ny-1; j++)
		{
			eps_r = Mat[Ind[i][j]][0];

			Fz_r = Fz[i][j];

			Fz[i][j] = K_Fz_a[j]*Fz[i][j] + K_Fz_b[j]*( (Hy[i][j] - Hy[i-1][j])*invdx -
														(Hx[i][j] - Hx[i][j-1])*invdy );

			Ez[i][j] = K_Ez_a[i]*Ez[i][j] + K_Ez_b[i]*(Fz[i][j] - Fz_r)/eps_r;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hx field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::calc_Hx()
{
	int i, j;
	double Gx_r, mu_r;

	/*for (i = 0; i<nx; i++)
	{
		for (j = 0; j<ny-1; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gx_r = Gx[i][j];

			Gx[i][j] = K_Gx_a[j]*Gx[i][j] - K_Gx_b[j]*( Ez[i][j+1] - Ez[i][j] );

			Hx[i][j] = Hx[i][j] + (K_Hx_a[i]*Gx[i][j] - K_Hx_b[i]*Gx_r)/mu_r;
		}}*/

	for (i = 0; i<nx; i++)
	{
		for (j = 0; j<n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gx_r = Gx[i][j];

			Gx[i][j] = K_Gx_a[j]*Gx[i][j] - K_Gx_b[j]*( Ez[i][j+1] - Ez[i][j] );

			Hx[i][j] = Hx[i][j] + (K_Hx_a[i]*Gx[i][j] - K_Hx_b[i]*Gx_r)/mu_r;
		}
	}

	for (i = 0; i<n_PML; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gx_r = Gx[i][j];

			Gx[i][j] = K_Gx_a[j]*Gx[i][j] - K_Gx_b[j]*( Ez[i][j+1] - Ez[i][j] );

			Hx[i][j] = Hx[i][j] + (K_Hx_a[i]*Gx[i][j] - K_Hx_b[i]*Gx_r)/mu_r;
		}
	}

	//outside PML
	for (i = n_PML; i<nx-n_PML; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Hx[i][j] = Hx[i][j] - K_Hx*( Ez[i][j+1] - Ez[i][j] )/mu_r;
		}
	}

	for (i = nx-n_PML; i<nx; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gx_r = Gx[i][j];

			Gx[i][j] = K_Gx_a[j]*Gx[i][j] - K_Gx_b[j]*( Ez[i][j+1] - Ez[i][j] );

			Hx[i][j] = Hx[i][j] + (K_Hx_a[i]*Gx[i][j] - K_Hx_b[i]*Gx_r)/mu_r;
		}
	}

	for (i = 0; i<nx; i++)
	{
		for (j = ny-n_PML; j<ny-1; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gx_r = Gx[i][j];

			Gx[i][j] = K_Gx_a[j]*Gx[i][j] - K_Gx_b[j]*( Ez[i][j+1] - Ez[i][j] );

			Hx[i][j] = Hx[i][j] + (K_Hx_a[i]*Gx[i][j] - K_Hx_b[i]*Gx_r)/mu_r;
		}
	}

}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::calc_Hy()
{
	int i, j;
	double Gy_r, mu_r;

	/*for (i = 0; i<nx-1; i++)
	{
		for (j = 0; j<ny; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gy_r = Gy[i][j];

			Gy[i][j] = K_Gy_a[i]*Gy[i][j] + K_Gy_b[i]*( Ez[i+1][j] - Ez[i][j] );

			Hy[i][j] = Hy[i][j] + (K_Hy_a[j]*Gy[i][j] - K_Hy_b[j]*Gy_r)/mu_r;
		}
	}*/
	
	for (i = 0; i<nx-1; i++)
	{
		for (j = 0; j<n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gy_r = Gy[i][j];

			Gy[i][j] = K_Gy_a[i]*Gy[i][j] + K_Gy_b[i]*( Ez[i+1][j] - Ez[i][j] );

			Hy[i][j] = Hy[i][j] + (K_Hy_a[j]*Gy[i][j] - K_Hy_b[j]*Gy_r)/mu_r;
		}
	}

	for (i = 0; i<n_PML; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gy_r = Gy[i][j];

			Gy[i][j] = K_Gy_a[i]*Gy[i][j] + K_Gy_b[i]*( Ez[i+1][j] - Ez[i][j] );

			Hy[i][j] = Hy[i][j] + (K_Hy_a[j]*Gy[i][j] - K_Hy_b[j]*Gy_r)/mu_r;
		}
	}

	//outside PML
	for (i = n_PML; i<nx-n_PML; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Hy[i][j] = Hy[i][j] + K_Hy*( Ez[i+1][j] - Ez[i][j] )/mu_r;
		}
	}

	for (i = nx-n_PML; i<nx-1; i++)
	{
		for (j = n_PML; j<ny-n_PML; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gy_r = Gy[i][j];

			Gy[i][j] = K_Gy_a[i]*Gy[i][j] + K_Gy_b[i]*( Ez[i+1][j] - Ez[i][j] );

			Hy[i][j] = Hy[i][j] + (K_Hy_a[j]*Gy[i][j] - K_Hy_b[j]*Gy_r)/mu_r;
		}
	}

	for (i = 0; i<nx-1; i++)
	{
		for (j = ny-n_PML; j<ny; j++)
		{		
			mu_r = Mat[Ind[i][j]][1];

			Gy_r = Gy[i][j];

			Gy[i][j] = K_Gy_a[i]*Gy[i][j] + K_Gy_b[i]*( Ez[i+1][j] - Ez[i][j] );

			Hy[i][j] = Hy[i][j] + (K_Hy_a[j]*Gy[i][j] - K_Hy_b[j]*Gy_r)/mu_r;
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::Init_ptSource(int **Coord, int nCoord, double **Par_Excit, int s_type)
{
	Coord_ptSource = Coord;
	n_ptSource = nCoord;
	Par_ptSource = Par_Excit;
	source_type = s_type;

	jel_plane_wave = 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init sinusoidal modulated Gaussian plane wave
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::Init_plWave_Gauss(double E0, double t0, double tw, double PHI,
								   int nxa, int nxb, int nya, int nyb)
{
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	double sigma = 0.0;
	//the incidence angle of the plane wave
	phi = PHI;

	//the size of the total field zone
    l_x = nxb - nxa + 1;
	l_y = nyb - nya + 1;
	cos_phi = cos(phi);
	sin_phi = sin(phi);
	
	l_1D = (int) (l_x*abs(cos_phi) + l_y*abs(sin_phi)) + 10 + 2*n_PML;
	//l_1D = nx;
	
	d = dx*cos_phi*cos_phi + dy*sin_phi*sin_phi;

	//fdtd_1D.Init(l_1D, dt, d, eps_r, mu_r);
	fdtd_1D.Init(l_1D, n_PML, dt, d, eps_r, mu_r, sigma);
	
	fdtd_1D.Init_PML_Param();

	fdtd_1D.Init_Gauss_Source(E0, t0, tw);
	
	fdtd_1D.Get_Data(Ez_1D,Hy_1D);
	
	n_TS_xa = nxa;
	n_TS_xb = nxb;
	n_TS_ya = nya;
	n_TS_yb = nyb;
	
	jel_plane_wave = 1;

	if (Init_TS_Param())
	{
		// memory allocation problems
		return 1;
	}
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init sinusoidal plane wave
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::Init_plWave_Sin(double E0, double omega, double PHI, int nxa,
								 int nxb, int nya, int nyb)
{
	double phase = 0.0;
	
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	double sigma = 0.0;
	//the incidence angle of the plane wave
	phi = PHI;

	//the size of the total field zone
    l_x = nxb - nxa + 1;
	l_y = nyb - nya + 1;
	cos_phi = cos(phi);
	sin_phi = sin(phi);
	
	l_1D = (int) (l_x*abs(cos_phi) + l_y*abs(sin_phi)) + 10 + 2*n_PML;
	//l_1D = nx;
	
	d = dx*cos_phi*cos_phi + dy*sin_phi*sin_phi;

	//correction term to shyncronize the 1D and 2D wave propagations
	double phase_vel, phase_vel_ref, corr_term;
	phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, phi, &phase_vel);
	if (abs(cos_phi) >= abs(sin_phi) ) //the x axis is closer
	{
		phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, 0.0, &phase_vel_ref);
	}
	else //the y axis is closer
	{
		phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, pi/2.0, &phase_vel_ref);
	}
	corr_term = phase_vel_ref/phase_vel;

    //fdtd_1D.Init(l_1D, dt, d, eps_r, mu_r);
	fdtd_1D.Init(l_1D, n_PML, dt, d, corr_term*eps_r, corr_term*mu_r, sigma);
	
	fdtd_1D.Init_PML_Param();
	
	fdtd_1D.Init_Sin_Source(E0, omega, phase);
	
	fdtd_1D.Get_Data(Ez_1D,Hy_1D);
	
	n_TS_xa = nxa;
	n_TS_xb = nxb;
	n_TS_ya = nya;
	n_TS_yb = nyb;
	
	jel_plane_wave = 1;

	if (Init_TS_Param())
	{
		// memory allocation problems
		return 1;
	}

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init sinusoidal modulated Gaussian plane wave
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_2D_TM::Init_plWave_GaussSin(double E0, double omega, double t0, double tw, 
									  double PHI, int nxa, int nxb, int nya, int nyb)
{
	double phase = 0.0;
	
	double eps_r = Mat[0][0];
	double mu_r  = Mat[0][1];
	double sigma = 0.0;
	//the incidence angle of the plane wave
	phi = PHI;

	//the size of the total field zone
    l_x = nxb - nxa + 1;
	l_y = nyb - nya + 1;
	cos_phi = cos(phi);
	sin_phi = sin(phi);
	
	l_1D = (int) (l_x*abs(cos_phi) + l_y*abs(sin_phi)) + 10 + 2*n_PML;
	//l_1D = nx;
	
	d = dx*cos_phi*cos_phi + dy*sin_phi*sin_phi;

	//correction term to shyncronize the 1D and 2D wave propagations
	double phase_vel, phase_vel_ref, corr_term;
	phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, phi, &phase_vel);
	if (abs(cos_phi) >= abs(sin_phi) )
	{
		phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, 0.0, &phase_vel_ref);
	}
	else
	{
		phase_velocity(dx, dy, dt, Mat[0][0], Mat[0][1], omega, pi/2, &phase_vel_ref);
	}
	corr_term = phase_vel_ref/phase_vel;

	//fdtd_1D.Init(l_1D, dt, d, eps_r, mu_r);
	fdtd_1D.Init(l_1D, n_PML, dt, d, corr_term*eps_r, corr_term*mu_r, sigma);
	
	fdtd_1D.Init_PML_Param();

	fdtd_1D.Init_GaussSin(E0, omega, phase, t0, tw);
	
	fdtd_1D.Get_Data(Ez_1D,Hy_1D);
	
	n_TS_xa = nxa;
	n_TS_xb = nxb;
	n_TS_ya = nya;
	n_TS_yb = nyb;
	
	jel_plane_wave = 1;

	if (Init_TS_Param())
	{
		// memory allocation problems
		return 1;
	}
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Curent source J
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::ptSource_J(double **&xx, double **&Par_ptS, double time)
{
	int i, i1, j1;
	double J0, omega, phi, t0, tw;
	
	for (i = 0; i <n_ptSource; i++ )
	{
		i1 = Coord_ptSource[i][0];
		j1 = Coord_ptSource[i][1];
		J0 = Par_ptS[i][0];//*dtDIVeps0/Mat[Ind[i1][j1][k1]][0];
		switch (source_type)
		{
			case 1: //Gaussian pulse
				tw = Par_ptS[i][1];
				t0 = Par_ptS[i][2];
				xx[i1][j1] += J0*exp( -pow( (time - t0)/tw ,2) );
				break;
			case 2: //Sinusoidal plane wave
				omega = Par_ptS[i][1];
				phi   = Par_ptS[i][2];
				xx[i1][j1] += J0*cos(omega*time + phi);
				break;
			case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
				tw    = Par_ptS[i][1];
				t0    = Par_ptS[i][2];
				omega = Par_ptS[i][3];
				phi   = Par_ptS[i][4];
				xx[i1][j1] += J0*cos( omega*(time-t0) + phi)*
					             exp( -pow( (time - t0)/tw ,2) );
		}	
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ez on total-field scattered field boundary - x direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::TS_Ez()
{
	int i, j;
	double eps_r;

	//the sign is changed due to the sign of Hy_inc
	for (i = n_TS_xa; i <= n_TS_xb; i++)
	{
		eps_r = Mat[Ind[i][n_TS_ya]][0];
		Ez[i][n_TS_ya] -= K_Ez_x*Hx_f_1D[i-n_TS_xa]/eps_r*sin_phi;

		eps_r = Mat[Ind[i][n_TS_yb]][0];
		Ez[i][n_TS_yb] += K_Ez_x*Hx_b_1D[i-n_TS_xa]/eps_r*sin_phi;
	}	

	for (j = n_TS_ya; j <= n_TS_yb; j++)
	{
		eps_r = Mat[Ind[n_TS_xa][j]][0];
		Ez[n_TS_xa][j] -= K_Ez_y*Hy_l_1D[j-n_TS_ya]/eps_r*cos_phi;
		
		eps_r = Mat[Ind[n_TS_xb][j]][0];
		Ez[n_TS_xb][j] += K_Ez_y*Hy_r_1D[j-n_TS_ya]/eps_r*cos_phi;
	}	
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hx on total-field scattered field boundary - x direction
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_2D_TM::TS_Hx()
{
	int i;
	double mu_r;
	
	for (i = n_TS_xa; i <= n_TS_xb; i++)