timeupdate.c

program FDTD for the photonic crystal structure

#include "./pFDTD.h"

void field_zero(float ***name);

static float tus;

void normal_E_field_update(int i, int j, int k);  //real field
void metal_E_field_update(int i, int j, int k);
void pml_E_field_update(int i, int j, int k);
void pcm_E_field_update(int i, int j, int k);
void normal_H_field_update(int i, int j, int k);
void pml_H_field_update(int i, int j, int k);
void pcm_H_field_update(int i, int j, int k);

void normal_iE_field_update(int i, int j, int k);  //imaginary field
void metal_iE_field_update(int i, int j, int k);
void pml_iE_field_update(int i, int j, int k);
void pcm_iE_field_update(int i, int j, int k);
void normal_iH_field_update(int i, int j, int k);
void pml_iH_field_update(int i, int j, int k);
void pcm_iH_field_update(int i, int j, int k);

void E_field_periodic_boundary_update_x();
void E_field_periodic_boundary_update_y();
void H_field_periodic_boundary_update_x();
void H_field_periodic_boundary_update_y();
void E_field_Gamma_boundary_update_x();
void E_field_Gamma_boundary_update_y();
void H_field_Gamma_boundary_update_x();
void H_field_Gamma_boundary_update_y();

void Ez_parity_boundary_update();
void Ez_parity_iboundary_update();
void Hz_parity_boundary_update();
void Hz_parity_iboundary_update();

void coefficient()
{
	int i,j,k;

	for(i=0;i<misize;i++)
	for(j=0;j<mjsize;j++)
	for(k=0;k<mksize;k++)
	{
		if(position[i][j][k])
		{  //Do nothing!
		}
		
		else
		{
			aax[i]=(2*eo*kx-dt*sigmax(i))/(2*eo*kx+dt*sigmax(i));
			aay[j]=(2*eo*ky-dt*sigmay(j))/(2*eo*ky+dt*sigmay(j));
			aaz[k]=(2*eo*kz-dt*sigmaz(k))/(2*eo*kz+dt*sigmaz(k));
			
			bbx[i]=(2*eo*dt)/(2*eo*kx+dt*sigmax(i));
			bby[j]=(2*eo*dt)/(2*eo*ky+dt*sigmay(j));
			bbz[k]=(2*eo*dt)/(2*eo*kz+dt*sigmaz(k));

			ccx[i]=(2*eo*kx-dt*sigmax(i))/(2*eo*kx+dt*sigmax(i));
			ccy[j]=(2*eo*ky-dt*sigmay(j))/(2*eo*ky+dt*sigmay(j));
			ccz[k]=(2*eo*kz-dt*sigmaz(k))/(2*eo*kz+dt*sigmaz(k));

			ddx[i][j][k]=2/(2*eo*kz+dt*sigmaz(k))/epsilonx[i][j][k];
			ddy[i][j][k]=2/(2*eo*kx+dt*sigmax(i))/epsilony[i][j][k];
			ddz[i][j][k]=2/(2*eo*ky+dt*sigmay(j))/epsilonz[i][j][k];

			eex[i]=kx+dt*sigmax(i+0.5)/2/eo;
			eey[j]=ky+dt*sigmay(j+0.5)/2/eo;
			eez[k]=kz+dt*sigmaz(k+0.5)/2/eo;

			ffx[i]=kx-dt*sigmax(i+0.5)/2/eo;
			ffy[j]=ky-dt*sigmay(j+0.5)/2/eo;
			ffz[k]=kz-dt*sigmaz(k+0.5)/2/eo;

			ggx[i]=(2*eo*kx-dt*sigmax(i+0.5))/(2*eo*kx+dt*sigmax(i+0.5));
			ggy[j]=(2*eo*ky-dt*sigmay(j+0.5))/(2*eo*ky+dt*sigmay(j+0.5));
			ggz[k]=(2*eo*kz-dt*sigmaz(k+0.5))/(2*eo*kz+dt*sigmaz(k+0.5));
			
			hhx[i]=(2*eo*dt)/(2*eo*kx+dt*sigmax(i+0.5));
			hhy[j]=(2*eo*dt)/(2*eo*ky+dt*sigmay(j+0.5));
			hhz[k]=(2*eo*dt)/(2*eo*kz+dt*sigmaz(k+0.5));

			iix[i]=(2*eo*kx-dt*sigmax(i+0.5))/(2*eo*kx+dt*sigmax(i+0.5));
			iiy[j]=(2*eo*ky-dt*sigmay(j+0.5))/(2*eo*ky+dt*sigmay(j+0.5));
			iiz[k]=(2*eo*kz-dt*sigmaz(k+0.5))/(2*eo*kz+dt*sigmaz(k+0.5));

			jjx[i]=(2*eo)/(2*eo*kx+dt*sigmax(i+0.5))/uo/ups;
			jjy[j]=(2*eo)/(2*eo*ky+dt*sigmay(j+0.5))/uo/ups;
			jjz[k]=(2*eo)/(2*eo*kz+dt*sigmaz(k+0.5))/uo/ups;

			kkx[i]=kx+dt*sigmax(i)/2/eo;
			kky[j]=ky+dt*sigmay(j)/2/eo;
			kkz[k]=kz+dt*sigmaz(k)/2/eo;

			llx[i]=kx-dt*sigmax(i)/2/eo;
			lly[j]=ky-dt*sigmay(j)/2/eo;
			llz[k]=kz-dt*sigmaz(k)/2/eo;
		}
	}
			tus=dt/uo/ups;
	
	printf("coefficient...ok\n");
}

