Project 291 - John Benjamin

Two dimensional simulation of electric field magnitude for dielectric device

In this project the propegation behaviour of electric fields were investigated in a system which has a dielectric and absorbing boundary. The investigated system was an infinite rectangular rod, using such a construct, the dimensionality of the governing equations (Maxwell's equations) can be reduced. The infinite rod is shown in the figure below:

1. The procedure and approach to implementation of Maxwell's equation

The approach in the finite difference time domain (FDTD) is determined by the style of materials used. Since perfectly matched layers (PML) is included in the model to absorb and remove energy from the system, the Maxwell equations has to be considered in the frequency domain. PML is strictly a frequency domain property, and has to be moved to the time domain using inverse Laplace transforms. Using inverse Laplace transforms the frequency domain is then moved to the time domain in the following manner:

$$j \omega [\,\,\,] = \frac{\partial}{\partial t}[\,\,\,]$$

and,

$$\frac{1}{j \omega} [\,\,\,] = \int_{0}^{t} [\,\,\,]dt'$$

The PML functions as an absorbing boundary placed at the edge of the problem space to simulate waves leaving the system. This is done by matching electromagnetic impedance with the problem space while gradually increasing the conductivity as the wave propegate deeper into the PML. The interal structure is shown below:

The dielectric has the following properties: Diagonally anisotropic permittivity \( [ \epsilon_{r} ] \), anisotropic permeability \( [ \mu_{r} ] \), and anisotropic conductivity \( [ \sigma ] \). The brackets represent that the quantities are forming a matrix when interacting with the field vectors of Maxwells equation.

1.1 Maxwell's equation from 3D to 2D

The Maxwell's curl equations (and constitutive relations) in frequency space with PML matrix \( [ S ] \) is given by:

$$\nabla \times E(\omega) = -j \omega \mu_{0} [\mu_{r}][S] H(\omega)$$ $$\nabla \times H(\omega) = [\sigma]E(\omega) + j \omega [S] D(\omega)$$ $$D(\omega) = \epsilon_{0} [\epsilon_{r}] E(\omega)$$

where the Electric field is given by E, the magnetic field by H, the displacement field by D, the vacuum permittivity \( \epsilon_{0} \), and the vacuum permeability \( \mu_{0} \). Before reducing the dimensions, the electric field is normalized: $$\tilde{E}=\sqrt{ \frac{\epsilon_{0}}{\mu_{0}} } E=\frac{1}{\eta_{0}} E,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, \tilde{D}=\sqrt{ \frac{1}{ \epsilon_{0} \mu_{0}} }=c_{0}D$$ Using a Cartesian coordinate system we note that the rod is extended infinitely far in the z-direction. This implies that all the functions are constant when displaced in the z-direction, i.e. \( \partial / \partial z = 0\). Expanding the Maxwell's equation in terms of its cartesian coordinate system using the vectors x,y,z:

$$C_{x}^{E}=\frac{\partial \tilde{E}_{z}}{\partial y}=-j\omega \big( \frac{\mu_{rx}}{c_{0}} S_{y} H_{x} \big)$$ $$C_{y}^{E}=-\frac{\partial \tilde{E}_{z}}{\partial x}=-j\omega \big( \frac{\mu_{ry}}{c_{0}} S_{y}^{-1} H_{y} \big)$$ $$C_{z}^{H}=\frac{\partial H_{y}}{\partial x}-\frac{\partial H_{x}}{\partial x}=(\eta_{0} \sigma_{z})\tilde{E}_{z} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{z} \big)$$

and,

$$C_{x}^{H}=\frac{\partial H_{z}}{\partial y}=(\eta_{0}\sigma_{x}) \tilde{E}_{x} + j\omega \big( \frac{1}{c_{0}} S_{y} \tilde{D}_{x} \big) $$ $$C_{y}^{H}=-\frac{\partial H_{z}}{\partial x}=(\eta_{0}\sigma_{y}) \tilde{E}_{y} + j\omega \big( \frac{1}{c_{0}} S_{y}^{-1} \tilde{D}_{y} \big) $$ $$C_{z}^{E}=\frac{\partial \tilde{E}_{y}}{\partial x}-\frac{\partial \tilde{E}_{x}}{\partial x}=-j\omega \big( \frac{\mu_{rz}}{c_{0}} S_{y} H_{z} \big)$$

for the PML oriented such that it absorbes radiation in the y-direction:

