A111Q2 blb3MM

NATL INST OF STANDARDS & TECH R.I.C.

A1 11 0261 6344

Randa, J/A lattice approach to volumes I QC100 .U5753 NO.1303 1986 C.1 NBS-PUB-C

*°*tAX> O* *• NBS TECHNICAL NOTE 1303

U.S. DEPARTMENT OF COMMERCE / National Bureau of Standards
Tm he National Bureau of Standards' was established by an act of Congress on March 3, 1901. The

m Bureau's overall goal is to strengthen and advance the nation's science and technology and facilitate their effective application for public benefit. To this end, the Bureau conducts research and provides: (1) a

basis for the nation's physical measurement system, (2) scientific and technological services for industry and

government, (3) a technical basis for equity in trade, and (4) technical services to promote public safely.

The Bureau's technical work is performed by the National Measurement Laboratory, the National

Engineering Laboratory, the Institute for Computer Sciences and Technology, and the Institute for Materials Science and Engineering.

The National Measurement Laboratory

Provides the national system of physical and chemical measurement;

coordinates the system with measurement systems of other nations and

furnishes essential services leading to accurate and uniform physical and

chemical measurement throughout the Nation's scientific community, in- dustry, and commerce; provides advisory and research services to other Government agencies; conducts physical and chemical research; develops,

produces, and distributes Standard Reference Materials; and provides

calibration services. The Laboratory consists of the following centers:

• Basic Standards2 • Radiation Research

• Chemical Physics

• Analytical Chemistry

The National Engineering Laboratory

Provides technology and technical services to the public and private sectors to address national needs and to solve national problems; conducts research in engineering and applied science in support of these efforts; builds and maintains competence in the necessary disciplines required to carry out this research and technical service; develops engineering data and measurement

capabilities; provides engineering measurement traceability services; develops

test methods and proposes engineering standards and code changes; develops

and proposes new engineering practices; and develops and improves

mechanisms to transfer results of its research to the ultimate user. The

Laboratory consists of the following centers:

Applied Mathematics

Electronics and Electrical Engineering2 Manufacturing Engineering

Building Technology

Fire Research

Chemical Engineering2

The Institute for Computer Sciences and Technology

Conducts research and provides scientific and technical services to aid Federal agencies in the selection, acquisition, application, and use of computer technology to improve effectiveness and economy in Government

operations in accordance with Public Law 89-306 (40 U.S.C. 759), relevant Executive Orders, and other directives; carries out this mission by managing

the Federal Information Processing Standards Program, developing Federal ADP standards guidelines, and managing Federal participation in ADP

voluntary standardization activities; provides scientific and technological advisory services and assistance to Federal agencies; and provides the technical foundation for computer-related policies of the Federal Government. The Institute consists of the following centers:

Programming Science and

Technology

Computer Systems

Engineering

The Institute for Materials Science and Engineering

Conducts research and provides measurements, data, standards, reference

materials, quantitative understanding and other technical information fundamental to the processing, structure, properties and performance of materials; addresses the scientific basis for new advanced materials technologies; plans

research around cross-country scientific themes such as nondestructive

evaluation and phase diagram development; oversees Bureau-wide technical programs in nuclear reactor radiation research and nondestructive evaluation; and broadly disseminates generic technical information resulting from

its programs. The Institute consists of the following Divisions:

• Ceramics

• Fracture and Deformation

• Polymers

Metallurgy

Reactor Radiation

'Headquarters and 1 aboraiories ai Gailhersburg, MD, unless otherwise noted; mailing address

Gailhersburg, MD 208W.

"Some divisions wilhin the center are located at Boulder, CO 80303.

lied ii Houldcr, CO, with some elements ai Gailhersburg, MD.
NBS

RESEARCH

INFORMATION CENTER A Lattice Approach to Volumes

Irradiated by Unknown Sources

J. Randa

M. Kanda

Electromagnetic Fields Division

Center for Electronics and Electrical Engineering

National Engineering Laboratory

National Bureau of Standards

Boulder, Colorado 80303

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director

Issued October 1986
National Bureau of Standards Technical Note 1303

Natl. Bur. Stand. (U.S.), Tech Note 1303, 64 pages (Oct. 1986)

CODEN: NBTNAE

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 1986

For sale by the Superintendent of Documents, U.S. Government Printing Office, Washington, DC 20402
CONTENTS

Page

Abstract 1

1. INTRODUCTION 2

2. STATIONARITY OF THE ACTION 5

2.1 Formulation 5

2.2 Nature of the Stationary Point and the Euclidean Action. . 9

2.3 Simple Example 12

2 A Comments 19

3. DIFFERENTIAL EQUATION SOLUTION 20

3.1 Method 20

3.2 Two-Dimensional Examples 24

3.3 Three-Dimensional Examples 28

4. DISCUSSION AND CONCLUSIONS 31

5. ACKNOWLEDGMENT 3^

6. REFERENCES 3^

in

A LATTICE APPROACH TO

VOLUMES IRRADIATED BY

UNKNOWN SOURCES

J. Randa and M. Kanda

Electromagnetic Fields Division

National Bureau of Standards

Boulder, CO 80303

We suggest an approach to the characterization of electromagnetic

environments irradiated by unknown sources. The approach is based on the

numerical solution of Maxwell's equations subject to the constraints imposed

by the measured values of the field at a small number of measurement points

and by boundary conditions. A thorough examination of two methods for the

numerical solution is presented. The examples attempted demonstrate the

approach but reveal that neither numerical technique is fully successful.

Possible future directions are suggested.

Key words: electromagnetic environment characterization; electromagnetic

environment effects; Hamilton's action principle; ill-posed problems;

numerical methods; successive over-relaxation method.
1 . INTRODUCTION

As the number of intentional and inadvertent sources of electromagnetic

(EM) radiation increases, the task of measuring and characterizing EM

environments becomes increasingly difficult. At the same time more sensitive

electronic devices are being used, and evidence is beginning to accumulate

(see e.g. [1]) that the human body may be more sensitive to low-level

nonionizing radiation than previously assumed. It is therefore increasingly

important to develop methods to measure and characterize the EM environments

in which devices and people are found.

The only systematic approach which has been developed and actually

implemented is the statistical approach [2-4]. It has been developed quite

thoroughly for applications to communications systems [5-7] and has been used

in a wide variety of applications, including EM noise in mines [8], EM

emissions from cars [9, 10] and trains [11], radio reception inside houses

[12, 13], and characterization of urban environments [14-16]. Statistical

methods are well suited to cataloguing general features of general

environments, such as the average field level in a typical room of a generic

suburban house. However, they require too many measurements to be efficient

for specific characteristics and particular situations. (A directional

scanning method which requires less measurement effort has been suggested

[17], but it has not yet been fully developed.)

In this paper we address the question of how to efficiently extract

information about a specific small to intermediate (several wavelength) size

environment, e.g. the fields in one particular room rather than the average of

many rooms. We assume some volume of interest which is free of primary

sources, but which may contain conductors carrying induced currents. No

knowledge of the source (s) irradiating the volume is assumed. Instead we

2
assume knowledge of the field at some number of measurement points. In

principle the field measured can be the electric or magnetic field, or both;

but in this paper we assume it is the electric field. If both were measured,

it would require holding fixed spatial derivatives of the vector potential at

measurement points, as well as the vector potential itself. We also restrict

our attention to the single-frequency case. The extension to multiple

frequencies could be effected through superposition, with a concomitant

increase in measurement effort.

The basic idea of the method we consider is to attempt to numerically

find an approximate solution to Maxwell's equations which is consistent with

the measured values of the field at measurement points and which is also

consistent with the appropriate boundary conditions at the surfaces of any

conductors within or bounding the volume of interest. Although it is not the

way in which the computation proceeds, one can think of the problem as being

divided into two distinct steps. Since the geometry of a real situation may

well be quite awkward, the normal modes of the system would not be known and

would need to be determined. The second step would be to determine what

combination of modes was consistent with the measurements. This second point

merits some discussion. We are interested in obtaining information from as

few measurements as possible, and therefore in the typical case there will be

insufficient information to uniquely determine the field in the volume of

interest. This raises two complications. Of the many possible solutions

consistent with the measured values of the field, we must be able to always

find the "same" one, the solution with some identifying feature (preferably a

useful one). For example, we would want to always find the one with the

largest (smallest) average energy density, or the most probable one, etc. The

other complication is that, because there are many solutions, the matrix

3
equation obtained by discretizing the system involves a singular matrix.

Consequently, methods requiring the inversion of that matrix are inapplicable.

In this paper we develop and illustrate the approach in some simple cases.

The particular geometry chosen is a rectangular box with perfectly conducting

walls. This example was obviously chosen for its simplicity rather than its

practical relevance, but it is not so far removed from reality. Introduction

of a few holes for windows and doors makes it a model of a metal building or

room. Another simplification which we have already mentioned is the

assumption that only one frequency is present.

We consider two different numerical techniques for finding solutions to

Maxwell's equations subject to the constraint that the field assume specific

values at the measurement points. Both methods first impose a spatial grid or

lattice on the volume of interest and are formulated in terms of the field

values at the lattice points. The first method employs Hamilton's principle

and the action functional for classical electromagnetism. For field

configurations which are solutions of Maxwell's equations, the action is

stationary with respect to small variations of the field; and so the

appropriate components of the field are fixed at boundaries, the field is

