EP3008619A2 - Computer simulation of electromagnetic fields - Google Patents

Abstract

A method and system is provided for solving for electromagnetic fields by approximating an electromagnetic function as a sum of basis functions multiplied by coefficients to be determined. The set of equations used to determine the coefficients results from taking derivatives of the action integral with respect to the coefficients (and/or other parameters) and setting the derivative equal to zero, thereby extremizing the action integral.

Description

COMPUTER SIMULATION OF ELECTROMAGNETIC FIELDS

CROSS-REFERENCE TO RELATED APPLICATIONS

[1001] This application is a continuation-in-part of United States Patent Application Number 14/261,415, filed on April 24, 2014, by Terje Graham Void, entitled "Computer Simulation of Electromagnetic Fields;" which is a continuation-in-part of United States Patent Application Number 13/923,305 (Docket number AN-4) entitled "Computer Simulation of Electromagnetic Fields" filed June 20, 2013, by Terje Graham Void; which is a continuation-in-part of United States Patent Application Number 13/914,507 (Docket number AN-3) entitled "Computer Simulation of Electromagnetic Fields," filed June 10, 2013, by Terje Graham Void; which is a continuation-in-part of United States Patent Application Number 13/136,010 (Docket number AN-2), entitled "Computer simulation of electromagnetic fields," filed July 19, 2011, by Terje Graham Void; which claims priority to United States Provisional Patent Application Number 61/365,366 (Docket number AN-1), entitled "Computer Simulation of Electromagnetic Fields," filed July 19, 2010, by Terje Void, and all of the above applications are each incorporated herein by reference, in their entirety.

FIELD

[1002]This specification generally relates to computations of electromagnetic fields.

BACKGROUND

[1003]The subject matter discussed in the background section should not be assumed to be prior art merely as a result of its mention in the background section. Similarly, a problem mentioned in the background section or associated with the subject matter of the background section should not be assumed to have been previously recognized in the prior art. The subject matter in the background section merely represents different approaches, which in and of themselves may also be inventions. There are a variety of numerical techniques for computing numerical fields. Boundary element methods are often used in electromagnetic simulations. The general idea of boundary element methods is that space is modeled as consisting of regions with uniform electromagnetic response properties separated by negligibly thin boundaries. For each uniform region, a parameterized solution to Maxwell's equations is written (for example, a sum of parameters times plane waves, or a sum of parameters times fields due to fictitious "equivalent sources"). The parameterization and the parameters are chosen so that the error is small in some sense.

[1004] However, the best method for minimizing the errors is not at all obvious. A common method is the "method of moments." However, the method of moments is rather poorly motivated and ad-hoc. Different choices of parameters are chosen for different notions of "error" and "small" for different situations. In other words, the choice of parameterization and parameters is an art form that depends on the user's intuition.

DESCRIPTION OF THE FIGURES

[1005] In the following drawings like reference numbers are used to refer to like elements.

Although the following figures depict various examples of the invention, the invention is not limited to the examples depicted in the figures.

[1006]FIG. 1 shows a block diagram of a machine that may be used for carrying out the computations for determining electromagnetic fields.

[1007]FIG. 2 is a block diagram of an example of system, which may store the machine code for solving for electromagnetic fields.

[1008]FIG. 3 is a flowchart of an embodiment of a method that is implemented by processor system.

[1009]FIG. 4 is a flowchart of an embodiment of a method for solving the electromagnetic equations.

[1010]FIG. 5 shows an example of homogeneous regions for computing the electromagnetic fields.

[1011]FIG. 6 shows an example of computational regions corresponding to the homogeneous regions of example with a thin boundary region.

[1012JFIG. 7 shows is a diagram illustrating elements of an example to which the least action method is applied.

[1013]FIG. 8 shows another sample problem having a dielectric sphere and a point dipole. [1014JFIG. 9 shows an example of a placement of virtual sources inside the sphere of FIG. 8 for computing the basis potential functions outside of the sphere.

[1015]FIG. 10 shows an example of a placement of virtual sources outside the sphere of FIG. 8 for computing basis function within the sphere.

[1016]FIG. 11 shows an example of an arrangement of virtual point dipoles on the polyhedron of FIG. 9.

[1017]FIG. 12 shows an example of using a mesh of simplices for computing a field. DETAILED DESCRIPTION

[1019] Although various embodiments of the invention may have been motivated by various deficiencies with the prior art, which may be discussed or alluded to in one or more places in the specification, the embodiments of the invention do not necessarily address any of these deficiencies. In other words, different embodiments of the invention may address different deficiencies that may be discussed in the specification. Some embodiments may only partially address some deficiencies or just one deficiency that may be discussed in the specification, and some embodiments may not address any of these deficiencies.

[1020] In general, at the beginning of the discussion of each of FIGs. 1, 2, 5-11 is a brief description of each element, which may have no more than the name of each of the elements in the one of FIGs. 1, 2, 5-11 that is being discussed. After the brief description of each element, each element is further discussed in numerical order. In general, each of FIGs. 1-11 is discussed in numerical order and the elements within FIGs. 1-11 are also usually discussed in numerical order to facilitate easily locating the discussion of a particular element. Nonetheless, there is no one location where all of the information of any element of FIGs. 1-11 is necessarily located. Unique information about any particular element or any other aspect of any of FIGs. 1-11 may be found in, or implied by, any part of the specification.

SYSTEM

[1021] FIG. 1 shows a block diagram of a machine 100 that may be used for carrying out the computations for determining electromagnetic fields. Machine 100 may include output system 102, input system 104, memory system 106, processor system 108, communications system 112, and input/output device 114. In other embodiments, machine 100 may include additional components and/or may not include all of the components listed above. [1022] Machine 100 is an example of what may be used for carrying out the computations for determining electromagnetic fields. Machine 100 may be a computer and/or a special purpose machine computing electromagnetic computations.

[1023] Output system 102 may include any one of, some of, any combination of, or all of a monitor system, a handheld display system, a printer system, a speaker system, a connection or interface system to a sound system, an interface system to peripheral devices and/or a connection and/or interface system to a computer system, intranet, and/or internet, for example. Output system 102 may be used to display indication the types of input needed, whether there are any issues with the inputs received, and output results of a computation of an

electromagnetic field, for example.

[1024] Input system 104 may include any one of, some of, any combination of, or all of a keyboard system, a mouse system, a track ball system, a track pad system, buttons on a handheld system, a scanner system, a microphone system, a connection to a sound system, and/or a connection and/or interface system to a computer system, intranet, and/or internet (e.g., IrDA, USB), for example. Input system 104 may be used to enter the problem parameters, such as the grid on which to perform the computations, the material parameters of each region within the grid, choose basis functions, and initial state, if applicable.

[1025] Memory system 106 may include, for example, any one of, some of, any combination of, or all of a long term storage system, such as a hard drive; a short term storage system, such as random access memory; a removable storage system, such as a floppy drive or a removable drive; and/or flash memory. Memory system 106 may include one or more machine-readable mediums that may store a variety of different types of information. The term machine-readable medium is used to refer to any non-transient medium capable carrying information that is readable by a machine. One example of a machine -readable medium is a computer-readable medium. Memory system 106 may store the machine instructions (e.g., a computer program, which may be referred to as a code) that cause machine 100 to compute electromagnetic field computations. Memory 106 may also store material parameters corresponding to specific materials and basis-functions. The basis function may be parameterized solutions to the electromagnetic equations within a homogenous region.

[1026] Processor system 108 may include any one of, some of, any combination of, or all of multiple parallel processors, a single processor, a system of processors having one or more central processors and/or one or more specialized processors dedicated to specific tasks. Processor system 110 may carry out the machine instructions stored in memory system 108, and compute the electromagnetic field.

[1027] Communications system 112 communicatively links output system 102, input system 104, memory system 106, processor system 108, and/or input/output system 114 to each other. Communications system 112 may include any one of, some of, any combination of, or all of electrical cables, fiber optic cables, and/or means of sending signals through air or water (e.g. wireless communications), or the like. Some examples of means of sending signals through air and/or water include systems for transmitting electromagnetic waves such as infrared and/or radio waves and/or systems for sending sound waves.

[1028] Input/output system 114 may include devices that have the dual function as input and output devices. For example, input/output system 114 may include one or more touch sensitive screens, which display an image and therefore are an output device and accept input when the screens are pressed by a finger or stylus, for example. The touch sensitive screens may be sensitive to heat and/or pressure. One or more of the input/output devices may be sensitive to a voltage or current produced by a stylus, for example. Input/output system 114 is optional, and may be used in addition to or in place of output system 102 and/or input device 104.

[1029] FIG. 2 is a block diagram of an example of system 200, which may store the machine code for solving for electromagnetic fields. System 200 may include action extremizer 201 and user interface 202 having tolerance choices 204, grid choices 206, output choices 208, basis function choices 210, and material parameter choices 212. Memory 106 may also include basis functions 214, initial state generator 218, and output generator 220. In other embodiments, system 200 may include additional components and/or may not include all of the components listed above.

[1030] In an embodiment, system 200 is an embodiment of memory system 106, and the blocks of system 200 represent different functions of the code that solves for electromagnetic fields performs, which may be different modules, units and/or portions of the code. Alternatively, system 200 may be a portion of the processor system, such as an Application Specific Integrated Circuit (ASIC), and/or another piece of hardware in which the code of system 200 is hardwired into system 100.

NOTATION [1031] The notation of Doran and Lasenby in Geometric Algebra for Physicists (Cambridge University Press, 2003) in SI units is used in the current specification. For reference, some details of the notation of Doran and Lasenby are given in the notation section.

[1032] A frame vector of 4D space-time is represented by γ_α for a E {0,1,2,3}, with γ₀ equal to the time-like frame vector and with γ₀γ₀ = 1, and γ γ =— 1 for i E {1,2,3}, and γ_αγ_β = 0 for α≠ β. A reciprocal frame vector of 4D space-time is represented by γ^α for a E {0,1,2,3}, with Y^aY_a ⁼ +1 for a £ {0,1,2,3} and γ^αΥ_β = 0 for α≠ β. A frame vector of the 3D space of the inertial frame associated with γ₀ is represented by ¾≡ γ^Υο for i E {1,2,3}. The

electromagnetic field F is defined to equal the antisymmetric part of the space-time vector derivative of the space-time vector potential A,

F = V Λ A, where the space-time vector derivative is related to the time and vector derivative in 3D space by the space-time potential A is related to the scalar potential Φ and vector potential A in 3D space

A (ε₀Φ + A ) Yo

μ₀

Λ is the outer product of geometric algebra, γ₀ is the time-like space-time velocity vector defining the inertial frame associated with the 3D space vector potential A and scalar potential Φ, and quantities with an over-script arrow such as A are always vectors in the 3D space of the inertial frame defined by γ₀. In this specification, the word "potential" is generic to the space- time vector potential, to the scalar and vector potential in 3D space together or separately, to the "media potential" and/or to other simply related quantities including a discontinuous potential that is related to a continuous potential by a gauge transformation. The electromagnetic field F is related to the electric field E and "magnetic B -field" or "magnetic intensity" B by

1

F = ε₀Ε + III B where

III≡ ¾¾¾ = YoYiY₂Y3 is the pseudoscalar of the geometric algebra of both 3D space and 4D space-time. The 3D cross product is related to the wedge product via the pseudoscalar by III a x b = a Λ b. In general, the field G is a function of F and a polarization field P, or G=f(F,P). Although the methods described herein are not limited to linearly polarizable material, for linearly polarizable media the field G is the sum of F and a polarization field P,

G =F+P, where the 4D space-time polarization field P is composed of the 3D space vector electric polarization P and magnetic polarization vector M according to

P = P - III - M

c so that G may be written as G = D + III

c Ή, where the electric displacement field D is given by D= ε₀Ε +P, and the "magnetic H-field" or "magnetic field" H is given by ϊϊ=— B-M.

μο

[1033] In this specification, the word "field" is generic to F or G, or their parts E, B, D, or H, or the simply related media field F_m . [1034] The polarizations P and M are generally modeled as functions of E and B, often linear functions, but not necessarily. The space-time vector field J representing charge-current density in 4D space-time volumes is related to the charge density p and 3D vector current density J by the equation and the space-time vector field K representing charge-current density on 3D boundaries in 4D space-time is related by the equation

Κ = ( σ + ^ κ) γ₀

[1035] The convention of Mathematica from Wolfram Research that expressions are grouped by rounded brackets (parentheses) while function arguments are enclosed by square brackets [like these] is used. Function arguments may be implied or written explicitly, so that, for example, the action may be written as either S or S[A, VA].

[1036] Extremizer 201 extremizes the action integral by adjusting coefficients of a sum of basis. After the coefficients have been adjusted by the extremizer, the sum of basis functions is a solution or an approximation of a solution to the electromagnetic fields. User interface 202 may be used for receiving input for setting up the problem to solve and outputting results of computing electromagnetic fields. Tolerance choice 204 is an input for entering the tolerance for performing the computations or another form of input for determining the accuracy with which to perform the computations.

[1037] Grid choices 206 includes one or more input fields for defining the space in which the problem is solved, including a description of boundaries between regions with different polarizabilities. Grid choices 206 may also include input fields for defining the initial state of the system. The grid choices 206 divide the computation region into two or more smaller computation regions or cells.

[1038] In an embodiment, a grid may be chosen to divide space or space-time into computation regions that are simplices, or to divide the boundary between regions into simplices. In the current specification, a simplex (which in plural form is simplexes or simplices) refers to a generalization of the notion of a triangle or tetrahedron to arbitrary dimension. Specifically, an n-simplex is an n-dimensional polytope, which is the convex hull of n + 1 vertices. For example, a 2-simplex is a triangle, a 3-simplex is a tetrahedron, and a 4-simplex is a pentachoron.

Simplices may be desirable when the parameterized potential is defined to be a linear function of position within each simplex. Using simplices, a continuous potential may be parameterized by the values of the continuous potential at all simplex vertices, and the field F is a constant within each simplex because the field F is defined to be a linear combination of derivatives of A. Also, the field G is a constant within each simplex if the field G is a function of only the field F, whether G is a linear or nonlinear function of F. Simplices may be defined that cover all of the simulation space, or may be combined with regions of other or arbitrary shape. Including simplices with or between regions of arbitrary shape may make it easier to define good potential basis functions.

[1039] An example of a situation where it may be convenient to have some regions that are simplices and some that are not is to cover a region containing a nonlinear material with simplices, and to use a single space-time region for a homogenous linearly polarizable medium outside the nonlinear material. In this example, the boundary between the nonlinear and linearly polarizable media is represented by lower dimensional simplices.

[1040] Basis functions 210 may include a set of default basis functions that the user can choose, such as those specified by boundary values and calculated with the help of Green's theorem as an integral of the boundaries of each simplex. Basis functions 210 may also include an input charge-current source or associated basis field. Basis function choices 210 may include a list of default basis functions, which the user may choose. In an embodiment, the default basis functions, when multiplied by a series of quantities (which may be referred to as parameters) and summed together, may form parameterized functions, and solving for the parameters solves the electromagnetic equations. The basis functions 210 may include potential basis functions and field basis functions. In this specification, the term basis function is generic to potential basis functions, field basis functions, and other basis functions.

[1041] In an embodiment, each potential basis function may be specified by potential values on boundaries along with the requirement that at non-boundary points the Euler-Lagrange equation is satisfied by the potential for any set of parameter values.

[1042] In another embodiment, each potential basis function may be specified by values of the scalar potential and tangential values of the vector potential on boundaries and the requirement that the potential satisfies the Lorenz condition at region points next to a boundary, along with the requirement that at non-boundary points the potential satisfies the wave equation in the media for any set of parameter values. [1043] Basis functions may be defined in terms of fields that would result from fictitious "virtual" charge and current densities, such as point dipole moments, in a corresponding fictitious space filled with an unbounded uniform polarizable medium. In one embodiment, one or more virtual point dipole moments at each of a series of points are used to define each basis function within each region with uniform polarizability. The positions, directions, and type (electric and/or magnetic) of the virtual dipole moments may be chosen to form a good set of basis functions. A good set of basis functions is a set basis functions capable of representing any function that may represent the electromagnetic fields in the non-boundary regions. In an embodiment a good set of basis functions is a complete set of basis functions. In an

embodiment, the virtual dipole moments and/or other virtual charge and current densities are located outside of a computation cell and/or region in which the basis function is being defined.

[1044] Material parameters 212 may include the material parameters that the user can choose for each region of the space in which the problem is being solved, for example. In an embodiment, material parameters may include combinations of material parameters associated with common materials, so that the user can just specify the name of the material and the combination of parameters associated with that material may be assigned to the region chosen.

[1045] Basis function 214 may include an algorithm for choosing basis functions when the user has not chosen and/or may include an algorithm for converting the user's input into a function that is called by the code as the basis functions that are used by least action extremizer 201 to extremize the action integral and solve for the fields.

An Embodiment

[1046] Action extremizer 201 is a module that solves for the electromagnetic fields by extremizing the action. In general, action extremizer 201 computes numerical values for a set of parameters of basis functions (e.g., the vector and scalar potentials from which the electric and magnetic fields may be derived). In an embodiment, the potential functions are a parameterized continuous electromagnetic potential. The values of the parameters are computed within a predetermined tolerance to extremize an action integral. The action integral is the integral of the Lagrangian density. The Lagrangian density contains terms for the electromagnetic field in vacuum or in polarizable media and for the interaction of the electromagnetic field with electric charge and current. The terms of the electromagnetic fields may describe linear or nonlinear polarizable media. [1047] A method is described for accurately and efficiently computing values of electromagnetic fields, suitable for use in various fields of applied physics and engineering ranging from magnetostatics to RF engineering to integrated optics. An electromagnetic potential is defined that is parameterized by the electromagnetic potential's tangential components at boundaries between regions with uniform polarizability, and parameter values are chosen that extremize the classical electromagnetic action. Useful general results are derived using the geometric algebra of space-time. The method is compared with the widely used method of moments.

The Principle of Least Action

[1048] The current section contains a summary of the principle of least action as applied to electromagnetism in space-time, using the notation of and results from [1 ,2], and SI units. For simplicity the problem is restricted to linearly polarizable media. The electromagnetic action S with a space-time vector charge-current source density J is a scalar- valued functional S [A, VA] of a continuous space-time vector field A and derivatives VA of A, over all of space-time volume u, given by the integral where the electromagnetic field F is a space-time bivector field defined as a derivative of the space-time potential A, a space-time vector field, by

F≡ VA A (2)

[1049] The total source J≡ J_f + J_b is the sum of a "free" source density J_f and a "bound" source density J_b = V · P where P = P— IM is the electric vector plus magnetic bivector

polarization. For simplicity, P is assumed to be a linear function of F. However, taking P to be a linear function of F is not necessary and G may be a nonlinear function of F in this case. A · J can be integrated by parts to rewrite the action as

S = J_u (i _F - G - A - J_f) du (3) where the displacement-magnetic field G≡ F + P is a linear function of F since P is. In the simple case of isotropic polarizability, the equation for the displacement-magnetic field G can be written as °^{[Ρ] ≡} (έ) ^Ρ · ^γ"^Υ" ⁺ © ^ΡΛγ"^Υ" (4) where the permittivity ε and permeability μ are scalar functions of position in the rest frame of the medium, associated with the unit time-like vector γ₀. The principle of least action says that variation 5S of the action S with respect to any variation δΑ of the field A must equal zero for any physically realizable field A:

[1050] After writing F in terms of A and using linear dependence of G on F, the first term of equation (5) can be integrated by parts so that equation (5) can be written as

J_u δΑ (V G - J_f) du = 0 (6)

[1051] Equation (6) is true if and only if the Euler-Lagrange equation,

V G - J_f = 0 (7) holds throughout space-time u.

[1052] F and therefore G do not depend on V · A, resulting in a gauge freedom. The gauge freedom may be removed by requiring that A satisfies the Lorenz condition,

V · A = 0 (8)

Maxwell's Equations

[1053] The following section contains a summary of the expression of the Euler-Lagrange equations using vectors and scalars in 3D space.

[1054] The space-time bivector field F may be projected into two parts, one part containing and one part not containing a unit time-like vector γ₀ associated with any particular inertial reference frame, as

F = F^■ YoYo + FAYoYo (9)

[1055] The first term of equation (9) is a 3D space vector proportional to the electric field E and the second is a 3D space bivector proportional to the magnetic intensity B. Using SI units, the first term of equation (9) is ^{Ρ = ε}ο^{6 +}Θ(έ)^{ιδ (10)} where ε₀ and μ₀ are constants, c = l/-j ₀ μ₀ and I is the unit pseudoscalar. Following the convention [1,2], an arrow over the letter is used to represent only 3D vectors.

[1056] Similarly, the space-time bivector field G may be expressed as the sum of parts projected onto and rejected by y₀,

G = G · YoYo + GAYoYo (11) the two terms of equation (11) give the displacement D≡ ε₀Ε + P = eEand magnetic field

1

H≡— B— M = - B according to

μο μ

G = D + Ιίϊ (12)

[1057] The scalar and space-time bivector parts of the product of the Euler-Lagrange equation with Yogive the first and second of Maxwell's equations,

V D = p_f (13) V X H -≤ = C (14)

[1058] The space-time bivector and pseudoscalar parts of the product of the mathematical identity V Λ V Λ A = 0 with γ₀ give the third and fourth of Maxwell's equations,

V x E + ^ = 0 (15)

V B = 0 (16)

Discontinuities in Polarizability

[1059] Discontinuities in polarizability of the medium with change in position in space correspond to discontinuities in permittivity ε or permeability μ. Although the space-time potential A is continuous everywhere by assumption, discontinuity in polarizability generally corresponds to discontinuities in derivatives of A and therefore in F and in G. If space-time is filled with bounded regions of 4D volume having continuous polarizability with discontinuities only at the boundaries between the bounded regions, then such discontinuities in F and G occur only at the boundary points. Let the collection of boundary points be called β and the collection of region points be called u and let the corresponding surface charge-current density and spatial charge-current density be called K_f and J_f , respectively. Then the equation from which the Euler-Lagrange equation was deduced can, in the instant case, be written as

J_u δΑ (V G - J_f)du + Jp 5A (n AG - K_f)dp = 0 (17) for any variation δΑ, where n is a unit normal space-time vector at the boundary, AG is the change in G across such a boundary in the direction n, and du and dp are scalar-valued integration differentials. The integration differential du is an infinitesimal volume in u and the integration differential dp is an infinitesimal area on the boundary surface P separating regions of space-time volume u with continuous polarizability.

The Method of Least Action

[1060] In the general theoretical framework above, any variation δΑ in the field A is considered. But in numerical computations, typically the field A is parameterized with a finite number of parameters c_a that may range over all possible values or may be constrained to range over some specified set of values

[1061] . For simplicity, consider A to be a linear function of parameters c_a, A = a_ac_a (18) where the a_a are basis fields for the field A, the subscripted index ranges over a set of integers, and use the convention of summing over repeated dummy indices. Each basis field a_a is defined and continuous at every point of space-time, including boundary points. Some parameters may be known or specified, corresponding to "known" or "incident" fields, and others unknown; for simplicity, let the parameters c_a for a≠ 0 be unknown and that for a = 0 be known, and choose c₀ = 1. The convention that a Greek index such as a ranges over all values while a Roman index such as i ranges over i≠ 0 is used herein.

[1062] Substituting A = a_ac_a, then extremize the action with respect to variations of A that are described by variations 5c_j of the parameters Cj ⁼ 3j 5Cj . The result is

true for any set of values of the variations 5Cj including 5c_j = 5j _j where 5j _j is the Dirac delta function, for which equation (17) above becomes the set of equations indexed by i, J_u a (V - G - J_f)du + J_p a (n - AG - K_f)dp = 0 (19)

Equation (19) is a general result that may be applied by the method of least action to compute electromagnetic fields. Specifically, equation (19) applies to both linear and nonlinear media, although for simplicity the discussion is continued below for linear media.

[1063] Since F is a linear function of A and G is a linear function of F, F and G can each be written as corresponding sums of basis functions,

F = f_ac_a (20) where f_a = V Λ a_a, and

G = g„c_a (21) where g_a = -f_a - γ₀γ₀ + ^f«^AYoYo-

[1064] Consequently, the integral function above, the set of equations (19), can be written in terms of the basis functions f_a, and g_a, as

J_u *i · (V · g_ac_a - J_f)du + J_p a,^■ (n · Ag_ac_a - K_f)dp = 0 (22)

[1065] It is convenient to put all terms containing unknown parameters (a = j≠ 0) on the left, and all others (including a = 0 terms) on the right; the result can be written explicitly as a matrix equation, where m_M≡ J_{u 3i} · (V · _gj)du + Jp _ai · (n · A_gj)dp (24) and b|≡ - /_u ¾^■ (V^■ go - J_f)du - Jp ¾^■ (n^■ Ago - K_f)dp (25)

[1066] The preceding system of equations, equations (23), (24) and (25), can be solved for the parameters Cj , and then using the value iti_j j , and b_i? the value of any field at any point can be calculated. [1067] The solving of equations (23), (24) and (25) is especially useful when combined with the use of "media potentials" (which are potential functions that are dependent on the media parameters) in polarizable media and the use of potentials that are discontinuous across the boundaries between regions containing material with different polarizabilities, as will be described in later sections, and when the basis functions are chosen so that the volume integrals are zero for any set of parameter values. Despite the simplicity of solving equations (23), (24) and (25), solving equations (23), (24) and (25) has not been previously used in computational electromagnetism.

[1068] Since solving equation (19) or equations (23), (24) and (25) to compute electromagnetic fields is based on the principle of least action, in this specification the solution method is referred to as the method of least action.

Vector Algebra Expression

[1069] Each space-time vector quantity w is related to a relative scalar w₀ and relative 3D space vector w in a frame associated with the unit time-like vector γ₀ by w = WYOYO = (w · Yo + w Λ γ₀)γ₀ = (w₀ + w)y₀

For the space-time vector position r, boundary normal n, basis potential a_i? free charge-current volume density J_f , free charge-current boundary density K_f , and derivative V, the corresponding scalar and 3D vector quantities are related by r = (ct + f )Yo n = (n · yo + n)y₀

[1070] The expressions for r, n, a_t J_f, K_f, and V can be used to rewrite the matrix and array elements above in terms of scalar and 3D vector quantities. In simple cases such that all media at non-boundary volume points V and boundary points s in 3D space are stationary in one inertial reference frame, n · γ₀ = 0 at all boundary points and these elements may be rewritten as i_,j≡ j h_j - a_td_j)) dv dt - ai · (n x hj)) ds dt

