fdtd_3d_pml_ey.cpp

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

			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fy_r = Fy_15[ia][ja][ka];
				
				Fy_15[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_15[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] + 
					          K_Ey_b_2[ii]*TwoOREp0*( Fy_15[ia][ja][ka] - Fy_r);
				ka++;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ey_16()
{
	long  i, j, k, jj, ja, ka;
	double Fy_r;

	//16
	for (i = 1; i < nPML_x_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fy_r = Fy_16[i][ja][ka];
				
				Fy_16[i][ja][ka] = K_a[Ind[i][j][k]]*Fy_16[i][ja][ka] + 
					               K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                       (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Ey[i][j][k] =  K_Ey_a_1[i]*Ey[i][j][k] + 
					           K_Ey_b_1[i]*( K_Ey_c_2[jj]*Fy_16[i][ja][ka] + K_Ey_d_2[jj]*Fy_r);
				ka++;
			}
			ja++;
			jj--;
		}
	}
}

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

	//17
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fy_r = Fy_17[ia][ja][ka];
				
				Fy_17[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_17[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Ey[i][j][k] +=  inv_TwoOREp0*( K_Ey_c_2[jj]*Fy_17[ia][ja][ka] + 
					                           K_Ey_d_2[jj]*Fy_r);
				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ey_18()
{
	long  i, j, k, ii, jj, ia, ja, ka;
	double Fy_r;

	//18
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_nPML_z_2; k++)
			{
				Fy_r = Fy_18[ia][ja][ka];
				
				Fy_18[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_18[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] + 
					          K_Ey_b_2[ii]*( K_Ey_c_2[jj]*Fy_18[ia][ja][ka] + K_Ey_d_2[jj]*Fy_r);
				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ey_19()
{
	long  i, j, k, kk, ka;
	double Fy_r, Gy_r;

	//19
	for (i = 1; i < nPML_x_1; i++)
	{	
		for (j = 0; 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++)
			{
				Fy_r = Fy_19[i][j][ka];
				
				Fy_19[i][j][ka] = K_a[Ind[i][j][k]]*Fy_19[i][j][ka] + 
					              K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                      (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_19[i][j][ka];

				Gy_19[i][j][ka] = K_Gy_a_2[kk]*Gy_19[i][j][ka] + 
					              K_Gy_b_2[kk]*(Fy_19[i][j][ka] - Fy_r);

				Ey[i][j][k] = K_Ey_a_1[i]*Ey[i][j][k] + 
					          K_Ey_b_1[i]*( K_Ey_c_1[j]*Gy_19[i][j][ka] + K_Ey_d_1[j]*Gy_r);
				ka++;
				kk--;
			}
		}
	}
}

void CFDTD_3D_PML::calc_Ey_20()
{
	long  i, j, k, kk, ia, ka;
	double Fy_r, Gy_r;

	//20
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		for (j = 0; 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++)
			{
				Fy_r = Fy_20[ia][j][ka];
				
				Fy_20[ia][j][ka] = K_a[Ind[i][j][k]]*Fy_20[ia][j][ka] + 
					               K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                       (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_20[ia][j][ka];

				Gy_20[ia][j][ka] = K_Gy_a_2[kk]*Gy_20[ia][j][ka] + 
					               K_Gy_b_2[kk]*(Fy_20[ia][j][ka] - Fy_r);

				Ey[i][j][k] +=  inv_TwoOREp0*( K_Ey_c_1[j]*Gy_20[ia][j][ka] + 
					                           K_Ey_d_1[j]*Gy_r);
				ka++;
				kk--;
			}
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ey_21()
{
	long  i, j, k, ii, kk, ia, ka;
	double Fy_r, Gy_r;

	//21
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		for (j = 0; 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++)
			{
				Fy_r = Fy_21[ia][j][ka];
				
				Fy_21[ia][j][ka] = K_a[Ind[i][j][k]]*Fy_21[ia][j][ka] + 
					               K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                       (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_21[ia][j][ka];

				Gy_21[ia][j][ka] = K_Gy_a_2[kk]*Gy_21[ia][j][ka] + 
					               K_Gy_b_2[kk]*(Fy_21[ia][j][ka] - Fy_r);

				Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] + 
					          K_Ey_b_2[ii]*( K_Ey_c_1[j]*Gy_21[ia][j][ka] + K_Ey_d_1[j]*Gy_r);
				ka++;
				kk--;
			}
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ey_22()
{
	long  i, j, k, kk, ja, ka;
	double Fy_r, Gy_r;

	//22
	for (i = 1; i < nPML_x_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_1_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++)
			{
				Fy_r = Fy_22[i][ja][ka];
				
				Fy_22[i][ja][ka] = K_a[Ind[i][j][k]]*Fy_22[i][ja][ka] + 
					               K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                       (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_22[i][ja][ka];

				Gy_22[i][ja][ka] = K_Gy_a_2[kk]*Gy_22[i][ja][ka] + 
					               K_Gy_b_2[kk]*(Fy_22[i][ja][ka] - Fy_r);

				Ey[i][j][k] = K_Ey_a_1[i]*Ey[i][j][k] + 
					          K_Ey_b_1[i]*TwoOREp0*( Gy_22[i][ja][ka] - Gy_r);
				ka++;
				kk--;
			}
			ja++;
		}
	}
}

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

	//23
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_1_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++)
			{
				Fy_r = Fy_23[ia][ja][ka];
				
				Fy_23[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_23[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
						 		                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Ey[i][j][k] = K_Gy_a_2[kk]*Ey[i][j][k] + 
					          K_Gy_b_2[kk]*(Fy_23[ia][ja][ka] - Fy_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ey_24()
{
	long  i, j, k, ii, kk, ia, ja, ka;
	double Fy_r, Gy_r;
	
	//24
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		for (j = nPML_y_1; j < ny_MIN_1_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++)
			{
				Fy_r = Fy_24[ia][ja][ka];
				
				Fy_24[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_24[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_24[ia][ja][ka];

				Gy_24[ia][ja][ka] = K_Gy_a_2[kk]*Gy_24[ia][ja][ka] + 
					                K_Gy_b_2[kk]*(Fy_24[ia][ja][ka] - Fy_r);

				Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] + 
					          K_Ey_b_2[ii]*TwoOREp0*( Gy_24[ia][ja][ka] - Gy_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void CFDTD_3D_PML::calc_Ey_25()
{
	long  i, j, k, jj, kk, ja, ka;
	double Fy_r, Gy_r;

	//25
	for (i = 1; i < nPML_x_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_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++)
			{
				Fy_r = Fy_25[i][ja][ka];
				
				Fy_25[i][ja][ka] = K_a[Ind[i][j][k]]*Fy_25[i][ja][ka] + 
					               K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                       (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_25[i][ja][ka];

				Gy_25[i][ja][ka] = K_Gy_a_2[kk]*Gy_25[i][ja][ka] + 
					               K_Gy_b_2[kk]*(Fy_25[i][ja][ka] - Fy_r);

				Ey[i][j][k] =  K_Ey_a_1[i]*Ey[i][j][k] + 
					           K_Ey_b_1[i]*( K_Ey_c_2[jj]*Gy_25[i][ja][ka] + K_Ey_d_2[jj]*Gy_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
	}
}

void CFDTD_3D_PML::calc_Ey_26()
{
	long  i, j, k, jj, kk, ia, ja, ka;
	double Fy_r, Gy_r;

	//26
	ia = 0;
	for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_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++)
			{
				Fy_r = Fy_26[ia][ja][ka];
				
				Fy_26[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_26[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_26[ia][ja][ka];

				Gy_26[ia][ja][ka] = K_Gy_a_2[kk]*Gy_26[ia][ja][ka] + 
					                K_Gy_b_2[kk]*(Fy_26[ia][ja][ka] - Fy_r);

				Ey[i][j][k] +=  inv_TwoOREp0*( K_Ey_c_2[jj]*Gy_26[ia][ja][ka] + 
					                           K_Ey_d_2[jj]*Gy_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void CFDTD_3D_PML::calc_Ey_27()
{
	long  i, j, k, ii, jj, kk, ia, ja, ka;
	double Fy_r, Gy_r;

	//27
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
	{	
		ja = 0;
		jj = nPML_y_2_MIN_1;
		for (j = ny_MIN_1_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++)
			{
				Fy_r = Fy_27[ia][ja][ka];
				
				Fy_27[ia][ja][ka] = K_a[Ind[i][j][k]]*Fy_27[ia][ja][ka] + 
					                K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz - 
								                        (Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );

				Gy_r = Gy_27[ia][ja][ka];

				Gy_27[ia][ja][ka] = K_Gy_a_2[kk]*Gy_27[ia][ja][ka] + 
					                K_Gy_b_2[kk]*(Fy_27[ia][ja][ka] - Fy_r);

				Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] + 
					          K_Ey_b_2[ii]*( K_Ey_c_2[jj]*Gy_27[ia][ja][ka] + K_Ey_d_2[jj]*Gy_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}