fdtd_xypbc_zpml.cpp

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

}

//Init Point Source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_ptSource(long **Coord, long nCoord, double **Par_Excit, long s_type)
{
	Coord_ptSource = Coord;
	n_ptSource = nCoord;
	Par_ptSource = Par_Excit;
	source_type = s_type;
	jel_plane_wave = 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Curent source J
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::ptSource_J(double ***&xx, double **&Par_ptS, double time)
{
	long i, i1, j1, k1;
	double J0, omega, phi, t0, tw;
	double dtDIVeps0 = dt/eps_0;
	double alfa;
	//time = time - dt/2;

	for (i = 0; i <n_ptSource; i++ )
	{
		i1 = Coord_ptSource[i][0];
		j1 = Coord_ptSource[i][1];
		k1 = Coord_ptSource[i][2];
		J0 = Par_ptS[i][0];
		switch (source_type)
		{
			case 1: //Gaussian pulse
				tw = Par_ptS[i][1];
				t0 = Par_ptS[i][2];
				xx[i1][j1][k1] += J0*exp( -pow( (time - t0)/tw ,2) );
				break;
			case 2: //Sinusoidal plane wave
				omega = Par_ptS[i][1];
				phi   = Par_ptS[i][2];
				alfa = omega*3.0/(2.0*pi);
				xx[i1][j1][k1] += J0*(1-exp(-alfa*time))*cos(omega*time + phi);
				break;
			case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
				tw    = Par_ptS[i][1];
				t0    = Par_ptS[i][2];
				omega = Par_ptS[i][3];
				phi   = Par_ptS[i][4];
				xx[i1][j1][k1] += J0*cos( omega*(time-t0) + phi)*
					              exp( -pow( (time - t0)/tw ,2) );
				break;
		}	
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_xyPBC_zPML::Init_Followed(long **Ind_Followed, long n, long n_t)
{
	Ind_Foll = Ind_Followed;
	length_Ind_Foll = n;
	n_tot = n_t;

	Hx_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hx_Foll)	
	{
		Free_Mem();
		return false;
	}

	Hy_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hy_Foll)	
	{
		Free_Mem();
		return false;
	}

	Hz_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hz_Foll)	
	{
		Free_Mem();
		return false;
	}

	Ex_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ex_Foll)	
	{
		Free_Mem();
		return false;
	}

	Ey_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ey_Foll)	
	{
		Free_Mem();
		return false;
	}

	Ez_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ez_Foll)	
	{
		Free_Mem();
		return false;
	}
	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init the geometry of the Wigner cell
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Init_Fourier_Transf(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 Frec_1, double Frec_2,
	                                      long N_frec, long Start_fourier)
{
	long i;

	nx_a_f = Nx_a_f;
	nx_b_f = Nx_b_f; 
	ny_a_f = Ny_a_f;
	ny_b_f = Ny_b_f; 
	nz_a_f = Nz_a_f;
	nz_b_f = Nz_b_f; 
	
	frec_1 = Frec_1;
	frec_2 = Frec_2;
	n_frec = N_frec; 
	start_fourier = Start_fourier;

	n_f_Points = (nx_b_f - nx_a_f + 1)*(ny_b_f - ny_a_f + 1)*(nz_b_f - nz_a_f + 1);

	Ex_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Ex_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Ey_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Ey_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Ez_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Ez_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);

	Hx_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Hx_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Hy_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Hy_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Hz_fourier_TR_re = Init_Matrix_2D<double>(n_f_Points,n_frec);
	Hz_fourier_TR_im = Init_Matrix_2D<double>(n_f_Points,n_frec);
	fourier_omega =  (double *) calloc(n_frec,sizeof(double));

	if (!Ex_fourier_TR_re || !Ex_fourier_TR_im || !Ey_fourier_TR_re || !Ey_fourier_TR_im ||
		!Ez_fourier_TR_re || !Ez_fourier_TR_im || !Hx_fourier_TR_re || !Hx_fourier_TR_im ||
		!Hy_fourier_TR_re || !Hy_fourier_TR_im || !Hz_fourier_TR_re || !Hz_fourier_TR_im ||
		!fourier_omega)
	{
		cout << "Memory allocation problem in Init_Fourier_Transf" << endl;
		return 1;
	}

	double delta_omega = TwoORpi*(frec_2 - frec_1)/(n_frec-1);
	fourier_omega[0] = TwoORpi*frec_1;
	for (i = 1; i< n_frec; i++)
	{
		fourier_omega[i] = fourier_omega[i-1] + delta_omega;
	}

	jel_fourier_transf_vol = 1;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init the geometry of the Wigner cell
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Init_Mask_Wigner(long ***&mask, long nx_W, long ny_W, long nz_W)
{
	Mask_Wigner = mask;

	long i, j, k;
	Vol_av = 0;
	for ( i = 0; i < nx_W; i++)
	{
		for ( j = 0; j < ny_W; j++)
		{
			for ( k = 0; k < nz_W; k++)
			{
				if (Mask_Wigner[i][j][k] >0)
					Vol_av++;
			}
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init for the averaged field components
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Init_Aver_xyz(long Nx_a_av, long Nx_b_av, long Ny_a_av, 
		                            long Ny_b_av, long Nz_a_av, long Nz_b_av, 
									long N_iter)
{
	nx_a_av = Nx_a_av;
	nx_b_av = Nx_b_av;
	ny_a_av = Ny_a_av;
	ny_b_av = Ny_b_av;
	nz_a_av = Nz_a_av;
	nz_b_av = Nz_b_av; 

	//Vol_av = (nx_b_av - nx_a_av + 1)*(ny_b_av - ny_a_av + 1)*(nz_b_av - nz_a_av + 1);

	Ex_Foll_av = (double *)calloc(N_iter,sizeof(double));
	Ey_Foll_av = (double *)calloc(N_iter,sizeof(double));
	Ez_Foll_av = (double *)calloc(N_iter,sizeof(double));
	Hx_Foll_av = (double *)calloc(N_iter,sizeof(double));
	Hy_Foll_av = (double *)calloc(N_iter,sizeof(double));
	Hz_Foll_av = (double *)calloc(N_iter,sizeof(double));

	if (!Ex_Foll_av || !Ey_Foll_av || !Ez_Foll_av ||
		!Hx_Foll_av || !Hy_Foll_av || !Hz_Foll_av)
	{
		Free_Mem();
		return 1;
	}

	aver_field_volume = 1;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Bragg Dispersion
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Init_aver_Bragg_WignerC(long nz_av, double K_x, double K_y, 
		                                      double **&G_xy, long n_G_xy, long N_iter)
{
	long i, j;

	kx = K_x;
	ky = K_y;
	
	Gxy = G_xy;
	n_Gxy = n_G_xy;

	nz_Bragg = nz_av; 

	Surf_Bragg = 0;
	for (i = 0; i<nx; i++)
	{
		for (j = 0; j<ny; j++)
		{
			if ( Mask_Wigner[i][j][0] > 0 )
	          Surf_Bragg++;
		}
	}

	Ex_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Ex_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Ey_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Ey_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Ez_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Ez_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hx_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hx_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hy_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hy_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hz_Br_av_re = Init_Matrix_2D<double>(n_Gxy,N_iter);
	Hz_Br_av_im = Init_Matrix_2D<double>(n_Gxy,N_iter);
	
	if (!Ex_Br_av_re || !Ex_Br_av_im || !Ey_Br_av_re || !Ey_Br_av_im ||
		!Ez_Br_av_re || !Ez_Br_av_im || !Hx_Br_av_re || !Hx_Br_av_im ||
		!Hy_Br_av_re || !Hy_Br_av_im || !Hz_Br_av_re || !Hz_Br_av_im)
	{
		Free_Mem();
		return 1;
	}

	jel_aver_Bragg_WignerUC = 1;

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Set_Data_Followed(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];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Get_Data_Followed(double **&h_x, double **&h_y, double **&h_z, 
										 double **&e_x, double **&e_y, double **&e_z)
{
	h_x = Hx_Foll;
	h_y = Hy_Foll;
	h_z = Hz_Foll;
	e_x = Ex_Foll;
	e_y = Ey_Foll;
	e_z = Ez_Foll;
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Get_Data_Aver(double *&hx_av, double *&hy_av, double *&hz_av, 
									 double *&ex_av, double *&ey_av, double *&ez_av)
{
	hx_av = Hx_Foll_av;
	hy_av = Hy_Foll_av;
	hz_av = Hz_Foll_av;
	ex_av = Ex_Foll_av;
	ey_av = Ey_Foll_av;
	ez_av = Ez_Foll_av;
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Fourier data components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Get_Data_Fourier(double **&ex_fourier_TR_re, double **&ex_fourier_TR_im,
										double **&ey_fourier_TR_re, double **&ey_fourier_TR_im,
										double **&ez_fourier_TR_re, double **&ez_fourier_TR_im,
										double **&hx_fourier_TR_re, double **&hx_fourier_TR_im,
										double **&hy_fourier_TR_re, double **&hy_fourier_TR_im,
										double **&hz_fourier_TR_re, double **&hz_fourier_TR_im,
										double *&omega)
{
	ex_fourier_TR_re = Ex_fourier_TR_re;
	ex_fourier_TR_im = Ex_fourier_TR_im;
	ey_fourier_TR_re = Ey_fourier_TR_re;
	ey_fourier_TR_im = Ey_fourier_TR_im;
	ez_fourier_TR_re = Ez_fourier_TR_re;
	ez_fourier_TR_im = Ez_fourier_TR_im;

	hx_fourier_TR_re = Hx_fourier_TR_re;
	hx_fourier_TR_im = Hx_fourier_TR_im;
	hy_fourier_TR_re = Hy_fourier_TR_re;
	hy_fourier_TR_im = Hy_fourier_TR_im;
	hz_fourier_TR_re = Hz_fourier_TR_re;
	hz_fourier_TR_im = Hz_fourier_TR_im;

	omega = fourier_omega;

}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_xyPBC_zPML::Get_Data_Bragg_Aver(double **&ex_avBr_re, double **&ex_avBr_im, double **&ey_avBr_re, 
		                                   double **&ey_avBr_im, double **&ez_avBr_re, double **&ez_avBr_im,
							               double **&hx_avBr_re, double **&hx_avBr_im, double **&hy_avBr_re, 
							               double **&hy_avBr_im, double **&hz_avBr_re, double **&hz_avBr_im)
{
	ex_avBr_re = Ex_Br_av_re;  ex_avBr_im = Ex_Br_av_im; 
	ey_avBr_re = Ey_Br_av_re;  ey_avBr_im = Ey_Br_av_im;
	ez_avBr_re = Ez_Br_av_re;  ez_avBr_im = Ez_Br_av_im;
	hx_avBr_re = Hx_Br_av_re;  hx_avBr_im = Hx_Br_av_im; 
	hy_avBr_re = Hy_Br_av_re;  hy_avBr_im = Hy_Br_av_im;
	hz_avBr_re = Hz_Br_av_re;  hz_avBr_im = Hz_Br_av_im;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save the field components
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Save_Field(long iter)
{
	if (save_3D_binary(Ex,nx,ny,nz,iter,path_name_Ex))
		return -1;
	if (save_3D_binary(Ey,nx,ny,nz,iter,path_name_Ey))
		return -1;
	if (save_3D_binary(Ez,nx,ny,nz-1,iter,path_name_Ez))
		return -1;
	if (save_3D_binary(Hx,nx,ny,nz-1,iter,path_name_Hx))
		return -1;
	if (save_3D_binary(Hy,nx,ny,nz-1,iter,path_name_Hy))
		return -1;
	if (save_3D_binary(Hz,nx,ny,nz,iter,path_name_Hz))
		return -1;
	
	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init to save slices from the field components
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_xyPBC_zPML::Init_Save_FieldSlices(long nxa, long nxb, long nya, long nyb, 
										     long nza, long nzb, long nZ_xy, long nY_xz, 
											 long nX_yz, long t1, long t2, long nr_Save, 
											 char *path)
{
	n_x_a = nxa; 
	n_x_b = nxb;
	n_y_a = nya; 
	n_y_b = nyb;
	n_z_a = nza; 
	n_z_b = nzb;

	nz_xy = nZ_xy;
	ny_xz = nY_xz; 
	nx_yz = nX_yz;

	saveFROMinst = t1;
    saveTOinst   = t2;
	nr_save = nr_Save;
	jel_Save_Slice = 1;

	path_name_Ex_1D = (char *) calloc(512,sizeof(char)); 
	if (!path_name_Ex_1D)
	{
		Free_Mem();
		return false;
	}
	
	path_name_Ex = (char *) calloc(512,sizeof(char)); 
	if (!path_name_Ex)
	{
		Free_Mem();
		return false;
	}
	path_name_Ey = (char *) calloc(512,sizeof(char));
	if (!path_name_Ey)
	{
		Free_Mem();
		return false;
	}
	path_name_Ez = (char *) calloc(512,sizeof(char));
	if (!path_name_Ez)
	{
		Free_Mem();
		return false;
	}
	path_name_Hx = (char *) calloc(512,sizeof(char)); 
	if (!path_name_Hx)
	{
		Free_Mem();
		return false;
	}
	path_name_Hy = (char *) calloc(512,sizeof(char));
	if (!path_name_Hy)
	{
		Free_Mem();
		return false;
	}
	path_name_Hz = (char *) calloc(512,sizeof(char));
	if (!path_name_Hz)
	{
		Free_Mem();
		return false;
	}
	
	strcpy(path_name_Ex_1D,path);
    strcat(path_name_Ex_1D,"/data_Ex_1D/Ex_1D");
	
	strcpy(path_name_Ex,path);
    strcat(path_name_Ex,"/data_Ex/Ex");
	
	strcpy(path_name_Ey,path);
    strcat(path_name_Ey,"/data_Ey/Ey");
	
	strcpy(path_name_Ez,path);
    strcat(path_name_Ez,"/data_Ez/Ez");
	
	strcpy(path_name_Hx,path);
    strcat(path_name_Hx,"/data_Hx/Hx");
	
	strcpy(path_name_Hy,path);
    strcat(path_name_Hy,"/data_Hy/Hy");
	
	strcpy(path_name_Hz,path);
    strcat(path_name_Hz,"/data_Hz/Hz");
	
	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save slices from the field components
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_xyPBC_zPML::Save_FieldSlices(long iter)
{
	if (jel_plane_wave == 1)
	{
		if (save_1D(Ex_1D,n_z_a,n_1D,iter,path_name_Ex_1D))
			return -1;
	}

	if (save_3D_xy(Ex,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Ex))
		return -1;
	if (save_3D_xz(Ex,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Ex))
		return -1;
	if (save_3D_yz(Ex,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Ex))
		return -1;
	
	if (save_3D_xy(Ey,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Ey))
		return -1;
	if (save_3D_xz(Ey,n_x_a,n_x_b,ny_xz,n_z_a,n_z_b,iter,path_name_Ey))
		return -1;
	if (save_3D_yz(Ey,nx_yz,n_y_a,n_y_b,n_z_a,n_z_b,iter,path_name_Ey))
		return -1;

    if (save_3D_xy(Ez,n_x_a,n_x_b,n_y_a,n_y_b,nz_xy,iter,path_name_Ez))
		return -1;