double  **D2dmatrix(int row,int col){//申請二維數組：0～row,0～col
	int ii,jj;
	double **temp;
	temp=(double **)malloc(long(row)*sizeof(double *));
	if(!temp){
		cout<<"Allocation Failed!"<<endl;
		exit(1);
	}
	for(ii=0;ii<row;ii++){
		temp[ii]=(double *)malloc(long(col)*sizeof(double));
		if(!temp[ii]){
			cout<<"Allocation Failed!"<<endl;
			exit(1);
		}
	}
	for(ii=0;ii<row;ii++)
		for(jj=0;jj<col;jj++)
			temp[ii][jj]=0.0;
	return temp;
}

class PML_2D_Arr_ez{
private:
	int XM,YM,Nx1,Nx2,Ny1,Ny2;
	double *pi1,*pi2,*pj1,*pj2,pzero;
public:
	PML_2D_Arr_ez(int i,int I1,int I2,int j,int J1,int J2);
	~PML_2D_Arr_ez();
	double& operator()(int i,int j);
};

PML_2D_Arr_ez::PML_2D_Arr_ez(int i,int I1,int I2,int j,int J1,int J2){
	int index;
	XM=i;YM=j;Nx1=I1;Nx2=I2;Ny1=J1;Ny2=J2;
	pi1=new double[(Nx1+1)*(YM+1)];VERIFY(pi1!=NULL);
	pi2=new double[(Nx2+1)*(YM+1)];VERIFY(pi2!=NULL);
	pj1=new double[(XM-Nx1-Nx2-1)*(Ny1+1)];VERIFY(pj1!=NULL);
	pj2=new double[(XM-Nx1-Nx2-1)*(Ny2+1)];VERIFY(pj2!=NULL);
	
	for(index=0;index<((Nx1+1)*(YM+1));index++) pi1[index]=0.0;
	for(index=0;index<((Nx2+1)*(YM+1));index++) pi2[index]=0.0;
	for(index=0;index<((XM-Nx1-Nx2-1)*(Ny1+1));index++) pj1[index]=0.0;
	for(index=0;index<((XM-Nx1-Nx2-1)*(Ny2+1));index++) pj2[index]=0.0;
	pzero=0.0;
}

PML_2D_Arr_ez::~PML_2D_Arr_ez(){
	delete pi1;
	delete pi2;
	delete pj1;
	delete pj2;
}

double& PML_2D_Arr_ez::operator ()(int i,int j ){
		VERIFY((i>=0)&&(j>=0));
		VERIFY((i<=XM)&&(j<=YM));
		if(i<=Nx1) return pi1[i*(YM+1)+j];
		else if(i>=(XM-Nx2)) return pi2[(i-XM+Nx2)*(YM+1)+j];
		else if(j<=Ny1) return pj1[(i-Nx1-1)*(Ny1+1)+j];
		else if(j>=(YM-Ny2)) return pj2[(i-Nx1-1)*(Ny2+1)+j-YM+Ny2];
		else{
		  cout<<i<<" "<<j<<" "<<"ez"<<"illegal access!"<<endl;
		  return pzero;
		}
}

class PML_2D_Arr_hz{
private:
	int XM,YM,Nx1,Nx2,Ny1,Ny2;
	double *pi1,*pi2,*pj1,*pj2,pzero;
public:
	PML_2D_Arr_hz(int i,int I1,int I2,int j,int J1,int J2);
	~PML_2D_Arr_hz();
	double& operator()(int i,int j);
};

PML_2D_Arr_hz::PML_2D_Arr_hz(int i,int I1,int I2,int j,int J1,int J2){
	int index;
	XM=i;YM=j;Nx1=I1;Nx2=I2;Ny1=J1;Ny2=J2;
	pi1=new double[Nx1*YM];VERIFY(pi1!=NULL);
	pi2=new double[Nx2*YM];VERIFY(pi2!=NULL);
	pj1=new double[(XM-Nx1-Nx2)*Ny1];VERIFY(pj1!=NULL);
	pj2=new double[(XM-Nx1-Nx2)*Ny2];VERIFY(pj2!=NULL);
	
	for(index=0;index<(Nx1*YM);index++) pi1[index]=0.0;
	for(index=0;index<(Nx2*YM);index++) pi2[index]=0.0;
	for(index=0;index<((XM-Nx1-Nx2)*Ny1);index++) pj1[index]=0.0;
	for(index=0;index<((XM-Nx1-Nx2)*Ny2);index++) pj2[index]=0.0;
	pzero=0.0;
}

PML_2D_Arr_hz::~PML_2D_Arr_hz(){
	delete pi1;
	delete pi2;
	delete pj1;
	delete pj2;
}

double& PML_2D_Arr_hz::operator ()(int i,int j){
		VERIFY((i>=0)&&(j>=0));
		VERIFY((i<XM)&&(j<YM));
		if(i<Nx1) return pi1[i*YM+j];
		else if(i>=(XM-Nx2)) return pi2[(i-XM+Nx2)*YM+j];
		else if(j<Ny1) return pj1[(i-Nx1)*Ny1+j];
		else if(j>=(YM-Ny2)) return pj2[(i-Nx1)*Ny2+j-YM+Ny2];
		else{
		  cout<<i<<" "<<j<<" "<<"hz"<<"illegal access!"<<endl;
		  return pzero;
		}
}