componentupdate.c

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

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];
	}
  }