/* 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*//* all inner updates come here *//* Update the hx values: */void UpdateInnerHx(void){  int i,j,k;  PRECISION nextY, nextZ; /* the upper neighbour-values in y and z direction */    for(i=0; i<(nx); i++)     for(j=0; j<(ny); j++)      for(k=0; k<(nz); k++)	{	 	  if (j < ny-1)	    nextY = ez[i][j+1][k];	  else	    nextY = pmlEZXTopXZ[i][0][k] + pmlEZYTopXZ[i][0][k];	  	  if (k < nz-1)	    nextZ = ey[i][j][k+1];	  else	    nextZ = pmlEYXTopXY[i][j+pmlWidth][0] + pmlEYZTopXY[i][j+pmlWidth][0];	  hx[i][j][k] = hx[i][j][k] + (materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ey[i][j][k]) - materialConstants[material[i][j][k]][DTMUDY]*(nextY - ez[i][j][k]));	}} /* Update the hy values: */void UpdateInnerHy(void){  int i,j,k;  PRECISION nextX, nextZ; /*  the upper neighbour-values in y and z direction */  for(i=0; i<(nx); i++)    for(j=0; j<(ny); j++)      for(k=0; k<(nz); k++)	{	  if (i < nx-1)	    nextX = ez[i+1][j][k];	  else	    nextX = pmlEZXTopYZ[0][j+pmlWidth][k+pmlWidth] + pmlEZYTopYZ[0][j+pmlWidth][k+pmlWidth];	  	  if (k < nz-1)	    nextZ = ex[i][j][k+1];	  else	    nextZ = pmlEXYTopXY[i][j+pmlWidth][0] + pmlEXZTopXY[i][j+pmlWidth][0];	  	  hy[i][j][k] = hy[i][j][k] + (materialConstants[material[i][j][k]][DTMUDX]*(nextX - ez[i][j][k]) - materialConstants[material[i][j][k]][DTMUDZ]*(nextZ - ex[i][j][k]));	}}/* Update the hz values: */void UpdateInnerHz(void){  int i,j,k;  PRECISION nextX, nextY;    for(i=0; i<(nx); i++)    for(j=0; j<(ny); j++)      for(k=0; k<(nz); k++)	{	  if (i < nx-1)	    nextX = ey[i+1][j][k];	  else	    nextX = pmlEYXTopYZ[0][j+pmlWidth][k+pmlWidth] +  pmlEYZTopYZ[0][j+pmlWidth][k+pmlWidth];	  	  if (j < ny-1)	    nextY = ex[i][j+1][k];	  else	    nextY = pmlEXYTopXZ[i][0][k] + pmlEXZTopXZ[i][0][k];	  	  hz[i][j][k] = hz[i][j][k] + (materialConstants[material[i][j][k]][DTMUDY]*(nextY - ex[i][j][k]) - materialConstants[material[i][j][k]][DTMUDX]*(nextX - ey[i][j][k]));	}}/* Update the ex values: */void UpdateInnerEx(void){  int i,j,k;  PRECISION lowerY, lowerZ;    for(i=0; i<(nx); i++)    for(j=0; j<(ny); j++)       for(k=0; k<(nz); k++) 	{	  if (j > 0)	    lowerY = hz[i][j-1][k];	  else	    lowerY = pmlHZXBottomXZ[i][pmlWidth-1][k] + pmlHZYBottomXZ[i][pmlWidth-1][k];	  	  if (k > 0)	    lowerZ = hy[i][j][k-1];	  else	    lowerZ = pmlHYXBottomXY[i][j+pmlWidth][pmlWidth-1] + pmlHYZBottomXY[i][j+pmlWidth][pmlWidth-1];	  	  switch(material[i][j][k])	    {	    case 0: 	      ex[i][j][k] = ex[i][j][k] + materialConstants[0][DTEPSDY] * (hz[i][j][k] - lowerY) - materialConstants[0][DTEPSDZ] * (hy[i][j][k] - lowerZ);	      break;	    case 1:	      ex[i][j][k] = 0;	      break;	    default:	      ex[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ex[i][j][k] + materialConstants[material[i][j][k]][DTEPSDY]*(hz[i][j][k] - lowerY) - materialConstants[material[i][j][k]][DTEPSDZ]*(hy[i][j][k] - lowerZ);	    }	}}/* Update the ey values: */  void UpdateInnerEy(void){  int i,j,k;  PRECISION lowerX, lowerZ;    for(i=0; i<(nx); i++)     for(j=0; j<(ny); j++)      for(k=0; k<(nz); k++) 	{	  if (i > 0)	    lowerX = hz[i-1][j][k];	  else	    lowerX = pmlHZXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHZYBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];	  	  if (k > 0)	    lowerZ = hx[i][j][k-1]; 	  else	    lowerZ = pmlHXYBottomXY[i][j+pmlWidth][pmlWidth-1] +  pmlHXZBottomXY[i][j+pmlWidth][pmlWidth-1];	  	  switch(material[i][j][k])	    {	    case 0: 	      ey[i][j][k] = ey[i][j][k] + materialConstants[0][DTEPSDZ] * (hx[i][j][k] - lowerZ) - materialConstants[0][DTEPSDX] * (hz[i][j][k] - lowerX);	      break;	    case 1:	      ey[i][j][k] = 0;	      break;	    default:	      ey[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ey[i][j][k] + materialConstants[material[i][j][k]][DTEPSDZ]*(hx[i][j][k] - lowerZ) - materialConstants[material[i][j][k]][DTEPSDX]*(hz[i][j][k] - lowerX);	    }	}} /* Update the ez values: */void UpdateInnerEz(void){  int i,j,k;  PRECISION lowerX, lowerY;  for(i=0; i<(nx); i++)     for(j=0; j<(ny); j++)      for(k=0; k<(nz); k++)	{	  if (i > 0)	    lowerX = hy[i-1][j][k];	  else	    lowerX = pmlHYXBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth] + pmlHYZBottomYZ[pmlWidth-1][j+pmlWidth][k+pmlWidth];	  	  if (j > 0)	    lowerY = hx[i][j-1][k];	  else	    lowerY = pmlHXYBottomXZ[i][pmlWidth-1][k] + pmlHXZBottomXZ[i][pmlWidth-1][k];	  switch(material[i][j][k])	    {	    case 0: 	      ez[i][j][k] = ez[i][j][k] + materialConstants[0][DTEPSDX] * (hy[i][j][k] - lowerX) - materialConstants[0][DTEPSDY] * (hx[i][j][k] - lowerY);	      break;	    case 1:	      ez[i][j][k] = 0;	      break;	    default:	      ez[i][j][k] = materialConstants[material[i][j][k]][DTSIGEPS]*ez[i][j][k] + materialConstants[material[i][j][k]][DTEPSDX]*(hy[i][j][k] - lowerX) - materialConstants[material[i][j][k]][DTEPSDY]*(hx[i][j][k] - lowerY);	    }	}}/* THIS IS THE PML SECTION */void UpdatePMLHx(void){  int i,j,k;  PRECISION mudt = MU_0 / dt;  PRECISION topY1, topY2, topZ1, topZ2; /* components which may be in different memory blocks than the one worked on */  PRECISION commonFactor;    /* 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 nz_minus_1        = nz-1;    /* start with the YZ planes */  /* there are no problems at any boundary ! */  for(i=0; i<pmlWidth; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<nz_plus_2pmlWidth; k++)	{ /* bottom */	    	    commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomYZ[i][j][k];	    pmlHXYBottomYZ[i][j][k] = pmlHXYBottomYZ[i][j][k] * 	       (mudt - 0.5 * pmlSigmaYStarBottomYZ[i][j][k]) / commonFactor -	      (pmlEZXBottomYZ[i][j+1][k] - pmlEZXBottomYZ[i][j][k] + 	       pmlEZYBottomYZ[i][j+1][k] - pmlEZYBottomYZ[i][j][k]) / 	      (dy * pmlSYBottomYZ[i][j][k] * commonFactor);	    	    commonFactor = mudt + 0.5 * pmlSigmaZStarBottomYZ[i][j][k];	    pmlHXZBottomYZ[i][j][k] = pmlHXZBottomYZ[i][j][k] * 	       (mudt - 0.5 * pmlSigmaZStarBottomYZ[i][j][k]) / commonFactor +	      (pmlEYZBottomYZ[i][j][k+1] - pmlEYZBottomYZ[i][j][k] + 	       pmlEYXBottomYZ[i][j][k+1] - pmlEYXBottomYZ[i][j][k]) /	      (dz * pmlSZBottomYZ[i][j][k] * commonFactor);	}  /* top, there are no problems at any boundary */  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 * pmlSigmaYStarBottomYZ[i][j][k];	    pmlHXYTopYZ[i][j][k] = pmlHXYTopYZ[i][j][k] * 	       (mudt - 0.5 * pmlSigmaYStarTopYZ[i][j][k]) / commonFactor -	      (pmlEZXTopYZ[i][j+1][k] - pmlEZXTopYZ[i][j][k] + 	       pmlEZYTopYZ[i][j+1][k] - pmlEZYTopYZ[i][j][k]) / 	      (dy * pmlSYTopYZ[i][j][k] * commonFactor);	    	    commonFactor = mudt + 0.5 * pmlSigmaZStarTopYZ[i][j][k];	    pmlHXZTopYZ[i][j][k] = pmlHXZTopYZ[i][j][k] * 	       (mudt - 0.5 * pmlSigmaZStarTopYZ[i][j][k]) / commonFactor +	      (pmlEYZTopYZ[i][j][k+1] - pmlEYZTopYZ[i][j][k] + 	       pmlEYXTopYZ[i][j][k+1] - pmlEYXTopYZ[i][j][k]) /	      (dz * pmlSZTopYZ[i][j][k] * commonFactor);	}  /* now the XY planes */  for(i=0; i<nx; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<pmlWidth; k++)	{ /* bottom */	  /* there are three boundary-problems:	     0           <= y < pmlWidth:     XZ - plane (bottom)	     pmlWidth    <= y < pmlWidth+ny:  SimSpace	     pmlWidth+ny <= y:                XZ - plane (top)	  */	  if (k < pmlWidth_minus_1)	    {	      topZ1 = pmlEYZBottomXY[i][j][k+1];	      topZ2 = pmlEYXBottomXY[i][j][k+1];	    }	  else	    {	      if (j < pmlWidth)		{		  topZ1 = pmlEYZBottomXZ[i][j][0];		  topZ2 = pmlEYXBottomXZ[i][j][0];		}	      else if (j < ny_plus_pmlWidth)		{		  topZ1 = ey[i][j-pmlWidth][0];		  topZ2 = 0.0;		}	      else		{		  topZ1 = pmlEYZTopXZ[i][j-ny_plus_pmlWidth][0];		  topZ2 = pmlEYXTopXZ[i][j-ny_plus_pmlWidth][0];		}	    }	  commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomXY[i][j][k];	  pmlHXYBottomXY[i][j][k] = pmlHXYBottomXY[i][j][k] * 	    (mudt - 0.5 * pmlSigmaYStarBottomXY[i][j][k]) / commonFactor -	    (pmlEZXBottomXY[i][j+1][k] - pmlEZXBottomXY[i][j][k] + 	     pmlEZYBottomXY[i][j+1][k] - pmlEZYBottomXY[i][j][k]) / 	    (dy * pmlSYBottomXY[i][j][k] * commonFactor);	  	  commonFactor = mudt + 0.5 * pmlSigmaZStarBottomXY[i][j][k];	  pmlHXZBottomXY[i][j][k] = pmlHXZBottomXY[i][j][k] * 	    (mudt - 0.5 * pmlSigmaZStarBottomXY[i][j][k]) / commonFactor +	    (topZ1 - pmlEYZBottomXY[i][j][k] + 	     topZ2 - pmlEYXBottomXY[i][j][k]) /	    (dz * pmlSZBottomXY[i][j][k] * commonFactor);	  	}    for(i=0; i<nx; i++)    for(j=0; j<ny_plus_2pmlWidth; j++)      for(k=0; k<pmlWidth; k++)	{ 	  /* top */	  /* there are no problems at this boundary */	  commonFactor =  mudt + 0.5 * pmlSigmaYStarTopXY[i][j][k];	  pmlHXYTopXY[i][j][k] = pmlHXYTopXY[i][j][k] * 	    (mudt - 0.5 * pmlSigmaYStarTopXY[i][j][k]) / commonFactor -	    (pmlEZXTopXY[i][j+1][k] - pmlEZXTopXY[i][j][k] + 	     pmlEZYTopXY[i][j+1][k] - pmlEZYTopXY[i][j][k]) / 	    (dy * pmlSYTopXY[i][j][k] * commonFactor);	  	  commonFactor = mudt + 0.5 * pmlSigmaZStarTopXY[i][j][k];	  pmlHXZTopXY[i][j][k] = pmlHXZTopXY[i][j][k] * 	    (mudt - 0.5 * pmlSigmaZStarTopXY[i][j][k]) / commonFactor +	    (pmlEYZTopXY[i][j][k+1] - pmlEYZTopXY[i][j][k] + 	     pmlEYXTopXY[i][j][k+1] - pmlEYXTopXY[i][j][k]) /	    (dz * pmlSZTopXY[i][j][k] * commonFactor);	  	}     /* and finally the XZ planes */  for(i=0; i<nx; i++)    for(j=0; j<pmlWidth; j++)      for(k=0; k<nz; k++)	{ /* bottom */	  /* there are boundaries at:	     j = pmlWidth-1 for all i,k      SimSpace is neighbour	     k = ny_plus_pmlWidth-1 for all i,j   TopXY is neighbour 	  */	  if (j < pmlWidth_minus_1)	    {	  	topY1 = pmlEZXBottomXZ[i][j+1][k];		topY2 = pmlEZYBottomXZ[i][j+1][k];	    }	  else	    {	      topY1 = ez[i][0][k];	      topY2 = 0.0;	    }	  	  if (k < nz_minus_1)	    {	      topZ1 = pmlEYZBottomXZ[i][j][k+1];	      topZ2 = pmlEYXBottomXZ[i][j][k+1];	    }	  else	    {	      topZ1 = pmlEYZTopXY[i][j][0];	      topZ2 = pmlEYXTopXY[i][j][0];	    }	  	  commonFactor =  mudt + 0.5 * pmlSigmaYStarBottomXZ[i][j][k];	  pmlHXYBottomXZ[i][j][k] = pmlHXYBottomXZ[i][j][k] * 	    (mudt - 0.5 * pmlSigmaYStarBottomXZ[i][j][k]) / commonFactor -	    (topY1 - pmlEZXBottomXZ[i][j][k] + 	     topY2 - pmlEZYBottomXZ[i][j][k]) /