FDTD_2D_Maxwell/tfsftmz.c

/*
    2D FDTD simulator
    Copyright (C) 2019 Emilia Blåsten

    This program is free software: you can redistribute it and/or
    modify it under the terms of the GNU Affero General Public License
    as published by the Free Software Foundation, either version 3 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
    Affero General Public License for more details.

    You should have received a copy of the GNU Affero General Public
    License along with this program.  If not, see
    <http://www.gnu.org/licenses/>.
*/
/* tfsftmz.c: Source code to implement a TFSF boundary for a TMz grid. 
 * The incident field is assumed to propagate along the x-direction and 
 * is calculated using an auxiliary 1D simulation. */

#include <stdio.h>
#include <stdlib.h>
#include "tfsftmz.h"
#include "updatetmz.h"
#include "ezinc.h"
#include "grid1dez.h"
#include "gridtmz.h"
#include "incident.h"

void tfsfInit(Tfsf *tf, Grid *g) {
  // initialize 1d grid
  gridInit1d(tf->g1, g->time, g->maxTime, g->cdtds);

  tf->firstX = 10;
  tf->firstY = 10;
  tf->lastX = g->sizeX-10;
  tf->lastY = g->sizeY-10;

  ezIncInit(g);                    // initialize source function

  return;
}

void tfsfUpdate(Tfsf *tf, Grid *g) {
  int mm, nn;

  // check if tfsfInit() has been called
  if( tf->firstX <= 0 || tf->firstY <=0 || tf->lastX >= g->sizeX ||
      tf->lastY >= g->sizeY || tf->firstX >= tf->lastX || tf->firstY >=
      tf->lastY ) {
    fprintf(stderr,
      "tfsfUpdate: tfsfInit must be called before tfsfUpdate.\n"
      "            Boundary location must be set to reasonable values.\n");
    exit(-1);
  }

  // correct Hy along left edge
  mm = tf->firstX - 1;
  for(nn = tf->firstY; nn <= tf->lastY; nn++)
    g->hy[mm][nn] -= g->hyUpdateEcoeff[mm][nn] * tf->g1->ez[mm+1];

  // correct Hy along right edge
  mm = tf->lastX;
  for(nn = tf->firstY; nn <= tf->lastY; nn++)
    g->hy[mm][nn] += g->hyUpdateEcoeff[mm][nn] * tf->g1->ez[mm];

  // correct Hx along the bottom
  nn = tf->firstY - 1;
  for(mm = tf->firstX; mm <= tf->lastX; mm++)
    g->hx[mm][nn] += g->hxUpdateEcoeff[mm][nn] * tf->g1->ez[mm];

  // correct Hx along the top
  nn = tf->lastY;
  for(mm = tf->firstX; mm <= tf->lastX; mm++)
    g->hx[mm][nn] -= g->hxUpdateEcoeff[mm][nn] * tf->g1->ez[mm];

  updateH1d(tf->g1); // update 1D magnetic field
  updateE1d(tf->g1); // update 1D electric field

  // incident wave
  tf->g1->ez[0] = incidentWaveOnLeftSide(tf->g1->time, tf->g1->cdtds);
  tf->g1->time++;

  // correct Ez adjacent to TFSF boudary
  // correct Ez field along left edge
  mm = tf->firstX;
  for(nn = tf->firstY; nn <= tf->lastY; nn++)
    g->ez[mm][nn] -= g->ezUpdateHcoeff[mm][nn] * tf->g1->hy[mm-1];

  // correct Ez field along right edge
  mm = tf->lastX;
  for(nn = tf->firstY; nn <= tf->lastY; nn++)
    g->ez[mm][nn] += g->ezUpdateHcoeff[mm][nn] * tf->g1->hy[mm];

  // no need to correct Ez along top and bottom
  // since incident Hx is zero (so need to be done
  // if the incident field has nonzero Hx)!!

  return;
}

Tfsf *tfsfCreate(int sizeX) {
  Tfsf *tf;
  tf = (Tfsf *)calloc(1, sizeof(Tfsf));
  if(!tf) {
    perror("tfsfCreate()");
    fprintf(stderr, "Failed to allocate memory for Tfsf.\n");
    exit(-1);
  }

  // Allocate memory for 1D Grid
  tf->g1 = gridCreate1d(sizeX);

  return tf;
}

void tfsfDestroy(Tfsf *tf) {
  gridDestroy1d(tf->g1);
  free(tf);
  return;
}