fixed at the measured values at measurement points, and it is allowed to vary

everywhere else until a stationary point of the action is found. The second

numerical technique is to attempt a direct solution of the relevant

differential equations, which are the vector Helmholtz equation supplemented

by the condition that the divergence of the vector potential vanish. This can

be considered an alternate approach to locating the stationary point of the

action, since the action stationarity leads to Maxwell's equations and thence

to the Helmholtz equation and divergence condition in our case. Again the

field is set equal to measured values at measurement points, and the Helmholtz

4
equation and divergence condition are solved simultaneously using a successive

over-relaxation (SOR) method. Neither of the two numerical methods is

entirely successful, but the SOR method has enough success to demonstrate the

general approach and to point out avenues for development.

In the next section we present the first numerical technique, based on

the action functional. The SOR technique is applied to the Helmholtz equation

and divergence condition in Section 3. In the final section we summarize the

work and comment on future prospects.

2. STATIONARITY OF THE ACTION

2.1 Formulation

In terms of the usual vector and scalar potentials A(x*,t)

and cj>(x,t), the action functional for classical (lossless) electromagnet ism is

[18, 19]

,t 2 &= J t dt Jd 3 x {![(% + f^ A) • e(x) - (vcj, + §£ A)

- (V x A) • y _1 (x) • (V x A) + A • J(x,t) - <j> p(x,t)},

where "? and ti are the permittivity and permeability, and

J and p are the current and charge densities. The spatial integral extends

over the volume of interest, subject to the restrictions imposed by the

variations of the potentials vanishing on the surface (see below). We work

with the potentials rather than the electric and magnetic fields because

Hamilton's principle requires that the quantities being varied be independent

degrees of freedom, whereas E and H are related. The four potentials

(A,(})) cannot all be independent either, of course, but one can be easily

eliminated when the gauge choice is made.

5
In this formulation, the two homogeneous Maxwell's equations follow

-> immediately from the definition of E and B in terms of the potentials:

E = - (v^ + f- a) I ^xE+|-B = 3t 3t

(2.2)

B = ^xA ' V • B =

The two remaining Maxwell's equations follow from Hamilton's principle [10].

If we consider small variations of the potentials and requires that the

resulting variation of the action vanish, we obtain

(2.3)

6 A Cl = * V • D = p,

6 n = + 7 x H - |- D = J,

A

8t

-> provided the variations of A and <j> vanish at the end points. To be more

exact, what is required is that 6<j> and 6A vanish on the surface of the volume

V at all times and that 6k vanish throughout the volume at times t x and t 2

.

What this means in practice is that in order to use this in a calculation, we

must specify the fields on the boundary surface at all times and throughout

the volume at the initial and final times. In any real application, such a

superabundance of information will not be available for the volume of interest,

In order to exploit the stationarity of the action we must expand the volume

considered beyond the region of interest, out to distances where the fields

can be assumed to be negligible. The same applies to the time; t 1 is chosen

before the fields are turned on and t 2 after they are turned off.
For the calculations that follow, it will be convenient to work in a

definite gauge. A propitious choice is cj>(x*,t) = 0. The action then

assumes the form

d- Jdt jd 3 x |l[(f- t) • % • (|r t) - (? x 1) . ?-« • (v- x 1)]

2 9t 3t

(2>4)

+ A • J}.

This gauge has the convenient feature that for the single-frequency case a

measurement of the electric field directly determines A,

E(x,t) = a) A (x,t - ir/2w).

Having extended the volume under consideration to virtually all of spacetime, we must restrict the problem to manageable size. In particular, since

we do not know the sources, the volume of interest is chosen to exclude them.

We envision the division depicted in Fig. 2-1. The volume marked Vj is the

region of interest; it is assumed to be free of primary sources, but it may

contain conductors with induced currents. Volume V3 is a buffer zone

separating Vj from V5, which is the rest of space, wherein are located any

primary sources. The idea is to divide the action into one piece from the

integral over volumes Vj and Vg and another piece from V5. Because Vg

contains unknown sources, we do not try to determine the fields there, which

means that we also will not know the fields on the surface between V5 and Vg.

The fields on that surface will be allowed to vary or will be fixed by a

reasonable guess. Obviously, near the boundary between Vg and Vg the

solution obtained will be very sensitive to the choice for the fields on that

boundary, and therefore it will not be reliable. As we move away from the

outer boundary of Vg the values of the fields should be influenced more by the

measurement points and less by the values on the surface between Vg and Vg.

For points far enough away from that surface—i.e. within V T — it should be

7

x
possible to obtain reliable solutions, given enough measurement points. The

hope is that "far enough away" and "enough points" are not so large as to

render the method impractical for most applications. The positions of the

measurement points will clearly affect the size required for Vg; it may well

be advantageous to make a few measurements on the perimeter of Vg. It would

probably also be advisable to choose the boundary of Vg to coincide with

conducting walls when it is feasible, in order to constrain the fields on the

boundary as much as possible.

