fdtd_2d_te_pml_lorentz.cpp

用C++编写的处理色散媒质的FDTD程序

// fdtd_2D_TE_PML_Lorentz.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "fdtd_2D_TE_PML_Lorentz.h"
#include "Matrix.h"
#include "Math.h"
#include "FDTD_2D_TE_LORENTZ.h"
#include "FDTD_1D_HzEy_LORENTZ.h"
#include "Save_File_Data.h"

#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif

/////////////////////////////////////////////////////////////////////////////
// The one and only application object

CWinApp theApp;

using namespace std;

int _tmain(int argc, TCHAR* argv[], TCHAR* envp[])
{
	int nRetCode = 0;

	// initialize MFC and print and error on failure
	if (!AfxWinInit(::GetModuleHandle(NULL), NULL, ::GetCommandLine(), 0))
	{
		// TODO: change error code to suit your needs
		cerr << _T("Fatal Error: MFC initialization failed") << endl;
		nRetCode = 1;
	}
	else
	{
		double pi = 3.14159265358979;
		double eps_0 = 8.854*1e-12; //[F/m]
		double mu_0  = 4*pi*1e-7; //[H/m]
		double c = 1/sqrt(eps_0*mu_0); //the speed of the light in the vacuum

		//the parameters of the photonic crystal
		double dx =  1.5e-8;
		double dy = dx;
		
		int n_x = 1320; //total number of cells in x direction
		int n_y = 610;  //total number of cells in y direction
		
		double dt;
		if (dx < dy)
		  dt = dx/c/2; //the time step, stability criteria dt<=dx/c!!!
		else
		  dt = dy/c/2; //the time step, stability criteria dt<=dx/c!!!
		
		int n_PML = 20; //the size of PML
		
		//Total field - Scattered field interface
		int n_x_a = n_PML + 10;
		int n_x_b = n_x - n_PML - 10;
		int n_y_a = n_PML + 10;
		int n_y_b = n_y - n_PML - 10;

		int num_iter = 100000;

		//the excitation
		//Gaussian pulse
		double H0 = 1; //[A/m]
	    double tw = 20*dt;
		double t0 = 4*tw;
		double omega = 1.3e15;
		double phi = -pi/2;
		int source_Type = 1; // 1- Gauss; 2- Sin; 3- Gauss-Sin

		int ** Index = NULL;
		Index = Init_Matrix_2D<int>(n_x,n_y);
		if(!Index)
		{
			return FALSE;
		}


		/////////////////////////////////////////////////////////////////////////////
		//Reads from file the elements of the Geometry matrix 
		/////////////////////////////////////////////////////////////////////////////
		CString path_Name = "Geom_40sz10_30_triang.geom";
		switch ( load_2D_int(Index,n_x,n_y,path_Name) )
		{
			case -1:
					printf("faild to open the file\n");
					return FALSE;
			case -2:
					printf("wrong file content\n");
					return FALSE;
			case -3:
					printf("wrong dimensions\n");
					return FALSE;
			default:
			{
					//CString name_Index = "Index";
					//int a1 = 0, b1 =0, c1 = 0;
					//if ( !save_2D_int( Index,a1,n_x,b1,n_y,c1,name_Index) )
					//			return FALSE;
			}
		}
		
		/////////////////////////////////////////////////////////////////////////////
		
		int n_Mat = 2; //number of materials in the geometry
		int *n_Mat_Lorentz = NULL; //the number of second orders Lorentz poles
		
		n_Mat_Lorentz = (int *)calloc(n_Mat, sizeof(int));
        if(!n_Mat_Lorentz)
		{
			Index = Free_Matrix_2D<int>(Index);
			return FALSE;
		}

		n_Mat_Lorentz[0] = 2 + 10*3 + 1; n_Mat_Lorentz[1] = 2 + 1*3 + 1;	//n_Mat_Lorentz[2] = 2 + 1*3 + 1;
		
		double ** Mat_Param = NULL;
		Mat_Param = Init_Matrix_2D<double>(n_Mat,n_Mat_Lorentz);
		if(!Mat_Param)
		{
			Index = Free_Matrix_2D<int>(Index);
			free(n_Mat_Lorentz);
			return FALSE;
		}
		
		//InP
		Mat_Param[0][0] = 10;//Number of Lorentz terms
		//eps_inf
		Mat_Param[0][1] = 2.9728;
		//am                             0 < delta_0 < 2                omega_0
		Mat_Param[0][2]  =  2.1809;     Mat_Param[0][3]  = 0.12385;    Mat_Param[0][4]  = 7.1897e+015;
		Mat_Param[0][5]  = -2.1426;	    Mat_Param[0][6]  = 0.3912;     Mat_Param[0][7]  = 1.9157e+015;
		Mat_Param[0][8]  =  0.87198;    Mat_Param[0][9]  = 0.33261;    Mat_Param[0][10] = 8.4884e+015;
		Mat_Param[0][11] =  1.0966;	    Mat_Param[0][12] = 0.11644;    Mat_Param[0][13] = 4.8676e+015;
		Mat_Param[0][14] = -2.9501;		Mat_Param[0][15] = 1.1994;     Mat_Param[0][16] = 1.6229e+015;
		Mat_Param[0][17] =  3.2874;     Mat_Param[0][18] = 0.43415;    Mat_Param[0][19] = 5.6233e+015;  
		Mat_Param[0][20] =  1.5489;	    Mat_Param[0][21] = 1.0763;     Mat_Param[0][22] = 1.9818e+015;
		Mat_Param[0][23] =  0.93079;    Mat_Param[0][24] = 0.37739;    Mat_Param[0][25] = 1.8472e+015;
		Mat_Param[0][26] =  0.61149;    Mat_Param[0][27] = 0.30277;    Mat_Param[0][28] = 1.9326e+015;
		Mat_Param[0][29] =  1.4274;     Mat_Param[0][30] = 1.178;      Mat_Param[0][31] = 1.4736e+015;
		//mu_r
		Mat_Param[0][32] = 1;

    	//Vacuum
		Mat_Param[1][0] = 1;//Number of Lorentz terms
		//eps_inf
		Mat_Param[1][1] = 1;
		//am                             0 < delta_0 < 2                omega_0
		Mat_Param[1][2] = 0;	         Mat_Param[1][3]  = 0;          Mat_Param[1][4]  = 0;
		//mu_r
		Mat_Param[1][5] = 1;
		
		
		
		//the averaged material coefficients needed for PML
		double av_eps_r = 10;
		double av_mu_r = 1;


		CFDTD_2D_TE_LORENTZ fdtd_TE;

		/////////////////////////////////////////////////////////////////////////////
		//Initialize fdtd_TE
		/////////////////////////////////////////////////////////////////////////////
		fdtd_TE.Init_Main_Param(Index, n_x, n_y, Mat_Param, n_Mat, n_PML, dt, dx, dy);
		
		int elment_W = 0; //switch to save the energy in all points
		if (elment_W)
			fdtd_TE.Init_W(n_x_a,n_x_b,n_y_a,n_y_b);
		
		int jel_plane_wave = 0;
		if (jel_plane_wave == 0)
		{
		//plane wave
		   fdtd_TE.Init_TotFScatF(n_x_a, n_x_b, n_y_a, n_y_b, H0, t0, tw, omega, phi,
			                      source_Type, av_eps_r, av_mu_r);
		}
		else
		{
		//point source
		   //fdtd_TE.Init_Gauss_Point_Source(n_x/2, n_y/2, H0, t0, tw);
		   fdtd_TE.Init_Sinus_Point_Source(n_PML+1, n_y/2, H0, omega, phi);
		}
		
		
		fdtd_TE.Set_PML_Param(av_eps_r, av_mu_r);

		/////////////////////////////////////////////////////////////////////////////
		//Initializations for data save
		/////////////////////////////////////////////////////////////////////////////
		int elment = 0; //switch to save the field components in all points
		
		// Create the directory to save output data
	    char path[512];
		GetCurrentDirectory(512, path); 
		CString Path;
		Path.Format("%s\\data",path);
		BOOL fret = CreateDirectory(Path,0);
		if (!fret)
		{
			DWORD gerr = GetLastError();
		}

		CString path_name_Hz;
			path_name_Hz.Format("%s\\Hz",Path);
		CString path_name_Ex;
			path_name_Ex.Format("%s\\Ex",Path);
		CString path_name_Ey;
			path_name_Ey.Format("%s\\Ey",Path);

		CString path_name_W;
			path_name_W.Format("%s\\W",Path);

		CString nev_Hz_1D;
			nev_Hz_1D.Format("%s\\Hz_1D",Path);
		
		double **hz = NULL, **ex = NULL, **ey = NULL;
		double *hz_1D = NULL, *ey_1D = NULL;

		double **w = NULL;
		int n_w_x_a = 0;
		int n_w_x_b = n_x_b - n_x_a;
		int n_w_y_a = 0;
		int n_w_y_b = n_y_b - n_y_a;

		int n_Ind_F = 4;
		int **Ind_F = NULL;
	    Ind_F = Init_Matrix_2D<int>(n_Ind_F,2);
		if(!Ind_F)
		{
			Index = Free_Matrix_2D<int>(Index);
			return FALSE;
		}
		
		Ind_F[0][0] = n_x_a + 5;            Ind_F[0][1] = n_y/2;
		Ind_F[1][0] = n_x_b - 5;            Ind_F[1][1] = n_y/2;
		Ind_F[2][0] = n_x/2;                Ind_F[2][1] = n_y_a + 5;
		Ind_F[3][0] = n_x/2;                Ind_F[3][1] = n_y_b - 5;
		
		//the field components
		fdtd_TE.Init_Followed(Ind_F, n_Ind_F, num_iter);
		
		int **Ind_F_W = NULL;
		CString name_Ind_F_W = "Ind_F_W";

		int n_Ind_F_W = (n_x_b - n_x_a - 1)/4;
		int elment_W_Foll = 1;
		if (elment_W_Foll ==1)
		{
			Ind_F_W = Init_Matrix_2D<int>(n_Ind_F_W,2);
			if(!Ind_F_W)
			{
				Index = Free_Matrix_2D<int>(Index);
				Ind_F = Free_Matrix_2D<int>(Ind_F);
				return FALSE;
			}
			
			/*j = 0;
			for (i = 0; i<n_Ind_F_W;  i++)
			{
				Ind_F_W[i][0] = n_x_a + j + 1;
				Ind_F_W[i][1] = n_y/2;
				j = j + 4;
			}
			*/
			Ind_F_W[0][0] = n_x_b - 5;
			Ind_F_W[0][1] = n_y/2;
			
			//int a1 = 0, b1 =0, b2 = 2, c1 = 0;
		    //if ( !save_2D_int( Ind_F_W,a1,n_Ind_F_W,b1,b2,c1,name_Ind_F_W) )
			//		return FALSE;
		}

	    //the energy	
		fdtd_TE.Init_Followed_W(Ind_F_W, n_Ind_F_W, num_iter);

		fdtd_TE.Get_Data_1D(hz_1D, ey_1D);
		fdtd_TE.Get_Data(hz, ex, ey);
		if (elment_W)
		     fdtd_TE.Get_W(w);
  		
		int n1 = 0; // to save 1D
		n_x_a = 0; n_y_a = 0; n_x_b = n_x; n_y_b = n_y; //to data save
				
		/////////////////////////////////////////////////////////////////////////////
		//start of the main FDTD loop
		/////////////////////////////////////////////////////////////////////////////
		for (int n = 0; n < num_iter; n++)
		{
			if (jel_plane_wave == 0)
			{
				fdtd_TE.calc_Ey_1D();
				fdtd_TE.calc_Hz_1D(n*dt);
			}
			
			
			fdtd_TE.calc_Hz_TE();
			if (jel_plane_wave)
				fdtd_TE.Point_Source(n*dt);

			fdtd_TE.calc_Ex_TE();
			fdtd_TE.calc_Ey_TE();
	
	        //the followed field components
			fdtd_TE.Set_Data_Followed(n);
			//the energy 
			if (elment_W && n%10 == 1 && n>30000)
			{
			    fdtd_TE.calc_W();
				if ( !save_2D(w,n_w_x_a,n_w_x_b,n_w_y_a,n_w_y_b,n,path_name_W) )
					return FALSE;
			}

			if (elment_W_Foll ==1)
			    fdtd_TE.Set_Data_Followed_W(n);

			if (elment == 1 && n%10 == 1 && n<30000)
			{
				
				/*if (jel_plane_wave == 0)
					if ( !save_1D(hz_1D, n1, n_x_b, n, nev_Hz_1D) )
						return FALSE;
				if ( !save_2D(hz,n_x_a,n_x_b,n_y_a,n_y_b,n,path_name_Hz) )
					return FALSE;
				if ( !save_2D(ex,n_x_a,n_x_b,n_y_a,n_y_b,n,path_name_Ex) )
					return FALSE;
				if ( !save_2D(ey,n_x_a,n_x_b,n_y_a,n_y_b,n,path_name_Ey) )
					return FALSE;*/
				
				if ( !save_2D_binary(hz,n_x,n_y,n,path_name_Hz) )
				{  
					printf( "Problem writing the file\n" ); 
					return FALSE; 
				}
			}

			if (n%100 == 1)
				cout<< n << endl;
		}

		/////////////////////////////////////////////////////////////////////////////
		//end of the main FDTD loop
		/////////////////////////////////////////////////////////////////////////////
		
		/////////////////////////////////////////////////////////////////////////////
		//Save the followed data
		/////////////////////////////////////////////////////////////////////////////
		double **hz_foll = NULL, **ex_foll = NULL, **ey_foll = NULL; 
		fdtd_TE.Get_Data_Followed(hz_foll, ex_foll, ey_foll);
		
		CString path_name_hz_foll;
		path_name_hz_foll.Format("%s\\Hz_foll",Path);
		CString path_name_ex_foll;
		path_name_ex_foll.Format("%s\\Ex_foll",Path);
		CString path_name_ey_foll;
		path_name_ey_foll.Format("%s\\Ey_foll",Path);
		
		double **w_foll;
		CString path_name_w_foll;
			
		int x1 =0, y1 = 0;
		if ( !save_2D(hz_foll,x1,n_Ind_F,y1,num_iter,n,path_name_hz_foll) )
			return FALSE;
		if ( !save_2D(ex_foll,x1,n_Ind_F,y1,num_iter,n,path_name_ex_foll) )
			return FALSE;
		if ( !save_2D(ey_foll,x1,n_Ind_F,y1,num_iter,n,path_name_ey_foll) )
			return FALSE;
		
		if (elment_W_Foll == 1)
		{
			fdtd_TE.Get_Data_Followed_W(w_foll);
			path_name_w_foll.Format("%s\\W_foll",Path);
			if ( !save_2D(w_foll,x1,n_Ind_F_W,y1,num_iter,n,path_name_w_foll) )
				return FALSE;
		}
		
		/////////////////////////////////////////////////////////////////////////////
		//Free the allocated memory 
		/////////////////////////////////////////////////////////////////////////////
		Index = Free_Matrix_2D<int>(Index);
		if (Ind_F)
			Ind_F = Free_Matrix_2D<int>(Ind_F);
		if (Ind_F_W)
			Ind_F_W = Free_Matrix_2D<int>(Ind_F_W);
		free(n_Mat_Lorentz);
		Mat_Param = Free_Matrix_2D<double>(Mat_Param);		
	}

	return nRetCode;
}