fdtd_2d_te_pml_period.cpp

2DFDTD程序

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

#include "stdafx.h"
#include "fdtd_2D_TE_PML_period.h"
#include "Matrix.h"
#include "Math.h"
#include "FDTD_2D_TE_PERIODIC.h"
#include "FDTD_1D_HzEy.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;
		double mu_0  = 4*pi*1e-7;
		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-008;//2.5e-9;//6.25e-9;
		double dy = dx;
		
		int n_x = 1320; //total number of cells in x direction
		int n_y = 52;  //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 num_iter = 100000;

		//the excitation
		//Gaussian pulse
		double H0 = 1; //[A/m]
	    double tw = 20*dt;
		double t0 = 4*tw;
		//Sin
		double omega = 1.525e15;
		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_period.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:
			{//remove the comments to write out the Geometry matrix
					//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
		double ** Mat_Type = NULL;
		Mat_Type = Init_Matrix_2D<double>(n_Mat,2);
		if(!Mat_Type)
		{
			Index = Free_Matrix_2D<int>(Index);
			return FALSE;
		}

		Mat_Type[0][0] = 10.1089;   //eps_r 
		Mat_Type[0][1] = 1;  //mu_r 
	
		Mat_Type[1][0] = 1; //eps_r 
		Mat_Type[1][1] = 1;  //mu_r 

		
		FDTD_2D_TE_PERIODIC fdtd_TE;

		/////////////////////////////////////////////////////////////////////////////
		//Initialize fdtd_TE
		/////////////////////////////////////////////////////////////////////////////
		fdtd_TE.Init_Main_Param(Index, n_x, n_y, Mat_Type, n_Mat, n_PML, dt, dx, dy);
		
		int jel_plane_wave = 1;
		
		if (jel_plane_wave == 1)
		{
		//plane wave
			fdtd_TE.Init_TotFScatF(n_x_a, n_x_b, H0, t0, tw, omega, phi,
			                      source_Type);
		}
		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();

		/////////////////////////////////////////////////////////////////////////////
		//Initializations for data save
		/////////////////////////////////////////////////////////////////////////////
		int elment = 0;
		
		// 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 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;

		int n_Ind_F = 2;
		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 - n_x_a - 5;      Ind_F[1][1] = n_y/2;
		
	
		fdtd_TE.Init_Followed(Ind_F, n_Ind_F, num_iter);
			
		fdtd_TE.Get_Data_1D(hz_1D, ey_1D);
		fdtd_TE.Get_Data(hz, ex, ey);
  		
		int n1 = 0; // to save 1D
		n_x_a = 0; n_x_b = n_x; //save data

		/////////////////////////////////////////////////////////////////////////////
		//start of the main FDTD loop
		/////////////////////////////////////////////////////////////////////////////
		for (int n = 0; n < num_iter; n++)
		{
			if (jel_plane_wave == 1)
			{
				fdtd_TE.calc_Ey_1D();
				fdtd_TE.calc_Hz_1D(n*dt);
			}
			
			
			fdtd_TE.calc_Hz_TE();
			if (jel_plane_wave == 1)
				fdtd_TE.incident_Hz();
			else
				fdtd_TE.Point_Source(n*dt);

			fdtd_TE.periodic_Hz();

			fdtd_TE.calc_Ex_TE();
			fdtd_TE.calc_Ey_TE();
			if (jel_plane_wave == 1)
			{
				fdtd_TE.incident_Ex();
				fdtd_TE.incident_Ey();
			}
			
			fdtd_TE.periodic_Ex();
	
	        fdtd_TE.Set_Data_Followed(n);

			if (elment == 1 && n%10 == 1 && n<3000)
			{
				if (jel_plane_wave == 1)
					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,n1,n_y,n,path_name_Hz) )
					return FALSE;
				if ( !save_2D(ex,n_x_a,n_x_b,n1,n_y,n,path_name_Ex) )
					return FALSE;
				if ( !save_2D(ey,n_x_a,n_x_b,n1,n_y,n,path_name_Ey) )
					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);
		
		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;
		
		/////////////////////////////////////////////////////////////////////////////
		//Free allocated memory 
		/////////////////////////////////////////////////////////////////////////////
		Index = Free_Matrix_2D<int>(Index);
		Mat_Type = Free_Matrix_2D<double>(Mat_Type);

	}

	return nRetCode;
}