The quantity we shall consider then is a reduced action, (Z, which is

defined as in eq. (2.4) but with the integral restricted to Vj €> Vg. (In the

general multiple-frequency case the volumes are four-dimensional space-time

volumes.) The next step is to discretize the expression for the reduced

action, converting the volume integral into a summation which approximates it.

There are any number of ways to do so; we are not interested in their relative

(dis)advantages at this time. The discretized reduced action will have the

general form

y\

a dt I a,B,YeV

I@

V

{ AV

1 3 * aeriKft^'U-c,BY (ft*(t) U

(2.5)

- tf. t(t» • M- Y

• (*" tCt» I •

The discrete indices a,B,Y label the spatial points, on which are centered the

volume elements AVa gy We have left the time variable continuous for now,

anticipating the single-frequency example below. The quantities

(3A/8t) agy and (V * A)agY will in general depend on the values of the
field A at a number of points on (or within) the surface bounding AV a gy

A concrete realization of this discretization of & will be given in the

example below. Imposition of the constraints required by the known values of

E and/or H at measurement points can be rather complicated, but for only one

-> -> frequency and for our gauge choice E a gy <* Aa gy, and a measurement of

-> the electric field at a point fixes A a gy at that point. A similar comment

applies to boundary conditions at perfectly conducting walls; for the

single-frequency case they can be imposed with relative ease.

The calculation then proceeds as follows. A grid is defined within

Vg@Vj, and Aagy is fixed at measurement points and appropriate

-> components are set equal to zero at conducting walls. All other A^y's

are varied, and we search numerically for a stationary point of (X in

eq.(2.5). When (if) a stationary point is found, then that set of A a gy

constitutes an approximate solution to Maxwell's equations within VggVj

which is consistent with measured field values and boundary conditions. The

problem of finding the stationary point is nontrivial, however, and we now

turn our attention to that.

2.2 Nature of the Stationary Point and the Euclidean Action

The technique presented in this paper depends on finding field

configurations such that the action is stationary with respect to small

variations of the field. Such field configurations are solutions of Maxwell's

equations. The standard method of finding the stationary point would be to

discretize the functional as in (2.5), set 3#./3A^gy = to obtain a

(large) matrix equation, and solve that equation. That approach founders here

because the matrix would be (very) singular. There are in general many

possible solutions consistent with boundary conditions, measured fields, and

unknown sources. We need to locate a stationary point and to know which

9
stationary point it is —does it correspond to the solution with the largest

maximum electric field in the volume, or the solution with the minimum power

density, etc. To locate the stationary points numerically we need to know

whether they are maxima, minima, or neither. The phrase "Principle of Least

Action" may lead one to expect a minimum, but in general this need be true

only for infinitesimal time intervals At = t 2 - t x

. To consider our specific

case, we adapt the treatment of Whittaker [20]. We use the form of eq. (2.4)

and assume for simplicity that *£ and "ft -1 are scalar constants. If

the fields are allowed to vary

A (x,t) » A(x,t) + ct(x,t),

a(x ,t j

) = a(x ,t 2

) = 0,

where 3 is small, then the second-order variation of the action is

(2.6)

6

2 a = \ ,

t 2dt

t,

d

3 x [ea 2 - - ($ x a) 2 ]. (2.7)

V

The usual argument is that because 5 but not 3 vanishes at t x and

t 2 , for small t 2 -t x the first term dominates the integral in (2.7) and the

stationary point is a minimum. That argument fails in a field theory,

however, where the second term involves spatial derivatives of the dynamical

variables. Because the functional form of &(x*,t) is arbitrary, its

spatial dependence can oscillate wildly and 6

2 CL can be positive or negative

(depending on the choice of 3) for any fixed nonzero At. Therefore the

stationary points of are neither maxima nor minima.

This changes when the action is discretized by the introduction of a

spatial grid. Then the spatial derivatives of dt are bounded by

(V x &) 2 < 0(& 2 /a 2 ), where a is the lattice spacing. The stationary point

is a minimum provided At is less than some number on the order of /e\Pa. For

10
larger values of At, 6

2d can again have either sign, and consequently the

minimum is again a saddle point. That is the case of interest since we need

At + °° in order to set all the fields equal to zero at t x and t 2

.

This poses a serious calculational problem. A stationary point which is

neither a maximum nor a minimum will be very difficult to locate numerically.

In earlier work [21, 22] we introduced a method which succeeded in some simple

cases, but subsequent tests with more complex field configurations (i.e.

higher modes) have exposed inadequacies in that method. Faced with the

impracticality of numerically locating saddle points of functions of very many

(10 3 -10 6) variables, we need a function to maximize or minimize. Our efforts

to construct such a function have not been marked by complete success. One

interesting attempt which was partially successful is to borrow a trick from

quantum field theory and work in Euclidean space-time [23]. To do so we take

time to be imaginary, t = - it with t real, thereby changing the sign of the

(9A/8t) 2 term in (2.5) and (2.7). We solve the problem in Euclidean space and

analytically continues the result back to Minkowski space, x+it. The

Euclidean action is

V d

3 xfl (|_XE( ;, T)) .t.(f^E

dx ) Tl (2.8)

+ \

(v- x A

E

)

• y" 1 • (v- x A

E

) - A

E

• J} ,

and this is to be minimized to obtain the physical solutions in Euclidean

space. The analytic continuation back to Minkowski space can be a difficulty

in general, and in fact it will prove to be the limiting factor for this

technique. Nevertheless, it is possible to handle some simple examples in

this manner.

11
As noted above, there will be many possible solutions consistent with a

few field measurements if the sources are unknown. The Euclidean solution is

obtained by minimizing the integral of the energy density

(V2 e Eg + V2

\x Hg). As will be seen below, for the single frequency case

this also corresponds to minimizing the energy density in Minkowski space.

Therefore, the solution obtained will be the one with the minimum average

energy density in the volume under consideration. This is a welcome bonus of

working in Euclidean space--not only can we find a solution, but it is a

useful one as well.

The entire procedure, including the analytic continuation and numerical

search method, is greatly clarified by an example, to which the next subsection is devoted.

2.3 Simple Example

A. Problem and Calculation

Having presented the general ideas of this approach, we now attempt to

implement it in a simple example—a rectangular waveguide with perfectly

conducting walls. This is obviously not supposed to be a practical

application, and many of the difficulties and nuances of the general case are

absent. It is a practice problem to demonstrate the idea and to provide a

basis on which to build toward solution of real problems. The rectangular

waveguide is chosen because it has simple known modes, because the boundary

conditions are easily imposed, and because it can be reduced to a two-(or even

one-) dimensional problem, thereby reducing the computational exercise. As we

shall see below, even in this simple case, it is not entirely trivial to

obtain the correct continuum answer, with which to compare the results of our

numerical computation.

12
We first restrict the problem to two dimensions, obtaining the general

forms for Ag and A^ (the vector potentials in Euclidean and Minkowski space,

respectively) and the relationship between them. To facilitate the analytic

continuation, we assume a standing wave due to perfect reflection from the

(very distant) end of the waveguide. For a given TE mode, TEmn , A^ can then

be written

A (x,t) = A(xi) cos kz cos cot = A(]_) cos kz cos tot,

(2.9)

/to 2

2 rm

2 nS

where a and b are the transverse dimensions of the waveguide (aSb), c is the

speed of light in the waveguide medium, xl = (x,y) is the two-dimensional

transverse position vector, and A^(x,t) and A(J_) are real. Equation (2.9)

would allow a simple restriction to two dimensions for a single mode.

However, the solutions of interest are those with the minimum energy density

consistent with the measurements, and those solutions are superpositions of

many modes. Consequently the true solution will not have a z dependence which

can be factored out as in (2.9), and the problem will not reduce to two

dimensions. In order to force it to do so, we impose the form of eq (2.9)

with k taken to be that of the lowest mode consistent with the measurements.

That yields a two-dimensional problem, but it is no longer electromagnetism:

if the solutions A(J_) are substituted back into (2.9), the resulting

A(x,t) need not be a solution of Maxwell's equations. Nonetheless, the action

and the problem are quite similar to electromagnetism, and it is a suitable

problem for demonstration purposes since one can obtain the true solutions for

comparison. This arbitrarily imposed z dependence does not reflect a

13
shortcoming of the method, but rather the difficulty of finding a simple

two-dimensional application.

To obtain the general form for Ag, we require that it be real and that

it satisfy the Euclidean wave equation

(V 2 + U e |^r) A

e (x,t) = . (2.10)

It can then be written as

A_ (x,t) = a

R (_L) e + b„(J_) e (2.11)

> ,iv -i(kz + ojt) * , h -i(kz - wx) . + c E

(J_) e + d

E

(J_) e + compl. conj .

,

where u need not be (in fact is not) real. If (2.11) is then continued

back to Minkowski space (x = it) and compared to (2.9), the relation between

Ag and A^ follows. We can write

A (x,x) = A(J_) cos kz cosh ux

,

(2.12)

where A(J_) , k, and oj are the same as in (2.9). Equation (2.12) enables us to

use measurements of A^ to constrain Ag. Conversely once Ag is determined by

the computation, A(J_) and hence also A^ are known.

Substituting the form for Ag (2.12) into the Euclidean action (2.8) leads

to

a+A b + A a

E

= c dx

-A -A

U) 2 dy Kpr + k

2

) [a 2 (1) + a

2 (D] + ^a2 (1)

(2.12)

+ [(9 A (I)) 2 + (9 A (I)) 2 + (3 A (]_) - 3 A (I)] 2

] y z

-1-

' x z

-1- '

v x y

-1- y x

-1-

; J

-2y J(l) • A(]_)},

14
where J(j_) is defined in terms of the wall current J(x,t) in the

obvious way (9), and where the constant C is given by

T 2

2U . dx cosh ut dz cos 2 kz. (2.14)

In writing (2.13) we have assumed that the material in the waveguide is

isotropic and that the range of the z integration is either very long or an

integral number of cycles. The induced current and the limits of the

transverse integrations require explanation. In applying Hamilton's principle

the variations of A must be zero on the boundary, which in our two-dimensional

-> case here means that we must specify A on the transverse boundary. In order

-> to be able to specify all components of A we choose the boundary to lie a few

skin depths within the conducting walls where the field can safely be assumed

to vanish. In figure 2-2 the x integration goes from x = -A to x = a + A, and

the y integration from -A to b + A, where A = N6, some suitable number of skin

depths. Then, however, the currents induced in the walls are contained within

the integration volume. Fortunately, we can show that the contribution to the

action from the volume within the conductor is negligible, and we can write

a b

a

E

= c dx dy !(£ + k

2

) [a*<]_) * a

2 (1)] +^(1)

(2.15)

+ [(3 A (I)) 2 + (3 A (I)) 2 + (3 A (I) - 3 A (I)) 2 L^y ]}.

z J_; (. x z

J_ J V xy J_ y X i_;JI

We still need to impose E tan = and Bnorm = at the conductor walls, which

is accomplished by the requirements that A(]_)^an var>ish at the walls.

The action is then discretized by introducing a rectangular grid into the

waveguide as indicated in Fig. 2-3. The spacing between points is Ax = a/N x ,

Ay = b/Ny. Derivatives at a point (i,j) are defined by

15
[j_ A)ij . ian..i)-A(i„i)

t

except at the boundaries, where

X X ,

and similarly for y derivatives. The Euclidean action then takes the form

N N

/7 x y

"-„ = C X I AA.. {[A (i,j) 2 + A (i,j) 2

] c tt 2 E . . . _ lj lL x ' J y '° J mn

i=0 j=0

+ 4it

2 A (i.J>« + [(A (l.J + 1) - A (i,j))/Ay] 2 Z Z Z (2.18)

+ [(A

z (I+l,j) - A

z (i,j))/Ax] 2 + [(Ax

(i,j + 1) - A

x (i,j))/Ay

- (A

y

(i+l,j) - A

y

(i,j))/ Ax] 2 },

Cmn

=

^ 8

" f* " ^ '

where all lengths are in units of the free-space wavelength A. The area

elements are

AA = Ax Ay 0<k,KN

= V-jAx Ay k or 1 = or N, (2.19)

and the "out of bounds" fields are defined as

(2.20)

A(i,N

y

+ 1) = 2A(i,N

y

) - A(i,N

y

- 1),

A(N X

+ 1, j) = 2A(Nx ,j) - A(NX- 1,j).

Having discretized the problem, we next set the values of the tangential

fields to zero on the boundaries, set the fields equal to their measured

values at measurement points, and vary all other fields to minimize 6cg in

16
(2.18). The manner in which we minimize^ is to step through the grid, at

each point (i,j) setting the fields equal to the value required by

8 <3.g/3A a

(i ,j ) = 0, a = 1-3, given the current values of the field at

neighboring points. This process is repeated until the action reaches a

minimum and does not change significantly during further passes through the

grid. (This is just the Gauss-Seidel method; see e.g. [24].)

B. Demonstration

To further demonstrate the method and to clarify the sort of solution

obtained, we present a simple one-dimensional computation and compare it to

the continuum results. The field was fixed at three measurement points,

A

y (|,|) = 0.707, A

y (|,|) = 1.0, Ay

(-jp,|) = 0.707, (2.21)

with Ax and Az equal to zero at all three points. This field pattern

corresponds to the TE 10 mode, or any number of higher modes. The waveguide

dimensions were taken to be a = 2. 0, b = 1 . 0. The computation to find the

minimum is straightforward and proceeds as outlined at the end of the previous

subsection. Obtaining an analytical result for comparison is less trivial.

To do so for the three-measurement case, we assume the y dependence of all

fields is a constant and expand

00 Ay(x) = I a. sin(j~). (2.22)

J-1

J

Ax and Az are zero at the measurement points and are therefore zero everywhere