(27) and i≡- j (φ_i(V · do - p_f) - ^ί_i · (V x ho - ^_tdo - J_f)) dV dt

^~ S_t - · (" ^x ho - K_f)) ds dt

(28)

Harmonic Time Dependence

[1071] If all medium boundaries are stationary in one inertial frame associated with γ₀, the time coordinate can be defined as t = r · γ₀ and the 3D space position vector as f = r Λ γ₀ at any point r of space-time; if in addition all time dependence is harmonic, such as for one Fourier component after applying a Fourier transform, the complex vector representation of fields with an implied factor of e^_lwt may be applied. In the instant case the general expressions for the matrix and array elements for frequency ω are identical to those of (26) and (27), but with d_t replaced by - io) and integration done over only the volume V and boundary s, and not time t.

[1072] The formulation of the general expression for the matrix and array elements allows simulation of electrical conductivity by using complex-valued permittivity ε. Since such a

1

permittivity has frequency dependence of— for small frequencies ω, multiplying every row in the matrix equation (22) containing a factor of such a permittivity by the frequency ω and canceling these factors symbolically before numerical calculation of the matrix elements results in a set of equations that can then be solved without any numerical problems for any frequency, including ω = 0. Other diverging factors such as (<)— (<) () , representing a physical resonance with numerical problems at frequency ω₀, can be similarly eliminated. Parameterized Potential Satisfying the Euler-Lagrange Equation

[1073] It is convenient to choose A = a_ac_a to satisfy the Euler-Lagrange equation

V · G = J_f where G = g_ac_a with g_a = g[V Λ a_a], for any set of unknown parameters Cj at all non-boundary points. The equation V · G = J_f (where G = g_ac_a with g_a = g[V Λ a_a]) will be satisfited for any set of unknown parameters Cj if the basis potential functions satisfy

V - _gi = 0 (29) at all non-boundary points. In the instant case the volume integrals in the matrix elements given by (24) and (25), or equivalently by (27) and (28), are all identically zero.

[1074] Choosing A to satisfy the Euler-Lagrange equation can generally be done if space is modeled as made up of homogeneous polarizable regions. If space is modeled as made up of homogeneous polarizable regions, each potential basis function a_a may be defined by specifying the value of potential basis function a_a at all points on boundaries between homogeneous regions, and by defining the function within each homogenous region to satisfy the

corresponding Euler-Lagrange equation while having specified boundary values. For homogeneously polarizable regions, computationally practical expressions for such basis functions can be found as described later.

The Euler-Lagrange Equation in Polarizable Media and Media Quantities

[1075] The following section on the Euler-Lagrange Equation in Polarizable Media describes a method for easily solving the Euler-Lagrange equation in regions with uniform polarizability, by introducing a "media potential" with components that are linear transformations of the potential components and that satisfy the wave equation. Introducing a "media potential" allows symbolic expression and fast numerical computations of basis functions that satisfy the Euler-Lagrange equation at non-boundary points.

[1076] The fundamental physical space-time vector operator V and vector quantities A and J, when expressed as relative quantities in an inertial frame characterized by y₀, for example, A = Αγ₀γ₀ = (A · γ₀ + A Λ γ₀)γ₀, are related to associated scalar and 3D vector quantities in SI units by

(30) ^A = (^eo1> ⁺ 0(^) ^S) o (31) where the quantity J_frepresents all "free" electric charge source p_f , and electric current source J_f , defined as any source that is not described by derivatives of the electric polarization P or magnetic polarization M of the medium. J_f , is related to the total charge-current source J by

