fdtd_2d_te_periodic.cpp

2DFDTD程序

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

#include "stdafx.h"
#include "fdtd_2D_TE_PML_period.h"
#include "FDTD_2D_TE_PERIODIC.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_PERIODIC::FDTD_2D_TE_PERIODIC()
{
	Hz = NULL; Ex = NULL; Gx = NULL; Ey = NULL; Gy = NULL;  

	K_E2_a = NULL; K_E2_b = NULL;
	K_E4_a = NULL; K_E4_b = NULL;
	K_E5_a = NULL; K_E5_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_PERIODIC::~FDTD_2D_TE_PERIODIC()
{
	Free_Mem();
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Main Parameters
///////////////////////////////////////////////////////////////////////////////////////
BOOL FDTD_2D_TE_PERIODIC::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);
	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_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_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;
	}

	return TRUE;

}

///////////////////////////////////////////////////////////////////////////////////////
//Set the PML matrices
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Set_PML_Param()
{
	for (i = 0; i<nx; i++)
	{
		K_E2_a[i] = 1.0;
		K_E2_b[i] = dt/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;
		}			
	}

	//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_max, ka_x;
		
	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*dt/(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];
	}

}

///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hz field
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::calc_Hz_TE()
{
	double onepdx = 1.0/dx;
	double onepdy = 1.0/dy;
	for (i = 1; i<nx-1; i++)
	{
		double prod_1 = K_E2_b[i]*onepdy;
		double prod_2 = K_E2_b[i]*onepdx;
		for (j = 1; j<ny; j++)
		{
			Hz[i][j] = K_E2_a[i]*Hz[i][j] + ( prod_1*(Ex[i][j] - Ex[i][j-1]) -
  				                       prod_2*(Ey[i][j] - Ey[i-1][j])  )/Mater[Index[i][j]][1];
		}
	}
}

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

			Gx[i][j] = Gx[i][j] + dtpdy*( 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_PERIODIC::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] + (Gy[i][j]-Gy_r)/(eps_0*Mater[Index[i][j]][0]);
		}
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Gauss
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Init_Gauss_Point_Source(int n_x_p, int n_y_p, double h0, double t_0, 
										 double t_w)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   jel_source_type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Sinusoidal
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Init_Sinus_Point_Source(int n_x_p, int n_y_p, double h0, double om, 
										 double phase)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   omega = om;
   phi = phase;
   jel_source_type = 2;
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Point Source -- Sinusoidal Modulated Gauss Pulse 
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Init_PulseGS_Point_Source(int n_x_p, int n_y_p, double h0, double t_0, 
							 			   double t_w, double om, double phase)
{
   n_x_P = n_x_p; 
   n_y_P = n_y_p;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   omega = om;
   phi = phase;
   jel_source_type = 3;
}

///////////////////////////////////////////////////////////////////////////////////////
//Soft Point Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Point_Source(double time)
{
	switch (jel_source_type)
	{
		case 1: //Gaussian pulse
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] + H0*exp( -pow( (time-t0)/tw ,2) ); 
			break;
		case 2: //Sinusoidal plane wave
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] = H0*cos( omega*time + phi);
			break;
		case 3: //Wave packet(sinusoidal modulated Gaussian pulse)
			Hz[n_x_P][n_y_P] = Hz[n_x_P][n_y_P] = H0*cos( omega*time + phi)*
				                                     exp( -pow( (time-t0)/tw ,2) );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Line Gauss Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Init_Gauss_Line_Source(int **&Co, int n_Co, double h0, double t_0, 
										 double t_w)
{
   Coord = Co;
   n_Coord = n_Co;
   H0 = h0;
   t0 = t_0;
   tw = t_w; 
   jel_source_type = 1;
}