for the minimumU

E solution. In terms of the aj of (2.22) the action is

abTT 2 V C

T~ Ia

?

31 +

^ 2) - (2.23)

J

J

The constraints of (2.21) can be combined to yield

17
a i +

I < a 8n+1 " a 8n-l) = 1

» n=1

(2.24)

and two other constraints which lead to all a^ which do not contribute to

(2.24) being zero.

In order to minimize &e (2.23) subject to the constraint (2.24) we first

use the constraint to eliminate a x and then set derivatives of (Xe with respect

to all other a^ equal to zero. This yields the infinite dimensional matrix

equation

AT

1

1

1

At

1

1a;

i

1 .

1 ,

1 . At

bT

,(2.25)

where

bn = ± a 8n±1

A^ = 2 + 2n 2 ± n/2,

(2.26)

This system can be solved by subtracting each row from the first to get

the relation

± A, - 1

) = —i

A n

~ 1 b x

. (2.27)

The first row of (2.25) then determines b7, and a x can be determined from

(2.24). The result is

a 8n±1 = ± — 1 a x

, n>0, (2.28)

An " 1

18
.. - [1 I (-^-^rJ-)]

-

n=1 An

- 1 An

4tt = j^ tahn [1 /3T1

) = 0.43530,

Equation (2.28) and the fact that all other a^ = can then be used in (2.22)

to determine the true solution to which to compare the result of the numerical

computation. Figure 2-4 shows the comparison, and it is obvious that a) the

two sets of results agree with each other, and b) they agree with the measured

values at x = a/4, a/2, and 3a/4.

2.4 Comments

The simple application in the preceding subsection illustrates both the

general lattice approach to complex environment characterization and the

attempt to implement this approach using Hamilton's principle and the action

in Euclidean space-time. The example demonstrates that the lattice

computation does produce the correct results. It indicates that in one

dimension, for a system very similar to Maxwell's equations, we are able to

find the minimum-energy-density solution consistent with measured field values

at a few points, without knowledge of the external sources.

Unfortunately, a fundamental obstacle arises when less trivial examples

are attempted. If there is any spatial variation in the phase of the field,

or if the different components of the vector potential are not all in phase

with each other, it becomes impossible to perform the analytic continuation as

-> required. Since we are taking A to be real, phase variation appears as an

admixture of cos mt and sin wt time dependences. One cannot continue both cos

19
and sin to imaginary time and still maintain the reality of A (which is

necessary in order to keep the Euclidean action positive definite).

Consequently, it does not appear that the trick of analytically continuing to

imaginary time will be applicable in practical problems.

That leaves us with the problem of finding a stationary point of the

action which is neither a maximum nor a minimum. Until that numerical problem

can be solved, the action-based approach holds little promise.

3. DIFFERENTIAL EQUATION SOLUTION

3.1 Method

As in the preceding section we assume a single frequency, and we continue

to choose the gauge so that the scalar potential vanishes. We then have

<j>(x,t) = 0, A(x,t) = A(x) e"

lut

, E(x) = ia)A(x), kx) = V x fox). ( 3' 1 )

Equation (3.1) guarantees that E and B will satisfy the two homogeneous

Maxwell equations. In a sourceless region of free space the two remaining

Maxwell equations can be written

V • A(x) = 0, (3.2a)

V

2 A(x) + o) 2 ue A(x) = 0. (3-2b)

We continue to measure all distances in units of the free space wavelength

A = 2tt/(u> /pe), so that w

2 ye = 4tt2 in (3-2b).

The method presented in this section is to discretize the differential

equations of (3.2) and then attempt to numerically solve them simultaneously

subject to the constraint that the field is fixed at measurement points.

Using finite difference forms for the derivatives, eqs (3.2a) and (3.2b)

become

20
Ax

( i + 1 , j , k ) -Ax (i-1,j,k) +A

y (i,j+1,k) -A

y (i,j-1,k) +Az (i,j,k + 1)

-Az (i,j,k-1) = (3.3a)

[A(i+1,j,k) + A(i-1,j,k) + A(i,j+1,k) + A(i,j-1,k) + A(i,j,k + 1)

+ A(i,j,k-1) - 6 A(i ,j ,k) ]/A 2 + 4rr

2 A(i,j,k) = 0, (3. 3b)

where we have introduced a rectangular grid with equal spacing A in x,y, and

z, and the positions of grid points are given by (x,y,z) = (i,j,k)A. As in

the previous section the cavity will have dimensions a,b,c, and (Nx »N

y

,N

z

) =

(a,b,c)/A are the number of intervals in each direction. Points on the

boundaries require special treatment since they lack neighbors in one or more

directions. There we use the equation at the neighboring interior point. For

example, the field at a point (a,y,z) = (Nx ,j,k)A is determined from

Ax (Nx ,j,k) - Ax (Nx -2,j,k) + A

y

(Nx -1,j+1,k) - A

y

(N x

-1 , j -1 , k)

+ Az (N x -1,j,k+1) - Az (N x -1,j,k-1) = 0, (3.4a)