J = J_f + J_b, where (j_b = p_b + (^j J_b γ₀ is the "bound" source, defined as source that is described by the derivatives p_b = —V · P and J_b = V x M— d_tP of the electric polarization P and magnetic polarization M of the medium, or equivalently by J_b≡ V · P where P≡ P—

IM. Note that since J_f is defined as any source that is not described by derivatives of the electric polarization P or magnetic polarization M of the medium, the induced charge p_b and current source J_b in an electrical conductor, described by a complex permittivity, is a "bound" source J_b, not a "free source" J_f . Consequently, in this specification, the charged current source J_f may be referred to as the "specified" source instead of "free" source to avoid confusion.

[1077] Next, the "media" quantities are defined for media at rest in an inertial frame

represented by y₀ by replacing ε₀→ ε and μ₀→ μ everywhere ε₀ and μ₀ would otherwise occur in the equations for free space, including in c = 1/ ^ε₀ μ₀, where ε₀ and μ₀ are the permittivity and permeability of free space, respectively, while ε and μ are the corresponding parameters of a particular media. Thus, the media vector derivative, media potential, and media charge-current are given in terms of components of the corresponding non-media quantities by

V_m≡Yo (( ) d_t + v) (33) A_m = (_£ φ + (¾ (1) Α) γ₀ (34) where v = 1/-J£[X, which is the speed of a light in a media having the material parameters ε and μ , while the speed of light in a vacuum is given by c = l/-j ₀ μ₀. When computing the harmonic time dependence at frequency ω, and using complex amplitudes, d_t is replaced by - io) and k≡ ω/c in V is replaced by k_m = ω/ν to form V_m. [1078] The Euler-Lagrange equation may be rewritten in terms of these quantities as

V_m · G_m = J_m (36) where

G_m≡V_m A A_m = £₀E + (±) (±) IB (37)

[1079] The wedge product of any space-time vector quantity, including V_m, with itself is still zero, and consequently,

V_m A G_m = 0, (38) which may be combined with the rewritten Euler-Lagrange equation as

V_mG_m = J_m (39) [1080] It is useful to require that A_m satisfies the "media Lorenz condition,"

V_m · A_m = 0 (40)

[1081] Although the "media Lorenz condition" is not equivalent to the Lorenz condition - the media Lorenz condition may be written as (^j δ,-φ + V · A = 0 while the Lorenz condition,

V · A = 0, may be written as (^j d_t(|) + V - A = 0 - either are equally valid gauges to which the electromagnetic equations may be constrained to remove the gauge invariance, thereby simplifying solving the electromagnetic equations numerically. Using the media Lorenz condition, the Euler-Lagrange equation simplifies to the wave equation

V_m ²A_m = J_m (41)

[1082] With the expressions of equations (37), (38), (39), (40), and (41), solutions to the Euler- Lagrange equation can be found in linearly polarizable media by transforming any boundary conditions on the potential A to boundary conditions on the media potential A_m and any specified source J_f to a media source J_m, solving the media wave equation for A_m, and then transforming the solution from A_m back to A.

Basis functions defined by boundary values

[1083] A useful method of defining a space-time potential in a bounded region, such as a potential basis function

[1084] 1) identify the space-time vector value a of the space time potential and the space-time vector value of the normal derivative n- Va_a of the space time potential, representing 8 degrees of freedom, on the region side of all boundary points,

[1085] 2) require that, within the region, the potential satisfies the Euler-Lagrange equation,

[1086] 3) require that within the region the potential also satisfies a gauge condition, and

[1087]4) use the gauge condition and Green's Theorem on the side of each of the boundary points that is within the region (in which the current computation is being performed), as constraints to reduce the number of degrees of freedom at each boundary point from 8 to 3. It may be convenient to choose the 3 remaining degrees of freedom to be the components of the potential that are tangential to the boundary in space-time. Specifically, the 4 space-time components of the potential at every boundary point plus the 4 space-time components of the normal derivative of the potential at every boundary point gives 8 degrees of freedom. Next, Greens Theorem is applied, which gives 4 equations (or 4 degrees of constraint) for every boundary point that relates all 8 of the degrees of freedom reducing the number of degrees of freedom to 4. The Lorenz condition provides one more equation (or 1 more degree of constraint), which relates only derivatives and further reduces the number of degrees of freedom to 3. One may choose values for any 3 of the 8 quantities (e.g., based on boundary conditions), and compute the remaining 5 quantities from the 3 chosen quantities.

[1088] It can be convenient to compute the potential from the boundary values using the corresponding medium potential, medium vector derivative, and other corresponding quantities such as the medium gauge condition and medium Euler-Lagrange equation. It may be convenient to choose the 3 degrees of freedom to be the tangential components of the potential function.

[1089] The 3 chosen degrees of freedom (e.g., the tangential components of the potential function) at each boundary point can be given by a finite number of parameters in a model of the potential, from which all 8 degrees of freedom of the potential and its normal derivative, and therefore the fields f_a and g_a, can be quickly calculated at all points on the boundary. The potential and field value at any point in the region may then be calculated from the values on the boundary using Greens Theorem [3].

[10901 Basis functions defined by the values on the boundary are easy to define and calculate if every point of space is either part of a region with uniform polarizability, or on the boundary between such regions. Then within each region, a medium-Lorenz condition may be chosen such that the medium potential satisfies the medium wave equation, for which Green's function and theorem are simple and numerically tractable.

[1091] The boundary between regions may be modeled or approximated as a set of connected triangles, with the potential specified at each vertex by 3 parameters and defined to vary linearly with position at all other points on each triangle.

[1092] Each such basis function may be chosen to have nonzero tangential boundary components on only one localized part of the boundary, such as by choosing the parameters of the basis function to be nonzero at only one vertex. Choosing nonzero tangential boundary components on only one localized part of the boundary results in the integrands of the boundary integrals of (19) being zero everywhere except in the localized part of the boundary having the nonzero tangential boundary components, since the integrand is independent of the

perpendicular component of the potential. The computation time required to compute the value of each integral from previously computed potential and field values is then linear in the size of the problem (e.g., the number V of triangle vertices).

[1093] The application of Greens theorem to define the basis potentials, however, requires the computation of order p matrix elements and inversion of the resulting p by p matrix, where p is the number of vertex points. The matrix mi,j of equations (7,9,11) is also is of order p by p and also must be inverted. The time required to solve electromagnetic problems by the application of Green's theorem method of least action may be dominated by these two computations.

[1094] This method results in good basis functions for a wide variety of geometries, including thick or thin wires and thick or thin sheets, with any electric and magnetic polarizability and any electric conductivity.

Stress Tensor

Conservation Laws

[1095] Noether's theorem depends on the satisfying the Euler-Lagrange equations that follow from Hamilton's principle of least action (briefly, Noether's theorem states that

any differentiable symmetry of the action of a physical system - e.g., a smooth transformation that does not affect the action - has a corresponding conservation law). Solutions using the method of least action that satisfy V^■ G— J_f = 0 at all non-boundary points and

(n^■ AG— K_f)dp = 0 for all basis functions a_i? satisfy a modified Noether's theorem that gives a modified conservation law for the electromagnetic stress-energy tensor T, where matrix elements of T are given by Τ_μν = F γ_μΡ γ_ν. Τ_μν represents the flux through a hyperplane perpendicular to γ_μ of the v component of the energy-momentum density.

Noether's Theorem applied to the invariance of the electromagnetic action with respect to the symmetry of translation in time is V^■ T = 0 at every point of space-time. From V^■ T = 0 it follows that n^■ ΔΤ = 0 at every point of any physical or mathematical surface with normal vector n where ΔΤ is the change in T across the surface. The equation n^■ ΔΤ = 0 represents the conservation of energy-momentum through the surface at the point being considered. The modified Noether's theorem is V -T = 0 at all non-boundary points, and f ^■ (n^■ ΔΤ)άβ for all basis functions a;. V -T = 0 describes the conservation of energy and momentum at all non- boundary points, and J ^^■ (n^■ ΔΤ)άβ describes a potential- weighted boundary integral conservation law for energy and momentum at boundary points. For localized basis functions, as described in the section "An embodiment of the method of basis functions defined by boundary values," J ^^■ (n^■ ΔΤ)άβ ensures that the integral or average of energy momentum over a small patch identified by the localized basis function a_j - but not at every infinitesimal point - is exactly conserved across the boundaries between regions of computation.

Gauge-equivalent Potentials

[1096] Discontinuity of the potential across boundaries appears to violate an initial assumption of the principle of least action, that the potential is a continuous function at all points.

[1097] However, it can be shown that given any potential A' that satisfies the Euler-Lagrange equation and the media Lorenz condition, and is continuous across the boundaries except for the perpendicular component of the potential A', the potential A = A'— VX satisfies the Euler- Lagrange equation, is continuous at all points including boundary points, and results in the same field VAA = VAA^f , where λ is any continuous scalar field such that at the boundaries, λ = 0 and n · (Α'— \7λ) is continuous, where n is the boundary normal. The resulting potential A generally cannot satisfy the Lorenz condition 0 = V · A or media Lorenz condition, since doing so requires that λ satisfies a 2^nd order differential equation, V · A' = \7²λ for the Lorenz condition, and the necessary boundary conditions over-determine a solution to such an equation. A and A' have the same tangential boundary values, result in the same field F, and differ by only a gauge transformation; we may say such potentials are "gauge-equivalent".

An embodiment of the method of basis functions defined by boundary values

[1098] The current section discusses the calculation of potentials and fields at boundary points. Specifically, the current section describes a method for calculating all components of the potential, and calculating the normal directional derivative of the potential at all boundary points enclosing a region with uniform polarizability, given the boundary conditions of only the tangential components of the potential. The electromagnetic field at all boundary points can then be calculated from the components of the potential and the normal directional derivative of the potential, which allows evaluation of the boundary integrals needed to extremize the action.

[1099] The first step is to specify the value of each basis potential a_a and the normal derivative of the basis potential a_a at boundary points by a) requiring that each basis potential a_a satisfies the corresponding homogenous or inhomogeneous Euler-Lagrange equation and the media Lorenz condition, and therefore the media wave equation, at all points enclosed by the given boundary, b) giving the values of each basis potential a_a tangential components at all boundary points, and c) using the media Lorenz condition and Greens Theorem for the media wave equation as constraints to calculate the basis potential normal component and all components of the normal derivative of the basis potential at boundary points. For a choice of basis potentials, that meet requirements (a), (b), and (c), the volume integrals in the equation expressing extremization of the action are all identically zero, and the coefficients C_j can be varied in the total potential A = a_ac_a to extremize the boundary integrals of the action.

[1100] To evaluate the boundary integrals needed to extremize the action, the values of each electromagnetic field basis function g_a must be found on each side of every boundary at every boundary point. The method described in the current section enables the efficient calculation of the normal component of a space-time vector potential a and the normal directional derivatives of all components of a at all boundary points, given the tangential components of a, for a potential that satisfies the Euler- Lagrange equation and the Lorenz condition such as any one of the basis potentials a_a. The value of any field basis function g_a can be calculated from the boundary values and normal derivatives of the potential basis function a_a. A useful method of defining the value and normal derivative of each space-time basis potential function a_a on the boundary of a region may be to require that within the region in which the computation is being performed that the potential satisfies V^■ g_a— δ_{0 a}J_f = 0 (where δ_{0 a} = 1 when a = 0 and δ_{0 α} = 0 for all other values of a). The potential function a and the potential function's normal derivatives at all boundary points are defined in terms of a finite number of variables. The potential a may represent, for example, any one of the basis functions a_a. The following simple model illustrates the method, although the same method can be used to define smoother albeit more complex models, or less smooth models.

[1101] First, a general curved boundary is approximated by connected boundary simplices. For a 2D boundary in 3D space (applicable to cases in which time-varying fields are represented by Fourier components in frequency space, for example), the boundary simplexes are triangles and it is useful to visualize the instant case as representative. The unit vector normal to the surface of the simplex s, may be written as n_s. At each vertex v, a unit vector n_s is specified, which may be normal to the physical boundary that is being approximated by the connected simplices, which may be defined as the normalized sum of the area vectors of the simplices that share that vertex. The resulting potential solution depends only weakly on the choice of normal vectors n_v, so the exact definition is not critical.

[1102] At each vertex v, we define the potential parameters to be the value of the space-time vector potential a_v≡ a[r_v] which have 4 degrees of freedom, and the derivative in the direction n_v of the components of a that are perpendicular to n_v, which are

d_va_vv≡ (n_v · V) (n_vn_vAa[r_v]), and which have 3 degrees of freedom. The derivative in the direction n_v of the component of a that is parallel to n_v is determined from the values of a on the boundary by the medium Lorenz condition. At any one vertex, this derivative in the direction n_v of the component of a that is parallel to n_v is generally different for each triangle that shares that vertex. All 4 components of the potential and all 4 components of the normal derivative of the potential at all points on the boundary are in expressed in terms of only 7 parameters at each vertex. The tangential components of the potential are continuous across any boundary, and so have the same values as the boundary is approached from either side, but the normal component of the potential and all components of the directional derivative of the potential are generally different on the two sides of a boundary between media with different polarizabilities; each side must be treated independently.

[1103] The potential a at any point on a boundary simplex is defined to be the linear function of position on that simplex that matches the vertex values a_v. Recall that the media potential a_m is just equal to a linear transformation of the components of a at each point, so it is possible to use either the function a or a_m as the potential (the form of the resulting equations are unchanged) and transform from one to the other as needed. [1104] Now, as a result of defining the potential a to be the linear function of position on that simplex that matches the vertex values a_v, one can derive expressions for the potential a and for the normal derivative (n_s · V)a of the potential at every boundary simplex point on both sides of the boundary, in terms of potential parameters that have 7 degrees of freedom per vertex. 7 degrees of freedom is too many degrees for the boundary conditions of a space-time vector function that satisfies a wave equation in the region enclosed by the boundary, making the potential in this region over-constrained. Green's theorem is applied to eliminate 4 of the 7 degrees of freedom in the parameters a_v (with 4 degrees) and n_v · Vn_vii_vAa_v (with 3 degrees) at each vertex v, as follows.

