fdtd_1d.cpp

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

//////////////////////////////////////////////////////////////////////
// fdtd_1D.cpp: implementation of the fdtd_1D class.
//////////////////////////////////////////////////////////////////////
#include <math.h>
#include "fdtd_1D.h"
#include "Load_Save_File_Data.h"
#include <sys/stat.h> //for mkdir()
#include <cstring>

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
CFDTD_1D_EzHx::CFDTD_1D_EzHx()
{
	Ez_1D = NULL;
	Hx_1D = NULL;
	pi = 3.14159265358979;
	eps_0 = 8.854e-12;
	mu_0 = 4*pi*1e-7;
}

CFDTD_1D_EzHx::~CFDTD_1D_EzHx()
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters in case of sinusoidal excitation
///////////////////////////////////////////////////////////////////////////////////////
bool CFDTD_1D_EzHx::Init(long  n_y, double dt, double dy, double eps_r, double mu_r, double sigma)
{
	ny = n_y;
	nyMIN1 = ny - 1;
	nyMIN2 = ny - 2;
	
	Ez_1D = (double *) calloc(ny,sizeof(double));
	if (!Ez_1D)
		return false;
	
	Hx_1D = (double *) calloc(ny,sizeof(double));
	if (!Hx_1D)
	{
		free(Ez_1D);
		Ez_1D = NULL;
		return false;
	}

	// MUR boundary conditions and coefficients for calc_Ez_1D() and calc_Hx_1D() 
	double speed_c = 1/sqrt(eps_0*eps_r*mu_0*mu_r);

	K_M_1 = (speed_c*dt-dy)/(speed_c*dt+dy); 
	K_M_2 = 2*dy/(speed_c*dt+dy);
	
	Ca = (1-sigma*dt/2/eps_0/eps_r)/(1+sigma*dt/2/eps_0/eps_r);
	Cb = (dt/eps_0/eps_r/dy)/(1+sigma*dt/2/eps_0/eps_r);
	Db = dt/(mu_r*mu_0*dy);

	long  i;
	for (i = 0; i<2; i++)
	{	  
		Ez_n_1[i] = 0.0;
		Ez_n[i]   = 0.0;
	}

	return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Init_Sin_Source(double E_0, double om, double Phi)
{
	E0 = E_0;
	omega = om;
	phi = Phi;

	alfa = 3*omega/(2*pi);

	Source_Type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize  a Gaussian source
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Init_Gauss_Source(double E_0, double t_0, double t_w)
{
	E0 = E_0;
	t0 = t_0;
	tw = t_w;
	twORtw = tw*tw;
    	
	Source_Type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal modulated Gaussian pulse
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Init_GaussSin(double E_0, double om, double Phi, double t_0, double t_w)
{
	E0 = E_0;
	omega = om;
	phi = Phi;
	t0 = t_0;
	tw = t_w;
	twORtw = tw*tw;

	Source_Type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field 
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::calc_Ez_1D(double time)
{
	long  i, j;

	for (i = 0; i < 2; i++)
	{
	  Ez_n_1[i] = Ez_n[i];
	}

    Ez_n[0] = Ez_1D[nyMIN2];
	Ez_n[1] = Ez_1D[nyMIN1];
    
	for (j = 1; j < ny; j++) 
	{ 
		Ez_1D[j] = Ca*Ez_1D[j] + Cb*(Hx_1D[j-1]-Hx_1D[j]);
	}

	// MUR Absorbing boundary conditions
	//at z=h
	Ez_1D[nyMIN1] = -Ez_n_1[0] + K_M_1*(Ez_1D[nyMIN2]+Ez_n_1[1]) + K_M_2*(Ez_n[1]+Ez_n[0]);

	//Calculate the incident Ez_1D field
	switch (Source_Type)
	{
		case 1:
			Ez_1D[0] = (1-exp(-alfa*time))*E0*cos(omega*time + phi);
			break;
		case 2:
			Ez_1D[0] = E0*exp( - (time - t0)*(time - t0)/(twORtw) );
			break;
		case 3:
			Ez_1D[0] = E0*cos(omega*time + phi)*exp( - (time - t0)*(time - t0)/(twORtw) );
			break;
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hx field 
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::calc_Hx_1D()
{
	long  j;
	for (j = 0; j < ny-1; j++ ) 
	{ 
		Hx_1D[j] = Hx_1D[j] + Db*(Ez_1D[j]-Ez_1D[j+1]);
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the address of Ex and Hy
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Get_Data(double *&X, double *&Y)
{
	X = Ez_1D;
	Y = Hx_1D;
}

///////////////////////////////////////////////////////////////////////////////////////
//Reinitialize with zero all the field components
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Reset_Field_Comp(void)
{
	long  j;
	for (j = 0; j < ny; j++)
	{
		Ez_1D[j] = 0.0;
		Hx_1D[j] = 0.0;
	}

	Ez_n_1[0] = 0.0; 
	Ez_n_1[1] = 0.0;
	Ez_n[0]   = 0.0;
	Ez_n[1]   = 0.0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Save out the workspace data - the possibility to continue computations
///////////////////////////////////////////////////////////////////////////////////////
int CFDTD_1D_EzHx::Save_FDTD_1D_Workspace(char *path)
{	
	char *path_name = NULL, *file_name = NULL;
	path_name =(char *) calloc(256,sizeof(char));
	file_name =(char *) calloc(256,sizeof(char));
	if (!path_name && !file_name)
	{
		if (path_name)
			free(path_name);
		if (file_name)
			free(file_name);
		return 1;
	}
	strcpy(path_name,path);
	strcat(path_name,"/data_fdtd_1D");

	//make directory to save F G E
	mode_t Mode = S_IRWXU;
	mkdir(path_name,Mode);//for UNIX-AIX

	//save Ez_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/Ez_1D");
	if (save_1D_binary(Ez_1D, ny, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save Hx_1D
	strcpy(file_name,path_name);
	strcat(file_name,"/Hx_1D");
	if(save_1D_binary(Hx_1D, ny, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	//save Ez_n
	strcpy(file_name,path_name);
	strcat(file_name,"/Ez_n_1D");
	if(save_1D_binary(Ez_n, 2, 0, file_name))
	{
		free(path_name);
		free(file_name);
		return 2; //faild to save data
	}

	free(path_name);
	free(file_name);

	return 0;
}

int CFDTD_1D_EzHx::Load_FDTD_1D_Workspace(char *path)
{
	char *file_name = NULL;
	file_name =(char *) calloc(256,sizeof(char));
	if (!file_name)
	{
		return 1;
	}
	
	//load Ez_1D
	strcpy(file_name,path);
	strcat(file_name,"/Ez_1D_0.dat");
	if (load_1D_binary(Ez_1D, ny, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Hx_1D
	strcpy(file_name,path);
	strcat(file_name,"/Hx_1D_0.dat");
	if (load_1D_binary(Hx_1D, ny, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	//load Ez_n
	strcpy(file_name,path);
	strcat(file_name,"/Ez_n_1D_0.dat");
	if (load_1D_binary(Ez_n, 2, file_name))
	{
		free(file_name);
		return 2; //faild to load data
	}

	free(file_name);

	return 0;
}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void CFDTD_1D_EzHx::Free_Mem()
{
	if (Ez_1D)
	{
		free(Ez_1D);
		Ez_1D = NULL;
	}

	if (Hx_1D)
	{
		free(Hx_1D);
		Hx_1D = NULL;
	}
}