[A(Nx ,j,k) + A(Nx -2,j,k) + A(Nx -1,j+1,k) + A(NX

-1 , j-1 ,k) + A(NX-1 , j , k+1

+ A(Nx -1,j,k-1) - 6A(Nx -1,j,k)]/A 2 + 4tt

2 A(Nx -1,j,k) = 0. (3.4b)

In general there will be more than one possible solution. A unique

solution would require measurements at sampling theorem separations (< A/2 for

infinite free space), or knowledge of the sources, or knowledge of

A on the boundary surface of the volume. This restricts the range of

numerical methods available. To see this, consider just the vector Helmholtz

equation, (3.3b). By arranging all field components A x (i,j,k), A

y (i,j,k),

Az (i,j,k) into one long vector z, one can put (3.3b) into the form of a matrix

equation Bz = 0. Those elements of z corresponding to measurement or boundary

points are known and can be brought over to the right hand side to yield an

equation of the form Ax = b. If there is more than one solution, then A is

singular and any method requiring the existence of A

-1 will fail.

21
This leads us to adopt an iterative technique —the successive

over-relaxation (SOR) method [24]. For simplicity, consider first just the

vector Helmholtz equation, (3. 3b). To apply SOR to this we step through the

grid, and at each point (i,j,k) we first compute the residual, defined by

RH (i,j,k) = A(i,j,k) - [A(i+1,j,k) + A(i-1,j,k) + A(i,j + 1,k) ,

. > * +

(3 ' 5)

+ A(i,j-1,k) + A(i,j,k + 1) + A(i,j,k-1)]/(6 - 4u 2 A

2 ).

(This is actually the residual divided by (6 - 4ir

2 A2 ), but we shall call it

the residual.) We then change the vector potential at that point according to

A'(i,j,k) = A(i,j,k) - ft H

RH (i,j,k), (3.6)

where &# is a number between and 2, chosen to optimize the rate of

convergence. For Q,^ = 1 , this is the Gauss-Seidel method, in which the

-> function (A) is set equal to the value required to satisfy the difference

equation(3«3b) , given the current values of the function at neighboring points,

For fi H greater (less) than one, the method is one of over-(under-) relaxation.

The term SOR is used generically to apply to both over and under relaxation

(and Gauss-Seidel as well) [24],

Application of the SOR method to the divergence condition (3. 3a) is

somewhat different due to the fact that (3.3a) relates different components of

A at different points, and we would like to treat all three components the

same. This can be accomplished in the following manner. We define the

residual as

R

v (i,j,k) = -Ax

(+) + Ax (-) - A

y

(+) + A

y

(-) - Az

(+) + Az (-), (3.7)

where

Ax (±) = Ax

(i ± 1 ,j ,k),

Ay

(±) = A

y

(i,j ± 1,k), (3.8)

A z (±) = A z

(i ,j ,k ± 1 ) .

22
The relaxation is then performed by spreading the residual equally over the

Aa (±)'s,

A^(±) = A

x (±) ± \

ft y

R v (i,j,k),

A*(±) = A (±) ± 1 n 7

R

v (i,j,k), (3.9)

A^,(±) = A

z (±) ± \ Qv

R v (i,j,k),

where the relaxation parameter fty e(0,2) is again chosen to optimize

convergence.

We have discussed separately the SOR method for the vector Helmholtz

equation and for the divergence condition, but in practice we want to solve

the two simultaneously. This is done in the obvious manner: at each point

(i,j,k) we first perform the divergence relaxation (3.9), thereby changing the

field at neighboring points, and then perform the Helmholtz relaxation (3.6),

changing the field at (i,j,k).

One final point we need to address before proceeding to examples is the

question whether this procedure converges. It does not. Considering just the

Helmholtz equation, even when all the boundary values are specified so that

there is a unique solution, SOR can be shown not to converge [24,25]. The

reason can be seen heuristically from eqs (3.5,3.6). For simplicity consider

Qpj = 1 (Gauss-Seidel) . Then if the sign of 4tt

2 A

2 were positive in the

denominator, SOR would replace the value of the field at a point by a number

somewhat smaller than the average of the field at neighboring points, and the

process would converge. With the negative sign, however, the field at the

point is replaced by a value larger than the average of neighboring values,

and the process does not converge—the fields just keep growing.

23
Why then do we use it? There are two reasons. The first is that it is

still possible to obtain approximate solutions by truncating the procedure at

some point even if it does not converge. If we monitor the volume average of

the residual, we can use as our result the field configuration corresponding

to the minimum average residual. This truncated SOR will be used in the next

subsections, with some success. The other reason for using this rather

suspect procedure is that there is not an obviously better one available which

suits our needs—i.e. does not require inversion of a singular matrix in the

problems of interest. The (scalar or vector) Helmholtz equation just happens

to be awkward numerically, due to the fact that when it is written in the form

Ax = b (cf. the discussion preceding eq (3-5)) the matrix A is not positive

definite, even when full information is given on the boundaries. How to

overcome this difficulty is the subject of current research (see e.g. [26] and

references therein), and in the summary we shall comment upon possible future

improvements on our scheme. For the present examples we use the truncated SOR

method.

3.2 Two-Dimensional Examples

In this subsection we present two-dimensional examples which were used to

develop the procedure and test dependence on various parameters. All the

examples assume a cavity of rectangular cross section with perfectly

conducting walls with dimensions axbxc. The TE£mn mode of the cavity is

given by

2,ttx . rrmy . nuz

A = cos sin —~ sin , x a b c

S,b . &7TX m-ny . nuz , ir>\ A = sin cos —r~ sin , (3.10)

y ma a b c

A z

= 0,

24
where we have chosen unit amplitude and assumed m * 0. The dimensions will

always be chosen to correspond to resonance of the mode used to specify

"measurement" results, (t/a) 2 + (m/b) 2 + (n/c) 2 = 4. This is necessary in

order to have any solution at all in the enclosed cavity in the absence of any

source. In these preliminary examples we restrict ourselves to cases where

only one component (Ax ) of the vector potential is nonzero, and then in order

to restrict the problem to two dimensions we impose constancy of A x in the x

direction. The divergence condition (3.2a) is then automatically satisfied,

and we need only treat the Helmholtz equation (3.2b).

In the first computation we assume the field is measured at one point,

the center of the cavity, where Ax

= 1.0, Ay= Az= 0. The cavity dimensions

are taken to be a = 2b = 2c = /2~7 which corresponds to the resonant frequency

for the TE 011 mode. The tangential components of A are fixed at zero ateach

wall, and Ax is required to be constant in x. The relaxation method outlined

above is used, and the iteration process is terminated when the magnitude of

the Ax residual averaged over the volume, <|R X |>, reaches a minimum. This

first example is atypical in that <|RX |> decreases asymptotically to a

constant, and so our solution is the field configuration any time after <|R X |>

has reached its limit value, which is about 0.4 x 10" 1*. A value of Rjj = 1.75

resulted in the fastest convergence to the limit.

The y dependence of Ax for two different values of z and two different

grid sizes is shown in fig. 3 - 1 . The solid and dashed curves represent the

field configuration for the TE 011 mode, which is the only mode resonant at

this frequency which is consistent with the measurement. It is evident that

the 16x8x8 grid is sufficiently fine for our purposes. The z dependence of A x

for fixed y is virtually identical to the y dependence for fixed z and

25
therefore is not shown. The final results for A

y

and Az are zero everywhere,

as they should be.

At the beginning of the computation, initial values of the field at all

points must be assigned. Two different starting configurations were tried.

In one, A is initially zero at all points except measurement points, where of

course it assumes the measured values. In the other, a smooth initial

configuration is generated by an algorithm which interpolates between

measurement points and boundary points. For this example, the same results

are obtained for either initial configuration. In general the smooth initial

configuration tends to yield a lower minimum <|RX |>, anc^ we snaH use it in

the following examples unless otherwise noted. We expect that when multiple

solutions exist the smooth starting configuration will result in the smoothest

of the possible solutions. This should correspond to the lowest mode.

If the measurement point is moved from the center in this TE 011 example,

the quality of the solution deteriorates. For example, if the measurement is

made at (y,z) = (b/4, c/2) then the computed solution at the center is about

6% higher than the correct solution, as shown in fig. 3~2.

Turning to a more complex field configuration, we consider the case of

four measurements with Ax

= 1.0 at (b/4, c/4) and (3b/4, 3c/4) and Ax= -1.0

at (3b/4, c/4) and (b/4, 3c/4). The lowest mode consistent with those values

is the TE 022 mode, and we choose a=b=c= /^so that it is resonant. In this

case we encounter what proves to be the typical behavior of <|R X |>. It

decreases with successive iterations at first, but it levels off, begins to

increase, and eventually blows up. The value ti^ = 1.75 again results in the

quickest minimum, and other values of Q# do not result in a significantly

better solution. We use as our approximate solution the field configuration

for which <|R X |> is a minimum with fl^ = 1.75. The results of figure 3~3 then

26
follow. The solid and dashed curves represent the TE 022 mode, and again the

agreement is very good. The same is true for the z dependence, which is not

shown.

Equally good results were obtained using four measurements representing

the TE 021 mode, using a cavity with dimensions a=b=2c = /21

. The measurement

points were taken to be the same as in the TE 022 case above, with the measured

values different of course. Figure 3

-iJa compares the computed A x to the TE 021

form as a function of y at z=c/2, and fig. 3 - ^b does the same as a function of

z for y=b/4. A 16 x 16 x 8 grid was used, with Jfy = 1.75

The results thus far have been very encouraging. Instructive

difficulties arise when a=b=c= /0 . 8

,' however. If we use as measurements the

TE 021 field at two points, (y,z) = (b/4, c/2) and (3b/4, c/2) , the agreement

of the computed results with the TE 021 field shape is not satisfactory. This

is demonstrated in fig. 3~5. In this and succeeding examples, the largest

dimension is divided into 16 divisions, with other dimensions divided

proportionately. Thus, this grid is 16x16x16. Away from the measurement

points the computed field falls off too rapidly. The problem persists when

the number of measurement points is increased to four, as is evident from fig.

3 -6. The value of the minimum <|R X | > ^ s somewhat larger than the TE 022 and

rectangular TE 021 cases: 0.25 x 10~ 2 as opposed to 0.20 x 10~ 2 and 0.11 x 10~ 2

respectively.

It is not completely clear what cause underlies these difficulties. The

numerical method being used does not actually converge to the true solution;

it approaches it and then wanders off, with the approximate solution being

obtained by taking the fields when "closest" to the true solution (as

determined by <|R X | > )« I n the present case (TE 021 with square cavity cross

section), the iterative process just does not get that close to the correct

27
result. The most likely reason appears to be that there are two possible

resonant modes with A

y

= Az= at this frequency, TE 021 and TE 012 , and the

computation does not pick out one pure mode.

3.3 Three-Dimensional Examples

We now allow variation in the x direction and treat the fully threedimensional case. This requires simultaneous solution of both the vector

Helmholtz equation (3.3b) and the divergence condition (3.3a). Again we start

with the simplest possible field configuration, assigning values corresponding

to the TE 011 mode at measurement points. Six measurement points were used,

located at (a/8, b/8, 7c/8) , (a/4, 7b/8, c/4) , (3a/8, b/2, c/2) , (5a/8, 3b/4,

3c/4), (3a/4, b/4, c/4), and (7a/8, 7b/8, 7c/8). The points are scattered

irregularly to cover the entire volume and to avoid having them all at zeros

of some low order mode. The box dimensions are taken to be 2a=b=c= 1 //2] and

an 8x16x16 grid is used.

The computation proceeds in much the same manner as the preceding ones,

except that now the grid has three instead of two dimensions, and the

divergence condition relaxation (3.9) is also performed at each point. The

quality of the solution depends on the choice of the relaxation parameters, %

and fty. The minimum value of 
 decreases as ft^ increases, but very

slowly, whereas the minimum value of 
 is smallest for small Q^ and

increases rather rapidly as % increases. The minima of both residuals are

smallest when Qy is large. The optimal choices seem to be around Q^ = 0.2,

fty = 1.75, which are the values we used. Representative results of the

computation are shown in fig. 3~7. The y dependence is satisfactory, but only

satisfactory, and a similar comment applies to the z dependence (not shown).

For the most part, the x dependence is very good—the computed A x is very

28
nearly flat. The one exception occurs along the line through the center of

the yz plane. Significant deviations from constancy occur, particularly at

the measurement point x/a = 3/8 where Ax is constrained to be 1

.

The TE 011 results are not bad in themselves, but they do not bode well

for more complicated cases, which will have greater spatial variation and more

than one nonzero component of A. The forebodings are borne out when we

usemeasurement results corresponding to the TE 022 mode. With six measurement

points we were unable to obtain acceptable results. In fact for some values

of £2fj and Qy the average residuals 
 an d 
 increased steadily

from the start of the computation—the best field configuration was the

initial guess. Increasing the number of measurement points to eleven did not

improve matters significantly. Also, the difficulties persisted when we used

different indicators of solution quality, such as max(R) or <R>

.

We conclude that in general three-dimensional problems, given just a few

measurement points, solving the Helmholtz equation and divergence condition

using the truncated SOR method fails to obtain useful approximate solutions.

To investigate whether the numerical method is more successful when the field

is specified at a large number of points, we did computations with all

components of the field fixed at all boundary points. The boundary points are

chosen because they can be thought of as driving the solution; when the field

is known at all boundary points there is a unique solution. We wish to

investigate whether our truncated SOR method can obtain something like the

correct solution and, if so, whether it can do so when some of the boundary

points are not given.

We have obtained results for a number of low order modes. A good example

is the TE 123 mode, which was the highest mode considered. For convenience, a

composite residual was defined by

29
R(a) = a <|R

R

|> + (1 - a)<|R

y

|>, (3.11)

and the solution was chosen as the point at which R(a) was a minimum. A

number of values of a were used, and the results were relatively insensitive

to the exact value. In the computations we used a = 0.5. Figure 3~8 compares

the computed Ax and A

y

(Az is zero) to the correct TE 123 results as functions

of y at x = a/4 and z = c/8 and c/2. The agreement for Ax is seen to be very

good in both cases, but the A

y

results fail to reproduce the extremum at the

center. In order to check that this is not due to some programming clumsiness

which introduced a spurious asymmetry in the treatment of Ax and A

y

, we also

computed the TE 213 mode with dimensions a and b interchanged. We found that

the Ax and A

y

results were also just interchanged, as should be the case. The

difficulties apparent in fig. 3~8 can therefore be taken as indications of

deficiencies in the truncated SOR method rather than the consequence of

(obvious) programming blunders. The poor results at central y will also be

reflected in the x and z dependence at fixed y. This is seen in figure 3~9,

where the (x,y) = (a/4, b/8) results are very good but the (a/2, b/2) result

is less than half the correct answer.

The TE 123 mode constitutes a rather demanding test. There is different

spatial variation in each of the three directions, two components are nonzero,

and the components have different spatial dependence. The numerical solutions

are quite good over most (-7/8) of the volume but fail in the central region.

A finer grid led to only slight improvement. When the fields are completely

specified on the boundary, this truncated SOR method could be useful in

generating a starting configuration for some other method which converges (as

it was used in [26]), but it does not produce satisfactory results by itself.

30
Finally, if the fields are specified only at one-fourth of the boundary

points (tangential fields still vanishing at all boundary points), no minimum

of R is obtained for the TE 123 mode. From the start of the computation each

iteration increases R. In principle we could fix unmeasured boundary points

by interpolation and might then obtain a minimum for R, but at this time that

appears to be beating a dead horse.

4. DISCUSSION AND CONCLUSIONS

We have suggested an approach to the characterization of EM environments

which is based on the numerical solution of Maxwell's equations subject to the

constraints imposed by boundary conditions and results of measurements of the

field at a relatively small number of points. Some success was achieved in

the simpler examples used to demonstrate the approach, but results in more

complicated three-dimensional examples were less impressive. The comparisons

made were of the field configurations, which constitute a far more demanding

test than does a global quantity such as the average power density. The

results of a less stringent test could well be less disappointing.

Nevertheless, it is clear that the method is not yet ready for real

applications. The basic obstacle to practical use of the approach is the lack

of a fully successful method for obtaining numerical solutions of Maxwell's

equations under these conditions.

Two different numerical techniques were investigated. The one based on

the action functional and Hamilton's principle founders due to our inability

to numerically locate a stationary point which is not an extremum. Having

previously tried a Gauss-Seidel type of procedure [21,22] and having found it

lacking, we performed an analytic continuation to imaginary time. This has

the advantage that the action then has a minimum. Furthermore, the solution

31
obtained is that with the minimum energy density consistent with the

measurements. This technique, however, proved to be inapplicable to cases in

which there was any spatial variation of the phase or any phase difference

between different components of the field. We therefore abandoned the attempt

to develop a numerical technique based on the stationarity of the action. If

one did insist on pursuing this line further, one possibility would be to

construct and numerically minimize the function defined by

2 2 2

l U 3A

v (i,j,k) J l 9A (i,j,k) J l 3A (i,j,k) J J * ^* 1J

i,j,k Y

D

2 is a nonnegative function whose zeros correspond to stationary points of

the action and therefore to solutions of Maxwell's equations. Obvious

problems with that tactic are that the construction of D

2 is quite awkward and

that there is no guarantee that it can be easily minimized.

The second technique we studied in detail was use of successive

over-relaxation (SOR) to simultaneously solve the vector Helmholtz equation

and the divergence condition, which are equivalent to Maxwell's equations in

the cases considered (sourceless, free space, one frequency). Since the SOR

method does not converge for the Helmholtz equation, the sequence of

iterations was truncated when the volume average of the magnitude of the

residual was a minimum. This produced good results in a number of cases but

failed in others. Even when the field was completely specified on the

boundary, the results were not always satisfactory.

There are various other possible numerical procedures which come to mind.

Mittra and coworkers at Illinois have had some success using a modified GaussSeidel method on problems of the form Ax = b where A is not positive definite.

The modification they make is to introduce a convergence factor which prevents

32
the Gauss-Seidel iterations from diverging [27]. This is similar to our

truncated SOR method, with the truncation replaced by a gradual damping (and

with ft = 1). It does not appear likely that the gradual damping would cure

the ills of the truncation method. The divergence problem is obviated by the

truncation method too. The difficulty is that the iterative procedure does

not approach very close to the correct solution in some cases, and the

introduction of a convergence factor does not address that problem.

A more promising alternative would be to try a truncated conjugate

gradient procedure, as has been used in ill-conditioned deconvolution problems

[28]. Some work would be required to adapt the procedure to our case, in

particular to develop a criterion for when to truncate the iterations. If the

method can be applied, it could prove very useful since proponents of the

procedure claim it yields a minimum-norm solution. For our problem that

corresponds to the minimum energy density solution.

One seemingly relevant method which proves inapplicable is the maximum

entropy method (MEM), as reviewed for example in [29]. It is highly effective

for reconstructing images which have been blurred, contaminated by noise, or

from which some fraction of the information is missing. Our problem falls

into the missing information category. Unfortunately, for the MEM to work it

requires several pieces of information per unit cell, where the cell size

corresponds to the volume over which the information at a point has been

smeared. Since there is no blurring in our problem there is no information

except at the measurement points, and the MEM is powerless to fill in the

empty space.

A final possibility we mention is to put the Helmholtz equation in the

form Ax = b, as described in Section 3, and then multiply by the adjoint of A

(A t

) to obtain

33
A

f Ax = A+ b. (4.2)

If all boundary values are specified the SOR method will converge for this

equation since the matrix A A is then positive definite. Even for

incomplete information A A will not have negative eigenvalues (although some

will be zero), and so SOR should be less unstable. As with D

2 of eq (4.1),

however, A A is not so simple to construct. Furthermore, the rate of convergence

for iterative solution of eq (4.2) would be expected to be very slow. (This

could be improved by using SOR on Ax = b to precondition x.)

In summary, neither of the numerical methods we investigated was

successful in all cases. There do exist other methods which may warrant

further investigation. The general approach (constrained solution of

Maxwell's equations) sounds attractive for characterization of complex

environments, but its practical implementation awaits development of better

numerical methods.

5. ACKNOWLEDGMENT

We are grateful to John Gary of NBS for a number of patient conversations

and helpful comments.

6. REFERENCES

[1] Lerner , E.J. Biological effects of electromagnetic fields. IEEE

Spectrum, May 1984, 57.

[2] Beckman, P. Probability in communication engineering. New York, NY:

Harcourt, Brace, and World, 1967.

[3] Skomal, E. Man-made radio noise. New York, NY: Van Nostrand Reinhold,

1978.

34 /
[4] Middleton, D. Introduction to statistical communication theory. New

York, NY: McGraw-Hill, 1960.

[5] Middleton, D. Statistical models of electromagnetic interference, IEEE

Trans. EMC- 19, no. 3, 106-127, August 1977.

[6] Berry, L. Understanding Middleton's canonical formula for Class-A

noise, IEEE Trans. EMC-23, no. 4, 337-344, November 1981.

[7] Middleton, D. Canonical and quasi-canonical probability models of

Class-A interference, IEEE Trans. EMC-25, no. 2, 76-106, May 1983-

[8] Kanda, M. Time and amplitude statistics for electromagnetic noise in

mines, IEEE Trans. EMC-17, no. 3, 122-129, August 1975.

[9] Shepherd, R. Measurements of amplitude probability distributions and

power of automobile ignition noise at high frequency, IEEE Trans. VT-23,

no. 3, 72-83, August 1974.

[10] Nakai, T. ; Kawasaki, Z.I. Automotive noise from a motorway: I,

measurements, IEEE Trans. EMC-26, no. 4, 169-174, November 1984.

[11] Nakai, T. ; Kawasaki, Z.I. APD's and CRD's for noise from bullet trains

in Japan, IEEE Int. Symp. on EMC (Tokyo, Japan), 250-255, October 1984.

[12] Hoffman, H.H.; Cox, D.C. Attenuation of 900 MHz radio waves propagating

into a metal building, IEEE Trans, on Antennas and Propagation, Vol.

AP-30, No. 4, 808-811, July 1982.

[13] Cox, D.C; Murray, R.R.; Norris, A.W. Measurements of 800 MHz radio

transmission into buildings with metallic walls. The Bell System

Technical Journal, Vol. 62, No. 9, 2695-2717; November 1984.

[14] Haskard, M.; Miller, M. ; Johnson, A.; Marconi, P. Radio noise

measurements in an urban environment, J. Elec. Electron. Eng.,

Australia, Vol. 2, no. 2, 94-102, June 1982.

35
[15] Struzak, R. ; Rymarowiez, Z. ; Moron, W. Urban MF radio noise survey in

Poland, Proc. 5th Symp. and Tech. Exhibition on EMC, 1-6, Zurich, 1984.

[16] Parsons, J.; Sheikh, A. Statistical characterization of VHF man-made

radio noise, Radio Elec. Eng., Vol. 53, no. 3, 99-106, March 1983.

[17] Randa, J.; Kanda, M. Directional Scanning of complex electromagnetic

environments, IEEE Trans. AP-33, no. 12, 141 3-1 41 6, December 1985.

[18] Landau, L.D. ; Lifshitz, E.M. The classical theory of fields, Chapter 4,

Reading, MA: Addison-Wesley, 1962.

[19] Morse, P.M.; Feshbach, H. Methods of theoretical physics, Chapter 3,

New York, NY: McGraw-Hill, 1953.

[20] Whittaker, E.T. Treatise on the analytical dynamics of particles and

rigid bodies, Chapter 9, New York, NY: Dover, 1944.

[21] Kanda, M. ; Randa, J.; Nahman, N.S. Possible estimation methodologies

for electromagnetic field distributions in complex environments, NBS

Tech. Note 1081, March 1985.

[22] Randa, J.; Kanda, M. A finite-element action approach to the

characterization of complex electromagnetic environments, Proc. Int.

Symp. on Antennas and Propagation, p. 48-53, Beijing, August 1985.

[23] See, e.g. Itzykson, C; Zuber, J. Quantum field theory. New York, NY:

McGraw-Hill, 1980.

[24] Hageman , L.A.; Young, D.M. Applied iterative methods. New York, NY:

Academic Press, 1981.

[25] Varga, R.S. Matrix iterative analysis. Englewood Cliffs, NJ

:

Prentice-Hall, 1962.

[26] Bayliss, A.; Goldstein, C.I.; Turkel, E. An iterative method for the

Helmholtz equation, J. Comp. Phys .

, Vol. 49, 443-457, 1983.

36
[27] Mittra, R. Electromagnetic Communication Laboratory, Univ. of Illinois,

Urbana, IL 61 801 ; private communication.

[28] Sarkar, T.P. et al. Deconvolution of impulse response from timelimited input and output: Theory and experiment, IEEE Trans. IM-34, no.

4, 541-546, December 1985.

[29] Gull, S.F.; Skilling, J. Maximum entropy method in image processing,

IEE Preceedings, Vol. 131, Pt. F, No. 6, 646-659, October 1984.

37
Sources

Fields = on surface

Figure 2-1 Schematic division of total volume into volume of interest (Vj)

,

buffer volume (VB )

, and remaining volume containing all sources (Vs)

y ^

c

f (a+A, b+A)

(-A, -A)

Figure 2-2 Rectangular waveguide dimensions and axes. Also depicted is the

boundary of the volume over which the action integral extends.

38
¥-

x

CO

COXXX

X

X

XX

XX

-K *- COo o

COoo

CO

CD •H

CO

CDo

-p

-p

COo

bOc

CO

enI CM

CD3

bO

•r-t

39
: X E33 rH C fO

•-l c

4-> o c •r-H o CO o c

0)

1) s £ •1-1 4-J -o £ 1CD

4-> c

•<H o 3

CD /—

.c

CO •p

0)

