(dz * pmlSZBottomXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZBottomXZ[i][j][k] * pmlComplexEXZBottomXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* TopXZ */  for(i=0; i<nx; i++)    for(j=0; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{	  if (j > 0)	    {	      lowerY1 = pmlHZXTopXZ[i][j-1][k];	      lowerY2 = pmlHZYTopXZ[i][j-1][k];	    }	  else	    {	      lowerY1 = hz[i][ny_minus_1][k];	      lowerY2 = 0.0;	    }	  	  if (k > 0)	    {	      lowerZ1 = pmlHYZTopXZ[i][j][k-1];	      lowerZ2 = pmlHYXTopXZ[i][j][k-1];	    }	  else	    {	      lowerZ1 = pmlHYZBottomXY[i][j+ny_plus_pmlWidth][pmlWidth_minus_1];	      lowerZ2 = pmlHYXBottomXY[i][j+ny_plus_pmlWidth][pmlWidth_minus_1];	    }	  	  currentEpsilon = materialConstants[pmlMaterialTopXZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopXZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYTopXZ[i][j][k]); 	  pmlEXYTopXZ[i][j][k] = pmlEXYTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopXZ[i][j][k])) / commonFactor +	    (pmlHZXTopXZ[i][j][k] - lowerY1 +	     pmlHZYTopXZ[i][j][k] - lowerY2) /	    (dy * pmlSYTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaYTopXZ[i][j][k] * pmlComplexEXYTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZTopXZ[i][j][k]); 	  pmlEXZTopXZ[i][j][k] = pmlEXZTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopXZ[i][j][k])) / commonFactor -	    (pmlHYZTopXZ[i][j][k] - lowerZ1 +	     pmlHYXTopXZ[i][j][k] - lowerZ2) /	    (dz * pmlSZTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZTopXZ[i][j][k] * pmlComplexEXZTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}}void UpdatePMLEy(void){  int i,j,k;  PRECISION lowerX1, lowerX2, lowerZ1, lowerZ2;  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 nz_minus_1        = nz-1;    /* BottomYZ */  for(i=1; i<pmlWidth; i++) /* 1 ! */    for(j=1; j<ny_plus_2pmlWidth; j++)      for(k=1; k<nz_plus_2pmlWidth; k++) /* 1 ! */	{	  currentEpsilon = materialConstants[pmlMaterialBottomYZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialBottomYZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZBottomYZ[i][j][k]); 	  pmlEYZBottomYZ[i][j][k] = pmlEYZBottomYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomYZ[i][j][k])) / commonFactor +	    (pmlHXYBottomYZ[i][j][k] - pmlHXYBottomYZ[i][j][k-1] +	     pmlHXZBottomYZ[i][j][k] - pmlHXZBottomYZ[i][j][k-1]) /	    (dz * pmlSZBottomYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZBottomYZ[i][j][k] * pmlComplexEYZBottomYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomYZ[i][j][k]); 	  pmlEYXBottomYZ[i][j][k] = pmlEYXBottomYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomYZ[i][j][k])) / commonFactor -	    (pmlHZXBottomYZ[i][j][k] - pmlHZXBottomYZ[i-1][j][k] +	     pmlHZYBottomYZ[i][j][k] - pmlHZYBottomYZ[i-1][j][k]) /	    (dx * pmlSXBottomYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomYZ[i][j][k] * pmlComplexEYXBottomYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* TopYZ */  for(i=0; i<pmlWidth; i++)    for(j=1; j<ny_plus_2pmlWidth; j++)      for(k=1; k<nz_plus_2pmlWidth; k++) /* 1 ! */	{	  if (i > 0)	    {	      lowerX1 = pmlHZXTopYZ[i-1][j][k];	      lowerX2 = pmlHZYTopYZ[i-1][j][k];	    }	  else	    {	      if (k < pmlWidth)		{		  lowerX1 = pmlHZXBottomXY[ nx_minus_1][j][k];		  lowerX2 = pmlHZYBottomXY[ nx_minus_1][j][k];		}	      else if (k < nz_plus_pmlWidth)		{		  if (j < pmlWidth)		    {		      lowerX1 = pmlHZXBottomXZ[ nx_minus_1][j][k-pmlWidth];		      lowerX2 = pmlHZYBottomXZ[ nx_minus_1][j][k-pmlWidth];		    }		  else if (j < ny_plus_pmlWidth)		    {		      lowerX1 = hz[ nx_minus_1][j-pmlWidth][k-pmlWidth];		      lowerX2 = 0.0;		    }		  else		    {		      lowerX1 = pmlHZXTopXZ[ nx_minus_1][j-ny_plus_pmlWidth][k-pmlWidth];		      lowerX2 = pmlHZYTopXZ[ nx_minus_1][j-ny_plus_pmlWidth][k-pmlWidth];		    }		}	      else		{		  lowerX1 = pmlHZXTopXY[ nx_minus_1][j][k-nz_plus_pmlWidth];		  lowerX2 = pmlHZYTopXY[ 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 + pmlSigmaZTopYZ[i][j][k]); 	  pmlEYZTopYZ[i][j][k] = pmlEYZTopYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopYZ[i][j][k])) / commonFactor +	    (pmlHXYTopYZ[i][j][k] - pmlHXYTopYZ[i][j][k-1] +	     pmlHXZTopYZ[i][j][k] - pmlHXZTopYZ[i][j][k-1]) /	    (dz * pmlSZTopYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZTopYZ[i][j][k] * pmlComplexEYZTopYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopYZ[i][j][k]); 	  pmlEYXTopYZ[i][j][k] = pmlEYXTopYZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopYZ[i][j][k])) / commonFactor -	    (pmlHZXTopYZ[i][j][k] - lowerX1 +	     pmlHZYTopYZ[i][j][k] - lowerX2) /	    (dx * pmlSXTopYZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopYZ[i][j][k] * pmlComplexEYXTopYZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Bottom XY */  for(i=0; i<nx; i++)    for(j=1; j<ny_plus_2pmlWidth; j++)      for(k=1; k<pmlWidth; k++) /* 1 ! */	{	  if (i > 0)	    {	      lowerX1 = pmlHZXBottomXY[i-1][j][k];	      lowerX2 = pmlHZYBottomXY[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHZXBottomYZ[pmlWidth_minus_1][j][k];	      lowerX2 = pmlHZYBottomYZ[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 + pmlSigmaZBottomXY[i][j][k]); 	  pmlEYZBottomXY[i][j][k] = pmlEYZBottomXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomXY[i][j][k])) / commonFactor +	    (pmlHXYBottomXY[i][j][k] - pmlHXYBottomXY[i][j][k-1] +	     pmlHXZBottomXY[i][j][k] - pmlHXZBottomXY[i][j][k-1]) /	    (dz * pmlSZBottomXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZBottomXY[i][j][k] * pmlComplexEYZBottomXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomXY[i][j][k]); 	  pmlEYXBottomXY[i][j][k] = pmlEYXBottomXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomXY[i][j][k])) / commonFactor -	    (pmlHZXBottomXY[i][j][k] - lowerX1 +	     pmlHZYBottomXY[i][j][k] - lowerX2) /	    (dx * pmlSXBottomXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomXY[i][j][k] * pmlComplexEYXBottomXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Top XY */  for(i=0; i<nx; i++)    for(j=1; j<ny_plus_2pmlWidth; j++)      for(k=0; k<pmlWidth; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHZXTopXY[i-1][j][k];	      lowerX2 = pmlHZYTopXY[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHZXBottomYZ[pmlWidth_minus_1][j][k+nz_plus_pmlWidth];	      lowerX2 = pmlHZYBottomYZ[pmlWidth_minus_1][j][k+nz_plus_pmlWidth];	    }	  if (k > 0)	    {	      lowerZ1 = pmlHXYTopXY[i][j][k-1];	      lowerZ2 = pmlHXZTopXY[i][j][k-1];	    }	  else	    {	      if (j < pmlWidth)		{		  lowerZ1 = pmlHXYBottomXZ[i][j][nz_minus_1];		  lowerZ2 = pmlHXZBottomXZ[i][j][nz_minus_1];		}	      else if (j < ny_plus_pmlWidth)		{		  lowerZ1 = hx[i][j-pmlWidth][nz_minus_1];		  lowerZ2 = 0.0;		}	      else		{		  lowerZ1 = pmlHXYTopXZ[i][j-ny_plus_pmlWidth][nz_minus_1];		  lowerZ2 = pmlHXZTopXZ[i][j-ny_plus_pmlWidth][nz_minus_1];		}	    }	  	  currentEpsilon = materialConstants[pmlMaterialTopXY[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopXY[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZTopXY[i][j][k]); 	  pmlEYZTopXY[i][j][k] = pmlEYZTopXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopXY[i][j][k])) / commonFactor +	    (pmlHXYTopXY[i][j][k] - lowerZ1 +	     pmlHXZTopXY[i][j][k] - lowerZ2) /	    (dz * pmlSZTopXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZTopXY[i][j][k] * pmlComplexEYZTopXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopXY[i][j][k]); 	  pmlEYXTopXY[i][j][k] = pmlEYXTopXY[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopXY[i][j][k])) / commonFactor -	    (pmlHZXTopXY[i][j][k] - lowerX1 +	     pmlHZYTopXY[i][j][k] - lowerX2) /	    (dx * pmlSXTopXY[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopXY[i][j][k] * pmlComplexEYXTopXY[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	}  /* Bottom XZ */    for(i=0; i<nx; i++)    for(j=1; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{	  if (i > 0)	    {	      lowerX1 = pmlHZXBottomXZ[i-1][j][k];	      lowerX2 = pmlHZYBottomXZ[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHZXBottomYZ[pmlWidth_minus_1][j][k+pmlWidth];	      lowerX2 = pmlHZYBottomYZ[pmlWidth_minus_1][j][k+pmlWidth];	    }	  	  if (k > 0)	    {	      lowerZ1 = pmlHXYBottomXZ[i][j][k-1];	      lowerZ2 = pmlHXZBottomXZ[i][j][k-1];	    }	  else	    {	      lowerZ1 = pmlHXYBottomXY[i][j][pmlWidth_minus_1];	      lowerZ2 = pmlHXZBottomXY[i][j][pmlWidth_minus_1];	    }	  	  currentEpsilon = materialConstants[pmlMaterialBottomXZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialBottomXZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZBottomXZ[i][j][k]); 	  pmlEYZBottomXZ[i][j][k] = pmlEYZBottomXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomXZ[i][j][k])) / commonFactor +	    (pmlHXYBottomXZ[i][j][k] - lowerZ1 +	     pmlHXZBottomXZ[i][j][k] - lowerZ2) /	    (dz * pmlSZBottomXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZBottomXZ[i][j][k] * pmlComplexEYZBottomXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXBottomXZ[i][j][k]); 	  pmlEYXBottomXZ[i][j][k] = pmlEYXBottomXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXBottomXZ[i][j][k])) / commonFactor -	    (pmlHZXBottomXZ[i][j][k] - lowerX1 +	     pmlHZYBottomXZ[i][j][k] - lowerX2) /	    (dx * pmlSXBottomXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXBottomXZ[i][j][k] * pmlComplexEYXBottomXZ[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 = pmlHZXTopXZ[i-1][j][k];	      lowerX2 = pmlHZYTopXZ[i-1][j][k];	    }	  else	    {	      lowerX1 = pmlHZXBottomYZ[pmlWidth_minus_1][j+ny_plus_pmlWidth][k+pmlWidth];	      lowerX2 = pmlHZYBottomYZ[pmlWidth_minus_1][j+ny_plus_pmlWidth][k+pmlWidth];	    }	  if (k > 0)	    {	      lowerZ1 = pmlHXYTopXZ[i][j][k-1];	      lowerZ2 = pmlHXZTopXZ[i][j][k-1];	    }	  else	    {	      lowerZ1 = pmlHXYBottomXY[i][j+ny_plus_pmlWidth][pmlWidth_minus_1];	      lowerZ2 = pmlHXZBottomXY[i][j+ny_plus_pmlWidth][pmlWidth_minus_1];	    }	  	  currentEpsilon = materialConstants[pmlMaterialTopXZ[i][j][k]][7];	  currentSigma   = materialConstants[pmlMaterialTopXZ[i][j][k]][8];	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZTopXZ[i][j][k]); 	  pmlEYZTopXZ[i][j][k] = pmlEYZTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopXZ[i][j][k])) / commonFactor +	    (pmlHXYTopXZ[i][j][k] - lowerZ1 +	     pmlHXZTopXZ[i][j][k] - lowerZ2) /	    (dz * pmlSZTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaZTopXZ[i][j][k] * pmlComplexEYZTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaXTopXZ[i][j][k]); 	  pmlEYXTopXZ[i][j][k] = pmlEYXTopXZ[i][j][k] *	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaXTopXZ[i][j][k])) / commonFactor -	    (pmlHZXTopXZ[i][j][k] - lowerX1 +	     pmlHZYTopXZ[i][j][k] - lowerX2) /	    (dx * pmlSXTopXZ[i][j][k] * commonFactor) -	    currentSigma * pmlSigmaXTopXZ[i][j][k] * pmlComplexEYXTopXZ[i][j][k] /	    (currentEpsilon * EPSILON_0 * commonFactor);	  	}}