I am using a template from Prof John B Schneider's book on understanding FDTD, with a minor modification. Increase the size of the E-field array by 1. The end edge of the E-field is now set to 0. Don't recompute this. Compute the H-field all the way to the end, based on this outer E-field.
Reference:
- Understanding the Finite-Difference Time-Domain Method, John B. Schneider, www.eecs.wsu.edu/~schneidj/ufdtd, 2010.
The below code is a modification of Program 3.1. Program 3.1 had a PEC at x=0, a PMC at x=200. The code below has a PEC at x=0, a PEC at x=201.
#include <stdio.h>
#include <math.h>
#define SIZE 200
int main()
{
double ez[SIZE+1] = {0.}, hy[SIZE] = {0.}, imp0=377.0;
int qTime, maxTime = 1000, mm;
char basename[80] = "sim", filename[100];
int frame=0;
FILE* snapshot;
/* do time stepping */
for (qTime = 0; qTime < maxTime; qTime++) {
/* Magnetic field */
for(mm=0; mm < SIZE; mm++)
hy[mm] = hy[mm] + (ez[mm+1]-ez[mm])/imp0;
/* Electric field */
for(mm=1; mm < SIZE; mm++)
ez[mm] = ez[mm] + (hy[mm]-hy[mm-1]) * imp0;
/* Additive current source at node 50 */
ez[50] += exp(-(qTime -30.) * (qTime - 30.) / 100.);
// Explicit forcing function at node 0
// ez[0] = exp(-(qTime - 30.) * (qTime - 30.) / 100.);
if (qTime % 10 == 0) {
sprintf(filename, "%s.%d", basename, frame++);
snapshot = fopen(filename, "w");
for(mm=0; mm<SIZE+1; mm++)
fprintf(snapshot, "%g\n", ez[mm]);
fclose(snapshot);
}
}
return 0;
}