$$ [ S ] =\begin{bmatrix} S_{y} & 0 & 0\\0 & S_{y}^{-1} & 0\\0 &0 &S_{y}\end{bmatrix},\,\,\,\,\,\,\,\, S_{y}=1 + \frac{1}{j \omega} \frac{\sigma'(y)}{\epsilon_{0}},\,\,\,\,\,\,\,\, \sigma'(y)= \frac{\epsilon_{0}}{2 \Delta t} \big( \frac{y'}{L_{y}} \big)^{3}$$

where \( L_{y}\) is the length of the PML, and \( y'\) is the current depth within the PML. Outside of the PML, \( [S] =1 \). The choice of the max value of \( [S] \) is usually the one shown above, but will depend on the exact parameters (source intensity, geometry of the chamber, material distribution, etc...), and is usually tweaked untill it produces a stable absorption. Selecting the wrong value of the max value can make the system numerically unstable. From the expanded curl equations it is noticed system becomes seperated into two independent sub-system of differential equations. To investigate the wave properties one of these modes can then be selected (given that the source is polarized). We select the mode corresponding to \( \{E_{x} (x,y),E_{y}(x,y),H_{z}(x,y)\} \). Moving the above equations to the time domain using the inverse Laplace transform:

$$C_{z}^{E}|_{t}=-\big( \frac{\mu_{ry}}{c_{0}} \big)\frac{\partial H_{z}}{\partial t}|_{t} - \big( \frac{\mu_{ry}\sigma'}{c_{0}\epsilon_{0}} \big) H_{z}|_{t}$$ $$C_{x}^{H}|_{t}=(\eta_{0}\sigma_{x})\tilde{E}_{x}|_{t} + \big( \frac{1}{c_{0}} \big)\frac{\partial \tilde{D}_{x}}{\partial t}|_{t} + \big( \frac{\sigma'}{c_{0}\epsilon_{0}} \big)\tilde{D}_{x}|_{t}$$ $$C_{y}^{H}|_{t} + \big( \frac{\sigma'}{\epsilon_{0}} \big) \int_{0}^{t} C_{y}^{H} dt = (\eta_{0}\sigma_{y})\tilde{E}_{y}|_{t} + \big( \frac{\eta_{0}\sigma_{y}\sigma'}{\epsilon_0} \big) \int_{0}^{t}\tilde{E}_{y}|_{t} dt + \big( \frac{1}{c_{0}} \big) \frac{\partial \tilde{D}_{y}}{\partial t}|_{t}$$

1.2a From a continuous model to a discrete model: Temporal discretization

The electric field and the magnetic field is temporally displaced relative to eachother. Using the electric field at \( t_{n} \), the magnetic field defined at the next half step in time can then be calculated \( t_{n}+\frac{\Delta t}{2} \). Using this magnetic field the electric field can be at the time step \( t_{n}+\Delta t =t_{n+1} \). This is how the temporal displacements work and how the fields evolve with time when in discrete form. Using finite differences the differentials and the integrals can be substituted with differences and sums, and from the differentials the next time-steps can be evaluated. The update equations are then given by:

$$H_{z}|_{t+\frac{\Delta t}{2}}=-h_{z1}C_{z}^{E}|_{t} + h_{z2}H_{z}|_{t-\frac{\Delta t}{2}}$$ $$\tilde{D}_{x}|_{t+\Delta t}= d_{x1}C_{x}^{H}|_{t+\frac{\Delta t}{2}} + d_{x1}C_{x}^{H}|_{t-\frac{\Delta t}{2}} - d_{x2}\tilde{E}_{x}|_{t} + \tilde{D}_{x}|_{t-\Delta t} - d_{x3}\tilde{D}_{x}|_{t}$$ $$\tilde{D}_{y}|_{t+\Delta t}= d_{y1}C_{y}^{H}|_{t+\frac{\Delta t}{2}} + d_{y1}C_{y}^{H}|_{t-\frac{\Delta t}{2}} + d_{y2} \sum_{t_{i}=\frac{\Delta t}{2}}^{t-\frac{\Delta t}{2}}C_{y}^{H}|_{t_{i}} - d_{x3}\tilde{E}_{y}|_{t} -d_{y4} \sum_{t_{i}=0}^{t}\tilde{E}_{y}|_{t_{i}} + \tilde{D}_{y}|_{t-\Delta t} $$

Note that linear interpolation has to be used on the time steps that are not defined. The full derivation of these update equations are given here. The constant parameters shown in the update equation above are given by:

$$h_{z1}=\frac{2c_{0}\epsilon_{0}\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}, \,\,\,\, h_{z2}=\frac{2\mu_{ry}\epsilon_{0}-\mu_{ry}\sigma'\Delta t}{2\mu_{ry}\epsilon_{0}+\mu_{ry}\sigma'\Delta t}$$ $$d_{x1}=c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{x2}=2c_{0}\eta_{0}\sigma_{x}\Delta t, \,\,\,\,\,\,\,\,\,\,\, d_{3}=\frac{2\sigma'\Delta t}{\epsilon_{0}}$$ $$d_{y1}=2c_{0}\Delta t\big(\frac{1}{2} + \frac{\sigma' \Delta t}{4 \epsilon_{0}} \big), \,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y2}=2\sigma' \eta_{0}(c_{0}\Delta t)^{2}$$ $$d_{y3}=2\sigma_{y} \eta_{0}c_{0}\Delta t, \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\, d_{y4}=2\sigma_{y}\sigma' (\eta_{0}c_{0}\Delta t)^{2}$$

1.2b From a continuous model to a discrete model: Spatially discretization

To discreteze space a Yee grid was used. A Yee grid is a spacial discretization where the electric fields and the magnetic fields are not defined in the same position, and thus spatially displaced from eachother. The motivation for using a Yee grid is that (1) the divergence equations are intrinsically fulfilled, (2) the material boundary conditions are intrinsically fulfilled. The Yee grid for the 2D system is show below:

The curl terms, \( \{C_{z}^{E},C_{x}^{H},C_{y}^{H}\} \) can then be written as:

$$C_{z}^{E}(i_{x},i_{y})=\frac{\tilde{E}_{y}(i_{x}+1,i_{y})-\tilde{E}_{y}(i_{x},i_{y})}{\Delta x}-\frac{\tilde{E}_{x}(i_{x},i_{y}+1)-\tilde{E}_{x}(i_{x},i_{y})}{\Delta y}$$ $$=c_{1}\tilde{E}_{y}(i_{x}+1,i_{y})-c_{1}\tilde{E}_{y}(i_{x},i_{y})-c_{2}\tilde{E}_{x}(i_{x},i_{y}+1)+c_{2}\tilde{E}_{x}(i_{x},i_{y})$$ $$C_{x}^{H}(i_{x},i_{y})=\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x},i_{y}-1)}{\Delta y}=c_{2}H_{z}(i_{x},i_{y})-c_{2}H_{z}(i_{x},i_{y}-1)$$ $$C_{y}^{H}(i_{x},i_{y})=-\frac{H_{z}(i_{x},i_{y})-H_{z}(i_{x}-1,i_{y})}{\Delta x}=-c_{1}H_{z}(i_{x},i_{y})+c_{1}H_{z}(i_{x}-1,i_{y})$$