void propagate()
{
	int i,j,k;

	//******* real E-field ******//
	for(i=1;i<pisize;i++)
	for(j=1;j<pjsize;j++)
	for(k=1;k<pksize;k++)
	{
		if(position[i][j][k]==1 && mepsilon[i][j][k]==0.0) // non-metal
			normal_E_field_update(i, j, k);
		else if(position[i][j][k]==1 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal
			metal_E_field_update(i, j, k);
		else if(position[i][j][k]==0 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal
			metal_E_field_update(i, j, k);
		else if(position[i][j][k]==1 && mepsilon[i][j][k]==1000) // PCM
			pcm_E_field_update(i, j, k);
		else if(position[i][j][k]==0 && mepsilon[i][j][k]==0.0) // non-metal & PML region
			pml_E_field_update(i, j, k);
		else
			{}
	}

	//******* imag E-field ******//
	if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
	{
		for(i=1;i<pisize;i++)
		for(j=1;j<pjsize;j++)
		for(k=1;k<pksize;k++)
		{
		if(position[i][j][k]==1 && mepsilon[i][j][k]==0.0) // non-metal
			normal_iE_field_update(i, j, k);
		else if(position[i][j][k]==1 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal
			metal_iE_field_update(i, j, k);
		else if(position[i][j][k]==0 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal
			metal_iE_field_update(i, j, k);
		else if(position[i][j][k]==1 && mepsilon[i][j][k]==1000) // PCM
			pcm_iE_field_update(i, j, k);
		else if(position[i][j][k]==0 && mepsilon[i][j][k]==0.0) // non-metal & PML region
			pml_iE_field_update(i, j, k);
		else
			{}
		}
	}

	//******* E-field periodic condition ******//
	if(use_periodic_x == 1 && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		E_field_periodic_boundary_update_x();
	if(use_periodic_y == 1 && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		E_field_periodic_boundary_update_y();
	if(use_periodic_x == 1 && (wave_vector_x==0.0 && wave_vector_y==0.0))
		E_field_Gamma_boundary_update_x();
	if(use_periodic_y == 1 && (wave_vector_x==0.0 && wave_vector_y==0.0))
		E_field_Gamma_boundary_update_y();

	//******* Ez_parity real-field ******//
	Ez_parity_boundary_update();

	//******* Ez_parity imag-field ******//
	if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		Ez_parity_iboundary_update();

	//******* real H-field ******//
	for(i=1;i<pisize;i++)
	for(j=1;j<pjsize;j++)
	for(k=1;k<pksize;k++)
	{
		if(position[i][j][k]==1 && (mepsilon[i][j][k]!=1000))
			normal_H_field_update(i, j, k);
		else if(position[i][j][k]==0 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal & PML region
			normal_H_field_update(i, j, k);
		else if(position[i][j][k]==1 && mepsilon[i][j][k]==1000) // PCM
			pcm_H_field_update(i, j, k);
		else if(position[i][j][k]==0 && mepsilon[i][j][k]==0.0) // non-metal & PML
			pml_H_field_update(i, j, k);
		else
			{}
	}

	//******* imag H-field ******//
	if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
	{
		for(i=1;i<pisize;i++)
		for(j=1;j<pjsize;j++)
		for(k=1;k<pksize;k++)
		{
		if(position[i][j][k]==1 && (mepsilon[i][j][k]!=1000))
			normal_iH_field_update(i, j, k);
		else if(position[i][j][k]==0 && (mepsilon[i][j][k]!=0.0 && mepsilon[i][j][k]!=1000)) // metal & PML region
			normal_iH_field_update(i, j, k);
		else if(position[i][j][k]==1 && mepsilon[i][j][k]==1000) // PCM
			pcm_iH_field_update(i, j, k);
		else if(position[i][j][k]==0 && mepsilon[i][j][k]==0.0) // non-metal & PML
			pml_iH_field_update(i, j, k);
		else
			{}
		}
	}

	//******* H-field periodic condition ******//
	if(use_periodic_x == 1 && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		H_field_periodic_boundary_update_x();
	if(use_periodic_y == 1 && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		H_field_periodic_boundary_update_y();
	if(use_periodic_x == 1 && (wave_vector_x==0.0 && wave_vector_y==0.0))
		H_field_Gamma_boundary_update_x();
	if(use_periodic_y == 1 && (wave_vector_x==0.0 && wave_vector_y==0.0))
		H_field_Gamma_boundary_update_y();

	//******* Hz_parity real-field ******//
	Hz_parity_boundary_update();

	//******* Hz_parity imag-field ******//
	if((use_periodic_x == 1 || use_periodic_y == 1) && (wave_vector_x!=0.0 || wave_vector_y!=0.0))
		Hz_parity_iboundary_update();

	printf("%d\n",t);

}

void normal_E_field_update(int i, int j, int k)
{
	Ex[i][j][k] = Ex[i][j][k] + (dt/eo/epsilonx[i][j][k])*((Hz[i][j][k]-Hz[i][j-1][k])/ds_y - (Hy[i][j][k]-Hy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])));
	Ey[i][j][k] = Ey[i][j][k] + (dt/eo/epsilony[i][j][k])*((Hx[i][j][k]-Hx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (Hz[i][j][k]-Hz[i-1][j][k])/ds_x);
	Ez[i][j][k] = Ez[i][j][k] + (dt/eo/epsilonz[i][j][k])*((Hy[i][j][k]-Hy[i-1][j][k])/ds_x - (Hx[i][j][k]-Hx[i][j-1][k])/ds_y);
}

void metal_E_field_update(int i, int j, int k)
{
	float km, bm; 
	float Ex_temp, Ey_temp, Ez_temp; 

	km = (2-mgamma[i][j][k]*dt)/(2+mgamma[i][j][k]*dt);
	bm = momega[i][j][k]*momega[i][j][k]*eo*dt / (2+mgamma[i][j][k]*dt);

	Ex_temp = Ex[i][j][k];  // store E-field at FDTD time 'n' 
	Ey_temp = Ey[i][j][k];  //          "               
	Ez_temp = Ez[i][j][k];  //          "

	Ex[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*Ex[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((Hz[i][j][k]-Hz[i][j-1][k])/ds_y - (Hy[i][j][k]-Hy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - 0.5*(1+km)*Jx[i][j][k]);
	Ey[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*Ey[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((Hx[i][j][k]-Hx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (Hz[i][j][k]-Hz[i-1][j][k])/ds_x - 0.5*(1+km)*Jy[i][j][k]);
	Ez[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*Ez[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((Hy[i][j][k]-Hy[i-1][j][k])/ds_x - (Hx[i][j][k]-Hx[i][j-1][k])/ds_y - 0.5*(1+km)*Jz[i][j][k]);

	Jx[i][j][k] = km*Jx[i][j][k] + bm*(Ex[i][j][k] + Ex_temp);
	Jy[i][j][k] = km*Jy[i][j][k] + bm*(Ey[i][j][k] + Ey_temp);
	Jz[i][j][k] = km*Jz[i][j][k] + bm*(Ez[i][j][k] + Ez_temp);
}

void pml_E_field_update(int i, int j, int k)
{
	float temp;

	temp=Dx[i][j][k];
	Dx[i][j][k] = aay[j]*Dx[i][j][k] + bby[j]*((Hz[i][j][k]-Hz[i][j-1][k])/ds_y - (Hy[i][j][k]-Hy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])));
	Ex[i][j][k] = ccz[k]*Ex[i][j][k] + ddx[i][j][k]*(eex[i]*Dx[i][j][k]-ffx[i]*temp);
	temp=Dy[i][j][k];
	Dy[i][j][k] = aaz[k]*Dy[i][j][k] + bbz[k]*((Hx[i][j][k]-Hx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (Hz[i][j][k]-Hz[i-1][j][k])/ds_x);
	Ey[i][j][k] = ccx[i]*Ey[i][j][k] + ddy[i][j][k]*(eey[j]*Dy[i][j][k]-ffy[j]*temp);
	temp=Dz[i][j][k];
	Dz[i][j][k] = aax[i]*Dz[i][j][k] + bbx[i]*((Hy[i][j][k]-Hy[i-1][j][k])/ds_x - (Hx[i][j][k]-Hx[i][j-1][k])/ds_y);
	Ez[i][j][k] = ccy[j]*Ez[i][j][k] + ddz[i][j][k]*(eez[k]*Dz[i][j][k]-ffz[k]*temp);
}

void pcm_E_field_update(int i, int j, int k)
{
	Ex[i][j][k] = 0.0;
	Ey[i][j][k] = 0.0;
	Ez[i][j][k] = 0.0;
}

void normal_H_field_update(int i, int j, int k)
{
	Hx[i][j][k] = Hx[i][j][k] - tus*((Ez[i][j+1][k]-Ez[i][j][k])/ds_y - (Ey[i][j][k+1]-Ey[i][j][k])/ds_nz[k]);
	Hy[i][j][k] = Hy[i][j][k] - tus*((Ex[i][j][k+1]-Ex[i][j][k])/ds_nz[k] - (Ez[i+1][j][k]-Ez[i][j][k])/ds_x);
	Hz[i][j][k] = Hz[i][j][k] - tus*((Ey[i+1][j][k]-Ey[i][j][k])/ds_x - (Ex[i][j+1][k]-Ex[i][j][k])/ds_y);
}

void pml_H_field_update(int i, int j, int k)
{
	float temp;

	temp=Bx[i][j][k];
	Bx[i][j][k] = ggy[j]*Bx[i][j][k] - hhy[j]*((Ez[i][j+1][k]-Ez[i][j][k])/ds_y - (Ey[i][j][k+1]-Ey[i][j][k])/ds_nz[k]);
	Hx[i][j][k] = iiz[k]*Hx[i][j][k] + jjz[k]*(kkx[i]*Bx[i][j][k]-llx[i]*temp);
	temp=By[i][j][k];
	By[i][j][k] = ggz[k]*By[i][j][k] - hhz[k]*((Ex[i][j][k+1]-Ex[i][j][k])/ds_nz[k] - (Ez[i+1][j][k]-Ez[i][j][k])/ds_x);
	Hy[i][j][k] = iix[i]*Hy[i][j][k] + jjx[i]*(kky[j]*By[i][j][k]-lly[j]*temp);
	temp=Bz[i][j][k];
	Bz[i][j][k] = ggx[i]*Bz[i][j][k] - hhx[i]*((Ey[i+1][j][k]-Ey[i][j][k])/ds_x - (Ex[i][j+1][k]-Ex[i][j][k])/ds_y);
	Hz[i][j][k] = iiy[j]*Hz[i][j][k] + jjy[j]*(kkz[k]*Bz[i][j][k]-llz[k]*temp);
}

void pcm_H_field_update(int i, int j, int k)
{
	Hx[i][j][k] = 0.0;
	Hy[i][j][k] = 0.0;
	Hz[i][j][k] = 0.0;
}

void normal_iE_field_update(int i, int j, int k)
{
	iEx[i][j][k] = iEx[i][j][k] + (dt/eo/epsilonx[i][j][k])*((iHz[i][j][k]-iHz[i][j-1][k])/ds_y - (iHy[i][j][k]-iHy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])));
	iEy[i][j][k] = iEy[i][j][k] + (dt/eo/epsilony[i][j][k])*((iHx[i][j][k]-iHx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (iHz[i][j][k]-iHz[i-1][j][k])/ds_x);
	iEz[i][j][k] = iEz[i][j][k] + (dt/eo/epsilonz[i][j][k])*((iHy[i][j][k]-iHy[i-1][j][k])/ds_x - (iHx[i][j][k]-iHx[i][j-1][k])/ds_y);
}

void metal_iE_field_update(int i, int j, int k)
{
	float km, bm; 
	float iEx_temp, iEy_temp, iEz_temp; 

	km = (2-mgamma[i][j][k]*dt)/(2+mgamma[i][j][k]*dt);
	bm = momega[i][j][k]*momega[i][j][k]*eo*dt / (2+mgamma[i][j][k]*dt);

	iEx_temp = iEx[i][j][k];  // store E-field at FDTD time 'n' 
	iEy_temp = iEy[i][j][k];  //          "               
	iEz_temp = iEz[i][j][k];  //          "

	iEx[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*iEx[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((iHz[i][j][k]-iHz[i][j-1][k])/ds_y - (iHy[i][j][k]-iHy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - 0.5*(1+km)*iJx[i][j][k]);
	iEy[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*iEy[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((iHx[i][j][k]-iHx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (iHz[i][j][k]-iHz[i-1][j][k])/ds_x - 0.5*(1+km)*iJy[i][j][k]);
	iEz[i][j][k] = ((2*eo*mepsilon[i][j][k]-dt*bm)/(2*eo*mepsilon[i][j][k]+dt*bm))*iEz[i][j][k] + (2*dt/(2*eo*mepsilon[i][j][k]+dt*bm))*((iHy[i][j][k]-iHy[i-1][j][k])/ds_x - (iHx[i][j][k]-iHx[i][j-1][k])/ds_y - 0.5*(1+km)*iJz[i][j][k]);

	iJx[i][j][k] = km*iJx[i][j][k] + bm*(iEx[i][j][k] + iEx_temp);
	iJy[i][j][k] = km*iJy[i][j][k] + bm*(iEy[i][j][k] + iEy_temp);
	iJz[i][j][k] = km*iJz[i][j][k] + bm*(iEz[i][j][k] + iEz_temp);
}

void pml_iE_field_update(int i, int j, int k)
{
	float temp;

	temp=iDx[i][j][k];
	iDx[i][j][k] = aay[j]*iDx[i][j][k] + bby[j]*((iHz[i][j][k]-iHz[i][j-1][k])/ds_y - (iHy[i][j][k]-iHy[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])));
	iEx[i][j][k] = ccz[k]*iEx[i][j][k] + ddx[i][j][k]*(eex[i]*iDx[i][j][k]-ffx[i]*temp);
	temp=iDy[i][j][k];
	iDy[i][j][k] = aaz[k]*iDy[i][j][k] + bbz[k]*((iHx[i][j][k]-iHx[i][j][k-1])*(2/(ds_nz[k-1]+ds_nz[k])) - (iHz[i][j][k]-iHz[i-1][j][k])/ds_x);
	iEy[i][j][k] = ccx[i]*iEy[i][j][k] + ddy[i][j][k]*(eey[j]*iDy[i][j][k]-ffy[j]*temp);
	temp=iDz[i][j][k];
	iDz[i][j][k] = aax[i]*iDz[i][j][k] + bbx[i]*((iHy[i][j][k]-iHy[i-1][j][k])/ds_x - (iHx[i][j][k]-iHx[i][j-1][k])/ds_y);
	iEz[i][j][k] = ccy[j]*iEz[i][j][k] + ddz[i][j][k]*(eez[k]*iDz[i][j][k]-ffz[k]*temp);
}

void pcm_iE_field_update(int i, int j, int k)
{
	iEx[i][j][k] = 0.0;
	iEy[i][j][k] = 0.0;
	iEz[i][j][k] = 0.0;
}

void normal_iH_field_update(int i, int j, int k)
{
	iHx[i][j][k] = iHx[i][j][k] - tus*((iEz[i][j+1][k]-iEz[i][j][k])/ds_y - (iEy[i][j][k+1]-iEy[i][j][k])/ds_nz[k]);
	iHy[i][j][k] = iHy[i][j][k] - tus*((iEx[i][j][k+1]-iEx[i][j][k])/ds_nz[k] - (iEz[i+1][j][k]-iEz[i][j][k])/ds_x);
	iHz[i][j][k] = iHz[i][j][k] - tus*((iEy[i+1][j][k]-iEy[i][j][k])/ds_x - (iEx[i][j+1][k]-iEx[i][j][k])/ds_y);
}

void pml_iH_field_update(int i, int j, int k)
{
	float temp;

	temp=iBx[i][j][k];
	iBx[i][j][k] = ggy[j]*iBx[i][j][k] - hhy[j]*((iEz[i][j+1][k]-iEz[i][j][k])/ds_y - (iEy[i][j][k+1]-iEy[i][j][k])/ds_nz[k]);
	iHx[i][j][k] = iiz[k]*iHx[i][j][k] + jjz[k]*(kkx[i]*iBx[i][j][k]-llx[i]*temp);
	temp=iBy[i][j][k];
	iBy[i][j][k] = ggz[k]*iBy[i][j][k] - hhz[k]*((iEx[i][j][k+1]-iEx[i][j][k])/ds_nz[k] - (iEz[i+1][j][k]-iEz[i][j][k])/ds_x);
	iHy[i][j][k] = iix[i]*iHy[i][j][k] + jjx[i]*(kky[j]*iBy[i][j][k]-lly[j]*temp);
	temp=iBz[i][j][k];
	iBz[i][j][k] = ggx[i]*iBz[i][j][k] - hhx[i]*((iEy[i+1][j][k]-iEy[i][j][k])/ds_x - (iEx[i][j+1][k]-iEx[i][j][k])/ds_y);
	iHz[i][j][k] = iiy[j]*iHz[i][j][k] + jjy[j]*(kkz[k]*iBz[i][j][k]-llz[k]*temp);
}

void pcm_iH_field_update(int i, int j, int k)
{
	iHx[i][j][k] = 0.0;
	iHy[i][j][k] = 0.0;
	iHz[i][j][k] = 0.0;
}

void Ez_parity_boundary_update()
{
	int i,j,k;

	if(xparity==1)  //xparity is for Hz_parity
	{
		for(j=1;j<pjsize;j++)
		for(k=1;k<pksize;k++)