to

L

03

O o «H LO **—

^—^ CDC

4-> •rH rH rH 3

CO "O

<D •r-l S- t-H o C CO O »

* •rH -P s—^ o oo

n) CO i • aj CO o M

•rH (\J

4-3 CM •P • CO C\J rH "—

<tH CO O cr

CD

Co Cw

co

O

•-1L -t-> cd rH • a 3 CD E CO CO o CD cfl o S- o

^r1

CDu3M

b-,
Mo

CmO CD H -U

-tJ CO

fC

O

LO

Cm -o c cu O CO O 3

"O CO

<-\ 03

CD 5 M

Cm 4-> C ^ .M - O o (X

4->

4-J c

CDE

CD

s_3

co

03

X'

E3

CDCo

c\j o

CO

4-> rH3

COCDU

•a

CD3aEOo

CmOco

CO •M Cm £5

S_

cCQ.eooI

t,3

bO M

4-3 c

CDC

3) CM Cm

C\J

41
C\J \

CO P

CO -Pc

•rH O

Q.c

0)s

CD

S-3

COCO

CDS • CO

SZ

x

-P CO

•rH

X

d ra .o co5

T— o 1 CD on m3 • bO T3

•i-( tH

<i-i

L

bO

C

•H o£

CO

H

COc • o •—

CO CM

•rH \

SO

05iX » e ^r o \ o £1

CMCO

t,

=5bO

[X,

42
d

ood

d

to

o

d^;

d

d

CVi d

^

vp oq p cvj d d d fH «-.

i

I I

x: CO

4-> 3

•i-l 3 CDN

N •.-i 03 CO

4-> c •a

•.-1 •-( o l

o. bO

4^> 0) c .c

<D H

E

Q)

i- • 3 =r

CO \

CO o

0) m