[1105] For simplicity, consider only the important practical case of a space composed of stationary bounded regions of uniform polarizability. For convenience, the computation is performed in the 3D algebra of an inertial frame of the boundaries. Assume harmonic time dependence, as is the case for Fourier components, and parameterize the potential by the potential's complex valued components on boundaries. Given such a basis potential a on the boundary of a region V, the potential's value a at any point in the region enclosed by the boundary s = dV of V or on the boundary s, such as a vertex point r_v, is related to the potential value a and normal derivative V_na of the potential on the entire boundary by Green Theorem, a[i ] = /_5=θν(ψ[? - iv a[r] - V_ni|/[f - ivMf])ds + J_v j_f [f] ψ[? - i ]dV (42) where V_n = n · V and a suitable Green's function for the potential a is [^f] = ^ ^eik|?l <⁴³)

[1106] Equations (42) and (43) apply if the medium's polarizability is zero. If the medium's polarizability is not zero, then a similar expression but with a, j_f , V, and k replaced by the media quantities a_m , j_m, V_m, k_m, as described earlier, is used.

[1107] Writing one equation for each boundary simplex vertex v of the simulation with p vertex points, the potential a[f_v] on the left hand side equals a the potential value at the position i of vertex v, and the potential a[f] and the normal derivative V_na[f] of the portential in the integrals is a linear function of the potential parameters a_v and d_va_vv on each triangle of the boundary. Consequently, equations (42) and (43) express linear relationships between the potential parameter a_v for one particular vertex on the left, and all potential parameters, a_v and d_va_vv = n_v · Vn_vii_vAa[iv] , for all vertices on the right. [1108] These p space-time vector equations make 4p scalar equations of constraint (where p is the number of vertex points). These constraints can be used to write all 7p parameters as linear functions of the 3p parameters n_vii_vAa_v (i.e., the components of the potential a that are perpendicular to the vertex normal vector n_v) at any boundary vertex v. The Lorenz or media Lorenz condition can then be applied at any boundary point to also calculate the 8th degree of freedom (the normal derivative of the normal component of the potential) at any boundary point.

[1109] Each basis function aj is chosen to be described by setting one component of the vertex parameter n_vii_vAa_v at one vertex equal to unity, and all other vertex parameters and the source j_f to zero. The basis function a₀ may be chosen to be described by setting the 3p vertex parameters at each vertex to equal zero but setting the source j_f to equal the specified source J_f , or J_m if using media quantities, of the simulation. There are then 3p basis functions a_j and associated coefficients , one corresponding to each of the 3p parameters n_vii_vAa_v.

[1110] Quantities can be numerically evaluated to obtain a 7p x 3p matrix for each region, giving the 4p components of the potential A = a_ac_a and the 3p normal derivatives of the potential components perpendicular to n_v at the p vertices as linear functions of the 3p basis function coefficients C_j and of position on each boundary triangle. The normal derivative of the potential component perpendicular to the simplex at any point on any boundary simplex can be numerically evaluated from these values.

[1111] From the resulting expressions, expressions for the field F as a linear function of the 3p parameters n_vii_vAa_v and for the field G as a linear function of F are written on both sides of every boundary point. The next step is to numerically evaluate the boundary integrals needed to find the coefficients Cj that extremize the action, and numerically solve for the coefficients Cj.

[1112] With a choice of basis functions described by the vertex parameters and within the approximation of using simplices to represent the physical boundary, each integral needed to extremize the action has a non-zero integrand over only a localized region of the boundary for which the tangential components of the corresponding basis function are non-zero. These integrals using localized basis potentials can be much faster to compute than integrals using non- localized basis potentials.

Calculation of potentials and fields at non-boundary points

[1113] Given the values of a potential A and the potential's normal derivatives at all boundary points - such as either a basis potential or the potential that extremizes the action - the potential A[f] and the fields F = VAA and G[F] can be calculated at any point f using Green's theorem as applied to A and to G, using the media quantities A_m and G_m for example.

[1114] To find the potential A in a region with zero polarizability, Green's theorem expression from the previous section can be used with the final potential solution A in place of a.

[1115] To find the field F in a region with zero polarizability, Green's theorem expression for Maxwell's equation can be used,

F[f '] = - f ']nF[f]ds - J_v ψ[? - f ']_Yoj[f]dV (44) where an appropriate Green's function is

Ψ[Γ] = (-ik + )φ[?] where k≡— . Alternatively, F can be numerically calculated from values of A at nearby points.

For a region with non-zero polarizability (including complex electric polarizability, representing electrical conductivity), the same expressions, but with the media quantities G_m, k_m, and j_m in place of F, k, and j, are used.

Additional Constraints

[1116] Each potential and field may be parameterized as a sum of parameters C_j times corresponding basis functions. The method of least action identifies parameter values that extremize the action.

[1117] For example the charge Q = J_v p dV = 0 in a 3D volume V may be rewritten using the

Euler-Lagrange equation V · G = J dotted with γ₀, or V · D = p , to give J_v V · D dV = Q.

Applying the divergence theorem gives J_s=av n · Dds = Q. If the charge is known - typically Q = 0 for non-zero frequencies when using complex amplitudes of Fourier expansions - then this equation is a linear constraint on the basis coefficients C_j via D.

[1118] The action may be extremized subject to a constraint equation such as the preceding example by Lagrange's method. Such constraints may be enforced by defining an auxiliary action function S' given by

S' = S + /_s=av n - Dds (46) where λ is a Lagrange multiplier in the current example, and then finding parameters Cj and λ that simultaneously satisfy the variational equations 5_C.S' = 0 and the constraint equation

[1119] Doing so may improve numerical performance even if not required for solution. For example, numerical errors may result in a very small but nonzero effective variation in the net charge and an associated longitudinal electric field wave, and these errors may be smaller using such a constraint.

Calculating Potential Values Using Virtual Sources

[1120] The phrase "physical source" is used to refer to any charge-current source that exists in the situation being simulated, and the phrase "virtual source" is used to refer to any charge- current source that does not exist in the situation being simulated, but rather is a mathematical artifice used to help calculate potentials and fields in analogy with virtual images in

electrostatics. Virtual sources on boundaries are often called equivalent sources in

computational electromagnetism literature.

[1121] Virtual sources may be used to define a basis potential a_j within any one particular region R. The basic idea is that a fictitious infinite space with uniform polarizability is considered, with fictitious virtual charge-current source jf* in addition to any physical charge- current j_j (take j_j = J for i = 0 and j_f = 0 otherwise), from which the potential that satisfies the Euler-Lagrange equation and the associated field F is calculated. The virtual charge-current source ) is chosen to be nonzero only at points not in region R. That is, the virtual charge- current source jf* is located at points on the boundary or outside of the region. The potential may be chosen to satisfy the wave equation and the Lorenz condition. The virtual charge- current source is parameterized, and the parameters are chosen so that the resulting potential meets the boundary conditions (e.g., the specified scalar potential and tangential vector potential values on the boundary) within an acceptable error. Parameterizing the virtual charge-current source results in an expression for the basis potential a_j at any point f in region R that is an acceptably good approximation, within region R, for the original problem involving a space with boundaries between regions of uniform polarizability. Every region is parameterized, and then the parameterization process is repeated for every basis potential.

[1122] Although the basis potential aj is calculated using the idea of a fictitious space and a virtual source, aj satisfies the Euler-Lagrange equation in the region R and satisfies the boundary condition regardless of the existence of any real source. Consequently, the basis potential a_j is also an acceptable potential function for the original problem even though in that problem there is or may be different polarizability and sources outside of region R. The result is that for each basis potential, for each region with uniform polarizability, there are virtual sources located outside the region (plus the physical source for i=0, located in the region, in an embodiment) from which the potential for that basis function in that region can be calculated.

[1123] If the virtual sources are well chosen, the resulting potential can satisfy the boundary conditions to a good approximation and satisfy the wave equation and the Lorenz condition exactly. For example, the virtual source for each basis potential and each region may be a collection of point electric and magnetic dipoles with well-chosen locations, and with parameterized amplitudes chosen to satisfy the boundary conditions within a good

approximation. Then, both the potential and the electromagnetic field associated with a given virtual source can be quickly calculated using simple well-known expressions for the potential and field due to dipoles.

[1124] By choosing a large enough set of parameters and corresponding virtual sources, the basis potential a_j can be chosen to yield the correct field to an approximation that is as good as desired. Note that the parameters of the basis potential are chosen for each basis potential aj for each region R to make aj acceptably continuous across all boundaries in the problem, and are subsequently treated as constants; the parameters chosen to approximate each basis potential aj are not the parameters Cj of the parameterized potential A =∑; a; that are determined by extremizing the action.

[1125] A useful measure of error in meeting the boundary conditions is the integral over the boundary of | (ΔΑ Λ η)γ₀ |², where ΔΑ is the difference between the space-time potential defined by virtual sources and the specified boundary values, and n is a unit space-time vector field normal to the boundary. Parameters describing the virtual sources for any basis function in any region can be chosen to minimize the measure of error found from the integral over the boundary of | (ΔΑ Λ η)γ₀1².

Calculating Potential Values as Approximate Solutions to the Euler-Lagrange Equation

[1126] Continuous basis functions that exactly or approximately satisfy the Euler-Lagrange equation at non-boundary points and have specified values at boundary points, or that satisfy the wave equation at non-boundary points and the Lorenz condition and have specified Lorenz- gauge-equivalent boundary values, may be defined in various ways. [1127] For example, if a region is very thin compared with the radius of curvature of the region, such as a matching region, or if the wavelength in the region is small compared with the radius of curvature of the boundary of the region, such as a piece of formed sheet metal, it may be possible to quickly identify simple acceptable approximate expressions for the potentials and fields. For example, in the case of the region very thin compared with the radius of curvature of the region, the potential may be chosen as the product of the potential value on the boundary times a function of the perpendicular distance from the boundary chosen so that the resulting potential approximately or exactly satisfies the Euler-Lagrange equation. The function of the perpendicular distance might, for example, be found as or approximated by a power series, an exponential function, or the product of a power series and exponential function.

[1128] In a polarizable medium, the medium potential A_m is generally used and transforms to the potential A after finding a desired solution. Computing A_m can be a useful method, because the differential equation for A_m is simple. Assuming for simplicity that there is no specified charge or current J_m inside the region of interest, the differential equation for A_m is

V_m ²A_m = 0 .

In any particular orthonormal coordinate system {x,y,z}, the equation for A_mbecomes

[1129] The preceding equation for A_m might be used, for example, if a region is very thin compared with the radius of curvature of the region. For simplicity, choose z to be

perpendicular to the boundary in the region of interest. It might be assumed that no quantities, including A_m, change quickly tangentially to the boundary, so that the x and y derivatives may be neglected. Then the above equation for A_m may be approximated as,

[1130] If there are quantities that vary harmonically with time and if the harmonically varying quantities of time are represented with complex fields, the above equation may be replaced by or (k_m ² + d_z ²)A_m * 0 ,

[1131] where k_m = | a)|/v_m is the medium wavenumber. The solution to the preceding equation may be written as an exponential function with a generally complex argument. Other approximations and solutions may also be used. For example, the region may contain a boundary of discontinuity in polarizability, in which case the potential approximation may be continuous but piecewise linear across the discontinuity.

[1132] From expressions for A_m, the potential A may be computed, and then the fields F and then G may be computed on the boundary of the region. Assuming the potential A and field G have already been determined on the other side of the boundary, the change AG in the field G across the boundary can be calculated. Using the calculation of AG, the boundary integrals required for extremizing the action may also be calculated.

METHOD

[1133] FIG. 3 is a flowchart of an embodiment of a method 300 that is implemented by processor system 300. In step 302, a determination is made of how to set up the grid, such as the size of the different regions of the problem area. In step 304, material parameters are chosen for each region. In step 306, a determination is made as to what basis functions to use for solving the electromagnetic fields. In an embodiment, the basis functions are parameterized solution to the electromagnetic equations in a homogeneous region. In step 308, the action is extremized (e.g., minimized), which may involve adjusting the parameters of the parameterized basis functions to extremize the action integral at the boundary between homogenous regions. In step 310, the output of action extremizer 201 is formatted and presented to the user.

[1134] In an embodiment, each of the steps of method 300 is a distinct step. In another embodiment, although depicted as distinct steps in FIG. 3, step 302-310 may not be distinct steps. In other embodiments, method 300 may not have all of the above steps and/or may have other steps in addition to or instead of those listed above. The steps of method 300 may be performed in another order. Subsets of the steps listed above as part of method 300 may be used to form their own method.

[1135] FIG. 4 is a flowchart of an embodiment of a machine implemented method 400 for solving the electromagnetic equations. Method 400 may be an embodiment, of step 308 of method 300. In step 402, the region in which the problem is being solved is divided into a series of homogenous regions, each homogeneous region being a region in which electromagnetic properties of the material such as permittivity, permeability, and electrical conductivity are homogeneous throughout the region. For example, the machine may choose the dimensions of the matrices and the number of iterations in the loops in the computer code that correspond to the representing the fields and parameters of the equations that are appropriate for the chosen regions. In step 404, choose a parameterized potential function that satisfies the Euler-Lagrange electromagnetic equations in the homogeneous regions but not necessarily on the boundaries between regions, which results in the volume integrals that appear in the equation for extremizing the action to be identically zero, leaving only boundary integrals to be calculated, and choose the functions to satisfy a gauge condition to simplify subsequent calculations (The Euler-Lagrange electromagnetic equations are typically written using tensor algebra or geometric algebra of space-time and are equivalent to Maxwell's inhomogeneous equations - the two equations without charge density or current - which are typically written using vector algebra of 3D space. The other two of Maxwell's equations are equivalent to a mathematical identity when written using tensor or geometric algebra of space-time. A practical choice of such a parameterized potential function is a sum of terms, each term equal to a parameter times a basis function. In step 405, the potential is further characterized by choosing how the potential is parameterized (tangential components of the potential values on the boundary in space-time) and calculating an expression that gives other necessary values (the non-tangential component, and all components of the normal derivative) in terms of the parameters (the tangential components of the potential in space-time are equal to the tangential components of the 3D space vector potential and the scalar potential). For example, the user may enter a choice of basis functions and/or the machine may automatically choose the basis functions depending on the homogeneous regions chosen in step 402. In step 406, the equations to be satisfied by the parameters are identified by writing them as symbolic equations appropriate for choices made so far. The machine may automatically write/determine the equations based on the choices made so far. Alternatively, the user may be offered a choice of equations to solve and/or may be provided with a field for entering the equation that the user desires to solve. For the choices made here, each of these equations contains a term for each parameter, multiplied typically by one boundary integral. In step 407, for any one of these equations that has a potentially diverging factor of 1/f where f may be very small or zero, as will typically occur if the problem contains an electrical conductor that is represented by a complex -value permittivity, the equation is first multiplied by f and simplified symbolically to cancel any factors of 1/f before numerical evaluation of coefficients of parameters in the equation. Step 407 may be performed as a result of user input making entering the choice or the code may automatically identify the 1/f dependence and choose the equations resulting from the multiplication by f as the equation to solve. In the last step, 408, the equations are solved by the computer for numerical values of the parameters.

[1136] In an embodiment, each of the steps of method 400 is a distinct step. In another embodiment, although depicted as distinct steps in FIG. 4, step 402-406 may not be distinct steps. In other embodiments, method 400 may not have all of the above steps and/or may have other steps in addition to or instead of those listed above. The steps of method 400 may be performed in another order. Subsets of the steps listed above as part of method 400 may be used to form their own method.

REGIONS AND BOUNDARIES

[1137] FIG. 5 shows an example 500 of homogeneous regions. Example 500 includes region 502, region 504, and boundary 506. One region may have a finite extent and be completely surrounded by the other, and one region may extend to infinity in all directions and completely surround the other. In other embodiments, example 500 may include additional components and/or may not include all of the components listed above.

[1138] Region 502 is homogeneous region 1, which has a homogeneous set of material parameters. In other words, throughout region 502 each of the material parameters is assigned a uniform value. The material parameters may include the permittivity (or susceptibility), permeability, polarization, magnetization, resistivity (or conductivity), for example. Parameter may also be included that specify the charge density and current density. A specified electric polarization or magnetization may be represented by the corresponding specified bound charge or bound current densities. Region 504 is homogeneous region 2, which also has a

homogeneous set of material parameters, which may be different than those of region 502. Boundary 506 separates region 502 and region 504, and in an embodiment, is where the computations are primarily formed. Region 502 and region 502 may represent two physical regions that have distinctly different sets of material parameters with the same values as assigned during the computations. A region having a continuously varying set of material parameters may be modeled in various ways, such as by dividing the region into smaller regions of homogeneous material parameters.

[1139] FIG. 6 shows an example 600 of computational regions corresponding to the homogeneous regions of example 500. Example 600 includes region 602, region 604, and boundary 606. In other embodiments, system 600 may include additional components and/or may not include all of the components listed above.

[1140] Region 602 is an interior region of homogeneous region 1, and region 604 is an interior region of homogeneous region 2. In this specification, region points may generally be considered to exclude boundary points, but the phrase "interior region point" or "non-boundary point" may be used to make it explicitly clear that the point is in the region and not a boundary point if it is important in a discussion that the point is not a boundary point. In region 602 and region 604 no actual computations may need to be carried out. Each potential basis function extends over all of space and the potential basis function, or a potential related by a gauge transformation to the potential basis function, is continuous at all points of space, which in the example of FIG. 6 includes regions 602 and 604. Each such continuous potential function, though, typically has a "kink" across a boundary, corresponding to discontinuity in the derivative of the function. The values of any given basis function aj in either or both region 602 and region 604 may be defined with the help of virtual sources, although different virtual sources will be used for each region and the material parameters, such as polarizability, of each region may differ. The computational parameters, or simply the parameters (e.g., the coefficient C_j multiplying each basis function in the expressions for any of the potential or electromagnetic fields, such as A =∑j aj Cj ), are independent of region or position in space or time.

[1141] Boundary 606 includes points representing a boundary between two regions.

Boundaries are modeled as mathematical surfaces, which are usually smooth. A real physical boundary may be very close to a mathematical surface, such as the boundary between a glass and air, or may have some small, but negligible, thickness, such as the boundary between a piece of wood and adjacent concrete. The computation of the electromagnetic fields by extremizing the action integrals, are primarily performed using integrals over the points of boundary 606, allowing the computational parameters Cj to be chosen to extremize the action.

[1142] The word "region" is generally used to refer to a volume containing material with uniform polarizability since solutions to the wave equation are well known in the instant case, but regions with non-uniform polarizability may also be used, especially if solutions are known. Also, although the polarizability may be discontinuous across the boundary between regions, the polarizability is not necessarily discontinuous. For example, a volume with uniform

polarizability can be divided into two regions, joined by a boundary across which the polarizability does not change. Dividing the volumes into two regions can be convenient for various reasons. For example, dividing the volume into two regions allows the use of different methods of defining basis functions in these two adjacent regions for faster or more accurate computations.

[1143] Volumes in 4D space-time may be denoted by u, and 3D boundaries in 4D space-time may be denoted by β . Volumes in 3D space may be denoted by V, and 2D boundaries in 3D space may be denoted by s.

EXAMPLE USING GREEN'S THEOREM

[1144] In following section illustrating an example using Green's Theorem, the problem of calculating the fields at all points for a simple example of a planar sheet of oscillating current parallel to a planar boundary between two media with different polarizabilities is solved.

[1145] FIG.7 shows a diagram illustrating elements of an example 700 to which the method of least action is applied. Example 700 includes medium 702, medium 704, interface 706, z-axis 708, and origin 710. A person familiar with electromagnetism may also use the ideas of incident wave 712, travel direction 714, transmitted wave 716, travel direction 718, reflected wave 720, and travel direction 722. The ideas of Fig.7 are not needed with the method using Green's theorem, but the ideas of Fig.7 are identified here to help explain the method.

[1146] Consider example 700 in which an oscillating current sheet is the source of a plane wave that is incident on a flat boundary between two regions with homogenous electric and magnetic polarizabilities, with the direction of wave propagation perpendicular to the surface. In FIG.7, medium 702 is the medium on the right side of the diagram and the functions and quantities associated within medium 702 have the subscript R. Medium 702 has a permittivity 8R and a permeability of κ. The phase velocity of waves in medium 702 regions is v_R = 1/V -R^£R_? ^and the wave number in medium 702 is k_R = o>/v_R, where ω is the angular frequency of the wave. In other problems other factors may be considered.

[1147] Similarly, medium 704 is the medium on the left side of the diagram and the functions and quantities associated within medium 704 have the subscript L. Medium 704 has a permittivity £_L and a permeability of μι_^. The phase velocity of waves in medium 704 regions is v_L = 1/VM-L^£L_> ^and the wave number in medium 702 is k_L = w/v_L. It may be assumed that mediums 702 or 704 are idealized in the sense that neither is conductive and both transmit electromagnetic waves without loss, or a complex permittivity may be used to represent electrically conductive media with resistivity. [1148] Interface 706 is the interface between the two media. In the coordinate system of FIG. 7, z-axis 708 indicates the position of the waves traveling within mediums 702 and 704. Origin 710 is the point on z-axis 708 at which the value of the position coordinate is 0. Although the physics does not depend on the coordinate system, the choice of coordinates simplifies the computations, so that issues that do not relate to the demonstrating the method do not need to be discussed.

[1149] Incident wave 712 is a plane wave of light entering from the right hand side of FIG. 7. The source of the incident wave is modeled as a sheet of oscillating current perpendicular to the z axis. The source of incident wave 712 and the material parameters of media 702 and 704 are the inputs to the problem. The incident wave (that would be due to the source if the incident wave were not interacting with the system being simulated) may alternatively be an input instead of the source of the wave, but in the current example the simpler method of using the source as input is applied. Travel direction 714 is the direction in which incident wave 712 travels.

[1150] Similarly, reflected wave 720 is a wave of light that was reflected from interface 706 as a result of incident wave 712 hitting interface 706. Travel direction 722 is the direction in which reflected wave 720 travels, which is the opposite direction as the travel direction 714 of incident wave 712. Whereas, transmitted wave 716 is a wave of light that was transmitted through interface 706 as a result of incident wave 712 hitting interface 706. Travel direction 718 is the direction in which transmitted wave 716 travels, which is the same direction as incident wave

712 (the wave vector of transmitted wave, k_L, is therefore parallel to k_R, but has a different magnitude).

[1151] Using the method of least action to solve for the resulting fields in example 700, a set of basis potential functions must first be specified. Applying the method of least action requires specifying boundary values and the differential equation that is satisfied at non-boundary points. It is not yet necessary to actually find expressions for the basis functions, and it is simpler not to.

[1152] The current example is done using Green's theorem methods. With Green's theorem, the basis functions satisfy the wave equation and the Lorenz condition in media at all non- boundary points, so that the volume integrals in the action are zero. The basis function a₀ is chosen to satisfy V^■ g₀ = J for the specified current source J and a₀ = 0 on the boundary, and basis function a_j to satisfy V^■ = 0 and a_j = 1 on the boundary. There are no other independent boundary conditions that meet the given symmetries of the instant problem, so there are only the two basis functions a₀ and needed to solve the problem of example 700.

[1153] Applying the method of least action requires the knowing of the parameterized value of the potential A on all boundaries, and of the change AG in the electromagnetic field G across all boundaries. The potential A = a₀c₀ + ajCj is directly parameterized on the boundaries, so the potential's value at boundary points is known. The field G is a function of the field F, which is a function of derivatives of the potential A. Specifically, the field F is given by F = VA = n n^■ VA + n n Λ VA, where n is a unit space-time vector perpendicular, or normal, to the boundary. The derivative tangential to the boundary, which is n Λ VA, is zero by the assumed symmetry in the current problem, and the derivative perpendicular to the boundary, which is n^■ VA , can be found using Greens function methods as follows.

[1154] Let the permittivity be e_R and permeability be μ_β for 0<z (points to the right), and e_L and μι_^ for z<0 (points to the left), and let the physical source J[z] be a an oscillating uniform current sheet located at z = Z_j > 0 with current in the x direction, perpendicular to z. The current sheet radiates with wave number k_R .

[1155] All information necessary to solve the instant problem with only the given information and using Greens function methods can be found without explicitly finding any basis functions. The method of extreme action is done first, and then the basis functions using Greens function methods are found explicitly.

[1156] Let the current sheet be oscillating parallel to the x axis, which has unit vector x. The parameterized potential A to be used to extremize the action is the sum of the basis potential functions, which may be written as,

A[z] = c₀a₀ [z] + C_ja z]

