FDTD_2D_Maxwell/grid1dez.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/>.
*/
/* grid1dez.c: Initialization function for the 1D auxiliary grid used
 * by the TFSF function to calculate the incident field. */

#include <math.h>
#include <stdio.h>
#include <stdlib.h> //for calloc
#include "grid1dez.h"

#define NLOSS    20    // number of lossy cells at the end of 1D grid
#define MAX_LOSS 0.35  // maximum loss factor in lossy layer

void gridInit1d(Grid1D *g, int time, int maxTime, double cdtds) {
  double imp0 = 377.0, depthInLayer, lossFactor;
  g->time = time;
  g->maxTime = maxTime;
  g->cdtds = cdtds;

  int mm;
  // set the electric- and magnetic-field update coefficients
  for(mm = 0; mm < g->sizeX - 1; mm++) {
    if( mm < g->sizeX - 1 - NLOSS) {
      g->ezUpdateEcoeff[mm] = 1.0;
      g->ezUpdateHcoeff[mm] = g->cdtds * imp0;
      g->hyUpdateHcoeff[mm] = 1.0;
      g->hyUpdateEcoeff[mm] = g->cdtds / imp0;
    } else {
      depthInLayer = mm - (g->sizeX - 1 - NLOSS) + 0.5;
      lossFactor = MAX_LOSS * pow(depthInLayer / NLOSS, 2);
      g->ezUpdateEcoeff[mm] = (1.0 - lossFactor) / (1.0 + lossFactor);
      g->ezUpdateHcoeff[mm] = g->cdtds * imp0 / (1.0 + lossFactor);
      depthInLayer += 0.5;
      lossFactor = MAX_LOSS * pow(depthInLayer / NLOSS, 2);
      g->hyUpdateHcoeff[mm] = (1.0 - lossFactor) / (1.0 + lossFactor);
      g->hyUpdateEcoeff[mm] = g->cdtds / imp0 / (1.0 + lossFactor);
      }
    }

  return;
}


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

  g->sizeX = sizeX + NLOSS;  // size of domain

  // Allocate memory for the arrays
  g->hy             = (double *)calloc(g->sizeX-1, sizeof(double));
  g->hyUpdateEcoeff = (double *)calloc(g->sizeX-1, sizeof(double));
  g->hyUpdateHcoeff = (double *)calloc(g->sizeX-1, sizeof(double));
  g->ez             = (double *)calloc(g->sizeX, sizeof(double));
  g->ezUpdateEcoeff = (double *)calloc(g->sizeX, sizeof(double));
  g->ezUpdateHcoeff = (double *)calloc(g->sizeX, sizeof(double));

  // Check memory allocated successfully
  if(!g->hy || !g->hyUpdateEcoeff || !g->hyUpdateHcoeff ||
     !g->ez || !g->ezUpdateEcoeff || !g->ezUpdateHcoeff) {
    perror("gridCreate1d()");
    fprintf(stderr, "Failed to allocate memory to Grid1D arrays.\n");
    exit(-1);
  }

  return g;
}


void gridDestroy1d(Grid1D *g) {
  free(g->hy);
  free(g->hyUpdateEcoeff);
  free(g->hyUpdateHcoeff);
  free(g->ez);
  free(g->ezUpdateEcoeff);
  free(g->ezUpdateHcoeff);
  free(g);

  return;
}