fdtd_3d_pml_ex.cpp

利用c++语言写的三维FDTD

		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fx_r = Fx_15[ia][ja][ka];
				
				Fx_15[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_15[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] += inv_TwoOREp0*(K_Ex_c_2[ii]*Fx_15[ia][ja][ka] +
					                         K_Ex_d_2[ii]*Fx_r);
				ka++;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ex_16()
{
	double Fx_r, Gx_r;
	long i, j, k, jj, ja, ka;

	//16
	for (i = 0; i < nPML_x_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fx_r = Fx_16[i][ja][ka];
				
				Fx_16[i][ja][ka] = K_a[Ind[i][j][k]]*Fx_16[i][ja][ka] +
					               K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                       (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_16[i][ja][ka];

				Gx_16[i][ja][ka] = K_Gx_a_2[jj]*Gx_16[i][ja][ka] + 
					               K_Gx_b_2[jj]*(Fx_16[i][ja][ka] - Fx_r);

				Ex[i][j][k] += inv_TwoOREp0*( K_Ex_c_1[i]*Gx_16[i][ja][ka] +
					                          K_Ex_d_1[i]*Gx_r);
				ka++;
			}
			ja++;
			jj--;
		}
	}
}

void CFDTD_3D_PML::calc_Ex_17()
{
	double Fx_r;
	long i, j, k, jj, ia, ja, ka;

	//17
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fx_r = Fx_17[ia][ja][ka];
				
				Fx_17[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_17[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Gx_a_2[jj]*Ex[i][j][k] + 
					          K_Gx_b_2[jj]*(Fx_17[ia][ja][ka] - Fx_r);
				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ex_18()
{
	double Fx_r, Gx_r;
	long i, j, k, ii, jj, ia, ja, ka;

	//18
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fx_r = Fx_18[ia][ja][ka];
				
				Fx_18[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_18[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_18[ia][ja][ka];

				Gx_18[ia][ja][ka] = K_Gx_a_2[jj]*Gx_18[ia][ja][ka] + 
					                K_Gx_b_2[jj]*(Fx_18[ia][ja][ka] - Fx_r);

				Ex[i][j][k] += inv_TwoOREp0*( K_Ex_c_2[ii]*Gx_18[ia][ja][ka] +
					                          K_Ex_d_2[ii]*Gx_r);

				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ex_19()
{
	double Fx_r, Gx_r;
	long i, j, k, kk, ka;

	//19
	for (i = 0; i < nPML_x_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_19[i][j][ka];
				
				Fx_19[i][j][ka] = K_a[Ind[i][j][k]]*Fx_19[i][j][ka] +
					              K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                      (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_19[i][j][ka];

				Gx_19[i][j][ka] = K_Gx_a_1[j]*Gx_19[i][j][ka] + 
					              K_Gx_b_1[j]*(Fx_19[i][j][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*( K_Ex_c_1[i]*Gx_19[i][j][ka] + K_Ex_d_1[i]*Gx_r);
				ka++;
				kk--;
			}
		}
	}
}

void CFDTD_3D_PML::calc_Ex_20()
{
	double Fx_r, Gx_r;
	long i, j, k, kk, ia, ka;

	//20
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_20[ia][j][ka];
				
				Fx_20[ia][j][ka] = K_a[Ind[i][j][k]]*Fx_20[ia][j][ka] +
					               K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                       (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_20[ia][j][ka];

				Gx_20[ia][j][ka] = K_Gx_a_1[j]*Gx_20[ia][j][ka] + 
					               K_Gx_b_1[j]*(Fx_20[ia][j][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*TwoOREp0*( Gx_20[ia][j][ka] - Gx_r);
				ka++;
				kk--;
			}
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ex_21()
{
	double Fx_r, Gx_r;
	long i, j, k, ii, kk, ia, ka;

	//21
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		for (j = 1; j < nPML_y_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_21[ia][j][ka];
				
				Fx_21[ia][j][ka] = K_a[Ind[i][j][k]]*Fx_21[ia][j][ka] +
					               K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                       (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_21[ia][j][ka];

				Gx_21[ia][j][ka] = K_Gx_a_1[j]*Gx_21[ia][j][ka] + 
					               K_Gx_b_1[j]*(Fx_21[ia][j][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*( K_Ex_c_2[ii]*Gx_21[ia][j][ka] + K_Ex_d_2[ii]*Gx_r);
				ka++;
				kk--;
			}
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ex_22()
{
	double Fx_r;
	long i, j, k, kk, ja, ka;

	//22
	for (i = 0; i < nPML_x_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_22[i][ja][ka];
				
				Fx_22[i][ja][ka] = K_a[Ind[i][j][k]]*Fx_22[i][ja][ka] +
					               K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                       (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*( K_Ex_c_1[i]*Fx_22[i][ja][ka] + K_Ex_d_1[i]*Fx_r);
				ka++;
				kk--;
			}
			ja++;
		}
	}
}

void CFDTD_3D_PML::calc_Ex_23()
{
	double Fx_r;
	long i, j, k, kk, ia, ja, ka;

	//23
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_23[ia][ja][ka];
				
				Fx_23[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_23[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*TwoOREp0*( Fx_23[ia][ja][ka] - Fx_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ex_24()
{
	double Fx_r;
	long i, j, k, ii, kk, ia, ja, ka;

	//24
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_24[ia][ja][ka];
				
				Fx_24[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_24[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*( K_Ex_c_2[ii]*Fx_24[ia][ja][ka] + K_Ex_d_2[ii]*Fx_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ex_25()
{
	double Fx_r, Gx_r;
	long i, j, k, jj, kk, ja, ka;

	//25
	for (i = 0; i < nPML_x_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_25[i][ja][ka];
				
				Fx_25[i][ja][ka] = K_a[Ind[i][j][k]]*Fx_25[i][ja][ka] +
					               K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                       (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_25[i][ja][ka];

				Gx_25[i][ja][ka] = K_Gx_a_2[jj]*Gx_25[i][ja][ka] + 
					               K_Gx_b_2[jj]*(Fx_25[i][ja][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*( K_Ex_c_1[i]*Gx_25[i][ja][ka] + K_Ex_d_1[i]*Gx_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
	}
}

void CFDTD_3D_PML::calc_Ex_26()
{
	double Fx_r, Gx_r;
	long i, j, k, jj, kk, ia, ja, ka;

	//26
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_26[ia][ja][ka];
				
				Fx_26[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_26[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_26[ia][ja][ka];

				Gx_26[ia][ja][ka] = K_Gx_a_2[jj]*Gx_26[ia][ja][ka] + 
					                K_Gx_b_2[jj]*(Fx_26[ia][ja][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*TwoOREp0*( Gx_26[ia][ja][ka] - Gx_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ex_27()
{
	double Fx_r, Gx_r;
	long i, j, k, ii, jj, kk, ia, ja, ka;
   
	//27
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fx_r = Fx_27[ia][ja][ka];
				
				Fx_27[ia][ja][ka] = K_a[Ind[i][j][k]]*Fx_27[ia][ja][ka] +
					                K_b[Ind[i][j][k]]*( (Hz[i][j][k] - Hz[i][j-1][k])*inv_dy - 
								                        (Hy[i][j][k] - Hy[i][j][k-1])*inv_dz );

				Gx_r = Gx_27[ia][ja][ka];

				Gx_27[ia][ja][ka] = K_Gx_a_2[jj]*Gx_27[ia][ja][ka] + 
					                K_Gx_b_2[jj]*(Fx_27[ia][ja][ka] - Fx_r);

				Ex[i][j][k] = K_Ex_a_2[kk]*Ex[i][j][k] + 
					          K_Ex_b_2[kk]*(K_Ex_c_2[ii]*Gx_27[ia][ja][ka] + K_Ex_d_2[ii]*Gx_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}