Electrodynamica

a program for simulating electromagnetic fields

Terje Vold

Continuum Technology, Inc.

Why Write An Electrodynamics Simulation Program?

Personal Professional Need

Market Need

Beautiful Physics and Mathematics

An Unexploited Method Exists

What Electrodynamica can Simulate

Space filled with piecewise homogeneous volumes 
with linear electromagnetic response (permittivity ϵ and permeability μ)

Complicated irregularly shaped volumes, mixed types of materials, any frequency

Input = Shapes, perms ϵ and μ, and source charge and current
Fundamental output = fields
Secondary output = impedance, energy, power

The unexploited method:  Geometric Algebra

Geometric algebra is a Clifford algebra with a geometric interpretation.

It contains conventional (Gibbs') vector algebra as a subsystem.

Elements of the algebra, called multivectors, 
are sums of scalars, vectors, pseudovectors, and pseudoscalars, 
much as complex numbers are sums of scalars and imaginary number. 

The geometric product is associative but not commutative (like matrix multiplication).

The symmetric part of the geometric product of two vectors is the same as the conventional dot product. 

The antisymmetric part of the geometric product of two vectors corresponds to the cross product. 

The product's associativity allows important results from real analysis to apply.

The square of any pseudoscalar or any pseudovector equals  -1.

The algebra has an 8x8 matrix representation.

Maxwell's equations

\(\((\(∇\& ⇀ \)\ \(+\ \ ∂\_t\))\) F\  = \ S\)

where

F  = (Overscript[E,  ⇀])   + I  
                    (Overscript[B,  ⇀])

and

\(S\  = \(\(\(ρ\)\(\ \)\)\(-\)\(j\& ⇀ \)\(\ \)\)\)

Maxwell's equations expanded

Explicitly write all four grades

				\(\(∇\& ⇀ \)\(\(·\)\(E\& ⇀ \)\) = ρ\)
				\(\(\(\(\(∇\& ⇀ \)\(\(·\)\((I\ \(B\& ⇀ \))\)\)\)\(\ \)\)\(+\)\(\ \)\(∂\_t\( 
                      E\& ⇀ \)\)\(\ \)\) = \ \(-\(j\& ⇀ \)\)\)
				\(\(\(\(\(∇\& ⇀ \)\(⋀\(E\& ⇀ \)\)\)\(\ \)\)\(+\)\(\ \)\(∂\_t\((I\ \(B\& ⇀ \))\)\)\(\ \)\) = \ 0\)
				\(\(∇\& ⇀ \)\(⋀\((I\ \(B\& ⇀ \))\)\)\  = \ 0\)

Use the identities
				\(\(a\& ⇀ \)·\((I\ \(b\& ⇀ \))\)\  = \ 
              I\ \((\(a\& ⇀ \)⋀\(b\& ⇀ \))\)\)
				\(\(a\& ⇀ \)⋀\((I\ \(b\& ⇀ \))\)\  = \ 
                I\ \((\(a\& ⇀ \)·\(b\& ⇀ \))\)\)
and the relationship
				\(\(I\ \(a\& ⇀ \)\)×\(b\& ⇀ \)\  = \ \(a\& ⇀ \)⋀\(b\& ⇀ \)\)

to rewrite as
				\(\(\(∇\& ⇀ \)\(\(·\)\(E\& ⇀ \)\) = ρ\)
				\(\(\(\(∇\& ⇀ \)\(\(\(×\)\(\ \)\) \(B\& ⇀ \)\)\)\(\ \)\(-\)\(\ \)\(∂\_t\( 
                        E\& ⇀ \)\)\(\ \)\) = \ \(j\& ⇀ \)\)
				\(\(\(\(∇\& ⇀ \)\(\(×\)\(E\& ⇀ \)\)\)\(\ \)\(+\)\(\ \)\(∂\_t\( 
                        B\& ⇀ \)\)\(\ \)\) = \ 
                0\)
				\(\(∇\& ⇀ \)\(\(\(·\)\(\ \)\) \(B\& ⇀ \)\)\  = \ 0\)\)

Maxwell's equations in homogeneous linear media

\(TraditionalForm\`\((\(∇\& ⇀ \)\ \(+\ \(1\/v\)\) ∂\_t)\)\ f = s\) 

where
			\(TraditionalForm\`f\  = \ \(e\& ⇀ \)\  + \ \(\(\(I\)\(\ \)\) \(b\& ⇀ \)\(\ \)\)\) 
with
			    \(\(e\& ⇀ \) = \(\@ϵ\) \(E\& ⇀ \)\)
and
			       \(\(b\& ⇀ \)\ \  = \(\@\(1\/μ\)\) \(B\& ⇀ \)\)  
and
		    \(s\ \  ≡ \ \ \(1\/\@ϵ\) ρ\ \ \  - \ \ \ \(\(\@μ\) \(j\& ⇀ \)\(\ \)\)\) .

Or, using
			 f[r, t] = 
        f[r] Exp[-i ω t]

we have
			\(TraditionalForm\`\((\(∇\& ⇀ \)\ \(-\ i\)\ κ)\)\ f = s\)

Maxwell's equations at boundaries

\(TraditionalForm\`\(\(\(ρ\_r\%R[\ \(f\^R\)[\(r\& ⇀ \)]]\)\(\ \)\)\(-\)\(\ \)\(ρ\_r\%L[\ \(f\^L\)[\(r\& ⇀ \)]]\)\(\ \)\) = \ \(1\/2\) \((\(u\&^\)\_r\%R\  - \ \(u\&^\)\_r\%L)\)\ θ\_r\) 

where
			\(TraditionalForm\`θ\_r\  ≡ \ σ\_r\  - \ \(\(κ\& ⇀ \)\_r\)\)

Statement of the Problem:

find field values that satisfy Maxwell's equation everywhere

Green's Theorem

Maxwell's equation in a volume V is equivalent to

		\(f[\(r\& ⇀ \)\^′] = ∮\_\(∂V\)g[\(r\& ⇀ \) - \(r\& ⇀ \)\^′] 
                  d \( σ\& ⇀ \)\ f[\(r\& ⇀ \)]\  + \ ∫\_V\ 
                g[\(r\& ⇀ \) - \(r\& ⇀ \)']\ \ dτ\ s[\(r\& ⇀ \)]\)

on the boundary ∂V of the volume,

where 
				g[(Overscript[r,  ⇀]) ] = ((Overscript[∇,  ⇀])  -  i κ)ϕ[(Overscript[r,  ⇀]) ]
				
			\(TraditionalForm\`ϕ[\(r\& ⇀ \)] = \ \(-\(1\/\(4\ π\ r\)\)\)\ Exp[i\ κ\ r]\)  
		
and s is the density of charge and current sources in volume V.

Green's Theorem Expressions Compared

Conventional (Gibbs') Vector Algebra		
	

		    (Overscript[E,  ⇀]) [(Overscript[r,  ⇀]) ']=(∮_(∂V))[-i ω μ \((\(n\& ⇀ \)×\(H\& ⇀ \))\)ϕ +\(\(\((\(n\& ⇀ \)×\(E\& ⇀ \))\)×\(∇\& ⇀ \)\) ϕ\) + \(\((\(n\& ⇀ \)·\(E\& ⇀ \))\) \(∇\& ⇀ \)ϕ\)] dσ 
	
		(Overscript[H,  ⇀]) [(Overscript[r,  ⇀]) ']=(∮_(∂V))[i ω ϵ \((\(n\& ⇀ \)×\(E\& ⇀ \))\)ϕ  + \(\(\((\(n\& ⇀ \)×\(H\& ⇀ \))\)×\(∇\& ⇀ \)\) ϕ\) + \(\((\(n\& ⇀ \)·\(H\& ⇀ \))\) \(∇\& ⇀ \)ϕ\)] dσ 


Geometric Algebra

				\(f[\(r\& ⇀ \)'] = ∮\_\(\(\ \)\(∂V\)\)\ g\ d \( σ\& ⇀ \)\ f\) .

Restatement of the Problem

Satisfy
		\(TraditionalForm\`\(\(\(\(\(f\^\(\(\(\ \)\(h\)\)\(\ \)\)\)[\(r\& ⇀ \)]\)\(\ \)\)\(\ \)\)\(-\)\(\ \)\(\ \)\(∮\_\(∂V\)\ \(g\^V\)[\(r\& ⇀ \)\^′ - \(r\& ⇀ \)]\ \ d \( σ\& ⇀ \)\^′\ \(f\^\(\(\(\ \)\(H[r\^′, 
                      V]\)\)\(\ \)\)\)[\(r\& ⇀ \)\^′]\)\(\ \)\(\ \)\) = \ ∫\_V\ \(\(\(g\^V\)[\(r\& ⇀ \)\^′ - \(r\& ⇀ \)]\) \(\(\ \)\(\ \)\) \(dτ\^′\)\(\ \)\(\(s\^V\)[\(r\& ⇀ \)\^′]\)\(\ \)\)\)
and
			\(TraditionalForm\`\(\(\(ρ\_r\%R[\ \(f\^R\)[\(r\& ⇀ \)]]\)\(\ \)\)\(-\)\(\ \)\(ρ\_r\%L[\ \(f\^L\)[\(r\& ⇀ \)]]\)\(\ \)\) = \ \(1\/2\) \((u\_r\%R\  - \ u\_r\%L)\)\ θ\_r\) 
			
on all boundaries.

Discretization

Find surface field values f that satisfy 


		\(TraditionalForm\`\(\(\(∑\_\(j = 1\)\%N\(∑\_\(g = R, L\)\(M\_\(\(\ \)\(i\ j\)\)\%\(\(\ \)\(h\ g\)\)\) 
              f\_j\%\(\(\ \)\(g\)\)\)\)\(\ \)\)\(=\)\(\ \)\(B\_i\%\(\(\ \)\(h\)\)\)\(\ \)\)\),			i ϵ {1,2,...,N},   h ϵ {R, L}

and

		(Subscript[N, i]^R) (Subscript[f, i]^R)-(Subscript[N, i]^L) 
            (Subscript[f, i]^L)  =  (Subscript[C, i]^R) - (Subscript[C, i]^L),			i ϵ {1,2,...,N}.

Solution

We have 3N multivector equations and 2N unknown multivectors;

either make these numbers equal and solve, or minimize total squared error

Example 1:  Some Simple Shapes

Description of Shapes and Sources for Example 1

Simple Shapes to test and demonstrate the simultaneous exercise of a number of capabilities of Electrodynamica:

	a sphere with a permittivity ϵ = 3 and permeability μ=10
	a cube with both permittivity ϵ=1 and permeability μ=1, 
		and with a dipole surface charge and current distribution.  
	a pair of ring-shapes surfaces with current (Helmholz coils) 
	a cylindrical standing wave generating a small oscillating electric field
	
	All sources are oscillating at a low frequency of 1/1000, 
	in units such that the sphere's radius is 1 and the speed of light is c=1.

Plots for Example 1

The sphere is a Mathematica Shape:

<< Graphics`Shapes`

testSphere =  Sphere[1, 6, 5] ;

Show[Graphics3D[testSphere]]

We set its permittivity to 3 and permeability to 10:

testSpherePerm = {3, 10};

Plots for Example 1

We also define a pair of rings and a cube; a plot of all surfaces in this example 1 is below.

Show[Graphics3D[ 
        surfaceList[
          testu] ],PlotRange -> {{-15, 15}, {-12, 
            12}, {-12, 12}},ViewPoint -> {1, -6, 
          1.5} ]

Plots for Example 1

Before solving, we look at the applied field due to the Helmholz coils.

Show[ helmholzPVFListMagneticPlot, Graphics3D[ surfaceList[testu] ], 
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> {{-30, 30}, {-30, 30}, {-30, 30}}, 
      ViewPoint -> {1, -5, 2}];

Plots for Example 1

For the electric field, we have

Show[Graphics3D[ surfaceList[testu] ], helmholzPVFListElectricPlot, 
      DisplayFunction -> $DisplayFunction, 
      PlotRange -> {{-30, 30}, {-30, 30}, {-30, 30}}, 
      ViewPoint -> {1, -5, 2}];

Plots for Example 1

Now we back way out, and Show a PlotRange of the field due to the Helmholz coils that is approximately 100 times larger:

Show[ helmholzPVFBigListMagneticPlot,Graphics3D[ surfaceList[testu] ],DisplayFunction -> $DisplayFunction,PlotRange -> {{-1000, 1000}, {-800, 
                800}, {-1000, 1000}},ViewPoint->{2,-10, 3}];

Plots for Example 1

This shows the beginnings of a radiated electromagnetic wave.   We can plot the magnetic and electric fields together:

Show[Graphics3D[ surfaceList[testu] ],helmholzPVFBigListMagneticPlot,helmholzPVFBigListElectricPlot,ViewPoint->{2,-10, 3},DisplayFunction -> $DisplayFunction,PlotRange -> {{-1000, 1000}, {-800, 
                800}, {-1000, 1000}}];

Plots for Example 1

We can also check just the nearly-uniform applied oscillating electric field:  
actually a cylinder wave, since uniform oscillating fields are prohibited by Maxwel's equations.   
With a PlotRange equal to about 100 times the width of the Helmholz coils, we have the following plot:

helmholzPVFBigListElectricPlot = 
      calcFieldListPlot[testuHelmPVFieldBigList, 1, 0, 5];

Show[Graphics3D[ surfaceList[testu] ],helmholzPVFBigListElectricPlot,ViewPoint->{2,-10, 5},DisplayFunction -> $DisplayFunction,PlotRange -> {{-1000, 1000}, {-1000, 
                1000}, {-1000, 1000}}];

Plots for Example 1

We back out still a bit further, to about 500 times the size of the Helmholz coils, 
and plot the field along one horizonal line passing through the Helmholz coils:

Show[Graphics3D[ surfaceList[testu] ],helmholzPvf3kListElectricPlot,helmholzPvf3kListMagneticPlot,ViewPoint->{2,-10, 4},DisplayFunction -> $DisplayFunction,PlotRange -> {{-500, 3000}, {-1000, 
                1000}, {-1000, 1000}}];

Plots for Example 1

We next solve for the field that results when all sources are active - 
Helmholz magnetic field, electric dipole cube, and the nearly-homogeneous small applied electric field - 
and both volumes are responding.

Below is a plot of the magnetic field on a square array of points:

Show[{Graphics3D[ surfaceList[testu] ], 
            testuPVFieldPlotPlane},ViewPoint->{0,-4, 1},DisplayFunction -> $DisplayFunction,PlotRange -> {{-15, 15}, {-12, 12}, {-10, 
                15}}];

This shows nothing surprising, but just a nearly uniform magnetic field, smoothly changing with no irregularities near either band of current, being applied to a weakly magnetic sphere - as it should.

Plots for Example 1

To look more closely at the calculated response of the magnetic sphere, 
we plot magnetic field values along a horizontal line that runs through the center of the sphere:

Show[Graphics3D[boundarySurfaceList[testu][[1]]], 
      plotxAxisBList, ViewPoint -> {0, -4, 0}];

Again, this shows nothing surprising.  The value inside the sphere larger than the exact result for a sphere by about 1% for algorithm c) above, and by about 3% for the algorithm used in this figure.

Plots for Example 1

Now we examine the electric field.  
Recall that a small, essentially homogeneous electric field is being applied, while the cube has a dipole surface charge:

Show[{Graphics3D[ surfaceList[testu] ], 
            testuPVElectricFieldPlotPlane},ViewPoint->{0,-4, 1},DisplayFunction -> $DisplayFunction,PlotRange -> {{-15, 15}, {-12, 12}, {-10, 
                15}}];

Plots for Example 1

In the plot above, the displayed field is too large near the cube and too small elsewhere, 
so we look close to the sphere, and then close to the cube, with different scales for the field.  

First we plot the electric field near the sphere:

Show[{Graphics3D[ surfaceList[testu] ], 
          testuPVElectricFieldPlotPlane},ViewPoint->{3,-4, 2},DisplayFunction -> $DisplayFunction,PlotRange -> {{-6, 6}, {-6, 6}, {-6, 6}}]

This is (at least qualitatively) as we expect:  the weakly dielectric sphere (ϵ = 3) is drawing the field in a bit, and the dipole field from the cube to the upper right of the sphere is weakly modifying the applied electric field.

Plots for Example 1

Near the cube, the field is dominated by the field due to the charge on the cube:

Show[{Graphics3D[ surfaceList[testu] ], 
            testuPVFieldPlanePointsPlot},ViewPoint->{3,-5, 2},DisplayFunction -> $DisplayFunction,PlotRange -> {{10, 16}, {-2, 2}, {0, 5}}];

Summary of Results for Example 1

No surprises:   

With some simple volumes and applied fields for which we know ahead of time what to expect, 

the numerical solution gives what we expect with small errors.

Example 2:  Shielding

Description of Shapes and Sources for Example 2

A cylindrical shell of material with permeability μ=100 and permittivity ϵ=1, 

with a uniform magnetic field applied by Helmholz coils.  
 
There are no other volumes or free field sources.  
We look at two orientations of the shell relative to the applied field.

Plots for Example 2

First we show a plot of solution field values on a square array of points for a cylinder with its axis parallel to the applied field.

Show[{Graphics3D[ surfaceList[testu] ], 
            testuPVFieldPlotPlane},ViewPoint->{0,-4, 1/3},DisplayFunction -> $DisplayFunction,PlotRange -> {{-15, 15}, {-12, 12}, {-10, 
                15}}];

Plots for Example 2

Looking more closely at the field near the cylinder, we have the following plot:

Show[{Graphics3D[ cs1 ], testuPVFieldPlotSmallPlane},ViewPoint->{0,-4, 1/3},DisplayFunction -> $DisplayFunction,PlotRange -> {{-4, 4}, {-4, 4}, {-4, 
                4}}];

This was a graphics cell with 14514 lines.

This shows both the known correct qualitative response of such a cylinder to an applied field, 
and also errors due to the large discrete surface segments 
and the approximation of a uniform field on each segment.

This was a graphics cell with 1631 lines.

Plots for Example 2

We plot the magnetic field with the shield axis perpendicular to the applied field. 
We display the "cut Shape" cs1 instead of the entire cylindrical surface:

Show[{Graphics3D[ {SurfaceColor[GrayLevel[0.6]], cs1} ], 
            testuPVFieldPlotSmallPlane},ViewPoint->{1,-4, 1/3},DisplayFunction -> $DisplayFunction,PlotRange -> {{-4, 4}, {-4, 4}, {-4, 
                4}}];

Errors due to the choice of large surface segments remain, though they are probably smaller.

Summary of Results for Example 2

The predicted response of a magnetic cylindrical shield 
to a uniform applied magnetic field 
is at least qualitatively as expected; 

but in this example, 
with surface segment dimensions chosen to be nearly as large as the radius of the cylinder, 
there are significant errors 
associated with the approximation of a constant field on each surface segment.

Example 3:  A Simple Circuit

Description of Shapes and Sources

A simple electrical circuit 
consisting of an oscillating current source connected to a small wire loop.  

The current source is modeled as a cube of material with ϵ=1 and μ=1, 
surface charge of opposite signs on two opposing faces, 
and surface current flowing between these faces on the four faces between them.  

A small wire loop with electrical conductivity σ 
is modeled with complex permittivity ϵ = 1 + iσ/ω , 
and its surface is modeled with a small number of square surface segments.

Plots for Example 3

The assembled simple circuit is plotted below:

Show[ Graphics3D[{SurfaceColor[RGBColor[1, 0, 0]], batterySurface}], 
      Graphics3D[wireSurface], ViewPoint -> {-2, -4, 3} , 
      PlotRange -> {{-1, 8}, {-3, 3}, {-1, 1}}]

Plots for Example 3

A plot of electric field values in the plane through the middle of the circuit:

planeFieldList = 
      calcPointVolumeFieldList[ ex2,  planeOfPoints];

planeFieldERePlot = 
      calcFieldListPlot[planeFieldList, 1, 0, 2.5];

Show[{Graphics3D[ cutWire ], planeFieldERePlot},ViewPoint->{1,-2, 5},DisplayFunction -> $DisplayFunction,PlotRange -> {{-3, 11}, {-5, 5}, {-10, 
              10}}]

Plots for Example 3

showAnimationWithViewOscillation[
      ex2, {electricFieldAnimation, 
        magneticFieldAnimation}, {{-3, 11}, {-6, 6}, {-5, 5}}, 5, 
      1]

Plots for Example 3

We plot the field values along one axis with a range of about  3 (10^6) (about (10^6) times the wire loop radius).  
We see a snapshot of radiated waves:

Plots for Example 3

These fields can be animated, too (again, click or double-click on the plot to animate it):

showAnimationWithViewOscillation[
      ex2, {farElectricFieldAnimation, 
        farMagneticFieldAnimation}, {{0, 3}, {-1, 1}, {-1, 1}}*10^6, 
      5, 1]

Summary of Results for Example 3

This simple circuit shows that expected features are predicted 
over at least 6 orders of magnitude of length scale, 
from details in and around the wire loop 
to the radiation field at a distance of about a million loop radii, 

although we expect quantitative errors due to the use of surface segments 
that are nearly as large as the dimensions of the structure being modeled.

Example 4:  A Dielectric Waveguide

Description of Shapes and Sources

A short segment of a dielectric waveguide with a square cross section, 
with an electric dipole point source on the axis of the waveguide, one wavelength from one end.  

Waveguide index of refraction equals 2. 
The waveguide surface is modeled with 
approximately 6 segments per interior wavelength in the axial direction, 
and one segment per side of the waveguide.  

With these dimensions, the waveguide allows only two longitudinal modes.

Plots for Example 4

First we show a plot of the waveguide segments 
and of a sphere that represents the point dipole source:

dipoleSource = Sphere[1/4, 8, 6];

waveguide = 
      TranslateShape[ 
        rectangularPrism[10, 1, 1, 60, 1, 1], {7, 0, 0}];

Show[ Graphics3D[ dipoleSource ], Graphics3D[waveguide], 
      PlotRange -> {{-1, 15}, {-2, 2}, {-2, 
            2}} ]

Plots for Example 3

Below is an animated plot of the solution fields, 
with an outline of the waveguide and again a sphere representing the dipole source. 
The fields are plotted on the axis of the waveguide on both sides of the dipole.

showAnimationWithViewOscillation[
      ex4, {xPosElectricAnimation, xPosMagneticAnimation, 
        xNegElectricAnimation, 
        xNegMagneticAnimation}, {{-15, 15}, {-6, 6}, {-8, 8}}, 5, 
      1/2]

We see poor coupling at the input face, 
a mix of a traveling and standing waves inside the waveguide, 
and a wave traveling away from the exit face that appears to be largely due to the source wave scattering around the outside of the waveguide.

Summary of Results for Example 4

As with the other examples, 
the surface segments are chosen to be about as large as is possible 
while still getting at least a qualitatively correct simulation.  

This simulation demonstrates that 
when surface segment dimensions are chosen to be comparable to the wavelength 
electromagnetic waves traveling between media with different indices of refraction 
are simulated at least qualitatively correctly. 

Better characterization will require using an example with a known analytic solution, 
and using a finer surface segmentation.

Summary Remarks

Everything works as desired.
Errors are large if surface segments are as large as the volumes.  

Planned development includes
	improved performance when using large surface segments
	better characterization of performance 
	packaging as a user-friendly Mathematica application