[1157] To extremize the action, the variation 5S of S with respect to variation of the parameters Cj for i≠ 0 must be equal to zero. In the current example the only non-zero value of i is i = 1, so the requirement that the variation of the action be equal to 0 may be written as

2L ^■ (n^■ Ag₀)dxdy c₀

[1158] By assumption in the instant problem all functions are independent on the position {x, y} on the boundary, which implies that the integrands are equal: aj^■ (n^■ Agi) Cj = -a ^■ (n^■ Agi) c₀

[1159] Rewriting the 4D space-time quantities in the preceding equation in terms of corresponding 3D quantities in an inertial frame of the boundary with frame vectors parallel to the coordinate axes, and using the fact that with only a 3D vector source (the current sheet) the scalar potential is zero in the current example and c₀ = 1, results in the equation a ^■ (z x Ahi) c = —a ^■ (z x Ah₀) so that cl satisfies

[1160] This represents the solution that extremizes the action. But at this point we do not yet have expressions for the changes Ahj in the magnetic field basis functions across the boundary.

To calculate a numerical solution an expression of the changes Ahj in the magnetic field basis functions across the boundary are needed. Expression of the changes Ahj in the magnetic field basis functions across the boundary are needed.

[1161] Changes Ahj in the magnetic field basis functions across the boundary are calculated using the field fj = v_3i = VyoYoaj = d_t - v) (ε₀φ_ί

Using complex vector notation for harmonic variation at frequency ω and the fact that in the current example the scalar potential φ_ί equals zero (since there is no charge density source) and since therefore V^■ a also equals zero (by the Lorenz condition) , fj is given by ϊω _Λ 1 _ \ 1 _ 1

f_j = V a. = -— iooa; H V Λ a.

V c A εμ₀ / ε² μ₀ εμ₀

→ 1 ^→

[1162] Comparing the above expression for fj with the expression F = ε₀θ; H III b;, results in and e, = iooa_j

[1163] In the current example with all fields independent of x and y, we have and therefore

L 1

h: =— b: =— z x d .a.

μ μ

[1164] Next the normal derivatives d_za^* _j of the potential on the boundary must be found.

Although in this simple example it is possible to correctly guess that a simple plane wave is a useful basis function and using plane waves the normal derivative and other quantities can be calculated, for the general case of arbitrary curved boundaries convenient basis functions are difficult to guess. Consequently, in this example the general method using Green's theorem methods is illustrated. The normal derivative values of the Green's theorem method will allow the calculation of the required changes Ahj in hj.

[1165] Greens' theorem can be inverted to find the normal derivative of the potential on the boundary as a function of the values of the potential on the boundary. The general 3D result as applied to a complex 3D vector potential a^* _j may be written

¼ [?'] = jf ("fi^■ ν_Γψ[? - f ']¾ [?] + ψ [Γ - Γ'] η - ν_Γ3_ί [Γ])ά₅ + JJJ i|/[f - f']J [f]dV s=av v

We can evaluate these integrals using polar coordinates {p, <p, z} to describe vectors f = p cos[<p] x + p sin[(p] y + z z with the z axis chosen to be perpendicular to the boundary. We take the current density J_j to be uniform and nonzero for i=0 at z = Z_j > 0, and equal to 0 for i=l . The integrals can be evaluated with the result a z'] = Sign[z']_X' [z - z']¾[z] - Sign[z']x[z - ζ'] ¾'[z] + χ[ζ, - ζ'] h where χ[ζ] = jj ψ[ρ, φ, ζ] ρ dp d<p = - ^ e^ik^*

and we've used the fact that ii is by definition a vector that points away from the region containing z' to write ii - V_r = -Sign[z']d_z

Although the basis functions and Greens function are continuous at all points, their derivatives are discontinuous at z=0, so we must be careful when evaluating quantities near z=0. We evaluate the boundary integral at z = +e/2 or z =— e/2 where e is a very small positive number, depending on which side of the boundary we're interested in. For functions that are continuous across the boundary, we can replace e with 0. The equation for a^*j [ζ'] then becomes

¾ [z'] = Sign[z']_X' [-z']ai [0] - Sign[z']x[-z'] 3_ζ¼ [±ε/2] + χ[ζ, - ζ'] Jj z' is now evaluated close to the boundary, at some small number with the same sign as z but with at least slightly larger magnitude than z - we use z' = ±e. Then

The source current amplitude Ji is chosen so that source current amplitude Ji creates a potential of unit amplitude at the boundary for i=0 of -7- e^lkzJjj = xS_{i 0}. Then d_zai [±e] = Sign[±e]ik_±(ai [±e] + 2^χδ _ο)

This is 4 equations, for +e and— e, for each value 0 or 1 of i. These 4 equations are d_za₀ [+e] = 2ik₊ d_za₀ [-e] = 0 d_zai [+e] = ^ik+ 9ζ¾ [-ε] = -ik+ [1166] Next, these derivatives are used to find the change across the boundary in the field:

1 1

A ₀ = h₀[+e] - h₀[-£] =— z x d_za₀[+e] z x d_za₀[-e] — 2yik_R-0

-2 i— ^kR and

1 1

Ahj = h +ε] - h -ε] =— z x d^ +e] z x d^ -s]

[1167] Then the transmission amplitude is

with transmission intensity

[1168] These results can be used to write either basis function at any point. For example, for i=l and the region L for which z' < 0, the i=l basis function is

¾[z'] = -a_nx[-z']a_![0] +χ[-ζ']β_η3₁[0] = +^Exp[-ik_Lz ] + ^Exp[-ik_Lz ]

= +Exp[— ik_Lz']

[1169] while the i=0 basis function equals 0 for z'<0, so the potential A for z' < 0 is [1170] A[z'] = C_!Exp[-ik_Lz']

[1171] which represents a leftward traveling wave on the left side of the boundary.

[1172] The expression for A[z'] might be used at the "start" of a solution by someone with some familiarity and intuition about the problem and using a different solution method. A key feature of the method disclosed here is that even for complex 3D problems, no familiarity or intuition is needed because the solution method is an unambiguous recipe that applies to all circumstances with no free parameters and can be easily calculated on a computer.

A 3D EXAMPLE

[1173] FIG. 8 shows a sample problem 800. Sample problem 800 includes dipole 802, polarizable sphere 804. Dipole 802 has a moment pointing in the direction of the arrow representing the dipole. Polarizable sphere 804 is made from a homogeneous polarizable material and the space outside of the sphere is nonpolarizable. The surface of the sphere 804 is approximated by connected triangles. The objective of problem 800 is to solve for the electric and magnetic fields induced by dipole 802 inside and outside of the sphere. The problem will not be completely solved here, but the calculation of basis functions using two methods that may be used to solve the problem, Green's theorem and virtual sources, will be discussed.

[1174] Space may be divided into two regions for the purposes of computing the

electromagnetic field in the example of FIG. 8. One region may be the interior of the sphere 804, the other region may be everything outside of sphere 804, and the surface of the sphere is the boundary between these two regions. Other methods of defining basis functions may additionally divide either or both of these two regions (the region including the inside of sphere 804 and the region including everywhere else) with homogeneous polarizability into two or more regions.

[1175] A general method for defining a set of basis functions using either Green's Theorem or virtual sources is to first i) for each basis function, specify the values of the tangential components of the space-time vector potential field at all boundary points between

homogeneous regions, then ii) specify expressions or a method for quickly calculating potential values and field values at any boundary point from the specified tangential potential

components, for a potential that satisfies the Euler-Lagrange equation at all non-boundary points. Specifying the values of the tangential components of the space-time vector potential field at all boundary points between homogeneous regions, and specifying expressions or a method for quickly calculating field values at any boundary point allows the action to be extremized. After extremizing the action, it is necessary to iii) quickly calculate at least field values and possibly potential values at any non-boundary point.

[1176] First consider step i) in FIG 8. For basis index i = 0, choose the potential be zero at all boundary points, and to satisfy the Euler-Lagrange equation, with J equal to the specified source dipole 802, at all non-boundary points. For each basis function i≠ 0, choose the potential to be nonzero at only one vertex on the boundary, to be a linear function of position on each triangular segment of the boundary, to be continuous at all points, and to satisfy the Euler-Lagrange equation with J = 0 at all non-boundary points. For each boundary vertex, define a vector ii that is normal to the physical boundary (such as, for example, an average of the normal vectors of the triangles that share that vertex), and choose three basis potential functions: one with nonzero scalar potential, and two with specified vector potential values perpendicular to the normal 3D space vector at that vertex (the normal vector component is one of the unknowns that is found).

[1177] Now consider step ii) in FIG 8: step ii) can be done using Green's theorem following the recipe previously outlined. Alternatively, if virtual sources are being used, for each basis function, choose a set of virtual sources inside the sphere to define the basis function outside the sphere such that the discrepancy between the potential calculated from the virtual sources and the specified boundary values for the basis potential is minimized, or more accurately reduced to within an acceptable tolerance. Similarly, a set of virtual sources may be chosen outside the sphere to define the basis function inside the sphere. Each virtual source may be a

parameterized collection of point charges, point electric dipoles and/or magnetic dipoles, or more complex sources. After the virtual sources are chosen to define a basis potential, the same virtual sources may be used to calculate the corresponding electric and magnetic field F or G. Next, for each of the basis potentials, which are derived from the physical and virtual sources, the potential inside the sphere may be defined as a linear combination of potentials due to virtual sources located outside the sphere. For convenience a linear combination of the same set of virtual sources may be chosen for each basis function, but for each basis function the

coefficients of the combination is chosen that gives the smallest integrated error relative to the specified values of the basis function on the boundary. Note that these coefficients are in addition to the parameters Cj used to extremize the action.

[1178] There is now a set of basis functions, aj , each one defined over all of space, that can be used to define a parameterized potential A =∑j aj Cj . The parameters Cj are chosen to extremize the action. The parameterized potential A results in an associated set of volume integrals that equal zero at all non-boundary points for any set of parameter values C_j. Although the potential resulting from using virtual sources may only approximately satisfy the Lorenz-gauge- equivalent boundary conditions, more virtual sources can always be chosen, so that the error is smaller than any required value.

[1179] The following paragraphs provide greater emphasis on technical details. [1180] FIG. 9 shows an example 900 of a placement of virtual sources inside the sphere.

