#include<stdio.h>
#include<math.h>
#include<fstream.h>


const double	pii=3.1415926535898;
const double	epsn0=1e-9/(36.*pii);	//Permit in free space
const double	miu0=1e-7*4.*pii;		//Permeab in free space
const double	Ltspd=3.0e+8;			//Light Speed(m/s)
const double	yta0=120*pii;			//wave impendance in free space(Omn)
const double	relpt=1.0;				//epsn/epsn0
const double	relpb=1.0;				//miu/miu0
const double	epsn=relpt*epsn0;
const double	miu=relpb*miu0;
//const double	fr=10.0e6;	//frequency(Hz)
//const double	dz=0.125;//2.0e-4;75e-6;//
//const double	dt=20.0*0.208e-9;//0.5e-12;0.125e-12;//


const double	dz=75e-6;
const double	dt=10.0*0.125e-12;
const int   MZ=1800;		//z axis maxium length(in deltas)		

const double xisu=(2.0*sqrt(epsn*miu)*dz/dt)*(2.0*sqrt(epsn*miu)*dz/dt);
const int   Nth=5;			//the depth of PML 
const int	K1=Nth;
const int	K2=MZ-Nth;

double etx1[Nth+1],ety1[Nth+1],etz1[Nth+1],etx2[Nth+1],ety2[Nth+1],etz2[Nth+1];
double mtx1[Nth+1],mty1[Nth+1],mtz1[Nth+1],mtx2[Nth+1],mty2[Nth+1],mtz2[Nth+1];

const double		tem1=dt/epsn;
const double		tem2=dt/miu;
const double		tem3=miu/epsn;
const double		g=1.44;	//the ratio of the geometric progression
const double		R0=0.00001;	//the reflect coefficent on the ObjectSpace-PML interface
const double		ct0=epsn*Ltspd*log(R0)/(2.0*(pow(g,(double)Nth)-1.0));
const double		ce0=ct0*(1.0-sqrt(g));
const double		cm0=tem3*ce0;

const int npml=3;
const double sigmmax=0.0-log(R0)*(npml+1)/2.0*epsn0*Ltspd/Nth;
const double sigma0=sigmmax/((npml+1)*pow(2,npml+1)*pow(Nth,npml));


//scattering.h
const int thick=50;			//the size of the cube
const int SN=3;
const int low=MZ/2-thick;
const int high=MZ/2+thick;
const int kk0=low-SN;		//
const int kk1=high+SN;		//total-field/scattered-field boundary
int index;
double EXi[MZ+1],HYi[MZ+1];	// plane wave incident along +z direction 


double inline conduct(int L){
	if (L<=0) return sigma0;
	else return sigma0*(pow(2.0*L+1.0,npml+1)-pow(2.0*L-1.0,npml+1));
}

double inline mconduct(int  L){
	if (L<=0) return tem3*sigma0;
	else return tem3*sigma0*(pow(2.0*(L+0.5)+1.0,npml+1)-pow(2.0*(L+0.5)-1.0,npml+1));
}

double sch(int t, const double aa=2.0){
	double fr=1e10;
	double enevop=1.0,w=2.0*pii*fr;
	int    T_max=(int)(0.5+aa/fr/dt);
	if (t<0) enevop=0.0;
	else if ((0<=t)&&(t<=T_max)) enevop=(0.5*(1.0-cos(w*t*dt*0.5/aa)));
		else enevop=1.0;		
	return(enevop*sin(w*t*dt));
}
double sins(int t){
	double fr=20.0e9;

	return(sin(2.0*pii*fr*t*dt));
}

double scp(int t){
	int tao=15;
	if(0<=t&&t<=10*tao) return (t-5.0*tao)*exp(0.0-(t-5.0*tao)*(t-5.0*tao)/(2.0*tao*tao))*100/11.9/1.8;
	else return(0.0);
}

void Zeroinit(void){
	int i,index;

	for(i=0;i<Nth+1;i++) etz1[i]=exp(0.0-tem1*conduct(i)/dz);
	for(i=0;i<Nth+1;i++) mtz1[i]=exp(0.0-tem2*mconduct(i)/dz);
	for(i=0;i<Nth+1;i++) etz2[i]=(1.0-exp(0.0-tem1*conduct(i)/dz))/conduct(i)*dz;
	for(i=0;i<Nth+1;i++) mtz2[i]=(1.0-exp(0.0-tem2*mconduct(i)/dz))/mconduct(i)*dz;
	for(index=0;index<=MZ;index++) EXi[index]=0.0;
	for(index=0;index<=MZ;index++) HYi[index]=0.0;

 }
int tri_matrix(int n, double *a, double *b, double *c, double *f)
{
	int kk;
	double w;
	f[0]=f[0]/b[0];
	w=b[0];
	for (kk=1; kk<=n-1;kk++)
	{
		b[kk-1]=c[kk-1]/w;
		w=b[kk]-a[kk-1]*b[kk-1];
		f[kk]=(f[kk]-a[kk-1]*f[kk-1])/w;
	}
	for(kk=n-1; kk>=1; kk--) f[kk-1]=f[kk-1]-b[kk-1]*f[kk];
	return(1);
}