calculatefdtd.h

基于电磁波方程

void ComputeE()
{
    int i,j,I,J;
    bool myBL;
    float rotEz;
    for(i=IaMin+1;i<=IaMax-1;i++)
    {
        for(j=JaMin+1;j<=JaMax-1;j++)
        {
            myBL=(i>=IsMin&&j>=JsMin&&i<=IsMax&&j<=JsMax);
            rotEz=Hx[i][j-1]-Hx[i][j]+Hy[i][j]-Hy[i-1][j];
            if(myBL==1)
            {
            J=j-JsMin;
            I=i-IsMin;
            Ez[i][j]=FE1[ob[I][J]]*Ez[i][j]+FE2[ob[I][J]]*rotEz;
            }
            else Ez[i][j]+=0.5*rotEz;
        }
    }
}

//--------------------------------------------------------------
void CalculateBoundaryEz()
{
    int i,j,k,L,N;
    L=0;
	for(i=IaMin+1;i<=IaMax-1;i++)
	{
		TSMurEz_x(i,JaMin,JaMin+1);
		TSMurEz_x(i,JaMax,JaMax-1);
	}
    for(j=JaMin+1;j<=JaMax-1;j++)
	{
		TSMurEz_y(j,IaMin,IaMin+1);
		TSMurEz_y(j,IaMax,IaMax-1);
	}
	TSMurEz_Point(0,IaMin,JaMin,1,1);
	TSMurEz_Point(1,IaMin,JaMax,1,-1);
	TSMurEz_Point(2,IaMax,JaMin,-1,1);
	TSMurEz_Point(3,IaMax,JaMax,-1,-1);
}

 
//-------------------------------------------------------------
void ComputeH(void)
{
    int i,j;
    float rotHx,rotHy;
    for(i=IaMin;i<=IaMax;i++)
    {
        for(j=JaMin;j<=JaMax-1;j++)
        {
            rotHx=Ez[i][j+1]-Ez[i][j];
            Hx[i][j]=FH1*Hx[i][j]-FH2*rotHx;
        }
    }
	for(i=IaMin;i<=IaMax-1;i++)
    {
        for(j=JaMin;j<=JaMax;j++)
        {
            rotHy=Ez[i+1][j]-Ez[i][j];
            Hy[i][j]=FH1*Hy[i][j]+FH2*rotHy;
        }
    }
}

void computeFDTD_0()
{
	char filename[256];
	init256BMP(IaMax+1,JaMax+1);
	printf("jlsdljsdlfjlsj:%s\n",myDirectory);
    for(int N=1;N<=TimeStop;N++)
    {
		AddSourceToFDTD(1,N,20,80);
		CalculateBoundaryEz();
        ComputeE();
		ComputeH();
		if(N%50==0) printf("%d\n",N);
        if((N%12==0)&&(N<=1200)) 
		{	
			sprintf(filename,"%s/Ez_%07d.bmp",myDirectory,N/12);
			New256BMPFile(filename,Ez,IaMax+1,JaMax+1);
		}
		printf("%f\n",Ez[100][100]);
	}
	
}

