fdtd_3d.cpp

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

			//k faces
			for (i = 0; i < l_x-1; i++)
			{
				for (j = 0; j < l_y; j++)
				{
					face_Ex_k0[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (-l_z+1)*dz*cos_teta; 
					face_Ex_k1[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi; 

					face_Hy_k0[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (-0.5-l_z+1)*dz*cos_teta; 
					face_Hy_k1[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + 0.5*dz*cos_teta;
				}
			}

			//k faces
			for (i = 0; i < l_x; i++)
			{
				for (j = 0; j < l_y-1; j++)
				{
					face_Ey_k0[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi +  (-l_z+1)*dz*cos_teta; 
					face_Ey_k1[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi;
										
					face_Hx_k0[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (-0.5-l_z+1)*dz*cos_teta; 
					face_Hx_k1[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + 0.5*dz*cos_teta; 
				}
			}			
		}

		//////////////////////////////////////////////////////////////////////////
		// 90 < teta < 180;  270 < phi < 360
		//////////////////////////////////////////////////////////////////////////
		if ( sin_phi < 0 && cos_phi > 0 )
		{
			//i faces
			for (j = 0; j < l_y; j++)
			{
				for (k = 0; k < l_z-1; k++)
				{
					face_Ez_i0[0][j][k] = (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta;  
					face_Ez_i1[0][j][k] = (l_x-1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta;  
					
					face_Hy_i0[0][j][k] = -0.5*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta;   
					face_Hy_i1[0][j][k] = (l_x-1+0.5)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta;    
				}
			}

			//i faces
			for (j = 0; j < l_y-1; j++)
			{
				for (k = 0; k < l_z; k++)
				{
					face_Ey_i0[0][j][k] = (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta;   
					face_Ey_i1[0][j][k] = (l_x-1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta;    

					face_Hz_i0[0][j][k] = -0.5*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta;  
					face_Hz_i1[0][j][k] = (l_x-1+0.5)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta;  
				}
			}
			
			
			//j faces
			for (i = 0; i < l_x; i++)
			{
				for (k = 0; k < l_z-1; k++)
				{
					face_Ez_j0[0][i][k] =  i*dx*sin_teta*cos_phi + (-l_y+1)*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta; 
					face_Ez_j1[0][i][k] =  i*dx*sin_teta*cos_phi + (k+0.5-l_z+1)*dz*cos_teta; 

					face_Hx_j0[0][i][k] =  i*dx*sin_teta*cos_phi + (-0.5-l_y+1) *dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta;  
					face_Hx_j1[0][i][k] =  i*dx*sin_teta*cos_phi +  0.5*dy*sin_teta*sin_phi + (k+0.5-l_z+1)*dz*cos_teta; 
				}
			}

			//j faces
			for (i = 0; i < l_x-1; i++)
			{
				for (k = 0; k < l_z; k++)
				{
					face_Ex_j0[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + (-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta;  
					face_Ex_j1[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + (k-l_z+1)*dz*cos_teta; 

					face_Hz_j0[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + (-0.5-l_y+1)*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta; 
					face_Hz_j1[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + 0.5*dy*sin_teta*sin_phi + (k-l_z+1)*dz*cos_teta; 
					
				}
			}
			
			//k faces
			for (i = 0; i < l_x-1; i++)
			{
				for (j = 0; j < l_y; j++)
				{
					face_Ex_k0[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (-l_z+1)*dz*cos_teta; 
					face_Ex_k1[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi;

					face_Hy_k0[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (-0.5-l_z+1)*dz*cos_teta; 
					face_Hy_k1[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + 0.5*dz*cos_teta;
				}
			}

			//k faces
			for (i = 0; i < l_x; i++)
			{
				for (j = 0; j < l_y-1; j++)
				{
					face_Ey_k0[0][i][j] = i*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (-l_z+1)*dz*cos_teta; 
					face_Ey_k1[0][i][j] = i*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi;
										
					face_Hx_k0[0][i][j] = i*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + (-0.5-l_z+1)*dz*cos_teta;
					face_Hx_k1[0][i][j] = i*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + 0.5*dz*cos_teta; 
				}
			}			
		}
	}
	
	double inv_dx = 1.0/dx, inv_dy = 1.0/dy, inv_dz = 1.0/dz;

	//some constants to reduce the computational time
	ct_Ex_1 = inv_dy*(-cos_gamma*sin_teta); //Hz_inc
	ct_Ex_2 = inv_dz*(-sin_gamma*cos_phi + cos_gamma*cos_teta*sin_phi); //Hy_inc
    ct_Ey_1 = inv_dx*(-cos_gamma*sin_teta); //Hz_inc
    ct_Ey_2 = inv_dz*(sin_gamma*sin_phi + cos_gamma*cos_teta*cos_phi); //Hx_inc
	ct_Ez_1 = inv_dx*(-sin_gamma*cos_phi + cos_gamma*cos_teta*sin_phi); //Hy_inc
	ct_Ez_2 = inv_dy*(sin_gamma*sin_phi + cos_gamma*cos_teta*cos_phi); //Hx_inc
	
	ct_Hx_1 = dt*inv_dy*(sin_gamma*sin_teta)/mu_0; //Ez_inc
	ct_Hx_2 = dt*inv_dz*(-cos_gamma*cos_phi - sin_gamma*cos_teta*sin_phi)/mu_0; //Ey_inc
	ct_Hy_1 = dt*inv_dx*(sin_gamma*sin_teta)/mu_0; //Ez_inc
	ct_Hy_2 = dt*inv_dz*(cos_gamma*sin_phi - sin_gamma*cos_teta*cos_phi)/mu_0; //Ex_inc
	ct_Hz_1 = dt*inv_dx*(-cos_gamma*cos_phi - sin_gamma*cos_teta*sin_phi)/mu_0; //Ey_inc
	ct_Hz_2 = dt*inv_dy*(cos_gamma*sin_phi - sin_gamma*cos_teta*cos_phi)/mu_0; //Ex_inc

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//calculates the phase velocity with Newton_Rhapson scheme
///////////////////////////////////////////////////////////////////////////////////////
int phase_velocity(double dx, double dy, double dz, double dt, 
				   double eps_r, double mu_r, double om, 
				   double teta, double phi, double *phase_vel)
{
	//permittivity of free space 
	double eps_0 = 8.8541878176203898505365630317107502606083701665994498081024171524053950954599821142852891607182008932e-12; // [F/m]
	//permeability of free space 
    double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]


	double vel_light = 1.0/sqrt(eps_0*eps_r*mu_0*mu_r);
	
	double dx2 = dx*dx;
	double dy2 = dy*dy;
	double dz2 = dz*dz;  

	double A = dx*sin(teta)*cos(phi)/2.0;
	double B = dy*sin(teta)*sin(phi)/2.0;
	double C = dz*cos(teta)/2.0;
	
	double D = sin(om*dt/2.0)/(vel_light*dt);
	D = D*D;
	double c_a, c_b, c_c;
	
	long cik = 0;
	
	double k_r = 0.0;
	double k = om/vel_light;
	while ( abs(k-k_r) > 1e-6 && cik < 1e5 )
	{
		k_r = k;
		c_a = sin(A*k)/dx;
		c_b = sin(B*k)/dy;
		c_c = sin(C*k)/dz;

		k -= ( c_a*c_a + c_b*c_b + c_c*c_c - D)/( A*sin(2.0*A*k)/dx2 + B*sin(2.0*B*k)/dy2 + C*sin(2.0*C*k)/dz2);
		cik++;
	}
  
	*phase_vel = om/k;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ex on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Ex(double ***Ex, long ***Ind, double *K_b, double **Hz_j0, double **Hz_j1,
		   double **Hy_k0, double **Hy_k1, double ct_Ex_1, double ct_Ex_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 i, j, k;
	long ind_i, ind_j, ind_k;
	
	//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;
			Ex[i][n_TS_ya][k] -= K_b[Ind[i][n_TS_ya][k]]*Hz_j0[ind_i][ind_k]*ct_Ex_1;
			Ex[i][n_TS_yb][k] += K_b[Ind[i][n_TS_yb][k]]*Hz_j1[ind_i][ind_k]*ct_Ex_1;
		}
	}

	//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;
			Ex[i][j][n_TS_za] += K_b[Ind[i][j][n_TS_za]]*Hy_k0[ind_i][ind_j]*ct_Ex_2;
			Ex[i][j][n_TS_zb] -= K_b[Ind[i][j][n_TS_zb]]*Hy_k1[ind_i][ind_j]*ct_Ex_2;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ey on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Ey(double ***Ey, long ***Ind, double *K_b, double **Hz_i0, double **Hz_i1,
		   double **Hx_k0, double **Hx_k1, double ct_Ey_1, double ct_Ey_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 i, j, k;
	long ind_i, ind_j, ind_k;

	//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;
			Ey[n_TS_xa][j][k] += K_b[Ind[n_TS_xa][j][k]]*Hz_i0[ind_j][ind_k]*ct_Ey_1;
			Ey[n_TS_xb][j][k] -= K_b[Ind[n_TS_xb][j][k]]*Hz_i1[ind_j][ind_k]*ct_Ey_1;
		}
	}

	//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;
			Ey[i][j][n_TS_za] -= K_b[Ind[i][j][n_TS_za]]*Hx_k0[ind_i][ind_j]*ct_Ey_2;
			Ey[i][j][n_TS_zb] += K_b[Ind[i][j][n_TS_zb]]*Hx_k1[ind_i][ind_j]*ct_Ey_2;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Ez on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Ez(double ***Ez, long  ***Ind, double *K_b, double **Hy_i0, double **Hy_i1,
		   double **Hx_j0, double **Hx_j1, double ct_Ez_1, double ct_Ez_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 i, j, k;
	long ind_i, ind_j, ind_k;

	//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;
			Ez[n_TS_xa][j][k] -= K_b[Ind[n_TS_xa][j][k]]*Hy_i0[ind_j][ind_k]*ct_Ez_1;
			Ez[n_TS_xb][j][k] += K_b[Ind[n_TS_xb][j][k]]*Hy_i1[ind_j][ind_k]*ct_Ez_1;
		}
	}

	//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;
			Ez[i][n_TS_ya][k] += K_b[Ind[i][n_TS_ya][k]]*Hx_j0[ind_i][ind_k]*ct_Ez_2;
			Ez[i][n_TS_yb][k] -= K_b[Ind[i][n_TS_yb][k]]*Hx_j1[ind_i][ind_k]*ct_Ez_2;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////
//Add the incident field to Hx on total-field scattered field boundary
///////////////////////////////////////////////////////////////////////////////////////
void TS_Hx(double ***Hx, long  ***Ind, double **Mat, double **Ez_j0, double **Ez_j1,
		   double **Ey_k0, double **Ey_k1, double ct_Hx_1, double ct_Hx_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_ya_MIN_1, long n_TS_za_MIN_1)
{
	long i, j, k;
	long ind_i, ind_j, ind_k;
	double 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];
			Hx[i][n_TS_ya_MIN_1][k] += Ez_j0[ind_i][ind_k]*ct_Hx_1/mu_r;
			mu_r = Mat[Ind[i][n_TS_yb][k]][1];
			Hx[i][n_TS_yb][k] -= Ez_j1[ind_i][ind_k]*ct_Hx_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];
			Hx[i][j][n_TS_za_MIN_1] -= Ey_k0[ind_i][ind_j]*ct_Hx_2/mu_r;
			mu_r = Mat[Ind[i][j][n_TS_zb]][1];
			Hx[i][j][n_TS_zb] += Ey_k1[ind_i][ind_j]*ct_Hx_2/mu_r;
		}
	}

	return;
}

///////////////////////////////////////////////////////////////////////////////////////