📄 pmlu.h

📁 二维FDTD算法
💻 H
📖 第 1 页 / 共 2 页
字号:
12 下一页
#ifndef PML_H
#define PML_H

#include "Def.h"


double PMLArrayE[NPML+1];		    //NPML数组(E->H)
double PMLArrayH[NPML+1];		    //NPML数组(H->E)

static double Dz[DOMAINX][DOMAINY];
static double iDz[DOMAINX][DOMAINY];

static double Bx[DOMAINX][DOMAINY];
static double iBx[DOMAINX][DOMAINY];

static double By[DOMAINX][DOMAINY];
static double iBy[DOMAINX][DOMAINY];


//起始中止点
int PMLBeginx=0;
int PMLEndx=0;
int PMLBeginy=0;
int PMLEndy=0;

double CA;
double CB;
double CP;
double CQ;
double C1;
double C2;
double C3;
double C4;

//初始化PML数组(PML区到散射场区，由大变小)
void InitialPMLArray()
{
	double m=4.0;
	double sigmaMax=(m+1)/(150.*PI*DELTAX);
	
	//电场渐变阻抗
	for (int i=0; i<NPML; i++)
	{
		double thick=(double)(NPML-i)/(double)NPML;		
		PMLArrayE[i]=sigmaMax*pow(thick,m);
	}
		
	//磁场渐变阻抗
	for (int j=0; j<NPML; j++)
	{
		double thick=(double)(NPML-j-0.5)/(double)NPML;
	    PMLArrayH[j]=sigmaMax*pow(thick,m);
	}
	PMLArrayE[NPML]=0.;
    PMLArrayH[NPML]=0.;

}



//PML初始化
void InitialPML()
{
	InitialPMLArray();  //阻抗	
	for (int i=0; i<DOMAINX;i++)
	{
	  for(int j=0;j<DOMAINY;j++)
	  {
	   Dz[i][j]=0.0;
	   Bx[i][j]=0.0;
	   By[i][j]=0.0;
	   iDz[i][j]=0.0;
	   iBx[i][j]=0.0;
	   iBy[i][j]=0.0;
	  }
	}
}
/********************************************************************/
//1.Hx,Hy推导场Dz，Dz推导Ez
void  PMLEz(double Coefficient1)
{
	int i,j;  //循环变量
	/**************************X平行方向******************************/
	
	PMLBeginx=NPML+1;
	PMLEndx=DOMAINX-NPML-1;
	//前--观察方向沿X轴正向
    PMLBeginy=1;
	PMLEndy=NPML+1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{	
		     C1=(2*EPSILON0-DELTAT*PMLArrayE[j])/(2*EPSILON0+DELTAT*PMLArrayE[j]);
			 C2=1/Z0*2/(2*EPSILON0+PMLArrayE[j]*DELTAT);
		     Dz[i][j]=Dz[i][j]+1/V*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));											

			 Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                          
             iDz[i][j]=Dz[i][j];
		}
	}

	//后--观察方向沿X轴正向
	PMLBeginy=DOMAINY-NPML-1;
	PMLEndy=DOMAINY-1;
	for (i=PMLBeginx; i<PMLEndx; i++)
	{
		for (j=PMLBeginy; j<PMLEndy; j++)
		{
            C1=(2*EPSILON0-DELTAT*PMLArrayE[DOMAINY-j-1])/(2*EPSILON0+DELTAT*PMLArrayE[DOMAINY-j-1]);
			C2=1/Z0*2/(2*EPSILON0+PMLArrayE[DOMAINY-j-1]*DELTAT);
			Dz[i][j]=Dz[i][j]+1/V*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
												
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                          
            iDz[i][j]=Dz[i][j];
		}
	}