Example 900 shows sphere 804 with a polyhedron 902 inside having virtual dipoles 904 at the vertices of the polyhedron that is inside the spherical boundary. Polyhedron 902 is used to aid in locating the dipoles 904, which are placed close to the inner surface of the sphere 802. Dipoles 904 are used to generate the basis potentials aj outside the sphere. The potential in the outer region from each dipole 904 in the inner region is used as a different basis function a_j having a unit dipole moment. For each basis potential aj defined outside the sphere with the help of just one virtual dipole, the basis potential aj is defined inside the sphere as a linear combination of potentials due to virtual dipoles located outside the sphere (which are discussed further in conjunction with FIG. 10, below), with the coefficients chosen to minimize (reduce to within a predetermined acceptable tolerance) the effective discontinuity of the basis potential aj across the surface of the sphere. A possible measure of effective discontinuity is described in another section. The end result is a set of basis functions, with each basis function defined outside the sphere by one virtual dipole (or one real dipole in the case of i = 0), and defined inside the sphere by a linear combination of virtual dipoles (a different combination for each value of the index i). The coefficients of the linear combination used to define the basis function for each basis index i are determined early in the solution by minimizing the discontinuity at boundary points, and are in addition to the parameters C_j that are later chosen to extremize the action.

[1181] FIG. 10 shows an example 1000 of a placement of virtual sources outside the sphere. Example 1000 shows sphere 802 with a polyhedron 1002 outside having virtual dipoles 1004 at the vertices of the polyhedron. Polyhedron 1002 is used to aid in locating the dipoles 1004, which are placed close to the outer surface of the sphere 802. Dipoles 1004 are used to generate the basis potentials a_j inside the sphere. For each basis function a_{i ?} already defined in the outer region and on the spherical boundary as the potential resulting from one dipole source, the potential a_j is defined inside the sphere by a linear combination of the potential resulting from dipoles 1004 outside the sphere, with the amplitudes of the dipoles chosen to minimize the discontinuity in the basis potential at the boundary. Each resulting basis potential a_j is defined at all points in space, and is used to define the parameterized potential A =∑j aj Cj .

[1182] FIG. 11 shows an example 1100 of an arrangement of dipoles on the polyhedron.

Example 1100 shows polyhedron 902 with virtual dipoles 904 at the vertices of the polyhedron and electric dipole 1102, magnetic dipole 1104, and magnetic dipole 1106. Electric dipole 1102 points towards the surface of the sphere and is perpendicular to the surface of the sphere. Magnetic dipoles 1104 and 1106 are perpendicular to the electric dipole 1102 and are therefore parallel to the surface of the sphere at the point on the sphere to which the electric dipole points.

[1183] It can be shown that only the scalar potential and the tangential part of the vector potential need to be defined and be continuous at boundary points, in order to specify a potential throughout all of space that uniquely determines the fields. All possible potentials can be specified, to some level of approximation, by dividing the boundary into boundary sections, and specifying the average scalar potential and tangential vector potential on each boundary section. One simple choice of dipoles to generate an approximation to such potential values on the boundary, and to generate the associated potential throughout a region, is one electric dipole near each boundary section with a moment perpendicular to the nearby boundary, and two magnetic dipoles near each boundary section with moments tangential to the nearby boundary. Three such dipole moments at each vertex are represented in Figure 11, which represents 36 dipole moments and corresponding degrees of freedom.

[1184] At a frequency equal to zero, one pair of virtual dipoles - one inside and one outside the sphere - at the correct locations can be used to define the exact solution for the fields of the simple ideal case of a real dipole located outside of a polarizable sphere. But in more general cases, such as non-zero frequencies, non-spherical objects, and applied fields from sources more complex than point dipoles, the exact solution method of solving for the field resulting from a charged object outside of the shape does not apply and may be difficult to implement even approximately, while the method of least action with basis functions defined with the help of virtual sources may be applied with good results.

[1185] Each basis function a_j (resulting from a dipoles near the surface the object inside and outside of the object) described above is approximately continuous everywhere, and a linear combination A =∑j aj Cj of these basis functions satisfies the Euler-Lagrange equation at every non-boundary point for any set of parameter values C_j for i≠ 0 and with c₀ = 1, but does not generally satisfy the Euler-Langrange equation at boundary points. As outlined earlier, values of the parameters Cj in the linear combination of basis functions that form the parametized potential extremize the action.

[1186] Although point dipole virtual sources are very useful for representing and visualizing geometry, and have very simple and easily calculated potentials and fields, for smoother fields, smoother virtual charge and current densities may be used instead of point dipoles.

EXAMPLE OF A 2D SIMPLEX MESH [1187] FIG. 12 shows an example 1200 of a mesh of simplices that may be used computing electromagnetic fields. Example 1200 includes wires 1204 and 1206, and simplices 1208 having nodes 1210 and boundaries 1212.

[1188] FIG. 12 shows a two-dimension problem, which was chosen because the situation of FIG. 12 can be graphically represented in 2D. FIG. 12 there is a metal object that varies in the x and y directions (in the plane of the paper), and extends very far in the z direction (perpendicular to the paper) so that the situation may approximate a 2D problem, with all quantities

independent of z and with harmonic dependence on time.

[1189] Wires 1204 have a current that travels out of the paper in the left (which are represented by large grey dots), wires 1206 have a current that travels into of the paper on the left (which are represented by large black dots). Wires 1204 and 1206 may represent a coil that is stretched in the z direction.

[1190] Computation of the resulting field may be performed on a mesh of simplices. Simplices 1208 are an example of a mesh of simplices that may be used for performing the field computations. Nodes 1210 are the vertexes of the simplices, and boundaries 1212 are the walls between the simplices (connecting the vertexes). The size of the simplices may be varied according to the amount of variation expected in a particular region. Areas with higher variation in the field may include a larger number of simplices, but the simplices will tend to be smaller than other areas, and areas with less variation in the fields may be have a smaller number of simplices, but the simplices will tend to be larger.

[1191] Such a problem can usefully model characteristics of a real 3D situation, but with much less computation. The equations are solved using integrals over the interior of each simplex 1208 and integrals over the boundaries 1212 of each simplex. As in FIG. 12, simplices need not be used throughout the entire region in which the field is being computed, but rather may be used, for example, for regions containing nonlinear media such that the field may be

approximated as being uniform within each simplex.

[1192] Other methods of defining computationally useful regions of space may also be used. For example, if the physical situation is symmetric with rotation about an axis, simplices in a 2D half-plane bounded by the axis may be defined, and then 3D volumes of rotation may be defined and the 3D volumes may be associated with virtual sources by rotating the simplices about the axis. Nonlinear Embodiments

[1193] The potential A is a function of the parameters c;. For simplicity take the

parameterization to be linear in the parameters. If a parameterization is chosen that is nonlinear in the parameters, a solution may be found that extremizes the action by any of various other methods, such as Newtonian iteration with linearized approximations to the potentials and fields.

[1194] Also, for simplicity it has been assumed that the fields G, D, and H are linear functions of the fields F, E, and B, respectively, but these may be nonlinear functions in which case the Lagrangian density may need additional corresponding terms and solutions may be found that extremize the action by any of various methods, such as using perturbation theory and choosing a set of basis functions that includes a fundamental frequency and harmonics.

[1195] An example of a nonlinear field G is an electrically and magnetically polarizable medium with an electric polarization P that saturates according to

E

P[E] = a_EE₀ ArcTan[— ] and magnetic polarization M that saturates according to

B

M[B] = o _BB₀ ArcTanf— ]

^Bo where a_E, E₀, a_B, and B₀ are constants that characterize the medium. The corresponding electric and magnetic interactions are

U_E [E] = ^ a_EE₀ ² Log[l + (E/E₀)²] U_B[B] = (a_B/c²)B₀ ² Log[l + (B/B₀)²] and the Lagrangian is

L = (1/2)F - F - A - J_f + U where U = U_E + U_B represents the interaction between the field and the medium. Setting to zero the variation of the action with respect to parameters gives the same set of equations (19) but with a_EE 1 / a_BB \

G≡ F + V_FU = E + _{- 7} + - III B - - " .„ _{- 7} )

l + (E/Eo)² c V 1 + (B/B₀)V where V_F is ^me multivector derivative with respect to F. Because of the denominators in this expression, G is a nonlinear function of E and B and therefore of F, so equation (19) is no longer linear in the solution parameters C_j and computing a solution may be more difficult.

[1196] In the limit of very large parameters E₀ and B₀, the example above reduces to expressions for linearly polarizable media. In particular, U_E→ - E^■ P where P = a_EE and

1 1

U_B→ - ^2 B^■ M where M = a_BB. We can combine this with F^■ F = E²— (l/c²)B² to rewrite the Lagrangian as L = F^■ G— A^■ J_f where G = (E + P) + (^j III(B— M), or alternatively and equivalently as

L = - F - F— A - J_f + U

2 " where U = ^ (E^■ P + ^ B^■ M) represents the interaction between the electromagnetic field and the linearly polarizable medium. CONTINUITY ERRORS

[1197] In an embodiment, least action method requires that the potential be continuous at all points, including boundary points, and so computation by the method of least action is most easily analyzed for errors and understood if the potential is continuous at all points. But the potential may be discontinuous in at least two acceptable ways and the method still be useful. This is detailed below.

[1198] The first way that the potential may be discontinuous is that it may be "effectively continuous" although not continuous, where effectively continuous means that a gauge transformation can remove any discontinuity. Specifically, even after requiring the Lorenz condition, 0 = V^■ A in terms of space-time quantities (or 0 =

(l/c)d_t + V^■ A in terms of 3D space quantities), for any solution A = γ₀(φ— A) to these equations, A' = A + \7ψ is also a solution if the scalar field ψ satisfies \7²ψ = 0 (or equivalently, in terms of 3D space quantities, at least for stationary media, A'

Φ = Φ H is also a solution if the scalar field ψ satisfies \7²ψ -— =- = 0). Since this is a 2nd

at at² '

order differential equation, the value of ψ is uniquely determined at every point of a 4D space- time volume, if either ψ or the directional derivative of ψ perpendicular to the boundary is specified on the 3D boundary of the 4D volume. We can therefore always choose such a field ψ to exactly cancel any discontinuity in the component of A perpendicular to any boundary in 4D space-time (or in terms of 3D space quantities, at least for stationary media, to exactly cancel any discontinuity in the component of A perpendicular to any boundary in 3D space). A potential that is continuous except for a discontinuous perpendicular vector part at a boundary is therefore effectively continuous and is not a problem.

