www.pudn.com > 2D-TM.rar > 2dpml_tm.cpp

 
 
#define source scp0(t)  //d_Guess pulse　第一步：加入激励源：微分高斯脉冲源 
                        //#define source scp(t) //Guess pulse　　高斯脉冲源 
                        //#define source sch(t)  // cos　　　　　升余弦脉冲源 
#include"parameters.h" 
#include"Myarray.h" 
#include"vector.h" 
#include"math.h" 
#include"2d_farfield.h" 
 
double phi=0;    //Incidence diretion　入射方向 
double *EXi=new double[dia+1];//Exi数组大小为124 
double *HYi=new double[dia];//Hyi数组为123 
 
double inwave(int which,int i,int j); 
double hxi(int ix,int jx)	{	return inwave(1,ix,jx);	}　　　//TM波 
double hyi(int ix,int jx)	{	return inwave(2,ix,jx);	} 
double ezi(int ix,int jx)	{	return inwave(3,ix,jx);	} 
 
void main(){//第三步 
 
	int i,j,t,T,flag;　 
	ofstream  fp; 
	fp.open("ez.dat",ios::trunc); 
 
	clock_t start, finish; 
	unsigned long int duration; 
 
	double	coef1[PML_Thickness+1],coef2[PML_Thickness+1],mcoef1[PML_Thickness+1],mcoef2[PML_Thickness+1];//PML层四个sigma参数 
　　　for(i=0;i<=PML_Thickness;i++){ 
		　coef1[i]=(1.0-0.5*dte*conduct(double(i)))/(1.0+0.5*dte*conduct(double(i))); 
		　coef2[i]=dte/(1.0+0.5*dte*conduct(double(i))); 
		　mcoef1[i]=(1.0-0.5*dtm*mconduct(double(i)))/(1.0+0.5*dtm*mconduct(double(i))); 
		　mcoef2[i]=dtm/(1.0+0.5*dtm*mconduct(double(i))); 
	　　　} 
 
	cout<<"Size="<>phi;　　　　　　　　　　　　　　　　　//之前已经定义过phi,即入射波的方向固定 
	cout<<"Input calculate time (N):"; 
	cin>>T; 
 
	double	**ez=D2dmatrix(MX+1,MY+1); 
	double	**psi=D2dmatrix(MX+1,MY+1); 
	double	**hx=D2dmatrix(MX+1,MY); 
	double  **hy=D2dmatrix(MX,MY+1); 
	PML_2D_Arr_ez	ezy(MX,PML_Thickness,PML_Thickness,MY,PML_Thickness,PML_Thickness); 
	for(i=0;i=ii0)&&(i<=ii1))	flag=0; 
				if((j==jj1)&&(i>=ii0)&&(i<=ii1))	flag=0; 
				if(flag==1)	 
					hx[i][j]+=-(dtm/dy)*(ez[i][j+1]-ez[i][j]); 
			}; 
 
		//boundary between PML zone and nonPML zone 
		for(i=I1+1;i=jj0)&&(j<=jj1))	flag=0; 
				if((i==ii1)&&(j>=jj0)&&(j<=jj1))	flag=0; 
				if(flag==1)	 
					hy[i][j]+=(dtm/dx)*(ez[i+1][j]-ez[i][j]); 
			}; 
 
		//boundary between PML zone and nonPML zone 
		for(j=J1+1;j=jj0)&&(j<=jj1))	flag=0; 
					if((i==ii1)&&(j>=jj0)&&(j<=jj1))	flag=0; 
					if((j==jj0)&&(i>=ii0)&&(i<=ii1))	flag=0; 
					if((j==jj1)&&(i>=ii0)&&(i<=ii1))	flag=0; 
					if(flag==1)								 
						ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j])-(dte/dy)*(hx[i][j]-hx[i][j-1]); 
				}; 
 
            ///target Ez:                        
			for(i=MX/2-Target_X;i<=MX/2+Target_X;i++) 
			for(j=MY/2-Target_Y;j<=MY/2+Target_Y;j++) 
			if(sqrt((i-MX/2)*(i-MX/2)+(j-MY/2)*(j-MY/2))<=Target_X)   //圆柱 
		    //	if((abs(i-MX/2)<=(Target_X))&&(abs(j-MY/2)<=(Target_X)))   //方柱 
				ez[i][j]=0.0; 
		 
		//PML zone 
		for(j=1;j=0&&phi<=90)//0<=phi<=90 
		rc=vector(ic-ii0,jc-jj0,0); 
	else 
		if(phi>90&&phi<=180)//90180&&phi<=270)//180270&&phi<360)//270