componentupdate.c

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

		  topX1 = pmlEYZTopXY[0][j][k-nz_plus_pmlWidth];
		  topX2 = pmlEYXTopXY[0][j][k-nz_plus_pmlWidth];
		}
	    }

	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomYZ[i][j][k];
	  pmlHZXBottomYZ[i][j][k] = pmlHZXBottomYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarBottomYZ[i][j][k]) / commonFactor -
	    (topX1 - pmlEYZBottomYZ[i][j][k] +
	     topX2 - pmlEYXBottomYZ[i][j][k]) /
	    (dx * pmlSXBottomYZ[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];
	  pmlHZYBottomYZ[i][j][k] = pmlHZYBottomYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarBottomYZ[i][j][k]) / commonFactor +
	    (pmlEXYBottomYZ[i][j+1][k] - pmlEXYBottomYZ[i][j][k] +
	     pmlEXZBottomYZ[i][j+1][k] - pmlEXZBottomYZ[i][j][k]) /
	    (dy * pmlSYBottomYZ[i][j][k] * commonFactor);
	}
  
  /* TopYZ - plane, no problem there */

  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<nz_plus_2pmlWidth; k++)
	{
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopYZ[i][j][k];
	  pmlHZXTopYZ[i][j][k] = pmlHZXTopYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarTopYZ[i][j][k]) / commonFactor -
	    (pmlEYZTopYZ[i+1][j][k] - pmlEYZTopYZ[i][j][k] +
	     pmlEYXTopYZ[i+1][j][k] - pmlEYXTopYZ[i][j][k]) /
	    (dx * pmlSXTopYZ[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarTopYZ[i][j][k];
	  pmlHZYTopYZ[i][j][k] = pmlHZYTopYZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarTopYZ[i][j][k]) / commonFactor +
	    (pmlEXYTopYZ[i][j+1][k] - pmlEXYTopYZ[i][j][k] +
	     pmlEXZTopYZ[i][j+1][k] - pmlEXZTopYZ[i][j][k]) /
	    (dy * pmlSYTopYZ[i][j][k] * commonFactor);
	  
	}
  
  /* Bottom XY */

  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEYZBottomXY[i+1][j][k];
	      topX2 = pmlEYXBottomXY[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEYZTopYZ[0][j][k];
	      topX2 = pmlEYXTopYZ[0][j][k];
	    }

	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomXY[i][j][k];
	  pmlHZXBottomXY[i][j][k] = pmlHZXBottomXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarBottomXY[i][j][k]) / commonFactor -
	    (topX1 - pmlEYZBottomXY[i][j][k] +
	     topX2 - pmlEYXBottomXY[i][j][k]) /
	    (dx * pmlSXBottomXY[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarBottomXY[i][j][k];
	  pmlHZYBottomXY[i][j][k] = pmlHZYBottomXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarBottomXY[i][j][k]) / commonFactor +
	    (pmlEXYBottomXY[i][j+1][k] - pmlEXYBottomXY[i][j][k] +
	     pmlEXZBottomXY[i][j+1][k] - pmlEXZBottomXY[i][j][k]) /
	    (dy * pmlSYBottomXY[i][j][k] * commonFactor);

	}

  /* Top XY */
  
  for(i=0; i<nx; i++)
    for(j=0; j<ny_plus_2pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEYZTopXY[i+1][j][k];
	      topX2 = pmlEYXTopXY[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEYZTopYZ[0][j][k+nz_plus_pmlWidth];
	      topX2 = pmlEYXTopYZ[0][j][k+nz_plus_pmlWidth];
	    }
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopXY[i][j][k];
	  pmlHZXTopXY[i][j][k] = pmlHZXTopXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarTopXY[i][j][k]) / commonFactor -
	    (topX1 - pmlEYZTopXY[i][j][k] +
	     topX2 - pmlEYXTopXY[i][j][k]) /
	    (dx * pmlSXTopXY[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarTopXY[i][j][k];
	  pmlHZYTopXY[i][j][k] = pmlHZYTopXY[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarTopXY[i][j][k]) / commonFactor +
	    (pmlEXYTopXY[i][j+1][k] - pmlEXYTopXY[i][j][k] +
	     pmlEXZTopXY[i][j+1][k] - pmlEXZTopXY[i][j][k]) /
	    (dy * pmlSYTopXY[i][j][k] * commonFactor);
	  
	}
  
  /* Bottom XZ */

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEYZBottomXZ[i+1][j][k];
	      topX2 = pmlEYXBottomXZ[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEYZTopYZ[0][j][k+pmlWidth];
	      topX2 = pmlEYXTopYZ[0][j][k+pmlWidth];
	    }
	  
	  if (j < pmlWidth_minus_1)
	    {
	      topY1 = pmlEXYBottomXZ[i][j+1][k];
	      topY2 = pmlEXZBottomXZ[i][j+1][k];
	    }
	  else
	    {
	      topY1 = ex[i][0][k];
	      topY2 = 0.0;
	    }
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarBottomXZ[i][j][k];
	  pmlHZXBottomXZ[i][j][k] = pmlHZXBottomXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarBottomXZ[i][j][k]) / commonFactor -
	    (topX1 - pmlEYZBottomXZ[i][j][k] +
	     topX2 - pmlEYXBottomXZ[i][j][k]) /
	    (dx * pmlSXBottomXZ[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarBottomXZ[i][j][k];
	  pmlHZYBottomXZ[i][j][k] = pmlHZYBottomXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarBottomXZ[i][j][k]) / commonFactor +
	    (topY1 - pmlEXYBottomXZ[i][j][k] +
	     topY2 - pmlEXZBottomXZ[i][j][k]) /
	    (dy * pmlSYBottomXZ[i][j][k] * commonFactor);

	}
  
  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  if (i <  nx_minus_1)
	    {
	      topX1 = pmlEYZTopXZ[i+1][j][k];
	      topX2 = pmlEYXTopXZ[i+1][j][k];
	    }
	  else
	    {
	      topX1 = pmlEYZTopYZ[0][j+ny_plus_pmlWidth][k+pmlWidth];
	      topX2 = pmlEYXTopYZ[0][j+ny_plus_pmlWidth][k+pmlWidth];
	    }
	  
	  commonFactor = mudt + 0.5 * pmlSigmaXStarTopXZ[i][j][k];
	  pmlHZXTopXZ[i][j][k] = pmlHZXTopXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaXStarTopXZ[i][j][k]) / commonFactor -
	    (topX1 - pmlEYZTopXZ[i][j][k] +
	     topX2 - pmlEYXTopXZ[i][j][k]) /
	    (dx * pmlSXTopXZ[i][j][k] * commonFactor);

	  commonFactor = mudt + 0.5 * pmlSigmaYStarTopXZ[i][j][k];
	  pmlHZYTopXZ[i][j][k] = pmlHZYTopXZ[i][j][k] *
	    (mudt - 0.5 * pmlSigmaYStarTopXZ[i][j][k]) / commonFactor +
	    (pmlEXYTopXZ[i][j+1][k] - pmlEXYTopXZ[i][j][k] +
	     pmlEXZTopXZ[i][j+1][k] - pmlEXZTopXZ[i][j][k]) /
	    (dy * pmlSYTopXZ[i][j][k] * commonFactor);

	}
}

