#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="<<MX<<"*"<<MY<<endl;
	cout<<"BandWide factor: "<<tao<<"\t"<<"Time step: "<<dt<<endl;//tao,dt
//	cout<<"Input wave direction (Phi):";
//	cin>>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<dia+1;i++) EXi[i]=0.0;//背景全部設為0
	for(i=0;i<dia;i++) HYi[i]=0.0;

	FarField FarEz(T,ii0-2,ii1+2,jj0-2,jj1+2);
	FarField FarHx(T,ii0-2,ii1+2,jj0-2,jj1+2);
	FarField FarHy(T,ii0-2,ii1+2,jj0-2,jj1+2);

	//loop T
	start=clock();
	int fg=0;
	for(t=1;t<T;t++){//開始疊代

		FarEz.NearToFarEz(t,ez);
		FarHx.NearToFarHx(t,hx);
		FarHy.NearToFarHy(t,hy);

		cout<<t<<"\t"<<ez[I1+5][I1+5]<<endl;

		// Incident H field iteration !	
		int K1=PML_Thickness;
		for(i=0;i<K1;i++) HYi[i]=mcoef1[K1-i-1]*HYi[i]+(mcoef2[K1-i-1]/dx)*(EXi[i]-EXi[i+1]);
		for(i=K1;i<dia-K1;i++)  HYi[i]+=(dtm/dx)*(EXi[i]-EXi[i+1]);
		for(i=dia-K1;i<dia;i++) HYi[i]=mcoef1[i-(dia-K1)]*HYi[i]+(mcoef2[i-(dia-K1)]/dx)*(EXi[i]-EXi[i+1]);
		// Incident E field iteration !		
		for(i=1;i<=K1;i++)  EXi[i]=coef1[K1-i]*EXi[i]+(coef2[K1-i]/dx)*(HYi[i-1]-HYi[i]);
		for(i=K1+1;i<dia-K1;i++)  EXi[i]+=(dte/dx)*(HYi[i-1]-HYi[i]);
		for(i=dia-K1;i<dia;i++)  EXi[i]=coef1[i-(dia-K1)]*EXi[i]+(coef2[i-(dia-K1)]/dx)*(HYi[i-1]-HYi[i]); 
		EXi[2*PML_Thickness]=source;//在連接邊界(10)處加入源

		// H field iteration in calculated zone on nth step begins!　計算連接邊界內的總場區

//hx
		//front HX: 
		j=jj0-1;
		for(i=ii0;i<=ii1;i++)
			hx[i][j]+=-(dtm/dy)*(ez[i][j+1]-ezi(i,j+1)-ez[i][j]);					
		//back HX: 
		j=jj1;
		for(i=ii0;i<=ii1;i++)
			hx[i][j]+=-(dtm/dy)*(ez[i][j+1]-ez[i][j]+ezi(i,j));
		//free space:
		for(i=I1+1;i<I2;i++)
			for(j=J1+1;j<J2-1;j++)
			{
				flag=1;
				if((j==(jj0-1))&&(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<I2;i++){
			hx[i][J1]=hx[i][J1]-dtm*(ez[i][J1+1]-ez[i][J1]-ezy(i,J1))/dy;
			hx[i][J2-1]=hx[i][J2-1]-dtm*(ez[i][J2]+ezy(i,J2)-ez[i][J2-1])/dy;
		}

		//PML zone
		for(i=0;i<=MX;i++){
			for(j=0;j<J1;j++)
				hx[i][j]=mcoef1[J1-j-1]*hx[i][j]-mcoef2[J1-j-1]*(ez[i][j+1]+ezy(i,j+1)-ez[i][j]-ezy(i,j))/dy;
			for(j=J2;j<MY;j++)
				hx[i][j]=mcoef1[j-J2]*hx[i][j]-mcoef2[j-J2]*(ez[i][j+1]+ezy(i,j+1)-ez[i][j]-ezy(i,j))/dy;
		}
		for(j=J1;j<J2;j++){
			for(i=0;i<=I1;i++)
				hx[i][j]=hx[i][j]-dtm*(ez[i][j+1]+ezy(i,j+1)-ez[i][j]-ezy(i,j))/dy;
			for(i=I2;i<=MX;i++)
				hx[i][j]=hx[i][j]-dtm*(ez[i][j+1]+ezy(i,j+1)-ez[i][j]-ezy(i,j))/dy;
		}
//hy
		//left HY: 
		i=ii0-1;
		for(j=jj0;j<=jj1;j++)
			hy[i][j]+=(dtm/dx)*(ez[i+1][j]-ezi(i+1,j)-ez[i][j]);				
		//right HY: 
		i=ii1;
		for(j=jj0;j<=jj1;j++)
			hy[i][j]+=(dtm/dx)*(ez[i+1][j]-ez[i][j]+ezi(i,j));		
		//free space:
		for(j=J1+1;j<J2;j++)
			for(i=I1+1;i<I2-1;i++)
			{
				flag=1;
				if((i==(ii0-1))&&(j>=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<J2;j++){
			hy[I1][j]=hy[I1][j]+dtm*(ez[I1+1][j]-ez[I1][j]-ezy(I1,j))/dx;
			hy[I2-1][j]=hy[I2-1][j]+dtm*(ez[I2][j]+ezy(I2,j)-ez[I2-1][j])/dx;
		}

		//PML zone
		for(j=0;j<=MY;j++){
			for(i=0;i<I1;i++)
				hy[i][j]=mcoef1[I1-i-1]*hy[i][j]+mcoef2[I1-i-1]*(ez[i+1][j]+ezy(i+1,j)-ez[i][j]-ezy(i,j))/dx;
			for(i=I2;i<MX;i++)
				hy[i][j]=mcoef1[i-I2]*hy[i][j]+mcoef2[i-I2]*(ez[i+1][j]+ezy(i+1,j)-ez[i][j]-ezy(i,j))/dx;
		}
		for(i=I1;i<I2;i++){
			for(j=0;j<=J1;j++)
				hy[i][j]=hy[i][j]+dtm*(ez[i+1][j]+ezy(i+1,j)-ez[i][j]-ezy(i,j))/dx;
			for(j=J2;j<=MY;j++)
				hy[i][j]=hy[i][j]+dtm*(ez[i+1][j]+ezy(i+1,j)-ez[i][j]-ezy(i,j))/dx;
		}
//ez
			//left EZ
			i=ii0;
			for(j=jj0+1;j<=jj1-1;j++)
				ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j]-hyi(i-1,j))-(dte/dy)*(hx[i][j]-hx[i][j-1]);
			i=ii0;	j=jj0;
			ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j]-hyi(i-1,j))-(dte/dy)*(hx[i][j]-hx[i][j-1]-hxi(i,j-1));
			i=ii0;	j=jj1;
			ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j]-hyi(i-1,j))-(dte/dy)*(hx[i][j]+hxi(i,j)-hx[i][j-1]);
			//right EZ:	
			i=ii1;
			for(j=jj0+1;j<=jj1-1;j++)
				ez[i][j]+=(dte/dx)*(hy[i][j]+hyi(i,j)-hy[i-1][j])-(dte/dy)*(hx[i][j]-hx[i][j-1]);
			i=ii1;	j=jj0;
			ez[i][j]+=(dte/dx)*(hy[i][j]+hyi(i,j)-hy[i-1][j])-(dte/dy)*(hx[i][j]-hx[i][j-1]-hxi(i,j-1));
			i=ii1;	j=jj1;
			ez[i][j]+=(dte/dx)*(hy[i][j]+hyi(i,j)-hy[i-1][j])-(dte/dy)*(hx[i][j]+hxi(i,j)-hx[i][j-1]);
			//front EZ:
			j=jj0;
			for(i=ii0+1;i<=ii1-1;i++)
				ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j])-(dte/dy)*(hx[i][j]-hx[i][j-1]-hxi(i,j-1));
			//back EZ:
			j=jj1;
			for(i=ii0+1;i<=ii1-1;i++)
				ez[i][j]+=(dte/dx)*(hy[i][j]-hy[i-1][j])-(dte/dy)*(hx[i][j]+hxi(i,j)-hx[i][j-1]);
	
            //free space
			for(i=I1+1;i<I2;i++)
				for(j=J1+1;j<J2;j++)
				{
					flag=1;
					if((i==ii0)&&(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<MY;j++){
			for(i=1;i<=I1;i++)
				ez[i][j]=coef1[I1-i]*ez[i][j]+coef2[I1-i]*(hy[i][j]-hy[i-1][j])/dx;
			for(i=I2;i<MX;i++)
				ez[i][j]=coef1[i-I2]*ez[i][j]+coef2[i-I2]*(hy[i][j]-hy[i-1][j])/dx;
		}
		for(i=I1+1;i<I2;i++){
			for(j=1;j<=J1;j++)
				ez[i][j]=ez[i][j]+dte*(hy[i][j]-hy[i-1][j])/dx;
			for(j=J2;j<MY;j++)
				ez[i][j]=ez[i][j]+dte*(hy[i][j]-hy[i-1][j])/dx;
		}
		for(i=1;i<MX;i++){
			for(j=1;j<=J1;j++)
				ezy(i,j)=coef1[J1-j]*ezy(i,j)-coef2[J1-j]*(hx[i][j]-hx[i][j-1])/dy;
			for(j=J2;j<MY;j++)
				ezy(i,j)=coef1[j-J2]*ezy(i,j)-coef2[j-J2]*(hx[i][j]-hx[i][j-1])/dy;
		}
		for(j=J1+1;j<J2;j++){
			for(i=1;i<=I1;i++)
				ezy(i,j)=ezy(i,j)-dte*(hx[i][j]-hx[i][j-1])/dy;
			for(i=I2;i<MX;i++)
				ezy(i,j)=ezy(i,j)-dte*(hx[i][j]-hx[i][j-1])/dy;
		}


	}// for time iteration

	//farfield:
	RCS(T,FarEz,FarHx,FarHy);


	finish=clock();
	cout<<"The  program takes"<<endl;
	duration=(finish-start)/CLOCKS_PER_SEC;
	cout<<duration/3600<<"\thours"<<endl<<(duration%3600)/60<<"\tminutes"<<endl<<(duration%3600)%60<<"\tseconds"<<endl;

	delete hx;
	delete hy;
	delete ez;
	delete EXi;
	delete HYi;
	fp.close();
}//main