[1199] The second way that a basis potential may be discontinuous is that the basis potential may be only approximately effectively continuous across a boundary, due to imperfect matching of basis potential values at points on the boundary between the two adjacent regions. The imperfect matching of basis potential values at points on the boundary, generally occurs with the method of virtual sources but not with the Green's Theorem method. In this case we may choose the functions defined on the two sides of the boundary for each basis potential so that some measure of total error due to discontinuity across the boundary is minimized. One useful measure of dissimilarity between the values of basis potentials on side I of the boundary, ^ai' ⁼ (εοΦϊ' + )YO_> ^and the values of basis potentials on side II of the boundary, a^¹ =

[1200] | (Δ , Λ ^η)γο Ι² = ε2 (φ_ι ^Ι - φ_ι ¹¹)² + ((a,¹ - a,¹¹)² - (fi^■ (a,¹ - a,"))²) [1201] where = a^ - a^¹

[1202] This expression for the local error in basis function i at a point on the boundary due to discontinuity across the boundary excludes the perpendicular vector component since the vector component perpendicular to the boundary does not need to be continuous as discussed above. The corresponding total error is the boundary integral [1203] This error expression can be used by choosing quantities that are used to define either or both (φ_ί', a ) and (φ_ί", a ¹), such as virtual source amplitudes, such that errors is minimized. This may be done before the action is extremized. This error in continuity of a basis potential function at boundary points is independent of the results of extremizing the action.

[1204] Defining potentials so that the potentials are identically effectively continuous at boundaries, such as by using potential basis functions parameterized by values of tangential components on boundaries as described previously, has advantages of easier theoretical analysis and more accurate results, but potentials defined so that they are only approximately continuous at boundaries as outlined above may be simpler and faster to compute.

APPLICATIONS

[1205] The system and method may be used to design the shape, placement, and materials for magnetic recording write heads, antenna, the core of a transformer, the core of electromagnets, permanent magnets, and/or electromagnets for generators, and/or transmission lines, electric conductors and electronic components in larger assemblies such as computers and cell phones, microwave devices, and optical devices, including optical and electro-optical integrated circuits.

ALTERNATIVES AND EXTENSIONS

[1206] Each embodiment disclosed herein may be used or otherwise combined with any of the other embodiments disclosed. Any element of any embodiment may be used in any

embodiment.

[1207] Although the invention has been described with reference to specific embodiments, it will be understood by those skilled in the art that various changes may be made and equivalents may be substituted for elements thereof without departing from the true spirit and scope of the invention. In addition, modifications may be made without departing from the essential teachings of the invention.

Claims

1. A machine implemented method comprising: computing, by a machine having a processor system including one or more processors and a memory system, numerical values for a set of parameters that extremize an action integral for a parameterized electromagnetic potential and an electromagnetic field associated with the parameterized electromagnetic potential,

the set of parameters, computed by the machine, extremizing the action integral, the action integral being an integral of an electromagnetic Lagrangian computing, by the machine, the electromagnetic field based on the set of parameters that extremize the action integral.

2. The method of claim 1, further comprising:

the parameterized electromagnetic potential being a linear combination of basis functions that satisfy an Euler-Lagrange equation for non-boundary points, and as a result of the parameterized electromagnetic potential being the linear combination of basis functions that satisfy the Euler-Lagrange equation for non-boundary points, the action integral including contributions only from integrals including boundary points.

3. The machine implemented method of any of claims 1 or 2, wherein time-dependent quantities vary harmonically with time at one frequency, and wherein if an equation that results from extremizing the action integral has a factor of

1/f, where / is a function of frequency that may be zero for at least one value of frequency, then this equation that results from extremizing the action integral is multiplied by a factor of f and simplified symbolically to cancel the factor of 1/f before a full set of equations is solved numerically for the set of parameters.

4. The method of any of claims 1-3, further comprising: constraining the potential by a gauge condition.

5. The method of any of claims 1-4, the gauge condition resulting in the Euler-Lagrange equation simplifying to a wave equation for a linear transformation of components of the electromagnetic potential.

6. The method of any of claims 1-5, wherein the gauge condition is a media Lorenz condition.

7. The method of any of claims 1-6, further comprising: parameterizing the electromagnetic potential to form the parameterized

electromagnetic potential by requiring that the electromagnetic potential agree with a parameterized potential function on boundaries between regions of computation, the potential function satisfying an Euler-Lagrange equation at non-boundary points.

8. The method of any of claims 1-7, further comprising: determining values of the electromagnetic potential at boundary points that are not given by the parameterized potential function, and values of normal derivatives of the electromagnetic potential at region points adjacent to the boundary of the region, by requiring that Green's theorem for the Euler-Lagrange equation and a gauge condition are satisfied at region points adjacent to the boundary of the region, and computing the electromagnetic potential and electromagnetic field in integrals over the boundaries that result when extremizing the action.

9. The method of any of claims 1-8, wherein the parameterized potential function at boundary points is the value of the potential's components tangential to the boundary in space-time.

The machine implemented method of any of claims 1-9, the computing of the set of parameters including at least solving a set of equations in a plurality of regions within which the polarizability is differentiable with respect to position and separated by boundaries, the potential being a linear combination of one or more basis potentials, with at least one basis function being defined in at least one region as a potential function of at least one or more virtual charge-current sources, the virtual charge-current source being located outside the region.

11. The machine implemented method of any of claims 1-10, the computing of the set of parameters including at least solving a set of equations in a plurality of regions within which the polarizability is differentiable with respect to position and separated by boundaries on which the polarizability is not differentiable, the potential being a linear combination of one or more basis potentials, with at least one basis function being defined in at least one region as a potential function of at least one or more virtual charge-current sources, the virtual charge-current source being located outside the region.

12. The machine implemented method of any of claims 1-11, further comprising: solving a set of equations in a plurality of regions, the potential being a linear

combination of one or more basis potentials, with at least one basis function being defined in at least one region as a potential function of:

at least one or more virtual charge-current sources, the virtual charge- current source being located outside the region, and

real charged current sources located inside the region.

13. The machine implemented method of any of claims 1-12, further comprising: solving a set of equations in a plurality of regions including one or more simplices, and the parameterized potential in one or more simplices satisfying the Euler-Lagrange

equations within a tolerance.

The machine implemented method of any of claims 1-13, wherein the

extremizing includes solving a set of equations

where

/ indicates to perform an integration,

u is a four dimensional volume having three spatial dimensions and one time dimension over which the integration is performed,

F is an electromagnetic field strength,

A is the electromagnetic potential,

J is free or total charge-current density,

U is the interaction between the field and polarizable media,

C₍ are the parameters of the potential A such that A ~ a_tCi and are parameters of the electromagnetic field F such that F ~∑_£ f_tCi where ¾ and are basis functions of A and F, respectively,

such that CLi is a function of space-time position, if is not a continuous function of space-time position, any discontinuities present are removable by a gauge

transformation,

and fi is a function of space-time position given by fi = V Λ

15. The machine implemented method of any of claims 1-14, wherein

U represents interaction between the field and linearly polarizable media, in which case the set of equations to be solved may be rewritten as

where

C ~∑i giCi is an electric displacement and magnetic field intensity where

i is a linear function of the field fi, and

if is free charge-current density.

16. The machine implemented method of any of claims 1-15, wherein the set of equations when written as a matrix equation are

M C = B

for an array C with elements equal to the unknown parameters c- for j≠ 0, where M is a matrix with elements m_i - for indices i≠ 0 and j≠ 0, and B is an array with elements bi for indices i≠ 0, given in terms of 4D space-time quantities by

mt_j = j a {V - g_j)dv + j a (ij - g_j)d

bi = - j a (V - g₀ - J) dv - j a (j] - Ag₀ - K) dp

in which integration over v is performed by a 4D volume integral and integration over β is performed by a 3D boundary integral in 4D space-time.

17. The machine implemented method of any of claims 1-16, wherein a media for which the computing is performed is

stationary in one inertial reference frame,

time dependence is harmonic at frequency a), and

the matrix and array elements are computed in terms of 3D space quantities by ^mi_j = j (φ_ϊ ^*(Υ ¹ d_j) - a ^■ (V x h_j + io>d )) dV ft, = - / (♦,·(? · 3o - P) - ¾^*■ (? x ¾ + _iw3o -/))

- / (♦,(. . * - .) - . (. * * - *)) * in which integration over s is a 2D boundary integral and integration over V is a 3D

volume integral; the 3D volume integral being neglectable if the 3D volume integral is equal to zero or less that a predetermined tolerance, in which case the matrix array elements are determined by the 2D boundary integrals in 3D space, by

™ij = / (φ_ϊ ^*(^η■ Ad_j) - ¾^*■ (n x ΔΛ,-)) ds bi = — / (φι (η^■ d₀— σ)— ¾^*■ (n x ΔΛ₀— K)^) ds

where

/ indicates to perform an integration,

Φι is an ith scalar potential basis function

n is a unit vector normal to a boundary of computation cell

Ad_j is a difference between the jth displacement basis function on a first side of the boundary and a jth displacement basis function on a second side of the boundary, the jth displacement basis function on the first side of the boundary being a displacement vector resulting from a jth scalar potential basis function and vector potential basis function on the first side of the boundary, and the jth displacement basis function on the second side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the second side of the boundary,

CLi is an ith vector potential basis function,

Ah_j is a difference between a jth magnetic intensity basis function on the first side of the boundary and a jth magnetic intensity basis function on the second side of the boundary,

d_j is a difference between a jth displacement basis function on one side of the

boundary and a jth displacement basis function on another side of the boundary, the jth magnetic intensity basis function on the first side of the boundary being a magnetic intensity vector resulting from the jth vector potential basis function on the first side of the boundary, and the jth magnetic intensity basis function on the second side of the boundary being a magnetic intensity vector resulting from the jth vector potential basis function on the second side of the boundary,

dV indicates that the integration is over the volume of the computation region, ds indicates that the integration is over the boundary of a computation region

Cj are unknown parameters that are being solved for,

Ad₀ is a difference between a 0th displacement basis function on one side of the

boundary and a 0th displacement basis function on another side of the boundary, the 0^th displacement basis function is a displacement function with a coefficient that has a known value,

Ah₀ is a difference between a 0th magnetic intensity basis function on one side of the boundary and a 0th magnetic intensity basis function on another side of the boundary, the 0th magnetic intensity basis function is a basis function with a known coefficient,

K, p,J, σ are the known volume charge density, volume current density, surface charge density on a boundary, and surface current density on a boundary, respectively, and

c₀ is the coefficient that has a known value.

18. A system comprising :

a processor system having one or more processors; and

a memory system including a one or more non-transitory machine readable media; the non-transitory machine readable media storing one or more machine instructions, which when implemented by the processor system cause the processor system to perform a method including at least

computing, by the machine, numerical values for a set of parameters of a

parameterized electromagnetic potential, by the processor system choosing parameter values that extremize an action integral whose integrand includes at least a Lagrangian density with terms for at least an electromagnetic field.

19. A non-transitory computer readable medium storing one or more machine

instruction, which when implemented by a processor causes the processor to perform a method comprising: computing, by a machine having a processor system including one or more processors and a memory system, numerical values for a set of parameters of a parameterized electromagnetic potential by the processor system extremizing an action integral within a predetermined tolerance for where the action integral is an integral of an electromagnetic Lagrangian density, the extremizing being performed by the machine at least numerically solving a set of equations for the set of parameters, the set of equations being equations that result from taking partial derivatives of the action integral with respect to the set of parameters.

20. A machine implemented method comprising: computing, by a machine having a processor system including one or more processors and a memory system, numerical values for a set of parameters of a parameterized effectively continuous electromagnetic potential;

the computing being performed, by the processor system, by at least extremizing an action

integral within a predetermined tolerance for which a Lagrangian density contains terms for at least an electromagnetic field;

the extremizing being performed, by the machine, by at least numerically solving a set of

equations for the set of parameters; the set of equations being a set of equations that result from taking partial derivatives of the action integral with respect to the set of parameters; and

storing, at least temporarily, results of the extremizing in the memory system.

21. The machine implemented method of any of claims 1-17 and 20, the potential being a linear function of the parameters, the linear function including at least a sum of products, each product being one parameter of the set of parameters times a potential basis function of a set of basis functions.

22. The machine implemented method of any of claims 1-17, 20, and 21, the electromagnetic displacement field and magnetic field of the Lagrangian density are equal to or approximated by linear functions

of the electric field, which is defined as a derivative of the potential, and

of the magnetic intensity, which is defined as a different derivative of the potential,

respectively.

23. The machine implemented method of any of claims 1-17, and 20-22, in which the

electromagnetic displacement and magnetic fields may be nonlinear functions of the electric field and magnetic intensity, resulting in a set of nonlinear equations in the parameters; and the numerical solving of the set of equations for the set of parameters of the potential.

24. The machine implemented method of claim of any of claims 1-17 and 20-23, the extremizing of the action being performed by solving a set of equations, each equation containing integrals over non-boundary points of a set of regions and integrals over boundary points between regions of the set of regions.

25. The machine implemented method of any of claims 1-17 and 20-24, the action being an integral whose integrand is a sum of at least

(1) a product of at least

(a) an electromagnetic field basis function and

a displacement-magnetic intensity field function and (2) a product of at least

(a) a potential basis function and

(b) a charge-current density function.

26. The machine implemented method of any of claims 1-17 and 20-25, the action integral having a region of integration that includes at least a polarizable medium, with polarizability that is either electric or magnetic or both; basis potential functions of the polarizable medium being products of a harmonic function of time and a function of spatial coordinates.

27. The machine implemented method of any of claims 1-17 and 20-26, the action integral having a region of integration that includes at least a polarizable medium; the solving being performed in a stationary frame of the polarizable medium.

28. The machine implemented method of any of claims 1-17 and 20-27, the parameterized potential being chosen so that for any parameter values, volume integrals over an interior of at least one region being equal to zero within a predetermined acceptable tolerance.

29. The machine implemented method of any of claims 1-17 and 20-28, the solving of the set of equations including solving the set of equations in a plurality of regions, the potential being a combination of one or more basis potentials, with at least one basis function being defined in at least one region as a potential function of at least one or more virtual charge-current sources, the virtual charge-current source being located outside the region or of real charge-current sources located inside the region or of both.

30. The machine implemented method of any of claims 1, 3-17, and 20-29, the solving of the set of equations including solving the set of equations in a plurality of regions, such that one or more regions are simplices and the parameterized potential in one or more simplices exactly or approximately satisfies the Euler-Lagrange equations.

31. The machine implemented method of any of claims 1-17 and 20-30, one or more of the potential functions being a linear function of position.

32. The machine implemented method of claim of any of claims 1-17 and 20-31, with time variations of electromagnetic potentials, fields, and charge and current densities approximated as varying harmonically with time, and electrical conductivity represented by complex permittivity.

33. The machine implemented method of any of claims 1-17 and 20-32, the potential being

parameterized as a sum of parameters times basis functions, the basis potential functions including two types of basis functions, one type of basis function being zero within an acceptable tolerance at all boundary points for which an electrical conductor is on at least one side, and another type of basis function being zero within an acceptable tolerance at all boundary points for which an electrical conductor is not on any side.

34. The machine implemented method of any of claims 1-17 and 20-33,

indexed that equations are solved are equations that result from a variation of the action being equal to zero, those indexed equations whose index corresponds to the first type of basis function, are equations resulting from multiplying by a factor equal to the frequency, and then symbolically simplifying to eliminate factors that diverge for frequencies that go to zero, before computing a numerical solution to the set of equations.

35. The machine implemented method of any of claims 1-17 and 20-34, in which each basis potential function is chosen so that its effective discontinuity at boundary points is minimized though not necessarily zero.

36. The machine implemented method of any of claims 1-17 and 20-35, in which space is made up of regions, each having uniform polarizability, separated by boundaries, and the basis functions are chosen so that they are continuous at all nonboundary points and effectively continuous at all boundary points for any parameter values, and the Euler- Lagrange equations are satisfied at all non-boundary points.

37. The machine implemented method of any of claims 1-14 and 20-36,

the solving including solving the set of equations

where / indicates to perform an integration,

Ω is a four dimensional volume having three spatial dimensions and one time dimension over which the integration is performed, F is the electromagnetic field, G is the displacement- magnetic intensity field, A is the potential function, and J is the current-charge field and c_£ are the parameters of the potential A such that A ~∑_£ a_£c_£, the field F such that F ~ ∑_£ fi Ci, and the field G such that G ~∑_£ ,g_£c_£, where ¾, f, and g; are the basis functions of A, F, and G, respectively, such that a_£ is a continuous, but not necessarily

differentiable, function of space-time position, /_£ is a function of space-time position given by /_£ = V Λ a_£, and g is a linear function of the field value _£ .

38. The machine implemented method of any of claims 1-14 and 20-37

M - C = B

for an array C with elements equal to the unknowns Cj for j≠ 0, where M is a matrix with elements m_£ for indices i ≠ 0 and j≠ 0, and B is an array with elements b_£ for indices i ≠ 0, given in terms of 4D space-time quantities by

mi,j = / ^ai ^" (7 ^" 3j)^dv + / ^ai ^" ( ^{" A}gj)^dv

= - / a-i ^■ (V^■ g₀ -/) dv - / a_£ ^■ (rj^■ Ag₀ - K) άη

in which the 4D volume integrals may equal zero or be negligible for some or all elements, in which case the associated elements are given by only boundary integrals in 4D space-time,

39. The machine implemented method of any of claims 1-14 and 20-38, wherein equations solved are equations of that extremize the action integral that have been transformed to equations having time dependence that is harmonic at frequency ω

