void UpdatePMLEz(void){  int i,j,k;  PRECISION lowerX1, lowerX2, lowerY1, lowerY2;  PRECISION epsdt = EPSILON_0 / dt;  PRECISION commonFactor;  PRECISION currentEpsilon;  PRECISION currentSigma;  /* speed up variables */  int ny_plus_2pmlWidth = ny + 2*pmlWidth;  int nz_plus_2pmlWidth = nz + 2*pmlWidth;   int ny_plus_pmlWidth  = ny + pmlWidth;  int nz_plus_pmlWidth  = nz + pmlWidth;  int pmlWidth_minus_1  = pmlWidth-1;  int nx_minus_1        = nx-1;  int ny_minus_1        = ny-1;  /* BottomYZ */  for(i=1; i<pmlWidth; i++) /* 1 ! */    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */      for(k=1; k<nz_plus_2pmlWidth; k++)	{	  currentEpsilon = materialConstants[pmlMaterialBottomYZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialBottomYZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomYZ[i][j][k]); 	  pmlEZXBottomYZ[i][j][k] = pmlEZXBottomYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomYZ[i][j][k])) / commonFactor +	    (pmlHYZBottomYZ[i][j][k] - pmlHYZBottomYZ[i-1][j][k] +	     pmlHYXBottomYZ[i][j][k] - pmlHYXBottomYZ[i-1][j][k]) /	    (dx * pmlSXBottomYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomYZ[i][j][k] * pmlComplexEZXBottomYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYBottomYZ[i][j][k]); 	  pmlEZYBottomYZ[i][j][k] = pmlEZYBottomYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomYZ[i][j][k])) / commonFactor -	    (pmlHXYBottomYZ[i][j][k] - pmlHXYBottomYZ[i][j-1][k] +	     pmlHXZBottomYZ[i][j][k] - pmlHXZBottomYZ[i][j-1][k]) /	    (dy * pmlSYBottomYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYBottomYZ[i][j][k] * pmlComplexEZYBottomYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* TopYZ */  for(i=0; i<pmlWidth; i++)    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */      for(k=1; k<nz_plus_2pmlWidth; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHYZTopYZ[i-1][j][k];	      lowerX2 = pmlHYXTopYZ[i-1][j][k];	    }	  else	    {	      if (k < pmlWidth)		{		  lowerX1 = pmlHYZBottomXY[ nx_minus_1][j][k];		  lowerX2 = pmlHYXBottomXY[ nx_minus_1][j][k];		}	      else if (k < nz_plus_pmlWidth)		{		  if (j < pmlWidth)		    {		      lowerX1 = pmlHYZBottomXZ[ nx_minus_1][j][k-pmlWidth];		      lowerX2 = pmlHYXBottomXZ[ nx_minus_1][j][k-pmlWidth];		    }		  else if (j< ny_plus_pmlWidth)		    {		      lowerX1 = hy[ nx_minus_1][j-pmlWidth][k-pmlWidth];		      lowerX2 = 0.0;		    }		  else		    {		      lowerX1 = pmlHYZTopXZ[ nx_minus_1][j-ny_plus_pmlWidth][k-pmlWidth];		      lowerX2 = pmlHYXTopXZ[ nx_minus_1][j-ny_plus_pmlWidth][k-pmlWidth];		    }		}	      else		{		  lowerX1 = pmlHYZTopXY[ nx_minus_1][j][k-nz_plus_pmlWidth];		  lowerX2 = pmlHYXTopXY[ nx_minus_1][j][k-nz_plus_pmlWidth];		}	    }	  	  currentEpsilon = materialConstants[pmlMaterialTopYZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopYZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopYZ[i][j][k]); 	  pmlEZXTopYZ[i][j][k] = pmlEZXTopYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopYZ[i][j][k])) / commonFactor +	    (pmlHYZTopYZ[i][j][k] - lowerX1 +	     pmlHYXTopYZ[i][j][k] - lowerX2) /	    (dx * pmlSXTopYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopYZ[i][j][k] * pmlComplexEZXTopYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYTopYZ[i][j][k]); 	  pmlEZYTopYZ[i][j][k] = pmlEZYTopYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopYZ[i][j][k])) / commonFactor -	    (pmlHXYTopYZ[i][j][k] - pmlHXYTopYZ[i][j-1][k] +	     pmlHXZTopYZ[i][j][k] - pmlHXZTopYZ[i][j-1][k]) /	    (dy * pmlSYTopYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYTopYZ[i][j][k] * pmlComplexEZYTopYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Bottom XY */  for(i=0; i<nx; i++)    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */      for(k=1; k<pmlWidth; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHYZBottomXY[i-1][j][k];	      lowerX2 = pmlHYXBottomXY[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHYZBottomYZ[pmlWidth_minus_1][j][k];	      lowerX2 = pmlHYXBottomYZ[pmlWidth_minus_1][j][k];	    }	  	  currentEpsilon = materialConstants[pmlMaterialBottomXY[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialBottomXY[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomXY[i][j][k]); 	  pmlEZXBottomXY[i][j][k] = pmlEZXBottomXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomXY[i][j][k])) / commonFactor +	    (pmlHYZBottomXY[i][j][k] - lowerX1 +	     pmlHYXBottomXY[i][j][k] - lowerX2) /	    (dx * pmlSXBottomXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomXY[i][j][k] * pmlComplexEZXBottomXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYBottomXY[i][j][k]); 	  pmlEZYBottomXY[i][j][k] = pmlEZYBottomXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomXY[i][j][k])) / commonFactor -	    (pmlHXYBottomXY[i][j][k] - pmlHXYBottomXY[i][j-1][k] +	     pmlHXZBottomXY[i][j][k] - pmlHXZBottomXY[i][j-1][k]) /	    (dy * pmlSYBottomXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYBottomXY[i][j][k] * pmlComplexEZYBottomXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Top XY */  for(i=0; i<nx; i++)    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */      for(k=0; k<pmlWidth; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHYZTopXY[i-1][j][k];	      lowerX2 = pmlHYXTopXY[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHYZBottomYZ[pmlWidth_minus_1][j][k+nz_plus_pmlWidth];	      lowerX2 = pmlHYXBottomYZ[pmlWidth_minus_1][j][k+nz_plus_pmlWidth];	    }	  currentEpsilon = materialConstants[pmlMaterialTopXY[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopXY[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopXY[i][j][k]); 	  pmlEZXTopXY[i][j][k] = pmlEZXTopXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopXY[i][j][k])) / commonFactor +	    (pmlHYZTopXY[i][j][k] - lowerX1 +	     pmlHYXTopXY[i][j][k] - lowerX2) /	    (dx * pmlSXTopXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopXY[i][j][k] * pmlComplexEZXTopXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYTopXY[i][j][k]); 	  pmlEZYTopXY[i][j][k] = pmlEZYTopXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopXY[i][j][k])) / commonFactor -	    (pmlHXYTopXY[i][j][k] - pmlHXYTopXY[i][j-1][k] +	     pmlHXZTopXY[i][j][k] - pmlHXZTopXY[i][j-1][k]) /	    (dy * pmlSYTopXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYTopXY[i][j][k] * pmlComplexEZYTopXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Bottom XZ */    for(i=0; i<nx; i++)     for(j=1; j<pmlWidth; j++) /* 1 ! */      for(k=0; k<nz; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHYZBottomXZ[i-1][j][k];	      lowerX2 = pmlHYXBottomXZ[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHYZBottomYZ[pmlWidth_minus_1][j][k+pmlWidth];	      lowerX2 = pmlHYXBottomYZ[pmlWidth_minus_1][j][k+pmlWidth];	    }	  currentEpsilon = materialConstants[pmlMaterialBottomXZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialBottomXZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomXZ[i][j][k]); 	  pmlEZXBottomXZ[i][j][k] = pmlEZXBottomXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomXZ[i][j][k])) / commonFactor +	    (pmlHYZBottomXZ[i][j][k] - lowerX1 +	     pmlHYXBottomXZ[i][j][k] - lowerX2) /	    (dx * pmlSXBottomXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomXZ[i][j][k] * pmlComplexEZXBottomXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYBottomXZ[i][j][k]); 	  pmlEZYBottomXZ[i][j][k] = pmlEZYBottomXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomXZ[i][j][k])) / commonFactor -	    (pmlHXYBottomXZ[i][j][k] - pmlHXYBottomXZ[i][j-1][k] +	     pmlHXZBottomXZ[i][j][k] - pmlHXZBottomXZ[i][j-1][k]) /	    (dy * pmlSYBottomXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYBottomXZ[i][j][k] * pmlComplexEZYBottomXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}    /* Top XZ */  for(i=0; i<nx; i++)    for(j=0; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHYZTopXZ[i-1][j][k];	      lowerX2 = pmlHYXTopXZ[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHYZBottomYZ[pmlWidth_minus_1][j+ny_plus_pmlWidth][k+pmlWidth];	      lowerX2 = pmlHYXBottomYZ[pmlWidth_minus_1][j+ny_plus_pmlWidth][k+pmlWidth];	    }	  if (j > 0)	    {	      lowerY1 = pmlHXYTopXZ[i][j-1][k];	      lowerY2 = pmlHXZTopXZ[i][j-1][k];	    }	  else	    {	      lowerY1 = hx[i][ny_minus_1][k];	      lowerY2 = 0.0;	    }	  	  currentEpsilon = materialConstants[pmlMaterialTopXZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopXZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopXZ[i][j][k]); 	  pmlEZXTopXZ[i][j][k] = pmlEZXTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopXZ[i][j][k])) / commonFactor +	    (pmlHYZTopXZ[i][j][k] - lowerX1 +	     pmlHYXTopXZ[i][j][k] - lowerX2) /	    (dx * pmlSXTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopXZ[i][j][k] * pmlComplexEZXTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYTopXZ[i][j][k]); 	  pmlEZYTopXZ[i][j][k] = pmlEZYTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopXZ[i][j][k])) / commonFactor -	    (pmlHXYTopXZ[i][j][k] - lowerY1 +	     pmlHXZTopXZ[i][j][k] - lowerY2) /	    (dy * pmlSYTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYTopXZ[i][j][k] * pmlComplexEZYTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}}void UpdatePMLComplexE(void){   /* speed up variables */  int ny_plus_2pmlWidth = ny + 2*pmlWidth;  int nz_plus_2pmlWidth = nz + 2*pmlWidth;   /* Bottom YZ */    for(i=0; i<pmlWidth; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<nz_plus_2pmlWidth; k++)	{	  pmlComplexEXYBottomYZ[i][j][k] += dt * pmlEXYBottomYZ[i][j][k];	  pmlComplexEXZBottomYZ[i][j][k] += dt * pmlEXZBottomYZ[i][j][k];	  pmlComplexEYXBottomYZ[i][j][k] += dt * pmlEYXBottomYZ[i][j][k];	  pmlComplexEYZBottomYZ[i][j][k] += dt * pmlEYZBottomYZ[i][j][k];	  pmlComplexEZXBottomYZ[i][j][k] += dt * pmlEZXBottomYZ[i][j][k];	  pmlComplexEZYBottomYZ[i][j][k] += dt * pmlEZYBottomYZ[i][j][k];	}  /* Top YZ */  for(i=0; i<pmlWidth; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<nz_plus_2pmlWidth; k++)	{	  pmlComplexEXYTopYZ[i][j][k] += dt * pmlEXYTopYZ[i][j][k];	  pmlComplexEXZTopYZ[i][j][k] += dt * pmlEXZTopYZ[i][j][k];	  pmlComplexEYXTopYZ[i][j][k] += dt * pmlEYXTopYZ[i][j][k];	  pmlComplexEYZTopYZ[i][j][k] += dt * pmlEYZTopYZ[i][j][k];	  pmlComplexEZXTopYZ[i][j][k] += dt * pmlEZXTopYZ[i][j][k];	  pmlComplexEZYTopYZ[i][j][k] += dt * pmlEZYTopYZ[i][j][k];	}  /* Bottom XY */  for(i=0; i<nx; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<pmlWidth; k++)	{	  pmlComplexEXYBottomXY[i][j][k] += dt * pmlEXYBottomXY[i][j][k];	  pmlComplexEXZBottomXY[i][j][k] += dt * pmlEXZBottomXY[i][j][k];	  pmlComplexEYXBottomXY[i][j][k] += dt * pmlEYXBottomXY[i][j][k];	  pmlComplexEYZBottomXY[i][j][k] += dt * pmlEYZBottomXY[i][j][k];	  pmlComplexEZXBottomXY[i][j][k] += dt * pmlEZXBottomXY[i][j][k];	  pmlComplexEZYBottomXY[i][j][k] += dt * pmlEZYBottomXY[i][j][k];	}  /* Top XY */  for(i=0; i<nx; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<pmlWidth; k++)	{	  pmlComplexEXYTopXY[i][j][k] += dt * pmlEXYTopXY[i][j][k];	  pmlComplexEXZTopXY[i][j][k] += dt * pmlEXZTopXY[i][j][k];	  pmlComplexEYXTopXY[i][j][k] += dt * pmlEYXTopXY[i][j][k];	  pmlComplexEYZTopXY[i][j][k] += dt * pmlEYZTopXY[i][j][k];	  pmlComplexEZXTopXY[i][j][k] += dt * pmlEZXTopXY[i][j][k];	  pmlComplexEZYTopXY[i][j][k] += dt * pmlEZYTopXY[i][j][k];	}  /* BottomXZ */  for(i=0; i<nx; i++)    for(j=0; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{	  pmlComplexEXYBottomXZ[i][j][k] += dt * pmlEXYBottomXZ[i][j][k];	  pmlComplexEXZBottomXZ[i][j][k] += dt * pmlEXZBottomXZ[i][j][k];	  pmlComplexEYXBottomXZ[i][j][k] += dt * pmlEYXBottomXZ[i][j][k];	  pmlComplexEYZBottomXZ[i][j][k] += dt * pmlEYZBottomXZ[i][j][k];	  pmlComplexEZXBottomXZ[i][j][k] += dt * pmlEZXBottomXZ[i][j][k];	  pmlComplexEZYBottomXZ[i][j][k] += dt * pmlEZYBottomXZ[i][j][k];	}    /* TopXZ */  for(i=0; i<nx; i++)    for(j=0; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{	  pmlComplexEXYTopXZ[i][j][k] += dt * pmlEXYTopXZ[i][j][k];	  pmlComplexEXZTopXZ[i][j][k] += dt * pmlEXZTopXZ[i][j][k];	  pmlComplexEYXTopXZ[i][j][k] += dt * pmlEYXTopXZ[i][j][k];	  pmlComplexEYZTopXZ[i][j][k] += dt * pmlEYZTopXZ[i][j][k];	  pmlComplexEZXTopXZ[i][j][k] += dt * pmlEZXTopXZ[i][j][k];	  pmlComplexEZYTopXZ[i][j][k] += dt * pmlEZYTopXZ[i][j][k];	}  }