fdtd_2d_te.cpp

fdtd的2d算法

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

#include "stdafx.h"
#include "fdtd_2D_TE_PML_a.h"
#include "FDTD_2D_TE.h"
#include "Matrix.h"
#include "Math.h"
#include <stdlib.h>         /* For _MAX_PATH definition */
#include <malloc.h>

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

//////////////////////////////////////////////////////////////////////
// Construction/Destruction
//////////////////////////////////////////////////////////////////////
FDTD_2D_TE::FDTD_2D_TE()
{
	Fz = NULL; Hz = NULL; Ex = NULL; Gx = NULL; Ey = NULL; Gy = NULL;  

	K_E1_a = NULL; K_E1_b = NULL; K_E2_a = NULL; K_E2_b = NULL;
	K_E3_a = NULL; K_E3_b = NULL; K_E4_a = NULL; K_E4_b = NULL;
	K_E5_a = NULL; K_E5_b = NULL; K_E6_a = NULL; K_E6_b = NULL;

	Hz_Foll = NULL; Ex_Foll = NULL; Ey_Foll = NULL;

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

}

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

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
BOOL FDTD_2D_TE::Init_Main_Param(int **&index, int n_x, int n_y, double **&mater, int n_mat, 
	  						     int n_pml, double d_t, double d_x, double d_y)
{
	Index = index; //contains the indices of the material type. Dimension [nx][ny]

	//dimensions of the computational space
    nx = n_x;
	ny = n_y; 
	
	//contains the material parameters [eps_r mu_r]
	Mater = mater;
	n_Mat = n_mat;

	//dimension of the PML region
	n_PML = n_pml;
	
	dx = d_x;
	dy = d_y;
	dt = d_t; //the time step

	Fz = Init_Matrix_2D<double>(nx,ny);
	if(!Fz)	
	{
		Free_Mem();
		return FALSE;
	}

	Hz = Init_Matrix_2D<double>(nx,ny);
	if(!Hz)
	{
		Free_Mem();
		return FALSE;
	}

	Gx = Init_Matrix_2D<double>(nx,ny-1);
	if(!Gx)
	{
		Free_Mem();
		return FALSE;
	}

	Ex = Init_Matrix_2D<double>(nx,ny-1);
	if(!Ex)
	{
		Free_Mem();
		return FALSE;
	}

	Gy = Init_Matrix_2D<double>(nx-1,ny);
	if(!Gy)
	{
		Free_Mem();
		return FALSE;
	}

	Ey = Init_Matrix_2D<double>(nx-1,ny);
	if(!Ey)
	{
		Free_Mem();
		return FALSE;
	}

	K_E1_a = (double *) calloc(ny,sizeof(double));
	if(!K_E1_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E1_b = (double *) calloc(ny,sizeof(double));
	if(!K_E1_b)
	{
		Free_Mem();
		return FALSE;
	}

	K_E2_a = (double *) calloc(nx,sizeof(double));
	if(!K_E2_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E2_b = (double *) calloc(nx,sizeof(double));
	if(!K_E2_b)
	{
		Free_Mem();
		return FALSE;
	}

	K_E3_a = (double *) calloc(ny-1,sizeof(double));
	if(!K_E3_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E3_b = (double *) calloc(ny-1,sizeof(double));
	if(!K_E3_b)
	{
		Free_Mem();
		return FALSE;
	}

	K_E4_a = (double *) calloc(nx,sizeof(double));
	if(!K_E4_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E4_b = (double *) calloc(nx,sizeof(double));
	if(!K_E4_b)
	{
		Free_Mem();
		return FALSE;
	}

	K_E5_a = (double *) calloc(nx-1,sizeof(double));
	if(!K_E5_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E5_b = (double *) calloc(nx-1,sizeof(double));
	if(!K_E5_b)
	{
		Free_Mem();
		return FALSE;
	}

	K_E6_a = (double *) calloc(ny,sizeof(double));
	if(!K_E6_a)
	{
		Free_Mem();
		return FALSE;
	}

	K_E6_b = (double *) calloc(ny,sizeof(double));
	if(!K_E6_b)
	{
		Free_Mem();
		return FALSE;
	}

	return TRUE;

}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::Set_PML_Param()
{
	for (i = 0; i<nx; i++)
	{
		K_E2_a[i] = 1.0;
		K_E2_b[i] = 1.0/mu_0;

		K_E4_a[i] = 1.0/eps_0;
		K_E4_b[i] = 1.0/eps_0;

		if (i < nx-1)
		{
			K_E5_a[i] = 1.0;
			K_E5_b[i] = dt/dx;
		}			
	}

	for (j = 0; j<ny; j++)
	{
		K_E1_a[j] = 1.0;
		K_E1_b[j] = dt;

		K_E6_a[j] = 1.0/eps_0;
		K_E6_b[j] = 1.0/eps_0;

		if (j < ny-1)
		{
			K_E3_a[j] = 1.0;
			K_E3_b[j] = dt/dy;
		}			
	}
	
	//PML parameters
	double ka_max = 1;
	int exponent = 4;
	double R_err = 1e-16;
	
	eps_r = Mater[0][0]; //(Mater[0][0] + Mater[1][0])/2;
	mu_r  = Mater[0][1]; //(Mater[0][1] + Mater[1][1])/2;

    double eta = sqrt(mu_0*mu_r/eps_0/eps_r);

	double sigma_x, sigma_y, sigma_max, ka_x, ka_y;
		
	sigma_max= -(exponent+1)*log(R_err)/(2*eta*n_PML*dx);

	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_E2_a[i]       = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		K_E2_a[nx-i-1] = K_E2_a[i];

		K_E2_b[i]       = 2*eps_0/(2*eps_0*ka_x+sigma_x*dt)/mu_0;
		K_E2_b[nx-i-1] = K_E2_b[i];

		K_E4_a[i]       = (2*eps_0*ka_x + sigma_x*dt)/(2*eps_0*eps_0);
		K_E4_a[nx-i-1] = K_E4_a[i];

		K_E4_b[i]       = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*eps_0);
		K_E4_b[nx-i-1] = K_E4_b[i];

		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_E5_a[i]       = (2*eps_0*ka_x - sigma_x*dt)/(2*eps_0*ka_x + sigma_x*dt);
		K_E5_a[nx-i-2] = K_E5_a[i];

		K_E5_b[i]       = 2*eps_0*dt/(2*eps_0*ka_x + sigma_x*dt)/dx;
		K_E5_b[nx-i-2] = K_E5_b[i];
	}

	sigma_max = -(exponent+1)*log(R_err)/(2*eta*n_PML*dy);
	for (j = 0; j<n_PML; j++)
	{
		sigma_y         = sigma_max*pow( (n_PML - j)/((double) n_PML) ,exponent);
		ka_y            = 1 + (ka_max - 1)*pow( (n_PML-j)/((double) n_PML) ,exponent);
		K_E1_a[j]       = (2*eps_0*ka_y - sigma_y*dt)/(2*eps_0*ka_y + sigma_y*dt);
		K_E1_a[ny-j-1] = K_E1_a[j];

		K_E1_b[j]       = 2*eps_0*dt/(2*eps_0*ka_y + sigma_y*dt);
		K_E1_b[ny-j-1] = K_E1_b[j];

		K_E6_a[j]       = (2*eps_0*ka_y + sigma_y*dt)/(2*eps_0*eps_0);      
		K_E6_a[ny-j-1] = K_E6_a[j];

		K_E6_b[j]       = (2*eps_0*ka_y - sigma_y*dt)/(2*eps_0*eps_0);      
		K_E6_b[ny-j-1] = K_E6_b[j];

		sigma_y         = sigma_max*pow( (n_PML - j - 0.5)/n_PML, exponent);
		ka_y            = 1 + (ka_max - 1)*pow( (n_PML - j - 0.5)/n_PML, exponent);
		K_E3_a[j]       = (2*eps_0*ka_y - sigma_y*dt)/(2*eps_0*ka_y + sigma_y*dt);
		K_E3_a[ny-j-2] = K_E3_a[j];

		K_E3_b[j]       = 2*eps_0*dt/(2*eps_0*ka_y + sigma_y*dt)/dy;
		K_E3_b[ny-j-2] = K_E3_b[j];

	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hz field
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::calc_Hz_TE()
{
	for (i = 1; i<nx-1; i++)
	{
		for (j = 1; j<ny-1; j++)
		{
			Fz_r = Fz[i][j];

			Fz[i][j] = K_E1_a[j]*Fz[i][j] + K_E1_b[j]*( (Ex[i][j] - Ex[i][j-1])/dy -
  				                                        (Ey[i][j] - Ey[i-1][j])/dx );

			Hz[i][j] = K_E2_a[i]*Hz[i][j] + K_E2_b[i]*(Fz[i][j]-Fz_r)/Mater[Index[i][j]][1];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ex field
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::calc_Ex_TE()
{
	for (i = 0; i<nx; i++)
	{
		for (j = 0; j<ny-1; j++)
		{		
			Gx_r = Gx[i][j];

			Gx[i][j] = K_E3_a[j]*Gx[i][j] + K_E3_b[j]*( Hz[i][j+1] - Hz[i][j] );

			Ex[i][j] = Ex[i][j] + (K_E4_a[i]*Gx[i][j]-K_E4_b[i]*Gx_r)/Mater[Index[i][j]][0];
		}
	}

}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ey field
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE::calc_Ey_TE()
{
	for (i = 0; i<nx-1; i++)
	{
		for (j = 0; j<ny; j++)
		{		
			Gy_r = Gy[i][j];

			Gy[i][j] = K_E5_a[i]*Gy[i][j] - K_E5_b[i]*( Hz[i+1][j] - Hz[i][j] );

			Ey[i][j] = Ey[i][j] + (K_E6_a[j]*Gy[i][j]-K_E6_b[j]*Gy_r)/Mater[Index[i][j]][0];
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Gauss
///////////////////////////////////////////////////////////////////////////////////////