where, \(c_{1}=\frac{1}{\Delta x}\) and \(c_{2}=\frac{1}{\Delta y}\).

1.3 The spatial and the temporal resolution

The general spatial resolution \( ds \) is given by the following:

$$ds:=\frac{\lambda_{0}}{10}=dx=dy,$$ where \( \lambda_{0} \) is the wavelengh produced by the source while in the vacuum. The selection of the value for the divisor is the main component for a high resolution result, however, as seen later also influces some visual errors. The temporal resolution \( dt \) is then given by the time light uses to propegate (since the propagation is restricted to 2D the fastest propegation is diagonally, which means it has to be scaled by \( (\sqrt{2})^{-1} \)): $$dt=\frac{ds}{\sqrt{2}c_{0}}.$$

But from certain books is also written as,

$$dt=\frac{1}{c_{0}\sqrt{\big(\frac{1}{\Delta x}\big)^{2}+\big(\frac{1}{\Delta y}\big)^{2}}}$$

1.4 Summary of the iteration procedure

The iteration process is shown in the figure below, and uses the formulas shown above.

For the full derivation of the theory see the following PDF form: [Theory.pdf]

2. The code constructed by C++ and ROOT

The simulation code was constructed using C++ and Cern ROOT. This code is shown below:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
#include <iostream>
#include <cmath>
#include <string>

using namespace std;

double** allocMatrix(int dimX,int dimY);
void initializeArrayZero(double **arr,double dimX,double dimY);
void initializeArrayOne(double **arr,double dimX,double dimY);
double findMaxfromArray(double **arr,int dimX,int dimY);

