fdtd_1d_hzey.cpp

2DFDTD程序

// FDTD_1D_HzEy.cpp: implementation of the FDTD_1D_HzEy class.
//
//////////////////////////////////////////////////////////////////////

#include "stdafx.h"
#include "fdtd_2D_TE_PML_period.h"
#include "FDTD_1D_HzEy.h"
#include "Math.h"

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
FDTD_1D_HzEy::FDTD_1D_HzEy()
{
	Ey_1D = NULL; Hz_1D = NULL;
	K_E1_a = NULL; K_E1_b = NULL; K_E2_a = NULL; K_E2_b = NULL;
	pi = 3.1415926535897932384626433832795;
	eps_0 = 8.854e-12; // [F/m]
	mu_0 = 4*pi*1e-7; // [H/m]
}

FDTD_1D_HzEy::~FDTD_1D_HzEy()
{
	if(Hz_1D)
	    free(Hz_1D); 	
	if(Ey_1D)
	    free(Ey_1D); 
	if(K_E1_a)
		free(K_E1_a);
	if(K_E1_b)
		free(K_E1_b);
	if(K_E2_a)
		free(K_E2_a);
	if(K_E2_b)
		free(K_E2_b);
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
BOOL FDTD_1D_HzEy::Init_Main_Param_1D(int n_x, int n_pml, double Eps_r, double Mu_r, double d_t, 
									  double d_x)
{
	//dimension of the computational space
    nx = n_x;
	eps_r = Eps_r;
	mu_r = Mu_r;
	dt = d_t;
	dx = d_x;
	n_PML = n_pml;
	
	Hz_1D = (double *) calloc(nx,sizeof(double));
	if(!Hz_1D)
	{
		return FALSE;
	}
	
	Ey_1D = (double *) calloc(nx-1,sizeof(double));
	if(!Ey_1D)
	{
		free(Hz_1D);
		return FALSE;
	}
	
	K_E1_a = (double *) calloc(nx,sizeof(double));
	if(!K_E1_a)
	{
		free(Hz_1D);
		free(Ey_1D);
		return FALSE;
	}
	
	K_E1_b = (double *) calloc(nx,sizeof(double));
	if(!K_E1_b)
	{
		free(Hz_1D);
		free(Ey_1D);
		free(K_E1_a);
		return FALSE;
	}
	
	K_E2_a = (double *) calloc(nx,sizeof(double));
	if(!K_E2_a)
	{
		free(Hz_1D);
		free(Ey_1D);
		free(K_E1_a);
		free(K_E1_b);
		return FALSE;
	}
	
	K_E2_b = (double *) calloc(nx,sizeof(double));
	if(!K_E2_b)
	{
		free(Hz_1D);
		free(Ey_1D);
		free(K_E1_a);
		free(K_E1_b);
		free(K_E2_a);
		return FALSE;
	}
	
	return TRUE;

}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Gaussian pulse
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::Init_Gauss_1D(double H_0, double t_0, double t_w)
{
	source_type = 1;
	H0 = H_0;
	t0 = t_0;
	tw = t_w;
}

//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Sinusoidal plane wave
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::Init_Sinus_1D(double H_0, double om, double Phi)
{
	source_type = 2;
	H0 = H_0;
	omega = om;
	phi = Phi;
}


//////////////////////////////////////////////////////////////////////
//Initialize the parameters of a Wave Packet (sinusoidal modulated Gauss pulse)
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::Init_Wave_Packet_1D(double H_0, double t_0, double t_w, 
									   double om, double Phi)
{
	source_type = 3;
	H0 = H_0;
	t0 = t_0;
	tw = t_w;
	omega = om;
	phi = Phi;
}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::Set_PML_Param_1D()
{
	for (i = 0; i<nx; i++)
	{
		K_E1_a[i] = 1.0;
		K_E1_b[i] = dt/(mu_0*mu_r*dx);

		if (i < nx-1)
		{
			K_E2_a[i] = 1.0;
			K_E2_b[i] = dt/(eps_0*eps_r*dx);
		}			
	}

	//PML parameters
	double ka_max = 1;
	int exponent = 4;
	double R_err = 1e-16;
	
	double eta = sqrt(mu_0*mu_r/eps_0/eps_r);

	double sigma_x, sigma_max, ka_x;
		
	sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx);
	//sigma_max= -(exponent+1)*log(R_err)*eps_0/(2*n_PML*dx*sqrt(mu_0*mu_r*eps_0*eps_r));

	for (i = 0; i<n_PML; i++)
	{
		sigma_x         = sigma_max*pow( (n_PML - i)/((double) n_PML) ,exponent);
		ka_x            = 1 + (ka_max - 1)*pow( (n_PML-i)/((double) n_PML) ,exponent);
		K_E1_a[nx-i-1]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		//K_E1_a[i]       = K_E1_a[nx-i-1];
		K_E1_b[nx-i-1]  = 2*eps_0*dt/(2*eps_0*ka_x+sigma_x*dt)/(mu_0*mu_r*dx);
		//K_E1_b[i]       = K_E1_b[nx-i-1];
		
		sigma_x         = sigma_max*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		ka_x            = 1 + (ka_max - 1)*pow( (n_PML - i - 0.5)/n_PML ,exponent);
		K_E2_a[nx-i-2]  = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		//K_E2_a[i]       = K_E2_a[nx-i-2];
		K_E2_b[nx-i-2]  = 2*dt/(2*eps_0*ka_x + sigma_x*dt)/(eps_r*dx);
		//K_E2_b[i]       = K_E2_b[nx-i-2];
	}
}

//////////////////////////////////////////////////////////////////////
//Compute the Hz component
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::calc_Hz_1D(double t)
{
	switch (source_type)
	{
		case 1: //Gaussian pulse
			Hz_1D[0] = H0*exp( -pow( (t-t0)/tw ,2) ); 
			break;
		case 2: //Sinusoidal plane wave
			Hz_1D[0] = H0*cos( omega*t + phi);
			break;
		case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
			Hz_1D[0] = H0*cos( omega*t + phi)*exp( -pow( (t-t0)/tw ,2) );
	}
	 
	for (i = 1; i < nx - 1; i++)
	{
		Hz_1D[i] = K_E1_a[i]*Hz_1D[i] - K_E1_b[i]*(Ey_1D[i] - Ey_1D[i-1]);
	}

}

//////////////////////////////////////////////////////////////////////
//Compute the Ey component
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::calc_Ey_1D()
{
	for (i = 0; i < nx - 1; i++)
	{
		Ey_1D[i] = K_E2_a[i]*Ey_1D[i] - K_E2_b[i]*(Hz_1D[i+1] - Hz_1D[i]);
	}
}


//////////////////////////////////////////////////////////////////////
//Access Hz_1D and Ey_1D
//////////////////////////////////////////////////////////////////////
void FDTD_1D_HzEy::Get_Data_1D(double *&X, double *&Y)
{
	X = Hz_1D;
	Y = Ey_1D;
}