fdtd_1d_ezhy.cpp

2维时域有限差分模拟

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

#include "stdafx.h"
#include "Fdtd_1D_EzHy.h"

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////

CFdtd_1D_EzHy::CFdtd_1D_EzHy()
{
	Ez = NULL;
	Hy = NULL;

	pi = 3.1415926535897932384626433832795;
	eps_0 = 8.854e-12; // [F/m]
	mu_0 = 4*pi*1e-7; // [H/m]
}

CFdtd_1D_EzHy::~CFdtd_1D_EzHy()
{
	if (Ez)
	{
		free(Ez);
		Ez = NULL;
	}

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

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the class parameters
///////////////////////////////////////////////////////////////////////////////////////
bool CFdtd_1D_EzHy::Init(int n_x, double dt, double dx, double eps_r, double mu_r)
{
	nx = n_x;

	Ez = (double *)calloc(nx,sizeof(double));
	if (!Ez)
	{
		return false;
	}

	Hy = (double *)calloc(nx-1,sizeof(double));
	if (!Hy)
	{
		free(Ez);
		Ez = NULL;
		return false;
	}

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

	K_E = dt/(eps_0*eps_r*dx);
	K_H = dt/(mu_0*mu_r*dx);

	double speed = 1/sqrt(eps_0*eps_r*mu_0*mu_r); //the speed of the light in the vacuum

	K_MUR_1 = ( speed*dt - dx )/( speed*dt + dx );
	K_MUR_2 =  2*dx/( speed*dt + dx );

	Ez_nxm2_n = 0.0;

    return true;
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal source
///////////////////////////////////////////////////////////////////////////////////////
void CFdtd_1D_EzHy::Init_Sin_Source(double E_0, double om, double Phi)
{
	double pi = 3.14159265358979;

	E0 = E_0;
	omega = om;
	phi = Phi;

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

	Source_Type = 1;
}

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

///////////////////////////////////////////////////////////////////////////////////////
//Initialize a sinusoidal modulated Gaussian pulse
///////////////////////////////////////////////////////////////////////////////////////
void CFdtd_1D_EzHy::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;

	Source_Type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ez field
///////////////////////////////////////////////////////////////////////////////////////
void CFdtd_1D_EzHy::calc_Ez(double time)
{
	int i;

	switch (Source_Type)
	{
		case 1:
			Ez[0] = (1-exp(-alfa*time))*E0*cos(omega*time + phi);
			break;
		case 2:
			Ez[0] = E0*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
		case 3:
			Ez[0] = E0*cos(omega*(time-t0) + phi)*exp( - (time - t0)*(time - t0)/(tw*tw) );
			break;
	}
	
	Ez_n_1[0] = Ez_n[0]; 
	Ez_n_1[1] = Ez_n[1];
	Ez_n[0]   = Ez[nx-1]; 
	Ez_n[1]   = Ez[nx-2];

	for (i = 1; i < nx-1; i++)
	{
		Ez[i] += K_E*( Hy[i] - Hy[i-1] );
	}

	//Ez[nx-1] = Ez_nxm2_n + K_MUR_1*( Ez[nx-2] - Ez[nx-1] );
	//Ez_nxm2_n = Ez[nx-2];
	Ez[nx-1] = -Ez_n_1[1] + K_MUR_1*(Ez_n_1[0] + Ez[nx-2]) + K_MUR_2*(Ez_n[0] + Ez_n[1]);
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void CFdtd_1D_EzHy::calc_Hy()
{
	int i;

	for (i = 0; i < nx-1; i++)
	{
		Hy[i] += K_H*( Ez[i+1] - Ez[i] );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Get the address of Ez and Hy
///////////////////////////////////////////////////////////////////////////////////////
void CFdtd_1D_EzHy::Get_Data(double *&ez, double *&hy)
{
	ez = Ez;
	hy = Hy;
}