e

l =r 3 \

o o

<i-i v--' ii

co N

4-> .H 4-J a CO

CO

J) C

S- o

•l-l •o 43

CD CO -P C

3 D

ft bO s H

o *H o C

o

<M a

o a

c rH o OJ • CO •i-l v£>

•rH 0-1 r- L X

CO CMVO ft <M «— E O X O Cri MD o E- ,-

CO

CO

1 on

0)

J-3,

u,

43
(7> d

x: 00

4-> 5 d ^ cq

CO n3 P 3 C

•H CD O N

r^ CO CD 43c -a-H E L

CD bO

S.3 0)

CO .C ^o CO H

CD

d e

S_ CO 3 "N. o o

Cw LO ^

^-' II

CO

N

d >^

4-> iH 4-> 3 CO

COCD

C

JO

•i-i •O 4-> <vj- CD CO 43 L d 3 3 a ho E -i O "m o co

<i-l

O

00

O

a

C iH

d

O <D . CO -.-1 vO

tH <H t-

£.

x

CO «VD

Q.

N «— E ° x C\J O U^D O E- <- d ^-. -Omi on r-H

CD

d 3

CO

u.

4U
x)

CMO

IIN

X / XyS

X^

X^-

XyS

x/

G>O *3

•/H 2

^-v CO

CO CO 00 4-' 5 O •H 1) O N

O, iH

CO

4-2c -a

0) -.-i r^

J_ o 3 1) CO .c

coH

0)e • e\j L O \

CO 3N

O

O

CD C|_< II 11 jO N

CO

J-^ t_ 4-i H O fl 3 4-i LO £

CO >> L O <M o >^ •H O

T> J-3 0) CO c

JJ

L

O 3 3 -H

Q, b0 JJ S -H O

^JO <w c O C 3 o O CO •o «c C rH o a; CO -r-t • •rH Q-H CO CO C x

CO -<o

d S ° x o w vn O E- «-

^_, C\J

r0 ^r o I oo

CD

£-3

CO

fc*

00 CD ^CM

O C\J ^t CD CX) T- o O o O O1 o1 O1 o1 o1

45
-^ CO 00 Oj -P 5 C

•H 0) O N

Q, —I

CO PCs3

CO

cfl 0)

bO

CD 3o

II II p

05

C

L

O 3 th

bO POc3

(in

CO

CO p

rH

00

•a

a)

4-5 Q.soo

CmOCo

00 •") • •rl IHOO

03 «MD e ox O E-" r- I3

bO

•rH

(X4

H6
X \X

X \

CMO

II

X / N X

x^s

X/ X

X - _ X X

*wX9 (

G)o

00oO

CDo

to ^ d

CO

(5

CMo

COo

CDo CD

CMO

CM ^f CD 00

C

O

O

O

O i- I I I I I

P>

•M2

co ro C >>

•MO Vh

CX

O P C C O

<D -H E 4-> CD O

S_

C 3 3

CO 5h

CO

CD CO E < O2 • P O —

II JO

CO P S- rH O

CO

CD

C

J_ o

"O P

CD CO P f_ 3 3

Q, bO E

CmCoo

TDH

CD

OO

Cmoco

CO -M • •M <M C\J

COQ.EOO - o I cn

CD

i_3

hO •M

tx,

47
<4-H OCOH

4->

4-3

COc

•rH o

4Jc

0)E

0) o L c

Z3 =S CO <4-l cdCD CO E < o5 • -P o

CO

i-H

O 3 <<-H CO

J- •o 4-> 1) CO 4J Sh 3 3

D, bO E •H o <i-i o Co

Cm

O

O

•o

C -H O 0)

CO -H

•rH Cv-, ScCQ.Eoo

-X3

LOon

CD

S-bO

rH

X <

48
1 X

k\ 1 X \.

CM

X \ X ] O

IIN

X

x

y

X yS

X

jr

X / XX \ X

X 1 o

00o

CD

COo

CO *-- -m =r c \

•-I o o on c \ E on

0) ^-< L3 -

CO ^

ca =r

0) \ E o

s.o \

<4-, JDm

Cm ^o

CO =T

CO \ o o

on o - • Cm -=r C\J

LTi X>

O

on ii . /-^ to • M \ JJ Cm O CO o \ Cm

CO .o

O

c0

—

c C II o O -H

CO -^ -i-5 COd

•rHN

O

t_ » c

co >> 3

CX Cm EO 4-2 CO

cj> cO -=t

s~^ C\J CO o 1 m

d)

J-3

bO

fc

CO CD <sf CM

O CM ^

CD CO o o O o O O1 o1 O1 o1 T