and matrix and array elements of the equations solved are given in terms of 3D space quantities by mi ⁼ j (ψί* {^ν ' ^dj) - _£ ^{* ■} (y x hj + i⁽i>dj) dV

+ j (φί^* (η^■ Adj) - a_£ ^{* ■} (n x ΔΛ_;)) dn m = -j(≠,'(v d - p-) - a,' - (v x k ₊ i_aae -n)iv

- / (*,' (*^■ Ad₀ - a) - a - (n , Ah₀ - K))_dn in which the 3D volume integrals may equal zero or below a predetermined tolerance for some or all elements, elements for which in which case elements having a volume integral that equals zero or is below a predetermined tolerance are given by only boundary integrals in 3D space, ^mi,j ⁼ / Φί^* (η ^{' Δίί};) ^{~~ a}i * ^{' x ΔΛ} )) dn

b_t = - /(φ_£ ^* (η^■ Δίί₀ - σ) - a_t ^*■ (n x ΔΛ₀ - K))dn

where

∑j_≠0 indicates to perform a summation over all values of indices j except j = 0,

/ indicates to perform an integration,

<t>i is the ith scalar potential basis function

n is a unit vector normal to a boundary of computation cell

Adj is a difference between the jth displacement basis function on a first side of the boundary and a jth displacement basis function on a second side of the boundary, the jth displacement basis function on the first side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the first side of the boundary, and the jth displacement basis function on the second side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the second side of the boundary,

CLi is an ith vector potential basis function,

Ahj is a difference between the jth magnetic intensity basis function on the first side of the boundary and the jth magnetic intensity basis function on the second side of the boundary,

Adj is a difference between the jth displacement basis function on one side of the boundary and the jth displacement basis function on another side of the boundary, the jth magnetic intensity basis function on the first side of the boundary being a magnetic intensity vector resulting from the jth vector potential basis function on the first side of the boundary, and the jth magnetic intensity basis function on the second side of the boundary being a magnetic intensity vector resulting from the jth vector potential basis function on the second side of the boundary,

dV indicates that the integration is over the boundary of the computation cell;

Cj are the unknown parameters that are being solved for,

Ad₀ is a difference between the 0th displacement basis function on one side of the boundary and the 0th displacement basis function on another side of the boundary, the 0^th displacement basis function is a displacement function with a coefficient that has a known value, Ah₀ is a difference between the Oth magnetic intensity basis function on one side of the boundary and the Oth magnetic intensity basis function on another side of the boundary, the Oth magnetic intensity basis function is a basis function with a known coefficient p,J, σ, K are the known volume charge density, volume current density, surface charge density on a boundary, and surface current density on a boundary, respectively

c₀ is the coefficient that has the known value and all time dependence is harmonic and

represented by complex value.

40. A system comprising:

a processor system including one or more processors; and

a memory system including one or more non-transient storage media coupled to the processor system;

the memory system storing one or more machine instructions, which when invoked cause the processor system to implement a method comprising:

computing, by the processor system, numerical values for a set of parameters of a parameterized electromagnetic potential;

the computing being performed, by the processor system, by at least extremizing an action

integral within a predetermined tolerance for which a Lagrangian density contains terms for at least an electromagnetic field;

the extremizing being performed, by the machine, by at least numerically solving a set of

equations for the set of parameters; the set of equations being a set of equations that result from taking partial derivatives of the action integral with respect to the set of parameters; and

storing, at least temporarily, results of the extremizing in the memory system.

41. The system of claims 18 or 40,

the parameterized electromagnetic potential being a linear function of the parameters, the linear function including at least a sum of products, each product being one parameter of the set of parameters times a potential basis function of a set of basis functions.

42. The system of claims 18, 40, or 41,

an electromagnetic displacement field and magnetic field of the Lagrangian density are linear functions or approximated by linear functions of the electric field and magnetic field, the electric field is computed by taking a first combination of one or more derivatives of the parameterized potential, and

the magnetic intensity is computed by taking second combination of one or more derivatives of the parameterized potential.

43. The system of any of claims 18 and 40-42,

presenting to a user an input interface for entering a functional relationship between an electromagnetic displacement and a magnetic field that are nonlinear functions of an electric field and a magnetic intensity, resulting in a set of nonlinear equations in the parameters.

44. The system of any of claims 18 and 40-43,

the extremizing of the action being performed by solving a set of equations, each equation containing integrals over non-boundary points of a set of regions and integrals over boundary points between regions of the set of regions.

45. The system of any of claims 18 and 40-43,, the action being an integral whose integrand is a sum of at least

(1) a product of at least

(a) an electromagnetic field basis function and

(b) a displacement-magnetic intensity field function and

(2) a product of at least

(a) a potential basis function and

(b) a charge-current density function.

46. The system of any of claims 18 and 40-45, the action integral having a region of integration that includes at least a polarizable medium, with a polarizability that is either electric or magnetic or both; basis potential functions of the polarizable medium being products of a harmonic function of time and a function of spatial coordinates.

47. The system of any of claims 18 and 40-46, the action integral having a region of integration that includes at least a polarizable medium; the solving being performed in a rest frame of the polarizable medium.

48. The system of any of claims 18 and 40-47, the parameterized potential being chosen so that for any parameter values, volume integrals over an interior of at least one region are equal to zero within a predetermined acceptable tolerance.

49. The system of any of claims 18 and 40-48, the solving of the set of equations including solving the set of equations in a plurality of regions, the potential being a combination of one or more basis potentials, with at least one basis function being defined in at least one region as a potential function of at least one or more virtual charge-current sources, the virtual charge- current source being located outside the region or of real charge-current sources located inside the region or of both.

50. The system of any of claims 18 and 40-49, the solving of the set of equations including solving the set of equations in a plurality of regions, such that one or more regions are simplices and the parameterized potential is chosen so that within an interior of one or more simplices the parameterized potential exactly or approximately satisfies the Euler-Lagrange equations regardless of the parameters chosen for the parameterized potential.

51. The system of any of claims 18 and 40-50, one or more of the potential functions being a linear function of position.

52. The system of any of claims 18 and 40-51, time variations of electromagnetic potentials, fields, and charge and current densities being approximated as varying harmonically with time, and electrical conductivity represented by complex permittivity.

53. The system of any of claims 18 and 40-52, with the potential parameterized as a sum of parameters times basis functions, the basis potential functions including at least two types of basis functions, one type of basis function being zero within an acceptable tolerance at all boundary points for which an electrical conductor is on at least one side of the boundary.

54. The system of claim 53,

the set of equations solved being equations that result from multiplying indexed equations that express variation of the action being equal to zero, whose index corresponds to the first type of basis function, by a factor equal to the frequency, and then symbolically simplify the indexed equations to eliminate factors that diverge for frequencies that go to zero.

55. The system of any of claims 18 and 40-54,

each basis potential function is chosen so that an effective discontinuity in the basis potential at boundary points is minimized though not necessarily zero.

56. The system of any of claims 18 and 40-55

space in which the computing is performed is made up of regions, each region having a uniform polarizability, the regions being separated by boundaries, and the basis functions are chosen so that the basis functions are continuous at all nonboundary points and effectively continuous at all boundary points for any parameter values, and the Euler-Lagrange equations are satisfied at all non-boundary points.

57. The system of any of claims 18 and 40-65 the solving including solving a set of equations including

-^ J ^j - G — A^■ άΩ = 0 , where / indicates to perform an integration,

Ω is a four dimensional volume having three spatial dimensions and one time dimension over which the integration is performed, F is the electromagnetic field, G is the displacement- magnetic intensity field, A is the potential function, and J is the current-charge field and Ci are the parameters of the potential A such that A ~∑_{£ t}Ci, the field F such that F ~∑_£ iC_{i ?} and the field G such that G ~∑_£ giC where ¾, f, and g; are the basis functions of A, F, and G, respectively, such that a_t is a continuous, but not necessarily differentiable, function of space- time position, f_t is a function of space-time position given by f_t = V A a and g_t is a linear function of the field value f_t.

58. The system of any of claims 18 and 40-57,

the set of equations solved include a matrix equation

M C = B

for an array C with elements equal to the unknowns c - for j≠ 0, where M is a matrix with elements m_i - for indices i≠ 0 and j≠ 0, and B is an array with elements bi for indices i≠ 0, given in terms of 4D space-time quantities by

mij = / «i ^■ (7^■ 9j)dv + Aflf )diy

bi = - j _t ^■ (V^■ g₀ - J) dv - j a_t ^■ (η^■ Ag₀ - K) άη

in which the 4D volume integrals may equal zero or be negligible for some or all elements, in which case the associated elements are given by only boundary integrals in 4D space-time, bi = - J di^■ (η^■ Ag₀

59. The system of any of claims 18 and 40-58, wherein time dependence of the

electromagnetic potential is harmonic, having frequency o)

and the solving including solving a matrix equations having array elements, which when given in terms of 3D space quantities are rrii _j = j (φ_ί ^*(ν^■ d_j) - a_t ^*■ (v x h_j + ioid_j†) dV

+ j (φ_ι ^*(η^■ Ad_j) - di^*■ (n x ΔΛ )) dn bi = - !{φ_ί ^*(ν - ά₀ - ρ) - α_ί ^*■ (Γ χ Λο + ϊωίίο - J))dV - j (φ^ (n^■ Ad₀ - σ) - ai nx ΔνδΟ Kdn. in which the 3D volume integrals are equal to zero or below a predetermined threshold for some or all elements, in which case the matrix elements associated with volume integrals that are equal to zero or below a predetermined threshold are given by only boundary integrals in 3D space,

mi,_j ⁼ / (φϊ^* (^η ' d_j) - ai^*■ (n x ΔΛ )) dn

b_j = - J(0_j ^*(n^■ Δίί₀ - σ) - a_t ^*■ (n x ΔΛ₀ - K))dn

where

∑j_≠o indicates to perform a summation over all values of indices j except j = 0,

/ indicates to perform an integration,

<2 i is an ith scalar potential basis function

n is a unit vector normal to a boundary of computation cell

Ad_j is a difference between a jth displacement basis function on a first side of a boundary and a jth displacement basis function on a second side of the boundary, the jth displacement basis function on the first side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the first side of the boundary, and the jth displacement basis function on the second side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the second side of the boundary,

CLi is an ith vector potential basis function,

Ah_j is a difference between a jth magnetic intensity basis function on the first side of the boundary and a jth magnetic intensity basis function on the second side of the boundary, dV indicates that the integration is over the boundary of the computation cell; Cj are the unknown parameters that are being solved for,

Ad₀ is a difference between a 0th displacement basis function on one side of the boundary and a 0th displacement basis function on another side of the boundary, the 0^th displacement basis function is a displacement function with a coefficient that has a known value,

Ah₀ is a difference between a 0th magnetic intensity basis function on one side of the boundary and a 0th magnetic intensity basis function on another side of the boundary, the 0th magnetic intensity basis function is a basis function with a known coefficient,

p,J, σ, K are a known volume charge density, volume current density, surface charge density on a boundary, and surface current density on the boundary, respectively

c₀ is the coefficient that has the known value and all time dependence is harmonic and represented by complex values.

60. A non-transitory computer readable medium storing one or more machine instruction, which when implemented by a processor causes the processor to perform a method comprising:

computing, by a machine having a processor system including one or more processors and a memory system, numerical values for a set of parameters of a parameterized effectively continuous electromagnetic potential;

the computing being performed, by the processor system, by at least extremizing an action

integral within a predetermined tolerance for which a Lagrangian density contains terms for at least an electromagnetic field;

the extremizing being performed, by the machine, by at least numerically solving a set of

equations for the set of parameters; the set of equations being a set of equations that result from taking partial derivatives of the action integral with respect to the set of parameters; and

storing, at least temporarily, results of the extremizing in the memory system.

61. A machine implemented method comprising:

iteratively solving, by a machine having a processor system including one or more processors and a memory system including a memory system having at least one or more nontransient-machine readable media, a set of equations, each equation being a result of taking a derivative of an action integral, which is an integral of a Langrangian with respect to time, the derivate of the action integral being a derivative with respect to a parameter of a Lagrangian density associated with the Lagrangian, the Lagrangian density has a dependence on the parameters, each parameter of the Lagrangian density being part of a set of parameters of a parameterized potential function, the dependence of the Lagrangian density on the set of parameters being that of a representation of the Lagrangian density in terms of the potential function;

prior to the iteratively solving, the set of parameters include parameters that have values that are not known;

the iteratively solving of the set of equations determines solution values for the parameters; computing, by the processor system, a value representative of an error associated with a current iteration of the iteratively solving; and

if the value representative of the error for the current iteration is on an opposite side of a

predetermined threshold compared to the value representative of the error for a prior iteration, ending the iteratively solving.

62. The machine implemented method of any of claims 1-17, 20-39, and 60-61, the potential being a linear function of the parameters, the linear function includes at least a sum of products, each product being one parameter of the set of parameters multiplied by a potential basis function of a set of basis functions.

63. The machine implemented method of any of claims 1-17, 20-39, and 60-62, the

electromagnetic displacement field and magnetic field of the Lagrangian density are equal to or approximated by a linear function of the potential function.

64. The machine implemented method of any of claims 1-17, 20-39, and 60-63, each equation is expressed as a sum of integrals over regions and integrals over boundaries of the regions.

65. The machine implemented method of any of claims 1-17, 20-39, and 60-64, a transformed version of equations representing the action integral is solved that has no express reference to time.

66. The machine implemented method of any of claims 1-17, 20-39, and 60-65, time variations of electromagnetic potentials, fields, and charge and current densities are approximated as varying harmonically with time; and

electrical conductivity is represented as having a complex permittivity.

67. The machine implemented method of claim 65, in the transformed equations the potential is parameterized as a linear function of the parameters and equal to a sum of parameters multiplied by basis functions plus one non-parameterized basis function, the set of basis functions including two types of basis functions

each of a first the two types of basis functions being chosen to be zero at all boundary points of any boundary including an electrical conductor where the basis function of the first type is defined, or

each of a second the two types of basis functions being chosen to be zero at all boundary points of a boundary not including an electrical conductor basis function of the second type is defined basis function of the first type are multiplied by the frequency before computing a numerical solution to the equations, to cancel factors of frequency in numerators.

68. The machine implemented method of any of claims 1-17, 20-39, and 60-67, the

parameterized potential is chosen so that for any parameter values, one or more of the volume integrals are equal to zero or have a magnitude that is less than a predetermined threshold, resulting in correspondingly reduced computational demands.

69. The machine implemented method of any of claims 1-17, 20-39, and 60-68, one or more basis functions being defined and calculated within one or more regions as a solution to Maxwell's equations in a fictitious auxiliary space using a fictitious virtual charge-current sources.

70. The machine implemented method of any of claims 1-17, 20-39, and 60-69, one or more regions are simplices and the parameterized potential in each simplex is a polynomial function of position.

71. The machine implemented method of any of claims 1-17, 20-39, and 60-70, the polynomial function of position is a linear function of position.

72. The machine implemented method of any of claims 1-17, 20-39, and 60-71, one of the two basis functions being a basis function of a field and another of the two basis functions being a basis function of potential, the integration including at least integrating a product of the basis function of the field and the basis function of the potential.

73. The machine implemented method of any of claims 1-17, 20-39, and 60-72, the basis function of the potential extremizing of the action, within a predetermined tolerance, in an interior portion of a computation cell regardless of what value is assigned to the unknown parameters.

74. (Original) The machine implemented method of any of claims 1-17, 20-39, and 60-73, the field being the displacement field, and the method comprising multiplying basis functions of the displacement field, that are associated with computation cells including a conductor, by a variable representing an angular frequency, while not multiplying basis functions of the displaced field, that are associated with computation cells not having a conductor, by the variable representing the angular frequency.

75. The machine implemented method of any of claims 1-17, 20-39, and 60-74, in a first

computation cell having a current, placing a fictitious charge on one or more boundaries of the computation cell sufficient to create the current; and

on each boundary having the fictitious charge, placing another fictitious charge of opposite polarity in an adjacent computation cell that cancels the charge of the boundary in the first computation cell.

76. The machine implemented method of any of claims 1-17, 20-39, and 60-75, choosing a first set of basis functions that are nonzero at boundaries between computation cells that have a conductor and are zero at boundaries not having conductors; and

choosing a second set of basis functions that are nonzero at boundaries between computation cells that do not having a conductor are zero at boundaries having conductors, therein decoupling computations of the unknown parameters that are coefficients of basis functions that are nonzero at boundaries that have conductors and basis functions that are nonzero at boundaries that do not have conductors.

77. The machine implemented method of any of claims 1-17, 20-39, and 60-76, further

comprising:

computing a first difference, which is a difference of at least scalar potential associated with an interior of a computation cell evaluated at a boundary between two computation cells and a scalar potential associated with the boundary evaluated at the boundary;

squaring the first different;

computing a second difference, which is a difference of at least vector potential associated with an interior of a computation cell evaluated at a boundary between two computation cells and a vector potential associated with the boundary evaluated at the boundary;

squaring the second difference;

computing a component of the second difference that is perpendicular to the boundary;

squaring the component;

subtracting results of the squaring of the component from results of the squaring of the second difference;

combining results of the subtracting and the results of squaring the first difference;

presenting a value resulting from the combining of the results as a measure of an accuracy of the solution values.

78. The machine implemented method of any of claims 1-17, 20-39, and 60-77, the solving being performed by iteratively computing new values of the unknown parameters that are expected to be values of the unknown parameters that result in a new value of the value representative of the error that represents less error than in a prior iteration, the method further comprising:

if the value resulting from the combining is less than a predetermined tolerance, terminating solving; and

if the value resulting from the combining is a value that is greater than the predetermined

tolerance, performing at least another iterative improvement of the unknown parameters.

79. The machine implemented method of any of claims 1-17, 20-39, and 60-78,

the solving including solving

∑_j≠0 / ( _ifi - Ad_j - a_i - fi x Ah_j)dVc_j = - / (Φ^ίϊ^■ Δά₀ - σ) - a_; ^■ (n x Ah₀ - K))dVc₀,

where

∑_j≠0 indicates to perform a summation over all values of indices j except j = 0,

/ indicates to perform an integration,

Φ; is an ith scalar potential basis function

n is a unit vector normal to a boundary of computation cell

Ad_j is a difference between a jth displacement basis function on a first side of the boundary and a jth displacement basis function on a second side of the boundary, the jth displacement basis function on the first side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the first side of the boundary, and the jth displacement basis function on the second side of the boundary being a displacement vector resulting from the jth scalar potential basis function and vector potential basis function on the second side of the boundary, a; is an ith vector potential basis function,

Ah_j is a difference between a jth magnetic intensity basis function on the first side of the

boundary and a jth magnetic intensity basis function on the second side of the boundary, dV indicates that the integration is over the volume of the computation cell;

Cj are the unknown parameters that are being solved for,

Δά₀ is a difference between a 0th displacement basis function on one side of the boundary and the 0th displacement basis function on another side of the boundary, the 0^th displacement basis function is a displacement function with a coefficient that has a known value,

A ₀ is a difference between a 0th magnetic intensity basis function on one side of the boundary and a 0th magnetic intensity basis function on another side of the boundary, the 0th magnetic intensity basis function is a basis function with a known coefficient,

K is the surface charge density,

c₀ is the coefficient that has the known value;

an elements of the array having indices i and j, being equal to / ( _ίΠ^■ Ad_j— a^* _j ^■ n x Ah_j)dV; and

an ith element of the result values being— / (Φ_ί(ίϊ^■ Δά₀— σ)— ¾^■ (n x Δη₀— K))dVc₀.

80. A computer readable medium storing thereon one or more machine instructions, which when implemented, causes the processor system to carry out the method of any of claims 1-17, 20-39, and 60-79.