//////////////////////////////////////////////////////////////////////////////////////
//																					//
//	This header definies the constants and functions in 1d_FDTD.                	//
//																					//
////////////////////////////////////////////////////////////////////////////////////// 

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

//general.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)
						
//FDTD.h
/*
const double	wlth=Ltspd/fr;						//wavelength(m)
const double	GridNumInaWavelenth=1000.0;
const double	delta=wlth/GridNumInaWavelenth;		//space step length(m)
const double	dz=delta;							//z axis spacial step
const double	CouNum=0.5;	 						//courant number equal to Ltspd*dt/delta
const double	dt=CouNum*delta/Ltspd;				//time step length(s)
*/
const double	dz=75e-6;
const double	dt=0.125e-12;
const int   MZ=800;		//z axis maxium length(in deltas)		
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=100;			//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];	// 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;
//	double dtt=10.0*dt;
	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;

 }