output.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
*/

/* returns field strength of a specific point */

/* convention: (0,0,0) -> (nx,ny,nz) is the simulation space, so pml-coordinates can be negative */

/* mode is one of the following:
   total H, H_x, H_y, H_z, E_x, E_y, E_z, total E
*/

/* currently there is no range checking ! */

PRECISION ReturnFieldStrength(int x, int y, int z, int mode)
{
  PRECISION result = 0.0;
  /* basically this is only a big 'if' */

  if (x < 0) /* use BottomYZ */
    switch (mode) {
      0: result = sqrt(pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth]);break;
      1: result = pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
      2: result = pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
      3: result = pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
      4: result = pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
      5: result = pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
      6: result = pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
    default: result = sqrt(pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth]);
    }
  else if (x<nx)
    if (z<0)
      switch (mode) { /* BottomXY */
	0: result = sqrt(pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth]);break;
	1: result = pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
	2: result = pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
	3: result = pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth];break;
	4: result = pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
	5: result = pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
	6: result = pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth];break;
      default: result = sqrt(pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth]);
      }
    else
      if (z < nz)
	if (y < 0)
	  switch (mode) { /* BottomXZ */
	    0: result = sqrt(pmlHXYBottomXZ[x][y+pmlWidth][z] * pmlHXYBottomXZ[x][y+pmlWidth][z] + pmlHXZBottomXZ[x][y+pmlWidth][z] * pmlHXZBottomXZ[x][y+pmlWidth][z] + pmlHYXBottomXZ[x][y+pmlWidth][z] * pmlHYXBottomXZ[x][y+pmlWidth][z] + pmlHYZBottomXZ[x][y+pmlWidth][z] * pmlHYZBottomXZ[x][y+pmlWidth][z] + pmlHZXBottomXZ[x][y+pmlWidth][z] * pmlHZXBottomXZ[x][y+pmlWidth][z] + pmlHZYBottomXZ[x][y+pmlWidth][z] * pmlHZYBottomXZ[x][y+pmlWidth][z]);break;
	    1: result = pmlHXYBottomXZ[x][y+pmlWidth][z] + pmlHXZBottomXZ[x][y+pmlWidth][z];break;
	    2: result = pmlHYXBottomXZ[x][y+pmlWidth][z] + pmlHYZBottomXZ[x][y+pmlWidth][z];break;
	    3: result = pmlHZXBottomXZ[x][y+pmlWidth][z] + pmlHZYBottomXZ[x][y+pmlWidth][z];break;
	    4: result = pmlEXYBottomXZ[x][y+pmlWidth][z] + pmlEXZBottomXZ[x][y+pmlWidth][z];break;
	    5: result = pmlEYXBottomXZ[x][y+pmlWidth][z] + pmlEYZBottomXZ[x][y+pmlWidth][z];break;
	    6: result = pmlEZXBottomXZ[x][y+pmlWidth][z] + pmlEZYBottomXZ[x][y+pmlWidth][z];break;
	  default: result = sqrt(pmlEXYBottomXZ[x][y+pmlWidth][z] * pmlEXYBottomXZ[x][y+pmlWidth][z] + pmlEXZBottomXZ[x][y+pmlWidth][z] * pmlEXZBottomXZ[x][y+pmlWidth][z] + pmlEYXBottomXZ[x][y+pmlWidth][z] * pmlEYXBottomXZ[x][y+pmlWidth][z] + pmlEYZBottomXZ[x][y+pmlWidth][z] * pmlEYZBottomXZ[x][y+pmlWidth][z] + pmlEZXBottomXZ[x][y+pmlWidth][z] * pmlEZXBottomXZ[x][y+pmlWidth][z] + pmlEZYBottomXZ[x][y+pmlWidth][z] * pmlEZYBottomXZ[x][y+pmlWidth][z]);
	  }
	else if(y < ny)
	  switch (mode){ /* Sim Space */
	    0: result = sqrt(hx[x][y][k] * hx[x][y][z] + hy[x][y][k] * hy[x][y][z] + hz[x][y][k] * hz[x][y][z]); break;
	    1: result = hx[x][y][k]; break;
	    2: result = hy[x][y][k]; break;
	    3: result = hz[x][y][k]; break;
	    4: result = ex[x][y][k]; break;
	    5: result = ey[x][y][k]; break;
	    6: result = ez[x][y][k]; break;
	  default: result = sqrt(ex[x][y][z] * ex[x][y][z] + ey[x][y][z] * ey[x][y][z] + ez[x][y][z] * ez[x][y][z]);
	  }
	else
	  switch (mode) { /* TopXZ */
	    0: result = sqrt(pmlHXYTopXZ[x][y-ny][z] * pmlHXYTopXZ[x][y-ny][z] + pmlHXZTopXZ[x][y-ny][z] * pmlHXZTopXZ[x][y-ny][z] + pmlHYXTopXZ[x][y-ny][z] * pmlHYXTopXZ[x][y-ny][z] + pmlHYZTopXZ[x][y-ny][z] * pmlHYZTopXZ[x][y-ny][z] + pmlHZXTopXZ[x][y-ny][z] * pmlHZXTopXZ[x][y-ny][z] + pmlHZYTopXZ[x][y-ny][z] * pmlHZYTopXZ[x][y-ny][z]);break;
	    1: result = pmlHXYTopXZ[x][y-ny][z] + pmlHXZTopXZ[x][y-ny][z];break;
	    2: result = pmlHYXTopXZ[x][y-ny][z] + pmlHYZTopXZ[x][y-ny][z];break;
	    3: result = pmlHZXTopXZ[x][y-ny][z] + pmlHZYTopXZ[x][y-ny][z];break;
	    4: result = pmlEXYTopXZ[x][y-ny][z] + pmlEXZTopXZ[x][y-ny][z];break;
	    5: result = pmlEYXTopXZ[x][y-ny][z] + pmlEYZTopXZ[x][y-ny][z];break;
	    6: result = pmlEZXTopXZ[x][y-ny][z] + pmlEZYTopXZ[x][y-ny][z];break;
	  default: result = sqrt(pmlEXYTopXZ[x][y-ny][z] * pmlEXYTopXZ[x][y-ny][z] + pmlEXZTopXZ[x][y-ny][z] * pmlEXZTopXZ[x][y-ny][z] + pmlEYXTopXZ[x][y-ny][z] * pmlEYXTopXZ[x][y-ny][z] + pmlEYZTopXZ[x][y-ny][z] * pmlEYZTopXZ[x][y-ny][z] + pmlEZXTopXZ[x][y-ny][z] * pmlEZXTopXZ[x][y-ny][z] + pmlEZYTopXZ[x][y-ny][z] * pmlEZYTopXZ[x][y-ny][z]);
	  }
      else
	switch (mode) { /* TopXY */
	0: result = sqrt(pmlHXYTopXY[x][y+pmlWidth][z-nz] * pmlHXYTopXY[x][y+pmlWidth][z-nz] + pmlHXZTopXY[x][y+pmlWidth][z-nz] * pmlHXZTopXY[x][y+pmlWidth][z-nz] + pmlHYXTopXY[x][y+pmlWidth][z-nz] * pmlHYXTopXY[x][y+pmlWidth][z-nz] + pmlHYZTopXY[x][y+pmlWidth][z-nz] * pmlHYZTopXY[x][y+pmlWidth][z-nz] + pmlHZXTopXY[x][y+pmlWidth][z-nz] * pmlHZXTopXY[x][y+pmlWidth][z-nz] + pmlHZYTopXY[x][y+pmlWidth][z-nz] * pmlHZYTopXY[x][y+pmlWidth][z-nz]);break;
	1: result = pmlHXYTopXY[x][y+pmlWidth][z-nz] + pmlHXZTopXY[x][y+pmlWidth][z-nz];break;
	2: result = pmlHYXTopXY[x][y+pmlWidth][z-nz] + pmlHYZTopXY[x][y+pmlWidth][z-nz];break;
	3: result = pmlHZXTopXY[x][y+pmlWidth][z-nz] + pmlHZYTopXY[x][y+pmlWidth][z-nz];break;
	4: result = pmlEXYTopXY[x][y+pmlWidth][z-nz] + pmlEXZTopXY[x][y+pmlWidth][z-nz];break;
	5: result = pmlEYXTopXY[x][y+pmlWidth][z-nz] + pmlEYZTopXY[x][y+pmlWidth][z-nz];break;
	6: result = pmlEZXTopXY[x][y+pmlWidth][z-nz] + pmlEZYTopXY[x][y+pmlWidth][z-nz];break;
      default: result = sqrt(pmlEXYTopXY[x][y+pmlWidth][z-nz] * pmlEXYTopXY[x][y+pmlWidth][z-nz] + pmlEXZTopXY[x][y+pmlWidth][z-nz] * pmlEXZTopXY[x][y+pmlWidth][z-nz] + pmlEYXTopXY[x][y+pmlWidth][z-nz] * pmlEYXTopXY[x][y+pmlWidth][z-nz] + pmlEYZTopXY[x][y+pmlWidth][z-nz] * pmlEYZTopXY[x][y+pmlWidth][z-nz] + pmlEZXTopXY[x][y+pmlWidth][z-nz] * pmlEZXTopXY[x][y+pmlWidth][z-nz] + pmlEZYTopXY[x][y+pmlWidth][z-nz] * pmlEZYTopXY[x][y+pmlWidth][z-nz]);
      }
  else
    switch (mode) { /* TopYZ */
      0: result = sqrt(pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth]);break;
      1: result = pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
      2: result = pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
      3: result = pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
      4: result = pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
      5: result = pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
      6: result = pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
    default: result = sqrt(pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth]);
    }
  return result;
}