49
CO ^c \ O o~> ac \

CD .O S m

0) ^

CO '-v

cd ^r E o

£- « O \

on o

CD '--

co -=r 03 \ o om

i«o - • <m =r -=r

LTi .O XI CO II

b0=r

•H \ .P

(h

O (S

C •> tSI \ <M

CO X3

O

CO ^ C C II o O -H

CO ^ -P

t- » c

as >, 3 sO J-3 CO I m

CD o G> CO 1^ CD LO ^t CO CM i—

1— o o o O o o o O o

50
Cvj C\l o CO

G>o

0)

s-<D5

00 CO

4-> CD •HO L O. O

<U( 4-> c

1^. TD CD H E O CD CD •H fc. Cw 3

CO - CO - CD ° E W

CD

H

X

•i-H

d o w

CO • 4-> C

rH

O 3 -H

m -Q CO -I-5 CD CO

*"• t, -i-> • o >> 3 ^-^ •O Q. CM

I) EN 4-> O O 3 O

Q.EHW

<vt o «\ O C CO d o -^

<t-i -H O CO II c C 0) ^-v O E SI CO .i-l » CO

•H T3

X

CO 0) • c> Q. CD "O E t, CD O pC CO O 4-> 3 ^

CO CMO t-1 on

CD

J51
X + S/ T— 1 ^f 1 '+/ 1 // O)

~(D00 + f // o

CD

CD o „ /

CO

4-> in S!

*3, X* + ///

COo co o

Cm 4->

""x... < C

£1 / >

TD CD

cd cd -II / 1^ V

CO *-»

"n" / X c> - to oo - CD \ ° E o / / E- X - X / +/

•H C\|

X

CD

OCON

4-> CO / O CO • H O C

X 11 3 ^ id

CO -P

CD CO ^

SL -P ^O

X

-1

lOO

D T—

cd e o P O LTi D O

Q. 1 E H ^r O CO \ X \ <tf OfC

ision =

(a \ o C CD ^ X + \ O E N

co •-• *

•h -a

x \

S- 1 ^ CO CD . + \

Q, CD "O

X COd E S- CD O £ CO O P 3 \ ^

X.

X + \\ I— Sc + \\\\

C\J o

en

CD

£-bO + \\ .,_

Cn \ o . \

s+ \ \. \ ^\

OQ CD ^

CM cD cD cD c:>

52
A —

/ fh/O r/0\ t —

- (y, z) - { (b/8, c/2) - X

1.0

+

+ —

0.8

+

+ —

+ + H0.6 — —

0.4

' ""X" —x—-x—-x—x- —X-- -X~ _j:

i \ A A A A A A A

A

0.2

I I I

—

0.2 0.4 0.6 0.8 1.0

x/a

Figure 3-7 c) Comparison of computed results to TE 0ll field for

three-dimensional computation. Six measurement points were

used, As a function of x for various y and z.

53
00 CO ^1- CM o CM ^1" CD 00 o o o o O O1 o1 O1 o1 T c

cd sz

co o

.H

CD

CO

O

4^>

i.

m CO - C a 3 H O

.O

CO

CD x:pc 3 o

CO

CD -O

t- CDX

0)

4-> 3 CD E CO O

«M

x

Cmoco

CO •H

COaEOo

COon

CD

5-

Cx.

54
o

00

o

o

CD

o

lo -Q

O

CO

o

CM

CD

oqcD^cMo^^-cpoop o o o o cJOOOi-1

I

I

I I I <

c

0)

•2-

CO a

CO <

M >>

S_

<•> CO

CM "O ~ c w

o

4-5

CO rH

co

D

O

.Q

CD

coa

X

CM

-a

CDa l

E cd oOK

Cm X

<4-l oc

o

CO •H

JCOexEoo

CO

I m

CD

bO

55
>>x

o

o

COc

o

•l-l o

CO a:

CO

I CO

bO

c

05co

CO •H

rd • cx^-. E oo O \

I on

J-s

00 CD ^t CM O CM ^t CD 00

o O o O O

1 o

1 O

1 o1

56
O) „ o

CM

cd

ii 00 >> o Xo

1^. SI o o

COco

4-> CD o O

CO <

LO o • 00 d N

fad

•l-l "* c

•H o CO

COco

CO

•l-l CO 03 • o

b)

Comp b/2)

CM o o i m

J-a

cist
NBS-114A (REV. 2-8C) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBS/TN-1303

2. Performing Organ. Report No. 3. Publication Date

October 1986

4. TITLE AND SUBTITLE

A LATTICE APPROACH TO VOLUMES IRRADIATED BY UNKNOWN SOURCES

5. AUTHOR(S)

J. Randa and M. Kanda

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

9. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State, ZIP)

10. SUPPLEMENTARY NOTES

"2 Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here) We suggest an approach to the characterization of electromagnetic

environments irradiated by unknown sources. The approach is based on the

numerical solution of Maxwell's equations subject to the constraints imposed

by the measured values of the field at a small number of measurement points

and by boundary conditions. A thorough examination of two methods for the

numerical solution is presented. The examples attempted demonstrate the

approach but reveal that neither technique is fully successful. Possible

future directions are suggested.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

electromagnetic environment characterization; electromagnetic environment

effects; Hamilton's action principle; ill-posed problems; numerical methods;

successive over-relaxation methods

13. AVAILABILITY

[V] Unlimited

2 For Official Distribution. Do Not Release to NTIS

[X~] Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402. 2 Order From National Technical Information Service (NTIS), Springfield, VA. 22161

14. NO. OF

PRINTED PAGES

64

15. Price

>U.S. GOVERNMENT PRINTING OFFICE: 19 86-0-77 3-597/65 18 USCOMM-DC 6043-P80


NBSTechnical Publications

Periodical

Journal of Research—The Journal of Research of the National Bureau of Standards reports NBS research and development in those disciplines of the physical and engineering sciences in which the Bureau is active. These include physics, chemistry, engineering, mathematics, and computer sciences. Papers cover a broad

range of subjects, with major emphasis on measurement methodology and the basic technology underlying

standardization. Also included from time to time are survey articles on topics closely related to the Bureau's

technical and scientific programs. Issued six times a year.

Nonperiodicals

Monographs—Major contributions to the technical literature on various subjects related to the Bureau's scientific and technical activities. Handbooks—Recommended codes of engineering and industrial practice (including safety codes) developed in cooperation with interested industries, professional organizations, and regulatory bodies.

Special Publications—Include proceedings of conferences sponsored by NBS, NBS annual reports, and other

special publications appropriate to this grouping such as wall charts, pocket cards, and bibliographies.

Applied Mathematics Series—Mathematical tables, manuals, and studies of special interest to physicists,

engineers, chemists, biologists, mathematicians, computer programmers, and others engaged in scientific and

technical work.

National Standard Reference Data Series—Provides quantitative data on the physical and chemical properties

of materials, compiled from the world's literature and critically evaluated. Developed under a worldwide pro- gram coordinated by NBS under the authority of the National Standard Data Act (Public Law 90-396).

NOTE: The Journal of Physical and Chemical Reference Data (JPCRD) is published quarterly for NBS by

the American Chemical Society (ACS) and the American Institute of Physics (AIP). Subscriptions, reprints,

and supplements are available from ACS, 1155 Sixteenth St., NW, Washington, DC 20056.

Building Science Series—Disseminates technical information developed at the Bureau on building materials,

components, systems, and whole structures. The series presents research results, test methods, and perfor- mance criteria related to the structural and environmental functions and the durability and safety

characteristics of building elements and systems.

Technical Notes—Studies or reports which are complete in themselves but restrictive in them treatment of a

subject. Analogous to monographs but not so comprehensive in scope or definitive in treatment of the subject

area. Often serve as a vehicle for final reports of work performed at NBS under the sponsorship of other government agencies.

Voluntary Product Standards—Developed under procedures published by the Department of Commerce in Part 10, Title 15, of the Code of Federal Regulations. The standards establish nationally recognized re- quirements for products, and provide all concerned interests with a basis for common understanding of the

characteristics of the products. NBS administers this program as a supplement to the activities of the private sector standardizing organizations.

Consumer Information Series—Practical information, based on NBS research and experience, covering areas of interest to the consumer. Easily understandable language and illustrations provide useful background

knowledge for shopping in today's technological marketplace.

Order the above NBS publications from: Superintendent of Documents, Government Printing Office,

Washington, DC 20402. Order the following NBS publications—FIPS and NBSIR 's—from the National Technical Information Service, Springfield, VA 22161.

Federal Information Processing Standards Publications (FIPS PUB)—Publications in this series collectively constitute the Federal Information Processing Standards Register. The Register serves as the official source of

information in the Federal Government regarding standards issued by NBS pursuant to the Federal Property

and Administrative Services Act of 1949 as amended, Public Law 89-306 (79 Stat. 1127), and as implemented

by Executive Order 1 1717 (38 FR 12315, dated May 1 1, 1973) and Part 6 of Title 15 CFR (Code of Federal

Regulations).

NBS Interagency Reports (NBSIR)—A special series of interim or final reports on work performed by NBS

for outside sponsors (both government and non-government). In general, initial distribution is handled by the sponsor; public distribution is by the National Technical Information Service, Springfield, VA 22161, in paper copy or microfiche form.
U.S. Department of Commerce

National Bureau of Standards

Gaithersburg, MD 20899

Official Business

Penalty for Private Use $300