fdtd_3d_pml_ez.cpp

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

			{
				Fz_r = Fz_15[ia][ja][ka];

				Fz_15[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_15[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Ez[i][j][k] =  K_Gz_a_2[ii]*Ez[i][j][k] + 
					           K_Gz_b_2[ii]*(Fz_15[ia][ja][ka] - Fz_r);
				ka++;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void calc_Ez_16(double ***Ez, double ***Hx, double ***Hy, double ***Fz_16,
				double ***Gz_16, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_1,
				double *K_Gz_b_1, double *K_Ez_a_2, double *K_Ez_b_2, double inv_dx, 
				double inv_dy, long  nPML_x_1, long  ny_MIN_nPML_y_2, long  ny_MIN_1, 
				long  nPML_z_1, long  nz_MIN_1_MIN_nPML_z_2, long  nPML_y_2_MIN_1, 
				double TwoOREp0)
{
	long  i, j, k, jj, ja, ka;
	double Fz_r, Gz_r;
	
	//16
	for (i = 1; 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_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_16[i][ja][ka];

				Fz_16[i][ja][ka] = K_a[Ind[i][j][k]]*Fz_16[i][ja][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_16[i][ja][ka];

				Gz_16[i][ja][ka] = K_Gz_a_1[i]*Gz_16[i][ja][ka] + 
					               K_Gz_b_1[i]*(Fz_16[i][ja][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*TwoOREp0*( Gz_16[i][ja][ka] - Gz_r);
				ka++;
			}
			ja++;
			jj--;
		}
	}
}

void calc_Ez_17(double ***Ez, double ***Hx, double ***Hy, double ***Fz_17, long  ***Ind,
			    double *K_a, double *K_b, double *K_Ez_a_2, double *K_Ez_b_2, 
			    double inv_dx, double inv_dy, long  nPML_x_1, long  nx_MIN_nPML_x_2, 
				long  ny_MIN_nPML_y_2, long  ny_MIN_1, long  nPML_z_1, 
				long  nz_MIN_1_MIN_nPML_z_2, long  nPML_y_2_MIN_1, double TwoOREp0)
{
	long  i, j, k, jj, ia, ja, ka;
	double Fz_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_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_17[ia][ja][ka];

				Fz_17[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_17[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*TwoOREp0*( Fz_17[ia][ja][ka] - Fz_r);
				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void calc_Ez_18(double ***Ez, double ***Hx, double ***Hy, double ***Fz_18, 
				double ***Gz_18, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_2,
				double *K_Gz_b_2, double *K_Ez_a_2, double *K_Ez_b_2, double inv_dx, 
				double inv_dy, long  nx_MIN_nPML_x_2, long  nx_MIN_1, long  ny_MIN_nPML_y_2, 
				long  ny_MIN_1, long  nPML_z_1, long  nz_MIN_1_MIN_nPML_z_2, 
				long  nPML_x_2_MIN_1, long  nPML_y_2_MIN_1, double TwoOREp0)
{
	long  i, j, k, ii, jj, ia, ja, ka;
	double Fz_r, Gz_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_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
			{
				Fz_r = Fz_18[ia][ja][ka];

				Fz_18[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_18[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_18[ia][ja][ka];

				Gz_18[ia][ja][ka] = K_Gz_a_2[ii]*Gz_18[ia][ja][ka] + 
					                K_Gz_b_2[ii]*(Fz_18[ia][ja][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*TwoOREp0*( Gz_18[ia][ja][ka] - Gz_r);
				ka++;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}

void calc_Ez_19(double ***Ez, double ***Hx, double ***Hy, double ***Fz_19, 
				double ***Gz_19, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_1,
				double *K_Gz_b_1, double *K_Ez_a_1, double *K_Ez_b_1, double *K_Ez_c_2, 
				double *K_Ez_d_2, double inv_dx, double inv_dy, long  nPML_x_1,
				long  nPML_y_1, long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, 
				long  nPML_z_2_MIN_1)
{
	long  i, j, k, kk, ka;
	double Fz_r, Gz_r;

	//19
	for (i = 1; 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_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_19[i][j][ka];

				Fz_19[i][j][ka] = K_a[Ind[i][j][k]]*Fz_19[i][j][ka] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_19[i][j][ka];

				Gz_19[i][j][ka] = K_Gz_a_1[i]*Gz_19[i][j][ka] + 
					              K_Gz_b_1[i]*(Fz_19[i][j][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_2[kk]*Gz_19[i][j][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
		}
	}
}

void calc_Ez_20(double ***Ez, double ***Hx, double ***Hy, double ***Fz_20, long  ***Ind,
			   double *K_a, double *K_b, double *K_Ez_a_1, double *K_Ez_b_1, 
			   double *K_Ez_c_2, double *K_Ez_d_2, double inv_dx, double inv_dy, 
			   long  nPML_x_1, long  nx_MIN_nPML_x_2, long  nPML_y_1, 
			   long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, long  nPML_z_2_MIN_1)
{
	long  i, j, k, kk, ia, ka;
	double Fz_r;

	//20
	ia = 0;
	for (i = nPML_x_1; i < nx_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_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_20[ia][j][ka];

				Fz_20[ia][j][ka] = K_a[Ind[i][j][k]]*Fz_20[ia][j][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_2[kk]*Fz_20[ia][j][ka] + K_Ez_d_2[kk]*Fz_r);
				ka++;
				kk--;
			}
		}
		ia++;
	}
}

void calc_Ez_21(double ***Ez, double ***Hx, double ***Hy, double ***Fz_21,
				double ***Gz_21, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_2,
				double *K_Gz_b_2, double *K_Ez_a_1, double *K_Ez_b_1, double *K_Ez_c_2,
				double *K_Ez_d_2, double inv_dx, double inv_dy, long  nx_MIN_nPML_x_2, 
				long  nx_MIN_1, long  nPML_y_1, long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, 
				long  nPML_x_2_MIN_1, long  nPML_z_2_MIN_1)
{
	long  i, j, k, ii, kk, ia, ka;
	double Fz_r, Gz_r;

	//21
	ia = 0;
	ii = nPML_x_2_MIN_1;
	for (i = nx_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_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_21[ia][j][ka];

				Fz_21[ia][j][ka] = K_a[Ind[i][j][k]]*Fz_21[ia][j][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_21[ia][j][ka];

				Gz_21[ia][j][ka] = K_Gz_a_2[ii]*Gz_21[ia][j][ka] + 
					               K_Gz_b_2[ii]*(Fz_21[ia][j][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_1[j]*Ez[i][j][k] + 
					          K_Ez_b_1[j]*( K_Ez_c_2[kk]*Gz_21[ia][j][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
		}
		ia++;
		ii--;
	}
}

void calc_Ez_22(double ***Ez, double ***Hx, double ***Hy, double ***Fz_22, 
				double ***Gz_22, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_1,
				double *K_Gz_b_1, double *K_Ez_c_2, double *K_Ez_d_2, double inv_dx, 
				double inv_dy, long  nPML_x_1, long  nPML_y_1, long  ny_MIN_nPML_y_2, 
				long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, long  nPML_z_2_MIN_1,
				double inv_TwoOREp0)
{
	long  i, j, k, kk, ja, ka;
	double Fz_r, Gz_r;

	//22
	for (i = 1; 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_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_22[i][ja][ka];

				Fz_22[i][ja][ka] = K_a[Ind[i][j][k]]*Fz_22[i][ja][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_22[i][ja][ka];

				Gz_22[i][ja][ka] = K_Gz_a_1[i]*Gz_22[i][ja][ka] + 
					               K_Gz_b_1[i]*(Fz_22[i][ja][ka] - Fz_r);

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_2[kk]*Gz_22[i][ja][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
			ja++;
		}
	}
}

void calc_Ez_23(double ***Ez, double ***Hx, double ***Hy, double ***Fz_23, long  ***Ind,
			   double *K_a, double *K_b, double *K_Ez_c_2, double *K_Ez_d_2,
			   double inv_dx, double inv_dy, long  nPML_x_1, long  nx_MIN_nPML_x_2, 
			   long  nPML_y_1, long  ny_MIN_nPML_y_2, long  nz_MIN_1_MIN_nPML_z_2, 
			   long  nz_MIN_1, long  nPML_z_2_MIN_1, double inv_TwoOREp0)
{
	long  i, j, k, kk, ia, ja, ka;
	double Fz_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_nPML_y_2; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_23[ia][ja][ka];

				Fz_23[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_23[ia][ja][ka] +
				                  K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                      (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_2[kk]*Fz_23[ia][ja][ka] + K_Ez_d_2[kk]*Fz_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
	}
}

void calc_Ez_24(double ***Ez, double ***Hx, double ***Hy, double ***Fz_24, 
				double ***Gz_24, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_2,
				double *K_Gz_b_2, double *K_Ez_c_2, double *K_Ez_d_2, double inv_dx,
				double inv_dy,  long  nx_MIN_nPML_x_2, long  nx_MIN_1, long  nPML_y_1, 
				long  ny_MIN_nPML_y_2, long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1,
				long  nPML_x_2_MIN_1, long  nPML_z_2_MIN_1, double inv_TwoOREp0)
{
	long  i, j, k, ii, kk, ia, ja, ka;
	double Fz_r, Gz_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_nPML_y_2; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_24[ia][ja][ka];

				Fz_24[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_24[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_24[ia][ja][ka];

				Gz_24[ia][ja][ka] = K_Gz_a_2[ii]*Gz_24[ia][ja][ka] + 
					                K_Gz_b_2[ii]*(Fz_24[ia][ja][ka] - Fz_r);

				Ez[i][j][k] += inv_TwoOREp0*( K_Ez_c_2[kk]*Gz_24[ia][ja][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
			ja++;
		}
		ia++;
		ii--;
	}
}

void calc_Ez_25(double ***Ez, double ***Hx, double ***Hy, double ***Fz_25, 
				double ***Gz_25, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_1,
				double *K_Gz_b_1, double *K_Ez_a_2, double *K_Ez_b_2, double *K_Ez_c_2,
				double *K_Ez_d_2, double inv_dx, double inv_dy, long  nPML_x_1, 
				long  ny_MIN_nPML_y_2, long  ny_MIN_1, long  nz_MIN_1_MIN_nPML_z_2, 
				long  nz_MIN_1, long  nPML_y_2_MIN_1, long  nPML_z_2_MIN_1)
{
	long  i, j, k, jj, kk, ja, ka;
	double Fz_r, Gz_r;

	//25
	for (i = 1; 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_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_25[i][ja][ka];

				Fz_25[i][ja][ka] = K_a[Ind[i][j][k]]*Fz_25[i][ja][ka] +
				                   K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                       (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_25[i][ja][ka];

				Gz_25[i][ja][ka] = K_Gz_a_1[i]*Gz_25[i][ja][ka] + 
					               K_Gz_b_1[i]*(Fz_25[i][ja][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_2[kk]*Gz_25[i][ja][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
	}
}

void calc_Ez_26(double ***Ez, double ***Hx, double ***Hy, double ***Fz_26, long  ***Ind,
			    double *K_a, double *K_b, double *K_Ez_a_2, double *K_Ez_b_2, 
			    double *K_Ez_c_2, double *K_Ez_d_2, double inv_dx, double inv_dy, 
			    long  nPML_x_1, long  nx_MIN_nPML_x_2, long  ny_MIN_nPML_y_2, long  ny_MIN_1, 
			    long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, long  nPML_y_2_MIN_1, 
				long  nPML_z_2_MIN_1)
{
	long  i, j, k, jj, kk, ia, ja, ka;
	double Fz_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_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_26[ia][ja][ka];

				Fz_26[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_26[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_2[kk]*Fz_26[ia][ja][ka] + K_Ez_d_2[kk]*Fz_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
	}
}

void calc_Ez_27(double ***Ez, double ***Hx, double ***Hy, double ***Fz_27,
				double ***Gz_27, long  ***Ind, double *K_a, double *K_b, double *K_Gz_a_2,
				double *K_Gz_b_2, double *K_Ez_a_2, double *K_Ez_b_2, double *K_Ez_c_2,
				double *K_Ez_d_2, double inv_dx, double inv_dy, long  nx_MIN_nPML_x_2, 
				long  nx_MIN_1, long  ny_MIN_nPML_y_2, long  ny_MIN_1, 
				long  nz_MIN_1_MIN_nPML_z_2, long  nz_MIN_1, long  nPML_x_2_MIN_1,
			    long  nPML_y_2_MIN_1, long  nPML_z_2_MIN_1)
{
	long  i, j, k, ii, jj, kk, ia, ja, ka;
	double Fz_r, Gz_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_nPML_y_2; j < ny_MIN_1; j++)
		{	
			ka = 0;
			kk = nPML_z_2_MIN_1;
			for (k = nz_MIN_1_MIN_nPML_z_2; k < nz_MIN_1; k++)
			{
				Fz_r = Fz_27[ia][ja][ka];

				Fz_27[ia][ja][ka] = K_a[Ind[i][j][k]]*Fz_27[ia][ja][ka] +
				                    K_b[Ind[i][j][k]]*( (Hy[i][j][k] - Hy[i-1][j][k])*inv_dx -
								                        (Hx[i][j][k] - Hx[i][j-1][k])*inv_dy );

			    Gz_r = Gz_27[ia][ja][ka];

				Gz_27[ia][ja][ka] = K_Gz_a_2[ii]*Gz_27[ia][ja][ka] + 
					                K_Gz_b_2[ii]*(Fz_27[ia][ja][ka] - Fz_r);

				Ez[i][j][k] = K_Ez_a_2[jj]*Ez[i][j][k] + 
					          K_Ez_b_2[jj]*( K_Ez_c_2[kk]*Gz_27[ia][ja][ka] + K_Ez_d_2[kk]*Gz_r);
				ka++;
				kk--;
			}
			ja++;
			jj--;
		}
		ia++;
		ii--;
	}
}