fdtd_3d.cpp

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

				  double &ct_Ez_1, double &ct_Ez_2, double &ct_Hx_1, double &ct_Hx_2, 
				  double &ct_Hy_1, double &ct_Hy_2, double &ct_Hz_1, double &ct_Hz_2)
{
	double mu_0  = 1.2566370614359172953850573533118011536788677597500423283899778369231265625144835994512139301368468271e-6;  // [H/m]
	long i, j, k;

	double cos_teta, sin_teta, cos_phi, sin_phi, sin_gamma, cos_gamma;
  
	cos_teta = cos(teta);
	sin_teta = sin(teta);
	
	cos_phi = cos(phi);
	sin_phi = sin(phi);

	cos_gamma = cos(gamma);
	sin_gamma = sin(gamma);

	*ll_1D_E = (double *) calloc(l_1D,sizeof(double));
	*ll_1D_H = (double *) calloc(l_1D-1,sizeof(double));

	*Hz_i0 = Init_Matrix_2D<double>(l_y-1,l_z);
	*Hy_i0 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*Hz_i1 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*Hy_i1 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*Hz_j0 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*Hx_j0 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*Hz_j1 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*Hx_j1 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*Hy_k0 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*Hx_k0 = Init_Matrix_2D<double>(l_x,l_y-1); 
	*Hy_k1 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*Hx_k1 = Init_Matrix_2D<double>(l_x,l_y-1); 

	*face_Hz_i0 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*face_Hy_i0 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*face_Hz_i1 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*face_Hy_i1 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*face_Hz_j0 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*face_Hx_j0 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*face_Hz_j1 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*face_Hx_j1 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*face_Hy_k0 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*face_Hx_k0 = Init_Matrix_2D<double>(l_x,l_y-1); 
	*face_Hy_k1 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*face_Hx_k1 = Init_Matrix_2D<double>(l_x,l_y-1); 

	*Ey_i0 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*Ez_i0 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*Ey_i1 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*Ez_i1 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*Ex_j0 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*Ez_j0 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*Ex_j1 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*Ez_j1 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*Ex_k0 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*Ey_k0 = Init_Matrix_2D<double>(l_x,l_y-1); 
	*Ex_k1 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*Ey_k1 = Init_Matrix_2D<double>(l_x,l_y-1); 

	*face_Ey_i0 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*face_Ez_i0 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*face_Ey_i1 = Init_Matrix_2D<double>(l_y-1,l_z); 
	*face_Ez_i1 = Init_Matrix_2D<double>(l_y,l_z-1); 
	*face_Ex_j0 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*face_Ez_j0 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*face_Ex_j1 = Init_Matrix_2D<double>(l_x-1,l_z); 
	*face_Ez_j1 = Init_Matrix_2D<double>(l_x,l_z-1); 
	*face_Ex_k0 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*face_Ey_k0 = Init_Matrix_2D<double>(l_x,l_y-1); 
	*face_Ex_k1 = Init_Matrix_2D<double>(l_x-1,l_y); 
	*face_Ey_k1 = Init_Matrix_2D<double>(l_x,l_y-1);

	if (!*ll_1D_E || !*ll_1D_H || !*Hz_i0 || !*Hy_i0 || !*Hz_i1 || !*Hy_i1 || !*Hz_j0 || 
		!*Hx_j0 || !*Hz_j1 || !*Hx_j1 || !*Hy_k0 || !*Hx_k0 || !*Hy_k1 || !*Hx_k1 || 
		!*face_Hz_i0 || !*face_Hy_i0 || !*face_Hz_i1 || !*face_Hy_i1 || !*face_Hz_j0 || 
		!*face_Hx_j0 || !*face_Hz_j1 || !*face_Hx_j1 || !*face_Hy_k0 || !*face_Hx_k0 || 
		!*face_Hy_k1 || !*face_Hx_k1 || !*Ey_i0 || !*Ez_i0 || !*Ey_i1 || !*Ez_i1 || 
		!*Ex_j0 || !*Ez_j0 || !*Ex_j1 || !*Ez_j1 || !*Ex_k0 || !*Ey_k0 || !*Ex_k1 || 
		!*Ey_k1 || !*face_Ey_i0 || !*face_Ez_i0 || !*face_Ey_i1 || !*face_Ez_i1 || 
		!*face_Ex_j0 || !*face_Ez_j0 || !*face_Ex_j1 || !*face_Ez_j1 || !*face_Ex_k0 || 
		!*face_Ey_k0 || !*face_Ex_k1 || !*face_Ey_k1)
	{
		return 1; //return in case of unsuccesfull memory allocation
	}

	for (i = 0; i < l_1D-1; i++)
	{
		ll_1D_E[0][i] = (i - n_PML - 5)*d;
		ll_1D_H[0][i] = (i + 0.5 - n_PML - 5)*d;
	}
	ll_1D_E[0][l_1D-1] = (l_1D-1 - n_PML - 5)*d;

	//////////////////////////////////////////////////////////////////////////
	// 0 <= teta <= 90
	//////////////////////////////////////////////////////////////////////////
	if ( sin_teta >= 0 && cos_teta >= 0 )
	{
		//////////////////////////////////////////////////////////////////////////
		// 0 <= teta <= 90;  0 <= phi <= 90
		//////////////////////////////////////////////////////////////////////////
		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*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Ez_i1[0][j][k] = (l_x-1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					
					face_Hy_i0[0][j][k] = -0.5*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;   
					face_Hy_i1[0][j][k] = (l_x-1+0.5)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*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)*dy*sin_teta*sin_phi + k*dz*cos_teta;   
					face_Ey_i1[0][j][k] = (l_x-1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*dz*cos_teta;    

					face_Hz_i0[0][j][k] = -0.5*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*dz*cos_teta;  
					face_Hz_i1[0][j][k] = (l_x-1+0.5)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*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 + (k+0.5)*dz*cos_teta; 
					face_Ez_j1[0][i][k] =  i*dx*sin_teta*cos_phi + (l_y-1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta; 

					face_Hx_j0[0][i][k] =  i*dx*sin_teta*cos_phi - 0.5*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Hx_j1[0][i][k] =  i*dx*sin_teta*cos_phi + (l_y-1+0.5)*dy*sin_teta*sin_phi + (k+0.5)*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 + k*dz*cos_teta;  
					face_Ex_j1[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + (l_y-1)*dy*sin_teta*sin_phi + k*dz*cos_teta; 

					face_Hz_j0[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi - 0.5*dy*sin_teta*sin_phi + k*dz*cos_teta; 
					face_Hz_j1[0][i][k] = (i+0.5)*dx*sin_teta*cos_phi + (l_y-1+0.5)*dy*sin_teta*sin_phi + k*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*dy*sin_teta*sin_phi; 
					face_Ex_k1[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi +(l_z-1)*dz*cos_teta; 

					face_Hy_k0[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi - 0.5*dz*cos_teta; 
					face_Hy_k1[0][i][j] = (i+0.5)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi +(l_z-1+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)*dy*sin_teta*sin_phi; 
					face_Ey_k1[0][i][j] =  i*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi +(l_z-1)*dz*cos_teta;
										
					face_Hx_k0[0][i][j] =  i*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi - 0.5*dz*cos_teta; 
					face_Hx_k1[0][i][j] =  i*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi +(l_z-1+0.5)*dz*cos_teta; 
				}
			}			
		}
		
		//////////////////////////////////////////////////////////////////////////
		// 0 <= teta <= 90;  90 < phi <= 180
		//////////////////////////////////////////////////////////////////////////
		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] =  (-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Ez_i1[0][j][k] = j*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					
					face_Hy_i0[0][j][k] = (-0.5-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;   
					face_Hy_i1[0][j][k] = 0.5*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi + (k+0.5)*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] = (-l_x+1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*dz*cos_teta;   
					face_Ey_i1[0][j][k] = (j+0.5)*dy*sin_teta*sin_phi + k*dz*cos_teta;    

					face_Hz_i0[0][j][k] = (-0.5-l_x+1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*dz*cos_teta;  
					face_Hz_i1[0][j][k] =  0.5*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi + k*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-l_x+1)*dx*sin_teta*cos_phi + (k+0.5)*dz*cos_teta; 
					face_Ez_j1[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi + (l_y-1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta; 

					face_Hx_j0[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi - 0.5*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Hx_j1[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi + (l_y-1+0.5)*dy*sin_teta*sin_phi + (k+0.5)*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-l_x+1)*dx*sin_teta*cos_phi + k*dz*cos_teta;  
					face_Ex_j1[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (l_y-1)*dy*sin_teta*sin_phi + k*dz*cos_teta; 

					face_Hz_j0[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi - 0.5*dy*sin_teta*sin_phi + k*dz*cos_teta; 
					face_Hz_j1[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (l_y-1+0.5)*dy*sin_teta*sin_phi + k*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-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi; 
					face_Ex_k1[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi +(l_z-1)*dz*cos_teta; 

					face_Hy_k0[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi - 0.5*dz*cos_teta; 
					face_Hy_k1[0][i][j] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + j*dy*sin_teta*sin_phi +(l_z-1+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)*dy*sin_teta*sin_phi; 
					face_Ey_k1[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi +(l_z-1)*dz*cos_teta;
										
					face_Hx_k0[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi - 0.5*dz*cos_teta; 
					face_Hx_k1[0][i][j] =  (i-l_x+1)*dx*sin_teta*cos_phi + (j+0.5)*dy*sin_teta*sin_phi +(l_z-1+0.5)*dz*cos_teta; 
				}
			}			
		}

		//////////////////////////////////////////////////////////////////////////
		// 0 <= teta <= 90;  180 < phi <= 270
		//////////////////////////////////////////////////////////////////////////
		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] =  (-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Ez_i1[0][j][k] = (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					
					face_Hy_i0[0][j][k] = (-0.5-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;   
					face_Hy_i1[0][j][k] = 0.5*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi + (k+0.5)*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] = (-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + k*dz*cos_teta;   
					face_Ey_i1[0][j][k] = (j+0.5-l_y+1)*dy*sin_teta*sin_phi + k*dz*cos_teta;    

					face_Hz_i0[0][j][k] = (-0.5-l_x+1)*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + k*dz*cos_teta;  
					face_Hz_i1[0][j][k] =  0.5*dx*sin_teta*cos_phi + (j+0.5-l_y+1)*dy*sin_teta*sin_phi + k*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-l_x+1)*dx*sin_teta*cos_phi + (-l_y+1)*dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta; 
					face_Ez_j1[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi + (k+0.5)*dz*cos_teta; 

					face_Hx_j0[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi + (-0.5-l_y+1) *dy*sin_teta*sin_phi + (k+0.5)*dz*cos_teta;  
					face_Hx_j1[0][i][k] =  (i-l_x+1)*dx*sin_teta*cos_phi +  0.5*dy*sin_teta*sin_phi + (k+0.5)*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-l_x+1)*dx*sin_teta*cos_phi + (-l_y+1)*dy*sin_teta*sin_phi + k*dz*cos_teta;  
					face_Ex_j1[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + k*dz*cos_teta; 

					face_Hz_j0[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + (-0.5-l_y+1)*dy*sin_teta*sin_phi + k*dz*cos_teta; 
					face_Hz_j1[0][i][k] = (i+0.5-l_x+1)*dx*sin_teta*cos_phi + 0.5*dy*sin_teta*sin_phi + k*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-l_x+1)*dx*sin_teta*cos_phi + (j-l_y+1)*dy*sin_teta*sin_phi; 
					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 +(l_z-1)*dz*cos_teta; 

					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*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 +(l_z-1+0.5)*dz*cos_teta;
				}