void Simulation(){
	// Input parameters:
	double x_max = 1.0;                     // Size of the chamber, x-direction
	double y_max = 1.0;			// Size of the chamber, y-direction
	double E0 = 1.0;                        // (Norm.) Electric field amplitude
	double frequency = 2.4 * pow(10.0,9.0); // Frequency of the source
	double diameter_dielectric = 0.4;       // Diameter dielectric, ratio of chamber
	double source_diameter = 0.5;           // Diameter EM source, ratio of chamber
	double source_xpos = 0.5;               // Center position EM source, ratio of chamber
	double epsilonrx = 2;                 // Relative permittivity of dielectric, x-direction
	double epsilonry = 2;                 // Relative permittivity of dielectric, y-direction
	double mury = 1.0;                      // Relative permeability of dielectric, y-direction
	double sigmax = 0.00;                   // Conductivity of dielectric, x-direction
	double sigmay = 0.00;                   // Conductivity of dielectric, y-direction
	double ntime = 500;                     // Desired time steps
	double PML_length_frac = 0.4;           // Thickness of absorbing boundary relative to the heigh
	double PML_strength    = 1.0;           // Strength of absorbing boundary. 0 = No PML
						// The larger PML strength is, the more absorption occurs, however it is less numerically stable
	const TString kGifName = "PML.gif";	// Used for GIF. The name of the gif that shows the result. Is placed in the same directory as the script
	const TString kGifSleep = "3";     	// Used for GIF. Time between each frame in the GIF. Given in centiseconds
	 
	// (1.0) Constructing the constants:
	double epsilon0 = 8.8541878128 * pow(10.0,-12.0);   // Free space permittivity
	double mu0 = 1.25663706212 * pow(10.0,-6.0);        // Free space permeability
	double c0 = sqrt(1/(epsilon0*mu0));                 // Free space propegation velocity (used in the resolution and normlization)
	double eta0 = sqrt(mu0/epsilon0);                   // Free space impedence (used in the normalization)
	double cMat = 0;
	if (epsilonrx > epsilonry){
		cMat = c0 * sqrt(1/(epsilonrx*mury));
	}
	else{
		cMat = c0 * sqrt(1/(epsilonry*mury));
	}
	double lambdaMat = cMat/frequency; 		    // Material wavelength  (used in the resolution)
	double H0 = (eta0)/(c0)*E0;                         // For the given (Norm.) E-Field, the associated propagating magnetic field has magnitude H0

	// (1.1) Resolution of the model: dt,dx,dy
	double ds = lambdaMat/10;                           // Spatial resolution of the model
	// Note: the following is the round function. My version does not support round, so this is a workaround given the input is positive
	int nx = int(std::floor(x_max/ds + 0.5));           // Getting the most snug integer fit with the partitions given by the length ds.
	int ny = int(std::floor(y_max/ds + 0.5));           // ------------------\\--------------------
	x_max = nx*ds;                                      // Redefining the size of the chamber such that they form a snug integer fit for the grid
	y_max = ny*ds;					    // Redefining the size of the chamber such that they form a snug integer fit for the grid
	double dx = x_max/nx;
	double dy = y_max/ny;
	double dt = 1/(c0*sqrt(pow(dx,-2)+pow(dy,-2)));      // Courant's stability criterion

	// (1.2) The key parameters related to the (integer) indices of the PML, source, and the dielectric:
	// For the PML:
	int nPML = int(std::floor(PML_length_frac*ny + 0.5));
	// For the Dielectric:
	int ny_dielectric = int(std::floor(0.5*ny + 0.5)); 						// Setting the dielectric in the middle of the chamber, y-direction
	int nx_dielectric = int(std::floor(0.5*nx + 0.5)); 						// Setting the dielectric in the middle of the chamber, x-direction
	double n_dielectric_radius = 0;
	double r_dielectric_radius = 0;
	if (dx > dy){
		n_dielectric_radius = 0.5*diameter_dielectric*nx;
		r_dielectric_radius = dx*n_dielectric_radius; // Used for the figure	
	}
	else{
		n_dielectric_radius = 0.5*diameter_dielectric*ny;
		r_dielectric_radius = dy*n_dielectric_radius;  // Used for the figure
	}
	double n_dielectric_radius_square = pow(n_dielectric_radius,2);
	// For the Source:
	int iblocks_source_radius; double ix_source_pos; int ix_source_start; int ix_source_end;
	if (nx % 2 ==1){ // The if and the else statement center the source properly, such that if source_xpos is 0.5 it is symmetric regardless of if nx is odd or even
		iblocks_source_radius = int(std::floor(0.5*source_diameter*nx + 0.5));     // Approximating the radius by number of index steps
		ix_source_pos = int(std::floor(source_xpos*nx + 0.5));                     // Finding the index where the center of the source is located
		ix_source_start = ix_source_pos-iblocks_source_radius;                     // Start position of the source
		ix_source_end   = ix_source_pos+iblocks_source_radius;                     // End position of the source
	}
	else{
		iblocks_source_radius = int(std::floor(0.5*source_diameter*nx + 0.5));     	  // Approximating the radius by number of index steps
		ix_source_pos = source_xpos*(nx+1);                     			  // Finding the index where the center of the source is located
		ix_source_start = int(std::floor(ix_source_pos-iblocks_source_radius+0.5 +0.5))+1;  // Start position of the source
		ix_source_end   = int(std::floor(ix_source_pos+iblocks_source_radius-0.5 +0.5))+1;  // End position of the source
	}

	// (1.3) The iteration constants that are not needed in matrix form.
	double c1 = 1/dx;
	double c2 = 1/dy;
	double dx1 = c0*dt;


	// (2.0) Allocating the iteration constants (for those needed in matrix form).
	int dimX = nx+2;             			     // Size of the problem space in the x-dimension +2 Due to the boundary on both sides
	int dimY = ny+nPML+2;          			     // Size of the problem space in the y-dimension +2 Due to the boundary on both sides
	double** erxx         = allocMatrix(dimX,dimY);      // Relative permittivity in the x-direction for all spatial coordinates
	double** eryy         = allocMatrix(dimX,dimY);      // Relative permittivity in the y-direction for all spatial coordinates 
	double** muryy        = allocMatrix(dimX,dimY);      // Relative permeability in the y-direction for all spatial coordinates              
	double** sigmaxx      = allocMatrix(dimX,dimY);      // Conductivity of the medium in every spatial coordinate (anisotropy:x-direction)
	double** sigmayy      = allocMatrix(dimX,dimY);      // Conductivity of the medium in every spatial coordinate (anisotropy:y-direction)
	double** sigmaYYprime = allocMatrix(dimX,dimY);      // Fictitious conductivity causing loss (absorption) in the PML region in the top of the chamber
	

	// (2.1) Setting the standard values based on free space
	initializeArrayOne(erxx, dimX, dimY);			// Is the relative permittivity (x-dirr.) is 1 in vacuum
	initializeArrayOne(eryy, dimX, dimY);			// Is the relative permittivity (y-dirr.) is 1 in vacuum
	initializeArrayOne(muryy, dimX, dimY);			// Is the relative permeability is 1 in vacuum
	initializeArrayZero(sigmaxx, dimX, dimY);		// The conductivity of vacuum is equal to 0
	initializeArrayZero(sigmayy, dimX, dimY);		// The conductivity of vacuum is equal to 0
	initializeArrayZero(sigmaYYprime, dimX, dimY);		// The ficticious conductivity of the PML is zero except in the PML region
	

	// (2.2) Filling in the values when in the dielectric:
	// There are alternative and quicker ways to create a circular dielectric element, but we will do it the intuitive way.
	// Using the circle formula.
	for (int iy=1; iy<dimY-1;iy++){
		int square_distanceY = (iy-ny_dielectric)*(iy-ny_dielectric);
		for (int ix=0; ix<dimX-1;ix++){
			int square_distanceX = (ix-nx_dielectric)*(ix-nx_dielectric);	
			if (square_distanceX+square_distanceY < n_dielectric_radius_square){
				erxx[ix][iy] = 1/epsilonrx;					// Since we first find the D-field, we need the inverse transformation to find the E-field
				eryy[ix][iy] = 1/epsilonry;					// Since we first find the D-field, we need the inverse transformation to find the E-field
				muryy[ix][iy] = mury;
				sigmaxx[ix][iy] = sigmax;
				sigmayy[ix][iy] = sigmay;
			}
		}
	}


	// (2.3) Filling in the values when in the PML (absorbing material):
	for (int iy=ny+1; iy<dimY-1;iy++){
		float PML_Depth = (iy-(ny+1))/nPML;     // Fraction of depth into the PML (0 to 1)
		float PML_shape = pow(PML_Depth,3);     // Determines the shape. Tests found that (fraction)^3 is the best shape
		float sigmaYYprime_iy = PML_strength * ((epsilon0)/(2*dt)) * PML_shape; 	// The ficticious conductivity found in the PML region
		for (int ix=1; ix<dimX-1;ix++){
			sigmaYYprime[ix][iy]=sigmaYYprime_iy;					// The ficticious conductivity found in the PML region
		}
	}


	// (2.4) Constructing the iteration constants (matrix)
	double** hz1 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Hz
	double** hz2 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Hz
	double** dx2 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dx             
	double** dx3 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dx
	double** dy1 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dy
	double** dy2 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dy
	double** dy3 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dy
	double** dy4 = allocMatrix(dimX,dimY);             // Used in updating (iterating) Dy


	// (2.5) Inserting the values into the coefficient matrices
	for (int iy=1; iy<dimY-1;iy++){
		for (int ix=1; ix<dimX-1;ix++){
			hz1[ix][iy] = (2*c0*epsilon0*dt)/((2*(muryy[ix][iy])*epsilon0)+((muryy[ix][iy])*(sigmaYYprime[ix][iy])*dt));
			hz2[ix][iy] = ((2*(muryy[ix][iy])*epsilon0)-((muryy[ix][iy])*(sigmaYYprime[ix][iy])*dt))/((2*(muryy[ix][iy])*epsilon0)+((muryy[ix][iy])*(sigmaYYprime[ix][iy])*dt));
			dx2[ix][iy] = 2*(sigmaxx[ix][iy])*eta0*c0*dt;
			dx3[ix][iy] = 2*(sigmaYYprime[ix][iy])*dt/epsilon0;
			dy1[ix][iy] = 2*c0*dt*(1/2 + (sigmaYYprime[ix][iy])*dt/(4*epsilon0));
			dy2[ix][iy] = 2*(sigmaYYprime[ix][iy])*eta0*pow(c0*dt,2);
			dy3[ix][iy] = 2*(sigmayy[ix][iy])*eta0*c0*dt;
			dy4[ix][iy] = 2*(sigmayy[ix][iy])*(sigmaYYprime[ix][iy])*pow(eta0*c0*dt,2);
		}
	}


	// (3.0) Constructing the field matrices
	double** CEz = allocMatrix(dimX,dimY);             // Matrix representing the curl of the E-field for every point, required for Hz iteration
	double** CHx = allocMatrix(dimX,dimY);             // Matrix representing the curl of the H-field for every point, required for Ex iteration
	double** CHxPRE1 = allocMatrix(dimX,dimY);         // Matrix representing the curl of the E-field for every point at last time step, required for Ex iteration
	double** CHy = allocMatrix(dimX,dimY);             // Matrix representing the curl of the E-field for every point, required for Ey iteration          
	double** CHyPRE1 = allocMatrix(dimX,dimY);         // Matrix representing the curl of the E-field for every point at last time step, required for Ey iteration          
	double** Ex = allocMatrix(dimX,dimY);              // The electric field in the x direction at the current time step
	double** Ey = allocMatrix(dimX,dimY);              // The electric field in the y direction at the current time step
	double** Hz = allocMatrix(dimX,dimY);              // The magnetic field in the z direction at the current time step (half t.step)
	double** HzPRE1 = allocMatrix(dimX,dimY);          // The magnetic field in the z direction at the previous time step (-half t.step)
	double** Dx = allocMatrix(dimX,dimY);              // The displacement field in the x direction at the current time step
	double** DxPRE1 = allocMatrix(dimX,dimY);          // The displacement field in the x direction at previous current time step (-1 t.step)
	double** DxPRE2 = allocMatrix(dimX,dimY);          // The displacement field in the x direction at previous current time step (-2 t.step)
	double** Dy = allocMatrix(dimX,dimY);              // The displacement field in the y direction at the current time step
	double** DyPRE1 = allocMatrix(dimX,dimY);          // The displacement field in the y direction at previous current time step (-1 t.step)
	double** DyPRE2 = allocMatrix(dimX,dimY);          // The displacement field in the y direction at previous current time step (-2 t.step)
	double** SumCHy = allocMatrix(dimX,dimY);          // Time integrated curl
	double** SumEy = allocMatrix(dimX,dimY);           // Time integrated electric field
	double** magE = allocMatrix(dimX,dimY);            // Magnitude of the electric field components

	// (3.1) Setting all initial field values to zero and applying the dirichlet boundary condition
	initializeArrayZero(CEz, dimX, dimY);
	initializeArrayZero(CHx, dimX, dimY);
	initializeArrayZero(CHxPRE1, dimX, dimY);
	initializeArrayZero(CHy, dimX, dimY);
	initializeArrayZero(CHyPRE1, dimX, dimY);
	initializeArrayZero(Ex, dimX, dimY);
	initializeArrayZero(Ey, dimX, dimY);
	initializeArrayZero(Hz, dimX, dimY);
	initializeArrayZero(HzPRE1, dimX, dimY);
	initializeArrayZero(Dx, dimX, dimY);
	initializeArrayZero(DxPRE1, dimX, dimY);
	initializeArrayZero(DxPRE2, dimX, dimY);
	initializeArrayZero(Dy, dimX, dimY);
	initializeArrayZero(DyPRE1, dimX, dimY);
	initializeArrayZero(DyPRE2, dimX, dimY);
	initializeArrayZero(SumCHy, dimX, dimY);
	initializeArrayZero(SumEy, dimX, dimY);
	initializeArrayZero(magE, dimX, dimY);

	// The graphics (pre-plotting)
	// The canvas, 
	TCanvas *c = new TCanvas("c","2D simulation of electric fields",0,0,600,900);
	gSystem->Unlink(kGifName.Data());  								     // Used for GIF. Delete existing GIF
	TPaveText *t = new TPaveText(0.02,0.92,0.24,0.99,"brNDC");					     // Used for GIF. Textbox
	TH2D *hist = new TH2D("hist","Norm. E-field magnitude",dimX,0,x_max+2*dx,dimY,0,y_max+dy*(2+nPML));  // Used for plot. The main graphical object used for plotting the E-field magnitude
	hist->GetXaxis()->SetTitle("x-direction, [m]"); hist->GetXaxis()->SetTitleOffset(1.1);		     // Sets the label of the x-axis with a pleasant spacing from the axis
	hist->GetYaxis()->SetTitle("y-direction, [m]"); hist->GetYaxis()->SetTitleOffset(1.1);		     // Sets the label of the y-axis with a pleasant spacing from the axis
	// Drawing the PLM
	double dx2cube = dx*(nx/2); double dy2cube = dy*(nPML/2);					     // The half width of the PML rectangle, in both directions (xy)
	double forigin[2]; forigin[0] = dx*(nx+2)/2; forigin[1] = dy*(1+ny)+dy*(nPML/2);		     // The center point of the PML rectangle in the problemspace
	auto pml = new TBox(forigin[0]-dx2cube,forigin[1]-dy2cube,forigin[0]+dx2cube,forigin[1]+dy2cube);    // Creating the rectangle object representing the position of the PML
	pml->SetFillStyle(3001);								             // Pseudo transparrent due to the fill style
	//Drawing the dielectric
	double rmax = r_dielectric_radius;							             // The radius of a circle describing the dielectric
	double forigin2[2]; forigin2[0] = dx*(nx+2)/2; forigin2[1] = dy*((ny/2)+1);			     // Setting the center point of the dielectric
	auto dielectric = new TEllipse(forigin2[0],forigin2[1],rmax);					     // Creating the circular object representing the position of the dielectric
	dielectric->SetFillStyle(3001); 								     // Pseudo transparrent due to the fill style
	// Drawing the source
	double ysource_start = dy*(ny-11-0.5); double ysource_end = dy*(ny-11+0.5);			     // The start and the end position of the source
	auto source = new TBox(dx*(ix_source_start),ysource_start,dx*(ix_source_end),ysource_end);		     // Constructing the source object
	source->SetFillStyle(3001);									     // Pseudo transparrent due to the fill style

	// (4) The Iteration loop
	for (int it = 0; it < ntime+1; it++){
		string itString = Form("i_{t} = %d",it);						     // Used for plot. Will be inserted into the textbox created by TPaveText
		double timer = it * dt;									     // Current simulation time
		string timeString = Form("t_{i} = %e s",timer);					   	     // Used for plot. Will be inserted into the textbox created by TPaveText
		
		// (4.1) Implementing the source
		for (int ix = ix_source_start; ix <= ix_source_end; ix++){
			Ex[ix][ny-11]=E0*sin(2*(TMath::Pi())*frequency*(dt*it));
			Hz[ix][ny-11]=H0*sin(2*(TMath::Pi())*frequency*(dt*it));
		}
				
		// (4.2) Iterating the fields. First time zero is plotted
		// Calculating the curl of the electric field at every point
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				CEz[ix][iy]= (c1)*(Ey[ix+1][iy]) - (c1)*(Ey[ix][iy]) - (c2)*(Ex[ix][iy+1]) + (c2)*(Ex[ix][iy]);
			}			
		}


		// Calculating the magnetic field of the next time-step
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				HzPRE1[ix][iy] = (Hz[ix][iy]);                 // Saving the previous value
				Hz[ix][iy] = -(hz1[ix][iy])*(CEz[ix][iy]) + (hz2[ix][iy])*(HzPRE1[ix][iy]);
			}			
		}
		

		// Calculating the curl of the magnetic field (dHz/dy)
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				CHxPRE1[ix][iy] = CHx[ix][iy];                 // Saving the previous value
				CHx[ix][iy] = (c2)*(Hz[ix][iy]) - (c2)*(Hz[ix][iy-1]);
			}			
		}
		
		
		// Calculating the displacement field Dx for the next time-step
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				DxPRE2[ix][iy] = DxPRE1[ix][iy];                 // Saving the previous value
				DxPRE1[ix][iy] = Dx[ix][iy];                 // Saving the previous value
				Dx[ix][iy] = (dx1)*(CHx[ix][iy]) + (dx1)*(CHxPRE1[ix][iy]) - (dx2[ix][iy])*(Ex[ix][iy]) + DxPRE2[ix][iy] - (dx3[ix][iy])*(DxPRE1[ix][iy]);
			}			
		}		
		
		// Calculating the curl of the magnetic field (-dHz/dx)
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				CHyPRE1[ix][iy] = CHy[ix][iy];
				CHy[ix][iy] = -(c1)*(Hz[ix][iy]) + (c1)*(Hz[ix-1][iy]);
			}			
		}		

		// Calculating the displacement field Dy for the next time-step
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				DyPRE2[ix][iy] = DyPRE1[ix][iy];
				DyPRE1[ix][iy] = Dy[ix][iy];
				Dy[ix][iy] = (dy1[ix][iy])*(CHy[ix][iy]) + (dy1[ix][iy])*(CHxPRE1[ix][iy]) + (dy2[ix][iy])*(SumCHy[ix][iy]) - (dy3[ix][iy])*(Ey[ix][iy]) - (dy4[ix][iy])*(SumEy[ix][iy]) + DyPRE2[ix][iy];
			}			
		}		

		// Calculating the electric field Ex for the next time-step
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				Ex[ix][iy] = (erxx[ix][iy])*(Dx[ix][iy]);
			}			
		}		

		// Calculating the electric field Ey for the next time-step
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				Ey[ix][iy] = (eryy[ix][iy])*(Dy[ix][iy]);
			}			
		}	

		// Updating the sums
		for (int iy=1; iy<dimY-1; iy++){
			for (int ix=1; ix<dimX-1; ix++){
				SumEy[ix][iy] = SumEy[ix][iy] + Ey[ix][iy];
				SumCHy[ix][iy] = SumCHy[ix][iy] + CHxPRE1[ix][iy];
			}			
		}
		
		// (4.3) Calculating the magnitude of the electric field
		for (int iy=0; iy<dimY; iy++){
			for (int ix=0; ix<dimX; ix++){
				magE[ix][iy] = sqrt((Ex[ix][iy])*(Ex[ix][iy])+(Ey[ix][iy])*(Ey[ix][iy]));
			}
		}
		double maxmagE = findMaxfromArray(magE,dimX,dimY);
		string maxmagEString = Form("max(#tilde{E}) = %f",maxmagE);
	
		// (4.4) Plotting
		hist->Reset("ICESM"); // clearing the histogram
		hist->SetMaximum(1.0*E0);
		hist->SetMinimum(-0.01*E0);	
		for (int iy=0; iy<dimY; iy++){
			for (int ix=0; ix<dimX; ix++){
				double magEc = magE[ix][iy];
				double xc    = ix*dx; 
				double yc    = iy*dy;
				hist->Fill(xc,yc,magEc);
			}
		}
		if(it==0){
			hist->Draw("colz");					// Used for plot. 2D Color plot
			hist->SetStats(0);   					// Remove the histogram statbox
		}
		else{
			c->Update();						// Used for plot. Updates the canvas with the new values
			t->Clear();						// Used for plot. Clear the textbox from the previous iteration
		}
		t->AddText(itString.c_str());	      				// Used for plot. Adding the current iteration number to the textbox object
		t->AddText(timeString.c_str());	      				// Used for plot. Adding the current time to the textbox object
		t->AddText(maxmagEString.c_str());				// Used for plot. Adding the max E-field to the textbox object
		t->Draw();	 						// Used for plot. Drawing the textbox with the values above
		pml->Draw("SAME");						// Drawing the rectangular PML					       
		dielectric->Draw("SAME");					// Drawing the cicular dielectric
		source->Draw("SAME");						// Drawing the rectangular source
		c->Print(kGifName+"+"+kGifSleep);				// Used for GIF. Adds another frame with a hold time
	}
	c->Print(kGifName+"++"+kGifSleep);					//Used for GIF. Makes the last frame of the iteration loop to the beginning
}

