2d_farfield.h

该代码可以用来计算二维TM波的FDTD程序。 外加理想匹配层（PML）。 有一定的广泛性。

//See:IEEE TAP:vol39,no4,page 429
#define	LightSpeed	3.0e+8

const int	InitAngle=180;//180;
const int	Increament=1;
const int	AngleNum=1;
const double	RadPerDegree=pii/180.0;

class FarField{
private:
	int		i1,i2,j1,j2;
	vector	Origin;
public:
	vector	(*FieldValue)[AngleNum];
	double phi;
	int		index;	
	FarField(int nmax,int ii1,int ii2,int jj1,int jj2);
	void	NearToFarEz(int timestep,double** NearField);
	void	NearToFarHx(int timestep,double** NearField);
	void	NearToFarHy(int timestep,double** NearField);
	void	DataOut();
};

FarField::FarField(int nmax,int ii1,int ii2,int jj1,int jj2){
	int		in,an;

	i1=ii1;i2=ii2;j1=jj1;j2=jj2;
	Origin=vector((i1+i2)/2.0,(j1+j2)/2.0,0.0);//中心点
	index=nmax+(int)(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/(LightSpeed*dt))+1;//最大需存贮空间		
	FieldValue=new vector[index][AngleNum];

	for(in=0;in<index;in++)
		for(an=0;an<AngleNum;an++)
			FieldValue[in][an]=(0.0,0.0,0.0);
}

