#include"fdtd.h"
#include"plasma.h"


//#define source sch(n)//升余弦函數
#define source scp(n)//高斯函數

void main()
{	
	int k,n,N;
		/*double *psi1=new double[MZ+1];
		double *psi2=new double[MZ+1];
		double *ex=new double[MZ+1];
	    double *hy=new double[MZ];*/

	complex<double> *ex=new complex<double>[MZ+1];
	complex<double> *hy=new complex<double>[MZ];
	complex<double> *psi=new complex<double>[MZ+1];
//	complex<double> *psi2=new complex<double>[MZ+1];


	for(k=0;k<MZ+1;k++)	ex[k]=0.0;
	for(k=0;k<MZ;k++)	hy[k]=0.0;
	for(k=0;k<MZ+1;k++)	psi[k]=0.0;
//	for(k=0;k<MZ+1;k++)	psi2[k]=0.0;

	ofstream fref,fdef,fsta;
	fref.open("ref.dat");
	fdef.open("def.dat");
	fsta.open("sta1.dat");

	cout<<"Input N:"<<endl;
	cin>>N;

	Zeroinit();

	for(n=1;n<=N;n++)
	{
		cout<<n<<endl;
	 
		// Incident  H_field iteration at nth time step!		
		for(k=0;k<K1;k++) HYi[k]=mtz1[K1-k-1]*HYi[k]+(mtz2[K1-k-1]/dz)*(EXi[k]-EXi[k+1]);//PML
		for(k=K1;k<K2;k++)  HYi[k]+=(tem2/dz)*(EXi[k]-EXi[k+1]);
		for(k=K2;k<MZ;k++) HYi[k]=mtz1[k-K2]*HYi[k]+(mtz2[k-K2]/dz)*(EXi[k]-EXi[k+1]);//PML
//		HYi[K1+1]=sins(n)/yta0;

		//FDTD H_iteration
		for(k=0;k<K1;k++) hy[k]=mtz1[K1-k-1]*hy[k]+(mtz2[K1-k-1]/dz)*(ex[k]-ex[k+1]);//PML
		hy[kk0-1]+=-(tem2/dz)*(ex[kk0]-EXi[kk0]-ex[kk0-1]);//Source
		hy[kk1]+=-(tem2/dz)*(ex[kk1+1]+EXi[kk1]-ex[kk1]);//Source
		for(k=K1;k<MZ;k++)
			if(k!=kk0-1&&k!=kk1)
				hy[k]+=-(tem2/dz)*(ex[k+1]-ex[k]);
		for(k=K2;k<MZ;k++) hy[k]=mtz1[k-K2]*hy[k]+(mtz2[k-K2]/dz)*(ex[k]-ex[k+1]);//PML


		// Incident field E_iteration at nth time step!		
		for(k=1;k<=K1;k++)  EXi[k]=etz1[K1-k]*EXi[k]+(etz2[K1-k]/dz)*(HYi[k-1]-HYi[k]);//PML
		for(k=K1+1;k<K2;k++)  EXi[k]+=(tem1/dz)*(HYi[k-1]-HYi[k]);
		for(k=K2;k<MZ;k++)  EXi[k]=etz1[k-K2]*EXi[k]+(etz2[k-K2]/dz)*(HYi[k-1]-HYi[k]); //PML
		EXi[K1+1]=source;    //K1+1 or K1 

		//FDTD E_iteration
		for(k=1;k<=K1;k++)  ex[k]=etz1[K1-k]*ex[k]+(etz2[K1-k]/dz)*(hy[k-1]-hy[k]);//PML
		ex[kk0]+=-(tem1/dz)*(hy[kk0]-hy[kk0-1]-HYi[kk0-1]);//Source
		ex[kk1]+=-(tem1/dz)*(hy[kk1]+HYi[kk1]-hy[kk1-1]);//Source
		for(k=K1+1;k<MZ;k++)
			if(k!=kk0&&k!=kk1)   //and
				if(k<MZ/2-thick||k>MZ/2+thick)  //or
					ex[k]+=-(tem1/dz)*(hy[k]-hy[k-1]);
				else{
					complex<double> temp=ex[k];
					ex[k]=((1.0-xi(0,k))*ex[k]+psi[k]-(tem1/dz)*(hy[k]-hy[k-1]))/(1.0+chi(0,k)-xi(0,k));
					psi[k]=(dChi0(k)-dXi0(k))*ex[k]+dXi0(k)*temp+exp(-1.0*(CollisionalFreq-j*wb)*dt)*psi[k];
					//psi2[k]=(dChi20(k)-dXi20(k))*ex[k]+dXi10(k)*temp+exp(-CollisionalFreq*dt)*psi2[k];
				   	};
		///////////
    //	ex[MZ/2+thick]=0.0;							
	//	EXi[MZ/2+thick]=0.0;
		///////////
		for(k=K2;k<MZ;k++)  ex[k]=etz1[k-K2]*ex[k]+(etz2[k-K2]/dz)*(hy[k-1]-hy[k]); //PML
	   if(n<2500)
		fref<<real(ex[kk0-2])<<endl;//反射
	   else
		   fref<<0.0<<endl;
		fdef<<real(ex[kk1-2])<<endl;//折射
		if(n<2500)
			fsta<<EXi[K1+11]<<endl;
          //  fsta<<ex[K1+11]<<endl;
		else
			fsta<<0.0<<endl;

	/*	if(n==1300)
			for(k=K1;k<K2;k++)
				fdef<<real(ex[k])<<endl;*/

	};

	fref.close();
	fdef.close();
	fsta.close();
}