componentupdate.c

a full 3D simulation of electromagnetic waves with efficient absorbing boundary a full 3D simulation

	    (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);
	  
	}
}