///////////////////////////////////////////////////////////////////////////////////////
//Soft Point Source
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Line_Source(double time)
{
	if (jel_source_type == 1) //Gauss pulse
	{
		for (i = 0; i< n_Coord; i++)
			Hz[Coord[i][0]][Coord[i][1]] = Hz[Coord[i][0]][Coord[i][1]] +
			                                       H0*exp( -pow((time-t0)/tw,2) );
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Init Total field - Scattered field formulation
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Init_TotFScatF(int n_xa, int n_xb, double H0, double t0, double tw, 
								double omega, double phi, int jel)
{
	n_x_a = n_xa;
	n_x_b = n_xb;

	Init_Main_Param_1D(n_x_b+1, n_PML, Mater[0][0], Mater[0][1], dt, dx);
	switch (jel)
	{
		case 1:
			Init_Gauss_1D(H0, t0, tw);
			break;
		case 2:
			Init_Sinus_1D(H0, omega, phi);
			break;
		case 3:
			Init_Wave_Packet_1D(H0, t0, tw, omega, phi);
	}
	Set_PML_Param_1D();
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Hz
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::incident_Hz()
{
	for (j = 0; j <= ny; j++)
	{
		Hz[n_x_a][j] = Hz[n_x_a][j] + 
						dt/(mu_0*Mater[Index[n_x_a][j]][1]*dx)*Ey_1D[0];
		Hz[n_x_b][j] = Hz[n_x_b][j] - 
						dt/(mu_0*Mater[Index[n_x_b][j]][1]*dx)*Ey_1D[n_x_b-n_x_a+1];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Ex
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::incident_Ex()
{
/*	for (i = n_x_a; i <= n_x_b; i++)
	{
		Ex[i][n_y_a-1] = Ex[i][n_y_a-1] -
						  dt/(eps_0*Mater[Index[i][n_y_a-1]][0]*dy)*Hz_1D[i-n_x_a];
		Ex[i][n_y_b] =   Ex[i][n_y_b] +
						  dt/(eps_0*Mater[Index[i][n_y_b]][0]*dy)*Hz_1D[i-n_x_a];
	}*/
}

///////////////////////////////////////////////////////////////////////////////////////
//Total field - Scattered field formulation -- Incident Ey
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::incident_Ey()
{
	for (j = 0; j <= ny; j++)
	{
		Ey[n_x_a-1][j] = Ey[n_x_a-1][j] +
						  dt/(eps_0*Mater[Index[n_x_a-1][j]][0]*dx)*Hz_1D[1];
		Ey[n_x_b][j] = Ey[n_x_b][j] -
						  dt/(eps_0*Mater[Index[n_x_b][j]][0]*dx)*Hz_1D[n_x_b-n_x_a+1];
	}

}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Hz
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::periodic_Hz()
{
	for (i = 0; i < nx; i++)
	{
		Hz[i][0] = Hz[i][ny-1];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Periodic boundary condition for Ex
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::periodic_Ex()
{
	for (i = 0; i < nx; i++)
	{
		Ex[i][ny-1] = Ex[i][0];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Free the allocated memory
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Free_Mem()
{
	if(Hz)
		Hz = Free_Matrix_2D<double>(Hz);
	if(Ex)
		Ex = Free_Matrix_2D<double>(Ex);
	if(Gx)
		Gx = Free_Matrix_2D<double>(Gx);
	if(Ey)
		Ey = Free_Matrix_2D<double>(Ey);
    if(Gy)
		Gy = Free_Matrix_2D<double>(Gy);
	
	if(K_E2_a)
	{
		free(K_E2_a); 	
		K_E2_a = NULL;
	}
	if(K_E2_b)
	{
		free(K_E2_b);
		K_E2_b = NULL;
	}
	if(K_E4_a)
	{
		free(K_E4_a); 	
		K_E4_a = NULL;
	}
	if(K_E4_b)
	{
		free(K_E4_b);
		K_E4_b = NULL;
	}
	if(K_E5_a)
	{
		free(K_E5_a);
		K_E5_a = NULL;
	}
	if(K_E5_b)
	{
		free(K_E5_b);
		K_E5_b = NULL;
	}
	
	if(Hz_Foll)
		Hz_Foll = Free_Matrix_2D<double>(Hz_Foll);
	if(Ex_Foll)
		Ex_Foll = Free_Matrix_2D<double>(Ex_Foll);
	if(Ey_Foll)
		Ey_Foll = Free_Matrix_2D<double>(Ey_Foll);
}

///////////////////////////////////////////////////////////////////////////////////////
//Initialize the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
BOOL FDTD_2D_TE_PERIODIC::Init_Followed(int **Ind_Followed, int n, int n_t)
{
	Ind_Foll = Ind_Followed;
	length_Ind_Foll = n;
	n_tot = n_t;

	Hz_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Hz_Foll)	
	{
		Free_Mem();
		return FALSE;
	}

	Ex_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ex_Foll)	
	{
		Free_Mem();
		return FALSE;
	}

	Ey_Foll = Init_Matrix_2D<double>(n,n_t);
	if(!Ey_Foll)	
	{
		Free_Mem();
		return FALSE;
	}

	return TRUE;
}

///////////////////////////////////////////////////////////////////////////////////////
//Sets the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Set_Data_Followed(int n_t)
{
	for (i = 0; i<length_Ind_Foll; i++)
	{
		Hz_Foll[i][n_t] = Hz[Ind_Foll[i][0]][Ind_Foll[i][1]];
		Ex_Foll[i][n_t] = Ex[Ind_Foll[i][0]][Ind_Foll[i][1]];
		Ey_Foll[i][n_t] = Ey[Ind_Foll[i][0]][Ind_Foll[i][1]];
	}
}

///////////////////////////////////////////////////////////////////////////////////////
//Returns the Followed field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Get_Data_Followed(double **&A, double **&B, double **&C)
{
	A = Hz_Foll;
	B = Ex_Foll;
	C = Ey_Foll;
}
	

///////////////////////////////////////////////////////////////////////////////////////
//Returns the field components
///////////////////////////////////////////////////////////////////////////////////////
void FDTD_2D_TE_PERIODIC::Get_Data(double **&A, double **&B, double **&C)
{
	A = Hz;
	B = Ex;
	C = Ey;
}