void computeFDTD_1()
{
	char *file1="kk1.dat",*file2="kk1.dat";
	float **tEz,**tHx,**tHy;
    for(int N=1;N<=TimeStop+WL/2;N++)
    {
		AddSourceToFDTD(1,N,(IaMax-1)/2,(IaMax-1)/2);
		CalculateBoundaryEz();
        ComputeE();	
		ComputeH();
		if(N%50==0) printf("%d\n",N);
		if(N==TimeStop)	
		{
			tEz=ArrayCopy(Ez,IaMax+1,JaMax+1);
			tHx=ArrayCopy(Hx,IaMax+1,JaMax);
			tHy=ArrayCopy(Hy,IaMax,JaMax+1);
			outputEHn(file1,myDirectory);
		}
	}
	outputEHn(file2,myDirectory);
    ImaginaryRealCompute(tEz,Ez,IaMax+1,JaMax+1,"RealEz","ImaginaryEz");
	ImaginaryRealCompute(tHx,Hx,IaMax+1,JaMax,"RealHx","ImaginaryHx");
	ImaginaryRealCompute(tHy,Hy,IaMax,JaMax+1,"RealHy","ImaginaryHy");
}
void computeFDTD_2()
{
	char filename[256];
	init256BMP(IaMax-1,JaMax-1);
    for(int N=1;N<=TimeStop;N++)
    {
		AddSqareWave(2,N);
		CalculateBoundaryEz();
        ComputeE();
		IntroduceInE();
		ComputeH();
		IntroduceInH();
		if(N%40==0) printf("%d\n",N);
        if((N%12==0)&&(N<=1200)) 
		{
			sprintf(filename,"%s/Ez_%07d.bmp",myDirectory,N/12);
			New256BMPFile(filename,Ez,IaMax-1,JaMax-1);
		}
	}
}
void computeFDTD_3()
{
	float **tEz,**tHx,**tHy;
	char *file1="kk1.dat",*file2="kk1.dat";
    for(int N=1;N<=TimeStop+WL/2;N++)
    {
		AddSqareWave(1,N);
		CalculateBoundaryEz();
        ComputeE();
		IntroduceInE();
		ComputeH();
		IntroduceInH();
		if(N%40==0) printf("%d\n",N);
		if(N==TimeStop)	
		{
			tEz=ArrayCopy(Ez,IaMax+1,JaMax+1);
			tHx=ArrayCopy(Hx,IaMax+1,JaMax);
			tHy=ArrayCopy(Hy,IaMax,JaMax+1);
			outputEHn(file1,myDirectory);
		}
	}
	outputEHn(file2,myDirectory);
    ImaginaryRealCompute(tEz,Ez,IaMax+1,JaMax+1,"RealEz","ImaginaryEz");
	ImaginaryRealCompute(tHx,Hx,IaMax+1,JaMax,"RealHx","ImaginaryHx");
	ImaginaryRealCompute(tHy,Hy,IaMax,JaMax+1,"RealHy","ImaginaryHy");
}
void computeFDTD_5()
{
	float **Receivers;
	int offLength=1;//炮检距和偏移距
	int startRecord=IsMin+5;//发射开始点
	int transmitPot,receivePot;//发射点和接收点
	int RecordNo=200;//发射次数
	Receivers=DefineArray(RecordNo,TimeStop);
	transmitPot=startRecord;
	for(int i=0;i<RecordNo;i++)
	{	
		printf("第%d次扫描！\n",i);
		receivePot=transmitPot+offLength;
		for(int N=0;N<TimeStop;N++)
		{
			CalculateBoundaryEz();
			ComputeE();
			AddSourceToFDTD(1,N,transmitPot,20);//离地面高度为5个空间步长（25）
			ComputeH();
			Receivers[i][N]=Ez[receivePot][20];
		}
		transmitPot=startRecord+1;
	}
	WriteToGrd(Receivers,RecordNo,TimeStop,"receviers");
	FreeArray(Receivers,RecordNo);	
}
void computeFDTD_4()
{
	char filename[256];
	init256BMP(IaMax+1,JaMax+1);
	float **Receivers;
	int offLength=1,NN;//炮检距和偏移距
	int startRecord=JtMin;//接收开始点
	int RecordNo=JtMax-JtMin+1;//接收点数
	int endRecord=startRecord+RecordNo-1;
//	int i,j;
	Receivers=DefineArray(RecordNo,TimeStop/3);
	for(int N=0;N<TimeStop;N++)
	{
		AddSqareWave(1,N);
		CalculateBoundaryEz();
        ComputeE();
		IntroduceInE();
		ComputeH();
		IntroduceInH();
		if((N%3==0))
		{
			NN=N/3;
			printf("%d\n",NN);
			for(int j=0;j<RecordNo;j++)
			{
				Receivers[j][NN]=Ez[ItMax+5][JtMin+j];
			}
		}
		if((N%12==0)&&(N<=2000)) 
		{
			sprintf(filename,"%s/Ez_%07d.bmp",myDirectory,N/12);
			
			New256BMPFile(filename,Ez,IaMax+1,JaMax+1);
		}
	}
	WriteToGrd(Receivers,RecordNo,TimeStop/3,"receviers");
	FreeArray(Receivers,RecordNo);		
}