fdtd_3d_pml.h

采用UMPL边界的FDTD编程

	long nx_MIN_1_MIN_nPML_x_2, ny_MIN_1_MIN_nPML_y_2, nz_MIN_1_MIN_nPML_z_2;

	long nPML_x_2_MIN_1, nPML_y_2_MIN_1, nPML_z_2_MIN_1;

	long nxDIV2, nyDIV2, nzDIV2;


	//the coordinates of the followed field components
	long **Ind_Foll; 
	 //number of the followed points and the number of iterations
	long length_Ind_Foll, n_tot;
	//stores the followed field components
	double **Hx_Foll, **Hy_Foll, **Hz_Foll, **Ex_Foll, **Ey_Foll, **Ez_Foll;
	
	//to save the field components
	long n_x_a, n_x_b, n_y_a, n_y_b, n_z_a, n_z_b; //identify the saved volume
	long n_x_b1, n_y_b1, n_z_b1;
	long nx_yz, ny_xz, nz_xy;
	char *path_name_Ex;
	char *path_name_Ey;
	char *path_name_Ez;
	char *path_name_Hx;
	char *path_name_Hy;
	char *path_name_Hz;
	char *path_name_Ez_1D;
	int jel_Save_Slice;
	long nr_save;

public:
	CFDTD_3D_PML(void);
	~CFDTD_3D_PML(void);

	int Init(long ***&Index, long n_x, long n_y, long n_z, double **&Mater, long nMat, 
		     long n_pml_x_1, long n_pml_x_2, long n_pml_y_1, long n_pml_y_2, long n_pml_z_1,
		     long n_pml_z_2, double d_t, double d_x, double d_y, double d_z,
		     double PML_eps_r_x_1, double PML_mu_r_x_1, double PML_eps_r_x_2, 
		     double PML_mu_r_x_2, double PML_eps_r_y_1, double PML_mu_r_y_1, 
		     double PML_eps_r_y_2, double PML_mu_r_y_2, double PML_eps_r_z_1, 
		     double PML_mu_r_z_1, double PML_eps_r_z_2, double PML_mu_r_z_2);
	//Index - describe the geometry, each element indicating the material type in 
	//        the computational grid 
	//nx, ny, nz - the dimensions of the computational space
	//Mater - contains the permittivity and conductivity. Dimension [n_Mat][2] 
	//n_PML - indicate the dimension of the PML
	
	void Init_PML_Par_x(double *K_Ey_a, double *K_Ey_b, double *K_Gz_a, double *K_Gz_b, 
						double *K_Hx_c, double *K_Hx_d, double *K_Ex_c, double *K_Ex_d, 
						double *K_Hy_a, double *K_Hy_b, double *K_Bz_a, double *K_Bz_b,
						double eps_r_x, double mu_r_x, long  nPML_x, double dx, double dt);

	void Init_PML_Par_y(double *K_Gx_a, double *K_Gx_b, double *K_Ez_a, double *K_Ez_b, 
						double *K_Hy_c, double *K_Hy_d, double *K_Ey_c, double *K_Ey_d,
						double *K_Bx_a, double *K_Bx_b, double *K_Hz_a, double *K_Hz_b,
						double eps_r_y, double mu_r_y, long  nPML_y, double dy, double dt);

	void Init_PML_Par_z(double *K_Ex_a, double *K_Ex_b, double *K_Gy_a, double *K_Gy_b, 
						double *K_Hz_c, double *K_Hz_d, double *K_Ez_c, double *K_Ez_d,
						double *K_Hx_a, double *K_Hx_b, double *K_By_a, double *K_By_b,
						double eps_r_z, double mu_r_z, long  nPML_z, double dz, double dt);


	void Init_ptSource(long **Coord, long nCoord, double **Par_Excit, int s_type,
					   int Pt_Source_Ex, int Pt_Source_Ey, int Pt_Source_Ez,
	                   int Pt_Source_Hx, int Pt_Source_Hy, int Pt_Source_Hz);

	int Init_TS(double Teta, double Phi, double Gamma, long N_TS,
				long nPML_1D, int source_type, double X0, double t0,
				double tw, double omega, double phase, double PML_eps_r_1D,
				double PML_mu_r_1D, int *Jel_TS_planes);

	void Init_bc_Symmetric(int jel_Sym_ia, int jel_Sym_ib, int jel_Sym_ja, 
		                   int jel_Sym_jb, int jel_Sym_ka, int jel_Sym_kb);

	void Init_bc_AntiSymmetric(int jel_A_Sym_ia, int jel_A_Sym_ib, int jel_A_Sym_ja, 
		                       int jel_A_Sym_jb, int jel_A_Sym_ka, int jel_A_Sym_kb);
	
	//collect the functions to compute the electric and magnetic field components
	void calc_E(double time);
	void calc_H(double time, long iter);
		
	void Get_Electric_Fields(double ***&E_x, double ***&E_y, double ***&E_z);

	//functions to save output data
	int Init_Followed(long **Ind_Followed, long n, long n_t);
	int 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 nr_Save,  long SaveFROMinst,
							  long SaveTOinst, char *path);
	
	void Get_Data_Followed(double **&h_x, double **&h_y, double **&h_z, 
						   double **&e_x, double **&e_y, double **&e_z);

	void Reset_Field_Comp();

	//free up the allocated memory	
	void Free_Mem();

	//Save intermediate data in binary file to be able to continue the computations
	#ifndef run_enviroment //IBM AIX
		int save_FGE_BH(long iter, double time, char *Path);
		int load_FGE_BH(char *Path);
	#endif