void FarField::NearToFarEz(int timestep,double** NearField){//近远场外推
	int	i,j,an;
	double tc;
	int m;
	vector norm,ms;

	for(an=0;an<AngleNum;an++){
		phi=(InitAngle+an*Increament)*RadPerDegree;

		//front: j=j1,i1<i<i2
		j=j1;
		norm=vector(0.0,-1.0,0.0);
		for(i=i1;i<i2;i++){
			ms=vector(0,0,0)-norm/vector(0.0,0.0,(NearField[i][j]+NearField[i+1][j])/2.0);
			tc=timestep*dt;
			tc-=(vector(i+0.5,j,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt+0.5);

			FieldValue[m-1][an]=FieldValue[m-1][an] + ms*dx*(0.5-tc/dt+m);
			FieldValue[m][an]=FieldValue[m][an] + ms*dx*2.0*(tc/dt-m);
			FieldValue[m+1][an]=FieldValue[m+1][an] - ms*dx*(0.5+tc/dt-m);
		}//i

		//back: j=j2,i1<i<i2
		j=j2;
		norm=vector(0.0,1.0,0.0);
		for(i=i1;i<i2;i++){
			ms=vector(0,0,0)-norm/vector(0.0,0.0,(NearField[i][j]+NearField[i+1][j])/2.0);
			tc=timestep*dt;
			tc-=(vector(i+0.5,j,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt+0.5);

			FieldValue[m-1][an]=FieldValue[m-1][an] + dx*ms*(0.5-tc/dt+m);
			FieldValue[m][an]=FieldValue[m][an] + dx*ms*2.0*(tc/dt-m);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dx*ms*(0.5+tc/dt-m);
		}//i

		//left: i=i1,j1<j<j2
		i=i1;
		norm=vector(-1.0,0.0,0.0);
		for(j=j1;j<j2;j++){
			ms=vector(0,0,0)-norm/vector(0.0,0.0,(NearField[i][j]+NearField[i][j+1])/2.0);
			tc=timestep*dt;
			tc-=(vector(i,j+0.5,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt+0.5);

			FieldValue[m-1][an]=FieldValue[m-1][an] + dy*ms*(0.5-tc/dt+m);
			FieldValue[m][an]=FieldValue[m][an] + dy*ms*2.0*(tc/dt-m);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dy*ms*(0.5+tc/dt-m);
		}//j

		//right: i=i2,j1<j<j2
		i=i2;
		norm=vector(1.0,0.0,0.0);
		for(j=j1;j<j2;j++){
			ms=vector(0,0,0)-norm/vector(0.0,0.0,(NearField[i][j]+NearField[i][j+1])/2.0);
			tc=timestep*dt;
			tc-=(vector(i,j+0.5,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt+0.5);

			FieldValue[m-1][an]=FieldValue[m-1][an] + dy*ms*(0.5-tc/dt+m);
			FieldValue[m][an]=FieldValue[m][an] + dy*ms*2.0*(tc/dt-m);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dy*ms*(0.5+tc/dt-m);
		}//j

	}//an
}

void FarField::NearToFarHx(int timestep,double** NearField){
	int	i,j,an;
	double tc;
	int m;
	vector norm,ms;

	for(an=0;an<AngleNum;an++){
		phi=(InitAngle+an*Increament)*RadPerDegree;

		//front: j=j1,i1<i<i2
		j=j1;
		norm=vector(0.0,-1.0,0.0);
		for(i=i1;i<i2;i++){
			ms=norm/vector((NearField[i][j-1]+NearField[i][j]+NearField[i+1][j-1]+NearField[i+1][j])/4.0,0.0,0.0);
			tc=timestep*dt;
			tc-=(vector(i+0.5,j,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt-0.5);//?

			FieldValue[m-1][an]=FieldValue[m-1][an] + ms*dx*(0.5-tc/dt+m+0.5);
			FieldValue[m][an]=FieldValue[m][an] + ms*dx*2.0*(tc/dt-m-0.5);
			FieldValue[m+1][an]=FieldValue[m+1][an] - ms*dx*(0.5+tc/dt-m-0.5);
		}//i

		//back: j=j2,i1<i<i2
		j=j2;
		norm=vector(0.0,1.0,0.0);
		for(i=i1;i<i2;i++){
			ms=norm/vector((NearField[i][j-1]+NearField[i][j]+NearField[i+1][j-1]+NearField[i+1][j])/4.0,0.0,0.0);
			tc=timestep*dt;
			tc-=(vector(i+0.5,j,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
			m=(int)(tc/dt-0.5);//?

			FieldValue[m-1][an]=FieldValue[m-1][an] + dx*ms*(0.5-tc/dt+m+0.5);
			FieldValue[m][an]=FieldValue[m][an] + dx*ms*2.0*(tc/dt-m-0.5);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dx*ms*(0.5+tc/dt-m-0.5);
		}//i

		//left: i=i1,j1<j<j2----No contribute as for Hx
		//right: i=i2,j1<j<j2----No contribute as for Hx

	}//an
}
void FarField::NearToFarHy(int timestep,double** NearField){　　
	int	i,j,an;
	double tc;
	int m;
	vector norm,ms;

	for(an=0;an<AngleNum;an++){
		phi=(InitAngle+an*Increament)*RadPerDegree;

		//front: j=j1,i1<i<i2----No contribute as for Hy
		//back: j=j2,i1<i<i2----No contribute as for Hy

		//left: i=i1,j1<j<j2
		i=i1;
		norm=vector(-1.0,0.0,0.0);
		for(j=j1;j<j2;j++){
			ms=norm/vector(0.0,(NearField[i-1][j]+NearField[i][j]+NearField[i-1][j+1]+NearField[i][j+1])/4.0,0.0);
			tc=timestep*dt;
			tc-=(vector(i,j+0.5,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
	//		tc+=fabs((vector(i2,j2,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed);
			m=(int)(tc/dt-0.5);//?

			FieldValue[m-1][an]=FieldValue[m-1][an] + dy*ms*(0.5-tc/dt+m+0.5);
			FieldValue[m][an]=FieldValue[m][an] + dy*ms*2.0*(tc/dt-m-0.5);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dy*ms*(0.5+tc/dt-m-0.5);
		}//j

		//right: i=i2,j1<j<j2
		i=i2;
		norm=vector(1.0,0.0,0.0);
		for(j=j1;j<j2;j++){
			ms=norm/vector(0.0,(NearField[i-1][j]+NearField[i][j]+NearField[i-1][j+1]+NearField[i][j+1])/4.0,0.0);
			tc=timestep*dt;
			tc-=(vector(i,j+0.5,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed;
			tc+=(sqrt((i2-i1)*(i2-i1)*dx*dx+(j2-j1)*(j2-j1)*dy*dy)/2.0)/LightSpeed;
	//		tc+=fabs((vector(i2,j2,0)-Origin)*vector(cos(phi)*dx,sin(phi)*dy,0)/LightSpeed);
			m=(int)(tc/dt-0.5);

			FieldValue[m-1][an]=FieldValue[m-1][an] + dy*ms*(0.5-tc/dt+m+0.5);
			FieldValue[m][an]=FieldValue[m][an] + dy*ms*2.0*(tc/dt-m-0.5);
			FieldValue[m+1][an]=FieldValue[m+1][an] - dy*ms*(0.5+tc/dt-m-0.5);
		}//j

	}//an
}

void FarField::DataOut(){
	int	in,an;
	ofstream  fp2;
	fp2.open("farfield0.dat",ios::trunc);

	for(in=0;in<index;in++){
		fp2<<endl;
		for(an=0;an<AngleNum;an++)
			fp2<<FieldValue[in][an].xx<<"   "<<FieldValue[in][an].yy<<"   "<<FieldValue[in][an].zz<<"   ";
	};
	fp2.close();
}


void RCS(int T,FarField Ez,FarField Hx,FarField Hy){  //main函数中调用求解“雷达散射截面”
	ofstream  fp1;
	fp1.open("rcs.dat",ios::trunc);
	ofstream  fp2;
	fp2.open("source.dat",ios::trunc);
	ofstream  fp3;
	fp3.open("farfield.dat",ios::trunc);
	cout<<"DFT calculating...."<<endl;

	int i,an,max=Ez.index;//index__最大需存贮空间
	double phi;
	vector u,w;
	double* ef=new double[max];
	double* ex=new double[max];
	double* ey=new double[max];
	double* ez=new double[max];

	//Source DFT: See—《王长清》Page:134
	int m,n, Ndft=T,Nfreq=1000;
	double df=1.0/(Ndft*dt);       //迭代总步数—Ndft；步数和dt越大，采样越密。Nfreq应越大。
	double* rGsource=new double[Ndft];   //real 
	double* iGsource=new double[Ndft];   //image
	double (*rGtarget)[AngleNum];
	double (*iGtarget)[AngleNum];
	rGtarget=new double[Ndft][AngleNum];
	iGtarget=new double[Ndft][AngleNum];

	for(m=0;m<Nfreq;m++){///步数和dt越大，采样越密。Nfreq应越大(m*df*1.0e-9)
		rGsource[m]=0.0;
		iGsource[m]=0.0;
		for(n=0;n<Ndft;n++){
			int t=n;
            #define sour scp0(t)  
			rGsource[m]+=dt*sour*cos(-2.0*pii*m*n/Ndft);
			iGsource[m]+=dt*sour*sin(-2.0*pii*m*n/Ndft);
		};
		fp2<<m*df*1.0e-9<<"\t"<<sqrt(rGsource[m]*rGsource[m]+iGsource[m]*iGsource[m])<<endl;
	};

	//Scatter field:
	for(m=0;m<Nfreq;m++){
		for(an=0;an<AngleNum;an++){ rGtarget[m][an]=0.0;iGtarget[m][an]=0.0;};
		for(i=0;i<max;i++)  //max:最大需存贮空间
		{
			if(m==10) fp3<<endl;
			for(an=0;an<AngleNum;an++){

			//calculate time domain far zone E field:
				phi=(InitAngle+an*Increament)*RadPerDegree;
				u=Ez.FieldValue[i][an]*(1.0/(4*pii*LightSpeed*dt));
				w=(Hx.FieldValue[i][an]+Hy.FieldValue[i][an])*(1.0/(4*pii*LightSpeed*dt));

				ez[i]=sqrt(me)*(-w.zz);// theta\-z
				ez[i]+=-u.xx*sin(phi)+u.yy*cos(phi);// phi\-X*sin+Y*cos
				ef[i]=-sqrt(me)*(-w.xx*sin(phi)+w.yy*cos(phi));// phi\-X*sin+Y*cos
				ef[i]+=-u.zz;// theta\-z
				ex[i]=-ef[i]*sin(phi);
				ey[i]=ef[i]*cos(phi);
		
				if(m==10) fp3<<ex[i]<<"   "<<ey[i]<<"   "<<ez[i]<<"   ";

			//DFT:
			//	if(i<Ndft)	{
				if(i<1300)	{
					n=i;
					rGtarget[m][an]+=dt*ez[n]*cos(-2.0*pii*m*n/Ndft);//ex=ey==0
					iGtarget[m][an]+=dt*ez[n]*sin(-2.0*pii*m*n/Ndft);
				};
			};
		};
	};

	//Calculate RCS:
	for(m=0;m<Nfreq;m++){
		fp1<<endl;
		fp1<<m*df*1.0e-9<<"\t";
		for(an=0;an<AngleNum;an++){
			double rcs;
			rcs=(rGtarget[m][an]*rGtarget[m][an]+iGtarget[m][an]*iGtarget[m][an]);
			rcs=rcs/(rGsource[m]*rGsource[m]+iGsource[m]*iGsource[m]);
			rcs=10*log10(2.0*pii*rcs*LightSpeed/(m*df));
			fp1<<rcs;
		};
	};

	cout<<"df= "<<df<<endl;
	fp1.close();
	fp2.close();
	fp3.close();
}