//sub_rutine for incident wave:第二步：設置入射波
double inwave(int which,int ic,int jc)
{

	vector ki(cos(phi*pii/180),sin(phi*pii/180),0);
	vector rc;

	phi=int(phi)%360;//%為取余數
	if(phi>=0&&phi<=90)//0<=phi<=90
		rc=vector(ic-ii0,jc-jj0,0);
	else
		if(phi>90&&phi<=180)//90<phi<=180
			rc=vector(ic-ii1,jc-jj0,0);
		else
			if(phi>180&&phi<=270)//180<phi<=270
				rc=vector(ic-ii1,jc-jj1,0);	
			else
				if(phi>270&&phi<360)//270<phi<=360
					rc=vector(ic-ii0,jc-jj1,0);	
	if(which==1) rc.yy+=0.5;   //磁場Hx,Hy
	if(which==2) rc.xx+=0.5;	
				
	//三次插值近似
	double d=ki*rc;
	int x1=int(d);
	int x0=x1-1;
	int x2=x1+1;
	int x3=x1+2;
	double ei=EXi[2*PML_Thickness+10+x0]*(d-x1)*(d-x2)*(d-x3)/((x0-x1)*(x0-x2)*(x0-x3));
	ei+=EXi[2*PML_Thickness+10+x1]*(d-x0)*(d-x2)*(d-x3)/((x1-x0)*(x1-x2)*(x1-x3));
	ei+=EXi[2*PML_Thickness+10+x2]*(d-x1)*(d-x0)*(d-x3)/((x2-x1)*(x2-x0)*(x2-x3));
	ei+=EXi[2*PML_Thickness+10+x3]*(d-x1)*(d-x2)*(d-x0)/((x3-x1)*(x3-x2)*(x3-x0));
	
	double hi;
	d=d-0.5;//關鍵
	if((d-int(d))<0.5){
		x1=int(d);
		x0=x1-1;
		x2=x1+1;
		x3=x1+2;
		d=d+0.5;
		hi=HYi[2*PML_Thickness+10+x0-1]*(d-x1)*(d-x2)*(d-x3)/((x0-x1)*(x0-x2)*(x0-x3));
		hi+=HYi[2*PML_Thickness+10+x1-1]*(d-x0)*(d-x2)*(d-x3)/((x1-x0)*(x1-x2)*(x1-x3));
		hi+=HYi[2*PML_Thickness+10+x2-1]*(d-x1)*(d-x0)*(d-x3)/((x2-x1)*(x2-x0)*(x2-x3));
		hi+=HYi[2*PML_Thickness+10+x3-1]*(d-x1)*(d-x2)*(d-x0)/((x3-x1)*(x3-x2)*(x3-x0));
		d=d-0.5;
	}
	else{
		x1=int(d);
		d=d-0.5;
		x0=x1-1;
		x2=x1+1;
		x3=x1+2;
		hi=HYi[2*PML_Thickness+10+x0]*(d-x1)*(d-x2)*(d-x3)/((x0-x1)*(x0-x2)*(x0-x3));
		hi+=HYi[2*PML_Thickness+10+x1]*(d-x0)*(d-x2)*(d-x3)/((x1-x0)*(x1-x2)*(x1-x3));
		hi+=HYi[2*PML_Thickness+10+x2]*(d-x1)*(d-x0)*(d-x3)/((x2-x1)*(x2-x0)*(x2-x3));
		hi+=HYi[2*PML_Thickness+10+x3]*(d-x1)*(d-x2)*(d-x0)/((x3-x1)*(x3-x2)*(x3-x0));
	};

	switch(which) 
	{
		case 1://hx
			return hi*sin(phi*pii/180);
		case 2://hy
			return -hi*cos(phi*pii/180);
		case 3://ez
			return ei;
		default:
			return 0;
	};
}