/**************************Y平行方向*****************************/
	PMLBeginy=NPML+1;
	PMLEndy=DOMAINY-NPML-1;

	//前--观察方向沿Y轴正向
	PMLBeginx=2;
	PMLEndx=NPML+1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{
            CA=(1/DELTAT-PMLArrayE[i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0))/(V*DELTAT);
            C1=1.0;
			C2=1/Z0*1/EPSILON0;
			Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
										          
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                         
            iDz[i][j]=Dz[i][j];
		}
	}

	//后--观察方向沿Y轴正向
	PMLBeginx=DOMAINX-NPML-1;
	PMLEndx=DOMAINX-1;
    for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{	
			CA=(1/DELTAT-PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(V*DELTAT);
		    C1=1.0;
			C2=1/Z0*1/EPSILON0;
			Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
													
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                         
            iDz[i][j]=Dz[i][j];
		}
	}
	 //角区1
	PMLBeginx=1;
	PMLEndx=NPML+1;
    PMLBeginy=1;
	PMLEndy=NPML+1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{	
			CA=(1/DELTAT-PMLArrayE[i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0))/(V*DELTAT);
            C1=(2*EPSILON0-DELTAT*PMLArrayE[j])/(2*EPSILON0+DELTAT*PMLArrayE[j]);
			C2=1/Z0*2/(2*EPSILON0+PMLArrayE[j]*DELTAT);

			Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
													
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                       
            iDz[i][j]=Dz[i][j];
		}
	}
    //角区2
    PMLBeginx=1;
	PMLEndx=NPML+1;
	PMLBeginy=DOMAINY-NPML-1;
	PMLEndy=DOMAINY-1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{		
            CA=(1/DELTAT-PMLArrayE[i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[i]/(2*EPSILON0))/(V*DELTAT);
            C1=(2*EPSILON0-DELTAT*PMLArrayE[DOMAINY-j-1])/(2*EPSILON0+DELTAT*PMLArrayE[DOMAINY-j-1]);
			C2=1/Z0*2/(2*EPSILON0+PMLArrayE[DOMAINY-j-1]*DELTAT);

		    Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
													   
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                          
            iDz[i][j]=Dz[i][j];
		}
	}
	//角区3
	PMLBeginx=DOMAINX-NPML-1;
	PMLEndx=DOMAINX-1;
	PMLBeginy=DOMAINY-NPML-1;
	PMLEndy=DOMAINY-1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{		
            CA=(1/DELTAT-PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(V*DELTAT);
            C1=(2*EPSILON0-DELTAT*PMLArrayE[DOMAINY-j-1])/(2*EPSILON0+DELTAT*PMLArrayE[DOMAINY-j-1]);
			C2=1/Z0*2/(2*EPSILON0+PMLArrayE[DOMAINY-j-1]*DELTAT);
			Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
													   
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                          
            iDz[i][j]=Dz[i][j];
			
		}
	}
	//角区4
	PMLBeginx=DOMAINX-NPML-1;
	PMLEndx=DOMAINX-1;
    PMLBeginy=1;
	PMLEndy=NPML+1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{		
			CA=(1/DELTAT-PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0));
            CB=1/(1/DELTAT+PMLArrayE[DOMAINX-1-i]/(2*EPSILON0))/(V*DELTAT);
            C1=(2*EPSILON0-DELTAT*PMLArrayE[j])/(2*EPSILON0+DELTAT*PMLArrayE[j]);
			C2=1/Z0*2/(2*EPSILON0+PMLArrayE[j]*DELTAT);
			Dz[i][j]=CA*Dz[i][j]+CB*Coefficient1*((Hy[i][j]-Hy[i-1][j])-\
							              (Hx[i][j]-Hx[i][j-1]));
													
			Ez[i][j]=C1*Ez[i][j]+C2*(Dz[i][j]-iDz[i][j]);                              
            iDz[i][j]=Dz[i][j];
		}
	}
	
  for (i=0;i<DOMAINX;i++)
  {
	  Ez[i][0]=0.;
	  Ez[i][DOMAINY-1]=0.;

  }
  for (j=0;j<DOMAINY;j++)
  {
      Ez[0][j]=0.;
	  Ez[DOMAINX-1][j]=0.;
  }
}

/**********************************************************************/
/**********************************************************************/
//2.Ez推导场Bx;Bx推导Hx
void  PMLHx(double Coefficient1)
{
	int i,j;  //循环变量
	/**************************X平行方向******************************/
	
	PMLBeginx=NPML+1;
	PMLEndx=DOMAINX-NPML-1;
	//前--观察方向沿X轴正向
    PMLBeginy=1;
	PMLEndy=NPML;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
		{	 
			 CP=(1/DELTAT-PMLArrayH[j]/(2*EPSILON0))/(1/DELTAT+PMLArrayH[j]/(2*EPSILON0));
		     CQ=Z0*1/(1/DELTAT+PMLArrayH[j]/(2*EPSILON0))/(V*DELTAT);
			 C3=1/MIU0;
			 C4=1/MIU0;

		     Bx[i][j]=CP*Bx[i][j]+CQ*Coefficient1*(Ez[i][j]-Ez[i][j+1]);				                      

			 Hx[i][j]=Hx[i][j]+C3*Bx[i][j]-C4*iBx[i][j];                                
             iBx[i][j]=Bx[i][j];
			}
	}

	//后--观察方向沿X轴正向
	PMLBeginy=DOMAINY-NPML-1;
	PMLEndy=DOMAINY-1;
	for ( i=PMLBeginx; i<PMLEndx; i++)
	{
		for ( j=PMLBeginy; j<PMLEndy; j++)
12 下一页
💿 文件大小 485 K
👤 上传用户 snowpilce
📂 所属分类数值算法/人工智能
🏷️ 相关标签

#FDTD #二维 #算法
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -