computeboundarysigma.c

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

/* Copyright notice:
   radarFDTD - A free 3D-FDTD simulation for EM-waves with GPML-ABCs ;-)
   Copyright (C) 2000 Carsten Aulbert
   ( 
   I used the following program as a 'manual', so in my opinion this needs to be quoted:
   ToyFDTD, version 1.02 
   The if-I-can-do-it-you-can-do-it FDTD! 
   Copyright (C) 1998,1999 Laurie E. Miller, Paul Hayes, Matthew O'Keefe 
   )

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

FILE *foutput;
char filename[1024];

#ifdef DEBUG
int tmpMaterial;
#endif

void SetPMLMaterial(void)
{
  int i,j,k;
  int relPosX, relPosY, relPosZ;

  for(i=0; i<nx+2*pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  if (i<pmlWidth)
	    relPosX = 0;
	  else if (i<nx+pmlWidth)
	    relPosX = i-pmlWidth;
	  else
	    relPosX = nx-1;

	  if (j<pmlWidth)
	    relPosY = 0;
	  else if (j<ny+pmlWidth)
	    relPosY = j-pmlWidth;
	  else
	    relPosY = ny-1;

	  if (k<pmlWidth)
	    relPosZ = 0;
	  else if (k<nz+pmlWidth)
	    relPosZ = k-pmlWidth;
	  else
	    relPosZ = nz-1;

	  if (i<pmlWidth)
	    pmlMaterialBottomYZ[i][j][k] = material[relPosX][relPosY][relPosZ];
	  else if (i<nx+pmlWidth)
	    if (k<pmlWidth)
	      pmlMaterialBottomXY[i-pmlWidth][j][k] = material[relPosX][relPosY][relPosZ];
	    else if(k<nz+pmlWidth)
	      if (j<pmlWidth)
		pmlMaterialBottomXZ[i-pmlWidth][j][k-pmlWidth] = material[relPosX][relPosY][relPosZ];
	      else if (j<ny+pmlWidth)
		{}
	      else 
		pmlMaterialTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth] = material[relPosX][relPosY][relPosZ];
	    else
	      pmlMaterialTopXY[i-pmlWidth][j][k-nz-pmlWidth] = material[relPosX][relPosY][relPosZ];
	  else
	    pmlMaterialTopYZ[i-nx-pmlWidth][j][k] = material[relPosX][relPosY][relPosZ];
	}

#ifdef DEBUG
  sprintf(filename, "./material.dat");
  foutput = fopen(filename, "w");
  for(i=0; i<nx+2*pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  if (i<pmlWidth)
	    tmpMaterial = pmlMaterialBottomYZ[i][j][k];
	  else if (i<nx+pmlWidth)
	    if (k<pmlWidth)
	      tmpMaterial = pmlMaterialBottomXY[i-pmlWidth][j][k];
	    else if(k<nz+pmlWidth)
	      if (j<pmlWidth)
		tmpMaterial = pmlMaterialBottomXZ[i-pmlWidth][j][k-pmlWidth];
	      else if (j<ny+pmlWidth)
		tmpMaterial = -1;
	      else 
		tmpMaterial = pmlMaterialTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth];
	    else
	      tmpMaterial = pmlMaterialTopXY[i-pmlWidth][j][k-nz-pmlWidth];
	  else
	    tmpMaterial = pmlMaterialBottomYZ[i-nx-pmlWidth][j][k];
	  fprintf(foutput, "%d %d %d %d\n", i,j,k,tmpMaterial);
	}
  fclose(foutput);
#endif
}

void SetConductivities(void)
{
  /*****************************************************************************************/
  /************************************* BottomYZ ******************************************/
  /*****************************************************************************************/
  
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaXBottomYZ[i][j][k]     = ReturnConductivity((pmlWidth-i) * dx, pmlWidth * dx, maxSigma);
	  pmlSigmaXStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaXBottomYZ[i][j][k];
	  pmlSXBottomYZ[i][j][k]         = ReturnStretching((pmlWidth-i) * dx, pmlWidth * dx, maxStretching, stretchSteepness);
	}
  
  /* computing sigmaY within the edges */

  for(i=0; i<pmlWidth; i++)
    for(j=0; j<pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaYBottomYZ[i][j][k]     = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaYBottomYZ[i][j][k];
	  pmlSYBottomYZ[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++)
	{
	  pmlSigmaYBottomYZ[i][j][k]     = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
	  pmlSigmaYStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaYBottomYZ[i][j][k];
	  pmlSYBottomYZ[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++)
	{
	  pmlSigmaZBottomYZ[i][j][k]     = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaZBottomYZ[i][j][k];
	  pmlSZBottomYZ[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++)
	{
	  pmlSigmaZBottomYZ[i][j][k]     = ReturnConductivity((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxSigma);
	  pmlSigmaZStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaZBottomYZ[i][j][k];
	  pmlSZBottomYZ[i][j][k]         = ReturnStretching((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
	}
  
  
  /*****************************************************************************************/
  /*************************************** TopYZ *******************************************/
  /*****************************************************************************************/
  
  for(i=0; i<pmlWidth; i++)
    for(j=0; j<ny+2*pmlWidth; j++)
      for(k=0; k<nz+2*pmlWidth; k++)
	{
	  pmlSigmaXTopYZ[i][j][k]     = ReturnConductivity((i+1) * dx, pmlWidth * dx, maxSigma);
	  pmlSigmaXStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaXTopYZ[i][j][k];
	  pmlSXTopYZ[i][j][k]         = ReturnStretching((i+1) * dx, pmlWidth * dx, maxStretching, stretchSteepness);
	}
  
  /* computing sigmaY within the edges */

  for(i=0; i<pmlWidth; i++)
    for(j=0; j<pmlWidth; j++)