computeboundarysigma.c

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

      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaYTopYZ[i][j][k]     = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaYTopYZ[i][j][k];
	  pmlSYTopYZ[i][j][k]         = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
  	}
  
  for(i=0; i<pmlWidth; i++)
    for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaYTopYZ[i][j][k]     = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaYTopYZ[i][j][k];
	  pmlSYTopYZ[i][j][k]         = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
  	}

  /* computing sigmaZ within the edges */

  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaZTopYZ[i][j][k]     = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaZTopYZ[i][j][k];
	  pmlSZTopYZ[i][j][k]         = ReturnStretching((pmlWidth-k) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
	}
  
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=nz+pmlWidth; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaZTopYZ[i][j][k]     = ReturnConductivity((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaZTopYZ[i][j][k];
	  pmlSZTopYZ[i][j][k]         = ReturnStretching((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
	}
    
  /*****************************************************************************************/
  /************************************* BottomXY ******************************************/
  /*****************************************************************************************/

  for(i=0; i<nx; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaZBottomXY[i][j][k]     = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaZBottomXY[i][j][k];
	  pmlSZBottomXY[i][j][k]         = ReturnStretching((pmlWidth-k) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
	}

  /* computing SigmaY within the edges */

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaYBottomXY[i][j][k]     = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaYBottomXY[i][j][k];
	  pmlSYBottomXY[i][j][k]         = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  for(i=0; i<nx; i++)
    for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaYBottomXY[i][j][k]     = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaYBottomXY[i][j][k];
	  pmlSYBottomXY[i][j][k]         = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  /*****************************************************************************************/
  /************************************** TopXY ********************************************/
  /*****************************************************************************************/

  for(i=0; i<nx; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaZTopXY[i][j][k]     = ReturnConductivity((k+1) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaZTopXY[i][j][k];
	  pmlSZTopXY[i][j][k]         = ReturnStretching((k+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
	}

  /* computing SigmaY within the edges */

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaYTopXY[i][j][k]     = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaYTopXY[i][j][k];
	  pmlSYTopXY[i][j][k]         = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  for(i=0; i<nx; i++)
    for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
      for(k=0; k<pmlWidth; k++)
	{
	  pmlSigmaYTopXY[i][j][k]     = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaYTopXY[i][j][k];
	  pmlSYTopXY[i][j][k]         = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  /*****************************************************************************************/
  /************************************ BottomXZ *******************************************/
  /*****************************************************************************************/

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  pmlSigmaYBottomXZ[i][j][k]     = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarBottomXZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXZ[i][j][k]][EPSILON]) * pmlSigmaYBottomXZ[i][j][k];
	  pmlSYBottomXZ[i][j][k]         = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  /*****************************************************************************************/
  /************************************** TopXZ ********************************************/
  /*****************************************************************************************/

  for(i=0; i<nx; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz; k++)
	{
	  pmlSigmaYTopXZ[i][j][k]     = ReturnConductivity((j+1) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarTopXZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXZ[i][j][k]][EPSILON]) * pmlSigmaYTopXZ[i][j][k];
	  pmlSYTopXZ[i][j][k]         = ReturnStretching((j+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
	}

  

#ifdef DEBUG
for(i=0; i<nx+2*pmlWidth; i++)
    {
      sprintf(filename, "pmlSigma%05d.dat", (int)i);

      /*sprintf(filename, "current/pic%05d.dat", (int) iteration / timeToPlot);*/
      foutput = fopen(filename, "w");
      for(j=0; j<ny+2*pmlWidth; j++)
	for(k=0; k<nz+2*pmlWidth; k++)
	  {
	    if (i<pmlWidth)
	      {
		fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomYZ[i][j][k], pmlSigmaYBottomYZ[i][j][k], pmlSigmaZBottomYZ[i][j][k]);
	      }
	    else if(i<pmlWidth+nx)
	      {
		if (k<pmlWidth)
		  {
		    fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomXY[i-pmlWidth][j][k], pmlSigmaYBottomXY[i-pmlWidth][j][k], pmlSigmaZBottomXY[i-pmlWidth][j][k]);
		  }
		else if (k<pmlWidth+nz)
		  {
		    if (j<pmlWidth)
		      {
			fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomXZ[i-pmlWidth][j][k-pmlWidth], pmlSigmaYBottomXZ[i-pmlWidth][j][k-pmlWidth], pmlSigmaZBottomXZ[i-pmlWidth][j][k-pmlWidth]);
		      }
		    else if (j<pmlWidth+ny)
		      {
			fprintf(foutput, "%d\t%d\tNaN\tNan\tNaN\n", j,k);
		      }
		    else
		      {
			fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth], pmlSigmaYTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth], pmlSigmaZTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth]);
		      }
		  }
		else
		  {
		    fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopXY[i-pmlWidth][j][k-nz-pmlWidth], pmlSigmaYTopXY[i-pmlWidth][j][k-nz-pmlWidth], pmlSigmaZTopXY[i-pmlWidth][j][k-nz-pmlWidth]);
		  }
	      }
	    else
	      {
		fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopYZ[i-nx-pmlWidth][j][k], pmlSigmaYTopYZ[i-nx-pmlWidth][j][k], pmlSigmaZTopYZ[i-nx-pmlWidth][j][k]);
	      }
	  }
      fclose(foutput);
    }
  
#endif

}