// Allocates space in memory for a 2D matrix of sizeX, sizeY
double** allocMatrix(int sizeX, int sizeY){
    double** matrix = new double*[sizeX];
    for (int ix = 0; ix < sizeX; ix++){
        matrix[ix] = new double[sizeY];
    }
    return matrix;
}

// Fills the memory allocated 2D matrix with zeros
void initializeArrayZero(double **arr,double sizeX,double sizeY){
	for (int i = 0; i < sizeX; i++){
		for (int j = 0; j < sizeY; j++){
			arr[i][j] = 0;
		}
	}
}

// Fills the memory allocated 2D matrix with ones
void initializeArrayOne(double **arr,double sizeX,double sizeY){
	for (int i = 0; i < sizeX; i++){
		for (int j = 0; j < sizeY; j++){
			arr[i][j] = 1;
		}
	}
}

// Returns the maximum value of a 2D matrix
double findMaxfromArray(double **arr,int sizeX,int sizeY){
	double max = arr[0][0];
	for (int i = 0; i < sizeX; i++){
		for (int j = 0; j < sizeY; j++){
			if (arr[i][j]>max){
				max=arr[i][j];
			}
		}
	}
	return max;
}

2.1 The results

Here are some examples with different input parameters showing some of the output that works well, as well as some graphical bugs and other odd behaviors. The roots of these odd behaviors were not found. The simplest case is an isotropic system with no dielectric and no PML with \( \frac{\lambda}{10} \) defining the grid lengths dx and dy:

