fdtd_3d.cpp

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

//Add the incident field to Hy on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Hy(double ***Hy, long  ***Ind, double **Mat, double **Ez_i0, double **Ez_i1,
		   double **Ex_k0, double **Ex_k1, double ct_Hy_1, double ct_Hy_2, 
		   long n_TS_xa, long n_TS_xb, long n_TS_ya, long n_TS_yb, long n_TS_za, 
		   long n_TS_zb, long n_TS_xa_MIN_1, long n_TS_za_MIN_1)
{
	long i, j, k;
	long ind_i, ind_j, ind_k;
	double mu_r;

	//i faces
	for (j = n_TS_ya; j <= n_TS_yb; j++)
	{
		ind_j = j - n_TS_ya;
		for (k = n_TS_za; k < n_TS_zb; k++)
		{
			ind_k = k - n_TS_za;
			mu_r = Mat[Ind[n_TS_xa_MIN_1][j][k]][1];
			Hy[n_TS_xa_MIN_1][j][k] -= Ez_i0[ind_j][ind_k]*ct_Hy_1/mu_r;
			mu_r = Mat[Ind[n_TS_xb][j][k]][1];
			Hy[n_TS_xb][j][k] += Ez_i1[ind_j][ind_k]*ct_Hy_1/mu_r;
		}
	}

	//k faces
	for (i = n_TS_xa; i < n_TS_xb; i++)
	{
		ind_i = i - n_TS_xa;
		for (j = n_TS_ya; j <= n_TS_yb; j++)
		{
			ind_j = j - n_TS_ya;
			mu_r = Mat[Ind[i][j][n_TS_za_MIN_1]][1];
			Hy[i][j][n_TS_za_MIN_1] += Ex_k0[ind_i][ind_j]*ct_Hy_2/mu_r;
			mu_r = Mat[Ind[i][j][n_TS_zb]][1];
			Hy[i][j][n_TS_zb] -= Ex_k1[ind_i][ind_j]*ct_Hy_2/mu_r;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hz on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Hz(double ***Hz, long  ***Ind, double **Mat, double **Ey_i0, double **Ey_i1,
		   double **Ex_j0, double **Ex_j1, double ct_Hz_1, double ct_Hz_2, 
		   long n_TS_xa, long n_TS_xb, long n_TS_ya, long n_TS_yb, long n_TS_za, 
		   long n_TS_zb, long n_TS_xa_MIN_1, long n_TS_ya_MIN_1)
{
	long i, j, k;
	long ind_i, ind_j, ind_k;
	double mu_r;

	//i faces
	for (j = n_TS_ya; j < n_TS_yb; j++)
	{
		ind_j = j - n_TS_ya;
		for (k = n_TS_za; k <= n_TS_zb; k++)
		{
			ind_k = k-n_TS_za;
			mu_r = Mat[Ind[n_TS_xa_MIN_1][j][k]][1];
			Hz[n_TS_xa_MIN_1][j][k] += Ey_i0[ind_j][ind_k]*ct_Hz_1/mu_r;
			mu_r = Mat[Ind[n_TS_xb][j][k]][1];
			Hz[n_TS_xb][j][k] -= Ey_i1[ind_j][ind_k]*ct_Hz_1/mu_r;
		}
	}

	//j faces
	for (i = n_TS_xa; i < n_TS_xb; i++)
	{
		ind_i = i - n_TS_xa;
		for (k = n_TS_za; k <= n_TS_zb; k++)
		{
			ind_k = k - n_TS_za;
			
			mu_r = Mat[Ind[i][n_TS_ya_MIN_1][k]][1];
			Hz[i][n_TS_ya_MIN_1][k] -= Ex_j0[ind_i][ind_k]*ct_Hz_2/mu_r;

			mu_r = Mat[Ind[i][n_TS_yb][k]][1];
			Hz[i][n_TS_yb][k] += Ex_j1[ind_i][ind_k]*ct_Hz_2/mu_r;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void Set_Data_Followed(double ***Ex, double ***Ey, double ***Ez,
					   double ***Hx, double ***Hy, double ***Hz,
					   double **Ex_foll, double **Ey_foll, double **Ez_foll,
					   double **Hx_foll, double **Hy_foll, double **Hz_foll,
					   long  **Ind_Foll, long  length_Ind_Foll, long  n_t)
{
	long  i, i1, j1, k1;
	for (i = 0; i<length_Ind_Foll; i++)
	{
		i1 = Ind_Foll[i][0];
		j1 = Ind_Foll[i][1];
		k1 = Ind_Foll[i][2];

		Ex_foll[i][n_t] = Ex[i1][j1][k1];
		
		Ey_foll[i][n_t] = Ey[i1][j1][k1];
		
		Ez_foll[i][n_t] = Ez[i1][j1][k1];
		
		Hx_foll[i][n_t] = Hx[i1][j1][k1];
		
		Hy_foll[i][n_t] = Hy[i1][j1][k1];
		
		Hz_foll[i][n_t] = Hz[i1][j1][k1];
	}
}


///////////////////////////////////////////////////////////////////////////////////////
//Calculate the energy
///////////////////////////////////////////////////////////////////////////////////////
void calc_W(double ***W, double ***Ex, double ***Ey, double ***Ez, double ***Hx, 
			double ***Hy, double ***Hz, long ***Ind, double **Mat, long nx_W_a, 
			long ny_W_a, long nz_W_a, long nx_a, long nx_b, long ny_a, long ny_b,
			long nz_a, long nz_b, double eps_0_DIV_2, double mu_0_DIV_2, double *W_av)
{
	long i, j, k, ia, ja, ka;
	
	double WaW = 0.0;

	ia = nx_a - nx_W_a ;
	for (i = nx_a; i < nx_b; i++)
	{
		ja = ny_a - ny_W_a;
		for (j = ny_a; j < ny_b; j++)
		{
			ka = nz_a - nz_W_a;
			for (k = nz_a; k < nz_b; k++)
			{
				W[ia][ja][ka] = eps_0_DIV_2*Mat[Ind[i][j][k]][0]*( Ex[i][j][k]*Ex[i][j][k] +
					              Ey[i][j][k]*Ey[i][j][k] + Ez[i][j][k]*Ez[i][j][k] ) +
				                mu_0_DIV_2*Mat[Ind[i][j][k]][1]*( Hx[i][j][k]*Hx[i][j][k] + 
								  Hy[i][j][k]*Hy[i][j][k] + Hz[i][j][k]*Hz[i][j][k] );
				
				WaW += W[ia][ja][ka];
				
				ka++;
			}
			ja++;
		}
		ia++;
	}
	*W_av = WaW;
}

///////////////////////////////////////////////////////////////////////////////////////
//Average the energy along x and z directions and time
///////////////////////////////////////////////////////////////////////////////////////
void aver_W_xz_t(double ***W, long nx_W, long ny_W, long nz_W, double *W_aver_xz, 
				 double *W_aver_xzt, double Surf_xz)
{
	long i, j, k;
	
	for (j = 0; j < ny_W; j++)
	{
		W_aver_xz[j] = 0.0;
 	    for (i = 0; i < nx_W; i++)
		{
	  		for (k = 0; k < nz_W; k++)
			{
				W_aver_xz[j] += W[i][j][k];
	  		}
		}
		W_aver_xz[j] = W_aver_xz[j]/Surf_xz;
		W_aver_xzt[j] += W_aver_xz[j];
	}
	
}

///////////////////////////////////////////////////////////////////////////////////////
//Average the energy along x and z directions and time
///////////////////////////////////////////////////////////////////////////////////////
int aver_xyz(double ***X, long nx_a_av, long nx_b_av, long ny_a_av, long ny_b_av,
			 long nz_a_av, long nz_b_av, double *X_av, double Vol)
{
	long i, j, k;
	double x_aver = 0.0;

 	for (i = nx_a_av; i <= nx_b_av; i++)
	{
		for (j = ny_a_av; j <= ny_b_av; j++)
		{
	  		for (k = nz_a_av; k <= nz_b_av; k++)
			{
				x_aver += X[i][j][k];
	  		}
		}
	}
	*X_av = x_aver/Vol;

	return 0;	
}

///////////////////////////////////////////////////////////////////////////////////////
//Average in x, y and z 
///////////////////////////////////////////////////////////////////////////////////////
int aver_xyz_WignerC(double ***X, long **Mask_Wigner, long n_Mask_Wigner, double *X_av)
{
	long i, j, k;
	long i1, j1, k1;
	double x_aver = 0.0;

	for (i = 0; i < n_Mask_Wigner; i++)
	{
		x_aver += X[Mask_Wigner[i][0]][Mask_Wigner[i][1]][Mask_Wigner[i][2]];
	}
	*X_av = x_aver/n_Mask_Wigner;

	return 0;	
}

///////////////////////////////////////////////////////////////////////////////////////
//The Bragg reflection over the Wigner cell
///////////////////////////////////////////////////////////////////////////////////////
int aver_Bragg_WignerC(double ***X, long **Mask_Wigner, long n_Mask_Wigner,
					   double kx, double ky, double kz, double Gx, double Gy, 
					   double Gz, double *X_avBr_re, double *X_avBr_im, 
					   double dx, double dy, double dz)
{
	long i;
	double x_aver_re = 0.0, x_aver_im = 0.0;

	double ct_x = (kx+Gx)*dx;
	double ct_y = (ky+Gy)*dy;
	double ct_z = (kz+Gz)*dz;

	for (i = 0; i < n_Mask_Wigner; i++)
	{
	  	x_aver_re += X[Mask_Wigner[i][0]][Mask_Wigner[i][1]][Mask_Wigner[i][2]]*
				        cos( ct_x*Mask_Wigner[i][0] + ct_y*Mask_Wigner[i][1] + ct_z*Mask_Wigner[i][2]);
		x_aver_im += X[Mask_Wigner[i][0]][Mask_Wigner[i][1]][Mask_Wigner[i][2]]*
				        sin( ct_x*Mask_Wigner[i][0] + ct_y*Mask_Wigner[i][1] + ct_z*Mask_Wigner[i][2]);
	}
	*X_avBr_re = x_aver_re/n_Mask_Wigner;
	*X_avBr_im = x_aver_im/n_Mask_Wigner;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Average the energy along x and z directions and time
///////////////////////////////////////////////////////////////////////////////////////
int aver_absEorH(double ***Ux, double ***Uy, double ***Uz, long nx_a_av, long nx_b_av, long ny_a_av, long ny_b_av,
			     long nz_a_av, long nz_b_av, double *U_av, double Vol)
{
	long i, j, k;
	double U_aver = 0.0, U;

 	for (i = nx_a_av; i <= nx_b_av; i++)
	{
		for (j = ny_a_av; j <= ny_b_av; j++)
		{
	  		for (k = nz_a_av; k <= nz_b_av; k++)
			{
				U = sqrt(Ux[i][j][k]*Ux[i][j][k] + Uy[i][j][k]*Uy[i][j][k] + Uz[i][j][k]*Uz[i][j][k] );
				U_aver += U;
	  		}
		}
	}
	*U_av = U_aver/Vol;

	return 0;	
}

///////////////////////////////////////////////////////////////////////////////////////
//Fourier transform of a field component in the specified volume
///////////////////////////////////////////////////////////////////////////////////////
int fourier_transf_vol(double ***U, long nx_a_f, long nx_b_f, long ny_a_f, long ny_b_f, 
				       long nz_a_f, long nz_b_f, double **FrT_re, double **FrT_im, 
				       double *FrT_om, long n_FrT_om, double time)
{
	long i, j, k, ff, pp;
	double coef_1, coef_2;
		 
	for (ff = 0; ff < n_FrT_om; ff++)
	{
		coef_1 = cos(FrT_om[ff]*time);
		coef_2 = sin(FrT_om[ff]*time);
		pp = 0;
 		for (i = nx_a_f; i <= nx_b_f; i++)
		{
			for (j = ny_a_f; j <= ny_b_f; j++)
			{
	  			for (k = nz_a_f; k <= nz_b_f; k++)
				{
					FrT_re[pp][ff] += U[i][j][k]*coef_1; 
					FrT_im[pp][ff] += U[i][j][k]*coef_2;
				}
				pp++;
	  		}
		}
	}

	return 0;	
}


///////////////////////////////////////////////////////////////////////////////////////
//Linear interpolation
///////////////////////////////////////////////////////////////////////////////////////
int interp_1D(double *x, int length_x, double *f, double *x_i, int length_x_i, double *g)
{
    int i, j; 
	int n_x_i, n_x;
	n_x_i = length_x_i;
	n_x = length_x;

	for (i = 0; i < n_x_i; i++)
	{
		if( (x_i[i] >= x[0]) && (x_i[i] <= x[n_x-1]))
		{
			for(j = 1; j < n_x; j++)
			{
				if( x_i[i] <= x[j] )
				{
					g[i] = (x_i[i]-x[j-1]) / (x[j]-x[j-1]) * (f[j]-f[j-1]) + f[j-1];
					break;
				}
			}
		}
		else
		{
			//cannot interpolate -- extrapolation would be needed
			return 1;
		}
	}

    return 0;
}