private:
	void Init_Alloc_FxGxEx();
	void Init_Alloc_FyGyEy();
	void Init_Alloc_FzGzEz();
	void Init_Alloc_BxHx();
	void Init_Alloc_ByHy();
	void Init_Alloc_BzHz();

	//the field components from Maxwell equations
	void calc_Ex_1();   void calc_Ex_2();   void calc_Ex_3();
	void calc_Ex_4();   void calc_Ex_5();   void calc_Ex_6();
	void calc_Ex_7();   void calc_Ex_8();   void calc_Ex_9();
	void calc_Ex_10();  void calc_Ex_11();  void calc_Ex_12();
	void calc_Ex_13();  
	void calc_Ex_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b); 
	void calc_Ex_15();
	void calc_Ex_16();  void calc_Ex_17();  void calc_Ex_18();
	void calc_Ex_19();  void calc_Ex_20();  void calc_Ex_21();
	void calc_Ex_22();  void calc_Ex_23();  void calc_Ex_24();
	void calc_Ex_25();  void calc_Ex_26();  void calc_Ex_27();

	void calc_Ey_1();   void calc_Ey_2();   void calc_Ey_3();
	void calc_Ey_4();   void calc_Ey_5();   void calc_Ey_6();
	void calc_Ey_7();   void calc_Ey_8();   void calc_Ey_9();
	void calc_Ey_10();  void calc_Ey_11();  void calc_Ey_12();
	void calc_Ey_13();  
	void calc_Ey_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b);  
	void calc_Ey_15();
	void calc_Ey_16();  void calc_Ey_17();  void calc_Ey_18();
	void calc_Ey_19();  void calc_Ey_20();  void calc_Ey_21();
	void calc_Ey_22();  void calc_Ey_23();  void calc_Ey_24();
	void calc_Ey_25();  void calc_Ey_26();  void calc_Ey_27();

	void calc_Ez_1();   void calc_Ez_2();   void calc_Ez_3();
	void calc_Ez_4();   void calc_Ez_5();   void calc_Ez_6();
	void calc_Ez_7();   void calc_Ez_8();   void calc_Ez_9();
	void calc_Ez_10();  void calc_Ez_11();  void calc_Ez_12();
	void calc_Ez_13(); 
	void calc_Ez_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b);
	void calc_Ez_15();
	void calc_Ez_16();  void calc_Ez_17();  void calc_Ez_18();
	void calc_Ez_19();  void calc_Ez_20();  void calc_Ez_21();
	void calc_Ez_22();  void calc_Ez_23();  void calc_Ez_24();
	void calc_Ez_25();  void calc_Ez_26();  void calc_Ez_27();
	
    void calc_Hx_1();   void calc_Hx_2();   void calc_Hx_3();
	void calc_Hx_4();   void calc_Hx_5();   void calc_Hx_6();
	void calc_Hx_7();   void calc_Hx_8();   void calc_Hx_9();
	void calc_Hx_10();  void calc_Hx_11();  void calc_Hx_12();
	void calc_Hx_13();  
	void calc_Hx_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b);
	void calc_Hx_15();
	void calc_Hx_16();  void calc_Hx_17();  void calc_Hx_18();
	void calc_Hx_19();  void calc_Hx_20();  void calc_Hx_21();
	void calc_Hx_22();  void calc_Hx_23();  void calc_Hx_24();
	void calc_Hx_25();  void calc_Hx_26();  void calc_Hx_27();

	void calc_Hy_1();   void calc_Hy_2();   void calc_Hy_3();
	void calc_Hy_4();   void calc_Hy_5();   void calc_Hy_6();
	void calc_Hy_7();   void calc_Hy_8();   void calc_Hy_9();
	void calc_Hy_10();  void calc_Hy_11();  void calc_Hy_12();
	void calc_Hy_13();  
	void calc_Hy_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b);
	void calc_Hy_15();
	void calc_Hy_16();  void calc_Hy_17();  void calc_Hy_18();
	void calc_Hy_19();  void calc_Hy_20();  void calc_Hy_21();
	void calc_Hy_22();  void calc_Hy_23();  void calc_Hy_24();
	void calc_Hy_25();  void calc_Hy_26();  void calc_Hy_27();
	
	void calc_Hz_1();   void calc_Hz_2();   void calc_Hz_3();
	void calc_Hz_4();   void calc_Hz_5();   void calc_Hz_6();
	void calc_Hz_7();   void calc_Hz_8();   void calc_Hz_9();
	void calc_Hz_10();  void calc_Hz_11();  void calc_Hz_12();
	void calc_Hz_13();  
	void calc_Hz_14(long nx_a, long nx_b, long ny_a, long ny_b, long nz_a, long nz_b);  
	void calc_Hz_15();
	void calc_Hz_16();  void calc_Hz_17();  void calc_Hz_18();
	void calc_Hz_19();  void calc_Hz_20();  void calc_Hz_21();
	void calc_Hz_22();  void calc_Hz_23();  void calc_Hz_24();
	void calc_Hz_25();  void calc_Hz_26();  void calc_Hz_27();

	//Symmetric boundary conditions
	void bc_Symmetric_iDir_a();
	void bc_Symmetric_iDir_b();

	void bc_Symmetric_jDir_a();
	void bc_Symmetric_jDir_b();

	void bc_Symmetric_kDir_a();
	void bc_Symmetric_kDir_b();

	//Anti-symmetric boundary conditions
	void bc_AntiSymmetric_iDir_a();
	void bc_AntiSymmetric_iDir_b();

	void bc_AntiSymmetric_jDir_a();
	void bc_AntiSymmetric_jDir_b();

	void bc_AntiSymmetric_kDir_a();
	void bc_AntiSymmetric_kDir_b();

	//Functions referring to the excitation
	
	//Point Source
	void ptSource_J(double ***&xx, double **&Par_ptS, double time);
	
	//TF_SF Ex
	void TS_Ex_j0();
	void TS_Ex_j1();
	void TS_Ex_k0();
	void TS_Ex_k1();

	//TF_SF Ey
	void TS_Ey_i0(); 
	void TS_Ey_i1(); 
	void TS_Ey_k0(); 
	void TS_Ey_k1();

	//TF_SF Ez
	void TS_Ez_i0(); 
	void TS_Ez_i1(); 
	void TS_Ez_j0(); 
	void TS_Ez_j1();

	//TF_SF Hx
	void TS_Hx_j0(); 
	void TS_Hx_j1(); 
	void TS_Hx_k0(); 
	void TS_Hx_k1();

	//TF_SF Hy
	void TS_Hy_i0();
	void TS_Hy_i1(); 
	void TS_Hy_k0(); 
	void TS_Hy_k1();

	//TF_SF Hz
	void TS_Hz_i0();
	void TS_Hz_i1(); 
	void TS_Hz_j0(); 
	void TS_Hz_j1();

	int Interp_1D(double *x, int length_x, double *f, double *x_i, int length_x_i, double *g);
	int Interp_H_TS();
	int Interp_E_TS();

	//Compute the phase velocity
	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);

	//the coordinates of the followed field components
	void Set_Data_Followed(long n_t);
	int Save_FieldSlices(long iter);

	void ErrorMessage(int er, int myrank, char* text);
};