As seen above, there is one huge graphical bug present which is causing the weird patterns. This bug seem to be associated with the resolution of the system as the following sorted out the graphical bug: by adjusting the resolution governed by \( ds=\frac{\lambda}{10} \) into the following, \( ds=\frac{\lambda}{19} \), causes the pattern to almost completely go away:

To resolve this bug, different definitions of the grid resolution has been invoked and tested as seen above, however all of them produced some graphical issues (to a large or small extent), but most importantly they depend on the other parameters as well such as the dielectric constant, as seen later. However, we can notice further issues: (1) Larger than expected field amplitudes, (2) no divergence around the source. The latter actually arrising from the fact that the fields does not spread propperly from one cell to the next in the x-direction, and will cause issues in the continuity of the wave equation to be broken when we add the dielectric. However, wavebehavior is attained as predicted: (1) reflections from the boundary condition, (2) production of standing waves due these reflections. Since the Yee grid should set the governing divergence equations, the mistake had to be with the Yee Grid. Taking a second look at the derivation and other books on FDTD electrodynamics no mistake was found by the Author. It is still interesting to see how the system performs when adding dielectric elements as this might give some further clues. For the next test we add a dielectric element: \( \epsilon_{r}=2, \mu_{r}=1, \sigma=0, ds=\frac{\lambda}{19} \):

This as in the first picture, we attain a lot of junk. Just as previously, we change the grid resolution, \( ds=\frac{\lambda}{25} \) as a temporary solution:

As mentioned earlier, here we see that there is some discontinuity that we should not be attaining from the Maxwell's equation on the boundary of the dielectric, which seems to stem from the inability to propegate along the x-direction. Otherwise we also have no standing waves beneath the dielectric. There thus seems to be two mains bugs, one bug related to the divergence, and one bug related to the divergence and wave behavior in the x-direction that makes waves exhibit odd behavior. Since the PML was also included in the derivation, it would be a shame to not introduce it. Including the PML:

The PML does not seem to absorb the waves as it should be doing, even when manually increasing the strength of the PML to 20 times the original strength it does not abosorb the waves properly. However, the code that simulates free space has issues the field amplitudes, and these large field values might mask the effect of the PML. What is left is two definite bugs that has been characterized, but unsucessfully resolved: The error with the divergence, and the resolution error.

Regarding the author

Created by John Benjamin in association with the course PHYS291 at the University of Bergen.
Email: jlo016@uib.no