/* That磗 all to do with the H - field, now I have to do the same with the E -field */

void UpdatePMLEx(void)
{
  int i,j,k;
  PRECISION lowerY1, lowerY2, lowerZ1, lowerZ2; /* these components belong to another memory space */
  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 pmlWidth_minus_1  = pmlWidth-1;

  int ny_minus_1        = ny-1;
  int nz_minus_1        = nz-1;
   

  /* Start again with the bottom YZ plane */
  /* There are no problems with x, y and z, as long as we start at 1 with j and k */
  
  for(i=1; i<pmlWidth; i++)
    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */
      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 + pmlSigmaYBottomYZ[i][j][k]);
 	  pmlEXYBottomYZ[i][j][k] = pmlEXYBottomYZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomYZ[i][j][k])) / commonFactor +
	    (pmlHZXBottomYZ[i][j][k] - pmlHZXBottomYZ[i][j-1][k] +
	     pmlHZYBottomYZ[i][j][k] - pmlHZYBottomYZ[i][j-1][k]) /
	    (dy * pmlSYBottomYZ[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaYBottomYZ[i][j][k] * pmlComplexEXYBottomYZ[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZBottomYZ[i][j][k]);
 	  pmlEXZBottomYZ[i][j][k] = pmlEXZBottomYZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomYZ[i][j][k])) / commonFactor -
	    (pmlHYZBottomYZ[i][j][k] - pmlHYZBottomYZ[i][j][k-1] +
	     pmlHYXBottomYZ[i][j][k] - pmlHYXBottomYZ[i][j][k-1]) /
	    (dz * pmlSZBottomYZ[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaZBottomYZ[i][j][k] * pmlComplexEXZBottomYZ[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	}
  /* Now the upper part: TopYZ */

  for(i=0; i<pmlWidth; i++)
    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */
      for(k=1; k<nz_plus_2pmlWidth; k++) /* 1 ! */
	{
	  currentEpsilon = materialConstants[pmlMaterialTopYZ[i][j][k]][7];
	  currentSigma   = materialConstants[pmlMaterialTopYZ[i][j][k]][8];
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYTopYZ[i][j][k]);
 	  pmlEXYTopYZ[i][j][k] = pmlEXYTopYZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopYZ[i][j][k])) / commonFactor +
	    (pmlHZXTopYZ[i][j][k] - pmlHZXTopYZ[i][j-1][k] +
	     pmlHZYTopYZ[i][j][k] - pmlHZYTopYZ[i][j-1][k]) /
	    (dy * pmlSYTopYZ[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaYTopYZ[i][j][k] * pmlComplexEXYTopYZ[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZTopYZ[i][j][k]);
 	  pmlEXZTopYZ[i][j][k] = pmlEXZTopYZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopYZ[i][j][k])) / commonFactor -
	    (pmlHYZTopYZ[i][j][k] - pmlHYZTopYZ[i][j][k-1] +
	     pmlHYXTopYZ[i][j][k] - pmlHYXTopYZ[i][j][k-1]) /
	    (dz * pmlSZTopYZ[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaZTopYZ[i][j][k] * pmlComplexEXZTopYZ[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	}
  
  /* BottomXY */

  for(i=0; i<nx; i++)
    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */
      for(k=1; k<pmlWidth; k++)
	{
	  currentEpsilon = materialConstants[pmlMaterialBottomXY[i][j][k]][7];
	  currentSigma   = materialConstants[pmlMaterialBottomXY[i][j][k]][8];
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaYBottomXY[i][j][k]);
 	  pmlEXYBottomXY[i][j][k] = pmlEXYBottomXY[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomXY[i][j][k])) / commonFactor +
	    (pmlHZXBottomXY[i][j][k] - pmlHZXBottomXY[i][j-1][k] +
	     pmlHZYBottomXY[i][j][k] - pmlHZYBottomXY[i][j-1][k]) /
	    (dy * pmlSYBottomXY[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaYBottomXY[i][j][k] * pmlComplexEXYBottomXY[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZBottomXY[i][j][k]);
 	  pmlEXZBottomXY[i][j][k] = pmlEXZBottomXY[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomXY[i][j][k])) / commonFactor -
	    (pmlHYZBottomXY[i][j][k] - pmlHYZBottomXY[i][j][k-1] +
	     pmlHYXBottomXY[i][j][k] - pmlHYXBottomXY[i][j][k-1]) /
	    (dz * pmlSZBottomXY[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaZBottomXY[i][j][k] * pmlComplexEXZBottomXY[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	}

  /* TopXY */

  for(i=0; i<nx; i++)
    for(j=1; j<ny_plus_2pmlWidth; j++) /* 1 ! */
      for(k=0; k<pmlWidth; k++)
	{
	  if (k > 0)
	    {
	      lowerZ1 = pmlHYZTopXY[i][j][k-1];
	      lowerZ2 = pmlHYXTopXY[i][j][k-1];
	    }
	  else
	    {
	      if (j < pmlWidth)
		{
		  lowerZ1 = pmlHYZBottomXZ[i][j][nz_minus_1];
		  lowerZ2 = pmlHYXBottomXZ[i][j][nz_minus_1];
		}
	      else if (j < ny_plus_pmlWidth)
		{
		  lowerZ1 = hy[i][j-pmlWidth][nz_minus_1];
		  lowerZ2 = 0.0;
		}
	      else
		{
		  lowerZ1 = pmlHYZTopXZ[i][j-ny_plus_pmlWidth][nz_minus_1];
		  lowerZ2 = pmlHYXTopXZ[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 + pmlSigmaYTopXY[i][j][k]);
 	  pmlEXYTopXY[i][j][k] = pmlEXYTopXY[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYTopXY[i][j][k])) / commonFactor +
	    (pmlHZXTopXY[i][j][k] - pmlHZXTopXY[i][j-1][k] +
	     pmlHZYTopXY[i][j][k] - pmlHZYTopXY[i][j-1][k]) /
	    (dy * pmlSYTopXY[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaYTopXY[i][j][k] * pmlComplexEXYTopXY[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZTopXY[i][j][k]);
 	  pmlEXZTopXY[i][j][k] = pmlEXZTopXY[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZTopXY[i][j][k])) / commonFactor -
	    (pmlHYZTopXY[i][j][k] - lowerZ1 +
	     pmlHYXTopXY[i][j][k] - lowerZ2) /
	    (dz * pmlSZTopXY[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaZTopXY[i][j][k] * pmlComplexEXZTopXY[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 (k > 0)
	    {
	      lowerZ1 = pmlHYZBottomXZ[i][j][k-1];
	      lowerZ2 = pmlHYXBottomXZ[i][j][k-1];
	    }
	  else 
	    {
	      lowerZ1 = pmlHYZBottomXY[i][j][pmlWidth_minus_1];
	      lowerZ2 = pmlHYXBottomXY[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 + pmlSigmaYBottomXZ[i][j][k]);
 	  pmlEXYBottomXZ[i][j][k] = pmlEXYBottomXZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaYBottomXZ[i][j][k])) / commonFactor +
	    (pmlHZXBottomXZ[i][j][k] - pmlHZXBottomXZ[i][j-1][k] +
	     pmlHZYBottomXZ[i][j][k] - pmlHZYBottomXZ[i][j-1][k]) /
	    (dy * pmlSYBottomXZ[i][j][k] * commonFactor) -
	    currentSigma * pmlSigmaYBottomXZ[i][j][k] * pmlComplexEXYBottomXZ[i][j][k] /
	    (currentEpsilon * EPSILON_0 * commonFactor);
	  
	  commonFactor = currentEpsilon * epsdt + 0.5 * (currentSigma + pmlSigmaZBottomXZ[i][j][k]);
 	  pmlEXZBottomXZ[i][j][k] = pmlEXZBottomXZ[i][j][k] *
	    (currentEpsilon * epsdt - 0.5 * (currentSigma + pmlSigmaZBottomXZ[i][j][k])) / commonFactor -
	    (pmlHYZBottomXZ[i][j][k] - lowerZ1 +
	     pmlHYXBottomXZ[i][j][k] - lowerZ2) /