Spherical coordinate system

From Wikipedia, the free encyclopedia

https://en.wikipedia.org/wiki/Spherical_coordinate_system

The physics convention. Spherical coordinates (r, θ, φ) as commonly used: (ISO 80000-2:2019): radial distance r (slant distance to origin), polar angle θ (theta) (angle with respect to positive polar axis), and azimuthal angle φ (phi) (angle of rotation from the initial meridian plane). This is the convention followed in this article.

In mathematics, a spherical coordinate system is a coordinate system for three-dimensional space where the position of a given point in space is specified by three numbers, (r, θ, φ): the radial distance of the radial line r connecting the point to the fixed point of origin (which is located on a fixed polar axis, or zenith direction axis, or z-axis); the polar angle θ of the radial line r; and the azimuthal angle φ of the radial line r.

The polar angle θ is measured between the z-axis and the radial line r. The azimuthal angle φ is measured between the orthogonal projection of the radial line r onto the reference x-y-plane—which is orthogonal to the z-axis and passes through the fixed point of origin—and either of the fixed x-axis or y-axis, both of which are orthogonal to the z-axis and to each other. (See graphic re the "physics convention".)

Once the radius is fixed, the three coordinates (r, θ, φ), known as a 3-tuple, provide a coordinate system on a sphere, typically called the spherical polar coordinates. Nota bene: the physics convention is followed in this article; (See both graphics re "physics convention" and re "mathematics convention").

The radial distance from the fixed point of origin is also called the radius, or radial line, or radial coordinate. The polar angle may be called inclination angle, zenith angle, normal angle, or the colatitude. The user may choose to ignore the inclination angle and use the elevation angle instead, which is measured upward between the reference plane and the radial line—i.e., from the reference plane upward (towards to the positive z-axis) to the radial line. The depression angle is the negative of the elevation angle. (See graphic re the "physics convention"—not "mathematics convention".)

Both the use of symbols and the naming order of tuple coordinates differ among the several sources and disciplines. This article will use the ISO convention frequently encountered in physics, where the naming tuple gives the order as: radial distance, polar angle, azimuthal angle, or ${\displaystyle (r,\theta ,\phi )}$. (See graphic re the "physics convention".) In contrast, the conventions in many mathematics books and texts give the naming order differently as: radial distance, "azimuthal angle", "polar angle", and ${\displaystyle (\rho ,\theta ,\phi )}$ or ${\displaystyle (r,\theta ,\phi )}$—which switches the uses and meanings of symbols θ and φ. Other conventions may also be used, such as r for a radius from the z-axis that is not from the point of origin. Particular care must be taken to check the meaning of the symbols.

The mathematics convention. Spherical coordinates (r, θ, φ) as typically used: radial distance r, azimuthal angle θ, and polar angle φ. + The meanings of θ and φ have been swapped—compared to the physics convention. The 'south'-direction x-axis is depicted but the 'north'-direction x-axis is not. (As in physics, ρ (rho) is often used instead of r to avoid confusion with the value r in cylindrical and 2D polar coordinates.)

According to the conventions of geographical coordinate systems, positions are measured by latitude, longitude, and height (altitude). There are a number of celestial coordinate systems based on different fundamental planes and with different terms for the various coordinates. The spherical coordinate systems used in mathematics normally use radians rather than degrees; (note 90 degrees equals π_/2 radians). And these systems of the mathematics convention may measure the azimuthal angle counterclockwise (i.e., from the south direction x-axis, or 180°, towards the east direction y-axis, or +90°)—rather than measure clockwise (i.e., from the north direction x-axis, or 0°, towards the east direction y-axis, or +90°), as done in the horizontal coordinate system.

The spherical coordinate system of the physics convention can be seen as a generalization of the polar coordinate system in three-dimensional space. It can be further extended to higher-dimensional spaces, and is then referred to as a hyperspherical coordinate system.

Definition

To define a spherical coordinate system, one must designate an origin point in space, O, and two orthogonal directions: the zenith reference direction and the azimuth reference direction. These choices determine a reference plane that is typically defined as containing the point of origin and the x– and y–axes, either of which may be designated as the azimuth reference direction. The reference plane is perpendicular (orthogonal) to the zenith direction, and typically is desiginated "horizontal" to the zenith direction's "vertical". The spherical coordinates of a point P then are defined as follows:

• The radius or radial distance is the Euclidean distance from the origin O to P.
• The inclination (or polar angle) is the signed angle from the zenith reference direction to the line segment OP. (Elevation may be used as the polar angle instead of inclination; see below.)
• The azimuth (or azimuthal angle) is the signed angle measured from the azimuth reference direction to the orthogonal projection of the radial line segment OP on the reference plane.

The sign of the azimuth is determined by designating the rotation that is the positive sense of turning about the zenith. This choice is arbitrary, and is part of the coordinate system definition. (If the inclination is either zero or 180 degrees (= π radians), the azimuth is arbitrary. If the radius is zero, both azimuth and inclination are arbitrary.)

The elevation is the signed angle from the x-y reference plane to the radial line segment OP, where positive angles are designated as upward, towards the zenith reference. Elevation is 90 degrees (= π/2 radians) minus inclination. Thus, if the inclination is 60 degrees (= π/3 radians), then the elevation is 30 degrees (= π/6 radians).

In linear algebra, the vector from the origin O to the point P is often called the position vector of P.

Conventions

Several different conventions exist for representing spherical coordinates and prescribing the naming order of their symbols. The 3-tuple number set ${\displaystyle (r,\theta ,\phi )}$ denotes radial distance, the polar angle—"inclination", or as the alternative, "elevation"—and the azimuthal angle. It is the common practice within the physics convention, as specified by ISO standard 80000-2:2019, and earlier in ISO 31-11 (1992).

As stated above, this article describes the ISO "physics convention"—unless otherwise noted.

However, some authors (including mathematicians) use the symbol ρ (rho) for radius, or radial distance, φ for inclination (or elevation) and θ for azimuth—while others keep the use of r for the radius; all which "provides a logical extension of the usual polar coordinates notation". As to order, some authors list the azimuth before the inclination (or the elevation) angle. Some combinations of these choices result in a left-handed coordinate system. The standard "physics convention" 3-tuple set ${\displaystyle (r,\theta ,\phi )}$ conflicts with the usual notation for two-dimensional polar coordinates and three-dimensional cylindrical coordinates, where θ is often used for the azimuth.

Angles are typically measured in degrees (°) or in radians (rad), where 360° = 2π rad. The use of degrees is most common in geography, astronomy, and engineering, where radians are commonly used in mathematics and theoretical physics. The unit for radial distance is usually determined by the context, as occurs in applications of the 'unit sphere', see applications.

When the system is used to designate physical three-space, it is customary to assign positive to azimuth angles measured in the counterclockwise sense from the reference direction on the reference plane—as seen from the "zenith" side of the plane. This convention is used in particular for geographical coordinates, where the "zenith" direction is north and the positive azimuth (longitude) angles are measured eastwards from some prime meridian.

Major conventions

coordinates set order	corresponding local geographical directions (Z, X, Y)	right/left-handed
(r, θinc, φaz,right)	(U, S, E)	right
(r, φaz,right, θel)	(U, E, N)	right
(r, θel, φaz,right)	(U, N, E)	left

Note: Easting (E), Northing (N), Upwardness (U). In the case of (U, S, E) the local azimuth angle would be measured counterclockwise from S to E.

Unique coordinates

Any spherical coordinate triplet (or tuple) ${\displaystyle (r,\theta ,\phi )}$ specifies a single point of three-dimensional space. On the reverse view, any single point has infinitely many equivalent spherical coordinates. That is, the user can add or subtract any number of full turns to the angular measures without changing the angles themselves, and therefore without changing the point. It is convenient in many contexts to use negative radial distances, the convention being ${\displaystyle (-r,-\theta ,\phi )}$, which is equivalent to ${\displaystyle (r,\theta ,\phi )}$ for any r, θ, and φ. Moreover, ${\displaystyle (r,-\theta ,\phi )}$ is equivalent to ${\displaystyle (r,\theta ,\phi +{180}^{\circ })}$.

When necessary to define a unique set of spherical coordinates for each point, the user must restrict the range, aka interval, of each coordinate. A common choice is:

• radial distance: r ≥ 0,
• polar angle: 0° ≤ θ ≤ 180°, or 0 rad ≤ θ ≤ π rad,
• azimuth : 0° ≤ φ < 360°, or 0 rad ≤ φ < 2π rad.

But instead of the interval [0°, 360°), the azimuth φ is typically restricted to the half-open interval (−180°, +180°], or (−π, +π ] radians, which is the standard convention for geographic longitude.

For the polar angle θ, the range (interval) for inclination is [0°, 180°], which is equivalent to elevation range (interval) [−90°, +90°]. In geography, the latitude is the elevation.

Even with these restrictions, if the polar angle (inclination) is 0° or 180°—elevation is −90° or +90°—then the azimuth angle is arbitrary; and if r is zero, both azimuth and polar angles are arbitrary. To define the coordinates as unique, the user can assert the convention that (in these cases) the arbitrary coordinates are set to zero.

Plotting

To plot any dot from its spherical coordinates (r, θ, φ), where θ is inclination, the user would: move r units from the origin in the zenith reference direction (z-axis); then rotate by the amount of the azimuth angle (φ) about the origin from the designated azimuth reference direction, (i.e., either the x– or y–axis, see Definition, above); and then rotate from the z-axis by the amount of the θ angle.

Applications

In the mathematics convention: A globe showing a unit sphere, with tuple coordinates of point P (red): its radial distance r (red, not labeled); its azimuthal angle θ (not labeled); and its polar angle of inclination φ (not labeled). The radial distance upward along the zenith–axis from the point of origin to the surface of the sphere is assigned the value unity, or 1. + In this image, r appears to equal 4/6, or .67, (of unity); i.e., four of the six 'nested shells' to the surface. The azimuth angle θ appears to equal positive 90°, as rotated counterclockwise from the azimuth-reference x–axis; and the inclination φ appears to equal 30°, as rotated from the zenith–axis. (Note the 'full' rotation, or inclination, from the zenith–axis to the y–axis is 90°).

Just as the two-dimensional Cartesian coordinate system is useful—has a wide set of applications—on a planar surface, a two-dimensional spherical coordinate system is useful on the surface of a sphere. For example, one sphere that is described in Cartesian coordinates with the equation x² + y² + z² = c² can be described in spherical coordinates by the simple equation r = c. (In this system—shown here in the mathematics convention—the sphere is adapted as a unit sphere, where the radius is set to unity and then can generally be ignored, see graphic.)

This (unit sphere) simplification is also useful when dealing with objects such as rotational matrices. Spherical coordinates are also useful in analyzing systems that have some degree of symmetry about a point, including: volume integrals inside a sphere; the potential energy field surrounding a concentrated mass or charge; or global weather simulation in a planet's atmosphere.

The output pattern of the industrial loudspeaker shown here uses spherical polar plots taken at six frequencies

Three dimensional modeling of loudspeaker output patterns can be used to predict their performance. A number of polar plots are required, taken at a wide selection of frequencies, as the pattern changes greatly with frequency. Polar plots help to show that many loudspeakers tend toward omnidirectionality at lower frequencies.

An important application of spherical coordinates provides for the separation of variables in two partial differential equations—the Laplace and the Helmholtz equations—that arise in many physical problems.The angular portions of the solutions to such equations take the form of spherical harmonics. Another application is ergonomic design, where r is the arm length of a stationary person and the angles describe the direction of the arm as it reaches out. The spherical coordinate system is also commonly used in 3D game development to rotate the camera around the player's position.

In geography

Main article: Geographic coordinate system

See also: ECEF

Instead of inclination, the geographic coordinate system uses elevation angle (or latitude), in the range (aka domain) −90° ≤ φ ≤ 90° and rotated north from the equator plane. Latitude (i.e., the angle of latitude) may be either geocentric latitude, measured (rotated) from the Earth's center—and designated variously by ψ, q, φ′, φ_c, φ_g—or geodetic latitude, measured (rotated) from the observer's local vertical, and typically designated φ. The polar angle (inclination), which is 90° minus the latitude and ranges from 0 to 180°, is called colatitude in geography.

The azimuth angle (or longitude) of a given position on Earth, commonly denoted by λ, is measured in degrees east or west from some conventional reference meridian (most commonly the IERS Reference Meridian); thus its domain (or range) is −180° ≤ λ ≤ 180° and a given reading is typically designated "East" or "West". For positions on the Earth or other solid celestial body, the reference plane is usually taken to be the plane perpendicular to the axis of rotation.

Instead of the radial distance r geographers commonly use altitude above or below some local reference surface (vertical datum), which, for example, may be the mean sea level. When needed, the radial distance can be computed from the altitude by adding the radius of Earth, which is approximately 6,360 ± 11 km (3,952 ± 7 miles).

However, modern geographical coordinate systems are quite complex, and the positions implied by these simple formulae may be inaccurate by several kilometers. The precise standard meanings of latitude, longitude and altitude are currently defined by the World Geodetic System (WGS), and take into account the flattening of the Earth at the poles (about 21 km or 13 miles) and many other details.

Planetary coordinate systems use formulations analogous to the geographic coordinate system.

In astronomy

A series of astronomical coordinate systems are used to measure the elevation angle from several fundamental planes. These reference planes include: the observer's horizon, the galactic equator (defined by the rotation of the Milky Way), the celestial equator (defined by Earth's rotation), the plane of the ecliptic (defined by Earth's orbit around the Sun), and the plane of the earth terminator (normal to the instantaneous direction to the Sun).

Coordinate system conversions

See also: List of common coordinate transformations § To spherical coordinates

As the spherical coordinate system is only one of many three-dimensional coordinate systems, there exist equations for converting coordinates between the spherical coordinate system and others.

Cartesian coordinates

The spherical coordinates of a point in the ISO convention (i.e. for physics: radius r, inclination θ, azimuth φ) can be obtained from its Cartesian coordinates (x, y, z) by the formulae

$${\displaystyle \begin{array}{rl}r& =\sqrt{x^2+y^2+z^2}\\ \theta & =\arccos \frac{z}{\sqrt{x^2+y^2+z^2}}=\arccos \frac{z}{r}=\{\begin{array}{ll}\arctan \frac{\sqrt{x^2+y^2}}{z}& \text{if }z>0\\ \pi +\arctan \frac{\sqrt{x^2+y^2}}{z}& \text{if }z<0\\ +\frac{\pi }{2}& \text{if }z=0\text{ and }\sqrt{x^2+y^2}\ne 0\\ \text{undefined}& \text{if }x=y=z=0\end{array}\\ \phi & =\mathrm{sgn}(y)\arccos \frac{x}{\sqrt{x^2+y^2}}=\{\begin{array}{ll}\arctan (\frac{y}{x})& \text{if }x>0,\\ \arctan (\frac{y}{x})+\pi & \text{if }x<0\text{ and }y\ge 0,\\ \arctan (\frac{y}{x})-\pi & \text{if }x<0\text{ and }y<0,\\ +\frac{\pi }{2}& \text{if }x=0\text{ and }y>0,\\ -\frac{\pi }{2}& \text{if }x=0\text{ and }y<0,\\ \text{undefined}& \text{if }x=0\text{ and }y=0.\end{array}\end{array}}$$

The inverse tangent denoted in φ = arctan y/x must be suitably defined, taking into account the correct quadrant of (x, y). See the article on atan2.

Alternatively, the conversion can be considered as two sequential rectangular to polar conversions: the first in the Cartesian xy plane from (x, y) to (R, φ), where R is the projection of r onto the xy-plane, and the second in the Cartesian zR-plane from (z, R) to (r, θ). The correct quadrants for φ and θ are implied by the correctness of the planar rectangular to polar conversions.

These formulae assume that the two systems have the same origin, that the spherical reference plane is the Cartesian xy plane, that θ is inclination from the z direction, and that the azimuth angles are measured from the Cartesian x axis (so that the y axis has φ = +90°). If θ measures elevation from the reference plane instead of inclination from the zenith the arccos above becomes an arcsin, and the cos θ and sin θ below become switched.

Conversely, the Cartesian coordinates may be retrieved from the spherical coordinates (radius r, inclination θ, azimuth φ), where r ∈ [0, ∞), θ ∈ [0, π], φ ∈ [0, 2π), by

$${\displaystyle \begin{array}{rl}x& =r\sin \theta \cos \phi ,\\ y& =r\sin \theta \sin \phi ,\\ z& =r\cos \theta .\end{array}}$$

Cylindrical coordinates

Main article: Cylindrical coordinate system

Cylindrical coordinates (axial radius ρ, azimuth φ, elevation z) may be converted into spherical coordinates (central radius r, inclination θ, azimuth φ), by the formulas

$${\displaystyle \begin{array}{rl}r& =\sqrt{{\rho }^2+z^2},\\ \theta & =\arctan \frac{\rho }{z}=\arccos \frac{z}{\sqrt{{\rho }^2+z^2}},\\ \phi & =\phi .\end{array}}$$

Conversely, the spherical coordinates may be converted into cylindrical coordinates by the formulae

$${\displaystyle \begin{array}{rl}\rho & =r\sin \theta ,\\ \phi & =\phi ,\\ z& =r\cos \theta .\end{array}}$$

These formulae assume that the two systems have the same origin and same reference plane, measure the azimuth angle φ in the same senses from the same axis, and that the spherical angle θ is inclination from the cylindrical z axis.

Generalization

See also: Ellipsoidal coordinates

It is also possible to deal with ellipsoids in Cartesian coordinates by using a modified version of the spherical coordinates.

Let P be an ellipsoid specified by the level set

$${\displaystyle a{x}^2+b{y}^2+c{z}^2=d.}$$

The modified spherical coordinates of a point in P in the ISO convention (i.e. for physics: radius r, inclination θ, azimuth φ) can be obtained from its Cartesian coordinates (x, y, z) by the formulae

$${\displaystyle \begin{array}{rl}x& =\frac{1}{\sqrt{a}}r\sin \theta \cos \phi ,\\ y& =\frac{1}{\sqrt{b}}r\sin \theta \sin \phi ,\\ z& =\frac{1}{\sqrt{c}}r\cos \theta ,\\ r^2& =a{x}^2+b{y}^2+c{z}^2.\end{array}}$$

An infinitesimal volume element is given by

$${\displaystyle \mathrm{d}V=\left|\frac{\mathrm{\partial }(x,y,z)}{\mathrm{\partial }(r,\theta ,\phi )}\right|drd\theta d\phi =\frac{1}{\sqrt{abc}}{r}^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi =\frac{1}{\sqrt{abc}}{r}^2\mathrm{d}r\mathrm{d}\mathrm{\Omega }.}$$

The square-root factor comes from the property of the determinant that allows a constant to be pulled out from a column:

$${\displaystyle \left|\begin{array}{ccc}ka& b& c\\ kd& e& f\\ kg& h& i\end{array}\right|=k\left|\begin{array}{ccc}a& b& c\\ d& e& f\\ g& h& i\end{array}\right|.}$$

Integration and differentiation in spherical coordinates

Unit vectors in spherical coordinates

The following equations (Iyanaga 1977) assume that the colatitude θ is the inclination from the positive z axis, as in the physics convention discussed.

The line element for an infinitesimal displacement from (r, θ, φ) to (r + dr, θ + dθ, φ + dφ) is

$${\displaystyle \mathrm{d}\mathbf{r}=\mathrm{d}r\hat{\mathbf{r}}+r\mathrm{d}\theta \hat{\mathit{\theta }}+r\sin \theta \mathrm{d}\phi \hat{\mathit{\phi }},}$$

where

$${\displaystyle \begin{array}{rl}\hat{\mathbf{r}}& =\sin \theta \cos \phi \hat{\mathbf{x}}+\sin \theta \sin \phi \hat{\mathbf{y}}+\cos \theta \hat{\mathbf{z}},\\ \hat{\mathit{\theta }}& =\cos \theta \cos \phi \hat{\mathbf{x}}+\cos \theta \sin \phi \hat{\mathbf{y}}-\sin \theta \hat{\mathbf{z}},\\ \hat{\mathit{\phi }}& =-\sin \phi \hat{\mathbf{x}}+\cos \phi \hat{\mathbf{y}}\end{array}}$$

are the local orthogonal unit vectors in the directions of increasing r, θ, and φ, respectively, and x̂, ŷ, and ẑ are the unit vectors in Cartesian coordinates. The linear transformation to this right-handed coordinate triplet is a rotation matrix,

$${\displaystyle R=\left(\begin{array}{ccc}\sin \theta \cos \phi & \sin \theta \sin \phi & \phantom{}\cos \theta \\ \cos \theta \cos \phi & \cos \theta \sin \phi & -\sin \theta \\ -\sin \phi & \cos \phi & \phantom{}0\end{array}\right).}$$

This gives the transformation from the spherical to the cartesian, the other way around is given by its inverse. Note: the matrix is an orthogonal matrix, that is, its inverse is simply its transpose.

The Cartesian unit vectors are thus related to the spherical unit vectors by:

$${\displaystyle \left[\begin{array}{c}\hat{\mathbf{x}}\\ \hat{\mathbf{y}}\\ \hat{\mathbf{z}}\end{array}\right]=\left[\begin{array}{ccc}\sin \theta \cos \phi & \cos \theta \cos \phi & -\sin \phi \\ \sin \theta \sin \phi & \cos \theta \sin \phi & \phantom{}\cos \phi \\ \cos \theta & -\sin \theta & \phantom{}0\end{array}\right]\left[\begin{array}{c}\hat{\mathit{r}}\\ \hat{\mathit{\theta }}\\ \hat{\mathit{\phi }}\end{array}\right]}$$

The general form of the formula to prove the differential line element, is

$${\displaystyle \mathrm{d}\mathbf{r}=\sum _i\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }{x}_i}\mathrm{d}{x}_i=\sum _i\left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }{x}_i}\right|\frac{\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }{x}_i}}{\left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }{x}_i}\right|}\mathrm{d}{x}_i=\sum _i\left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }{x}_i}\right|\mathrm{d}{x}_i{\hat{\mathit{x}}}_i,}$$

that is, the change in ${\displaystyle \mathbf{r}}$ is decomposed into individual changes corresponding to changes in the individual coordinates.

To apply this to the present case, one needs to calculate how ${\displaystyle \mathbf{r}}$ changes with each of the coordinates. In the conventions used,

$${\displaystyle \mathbf{r}=\left[\begin{array}{c}r\sin \theta \cos \phi \\ r\sin \theta \sin \phi \\ r\cos \theta \end{array}\right].}$$

Thus,

$${\displaystyle \frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }r}=\left[\begin{array}{c}\sin \theta \cos \phi \\ \sin \theta \sin \phi \\ \cos \theta \end{array}\right]=\hat{\mathbf{r}},\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\theta }=\left[\begin{array}{c}r\cos \theta \cos \phi \\ r\cos \theta \sin \phi \\ -r\sin \theta \end{array}\right]=r\hat{\mathit{\theta }},\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\phi }=\left[\begin{array}{c}-r\sin \theta \sin \phi \\ \phantom{}r\sin \theta \cos \phi \\ 0\end{array}\right]=r\sin \theta \hat{\mathit{\phi }}.}$$

The desired coefficients are the magnitudes of these vectors:

$${\displaystyle \left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }r}\right|=1,\left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\theta }\right|=r,\left|\frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\phi }\right|=r\sin \theta .}$$

The surface element spanning from θ to θ + dθ and φ to φ + dφ on a spherical surface at (constant) radius r is then

$${\displaystyle \mathrm{d}{S}_r=\Vert \frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\theta }\times \frac{\mathrm{\partial }\mathbf{r}}{\mathrm{\partial }\phi }\Vert \mathrm{d}\theta \mathrm{d}\phi =\left|r\hat{\mathit{\theta }}\times r\sin \theta \hat{\mathit{\phi }}\right|\mathrm{d}\theta \mathrm{d}\phi =r^2\sin \theta \mathrm{d}\theta \mathrm{d}\phi .}$$

Thus the differential solid angle is

$${\displaystyle \mathrm{d}\mathrm{\Omega }=\frac{\mathrm{d}{S}_r}{r^2}=\sin \theta \mathrm{d}\theta \mathrm{d}\phi .}$$

The surface element in a surface of polar angle θ constant (a cone with vertex at the origin) is

$${\displaystyle \mathrm{d}{S}_{\theta }=r\sin \theta \mathrm{d}\phi \mathrm{d}r.}$$

The surface element in a surface of azimuth φ constant (a vertical half-plane) is

$${\displaystyle \mathrm{d}{S}_{\phi }=r\mathrm{d}r\mathrm{d}\theta .}$$

The volume element spanning from r to r + dr, θ to θ + dθ, and φ to φ + dφ is specified by the determinant of the Jacobian matrix of partial derivatives,

$${\displaystyle J=\frac{\mathrm{\partial }(x,y,z)}{\mathrm{\partial }(r,\theta ,\phi )}=\left(\begin{array}{ccc}\sin \theta \cos \phi & r\cos \theta \cos \phi & -r\sin \theta \sin \phi \\ \sin \theta \sin \phi & r\cos \theta \sin \phi & \phantom{}r\sin \theta \cos \phi \\ \cos \theta & -r\sin \theta & \phantom{}0\end{array}\right),}$$

namely

$${\displaystyle \mathrm{d}V=\left|\frac{\mathrm{\partial }(x,y,z)}{\mathrm{\partial }(r,\theta ,\phi )}\right|\mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi =r^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi =r^2\mathrm{d}r\mathrm{d}\mathrm{\Omega }.}$$

Thus, for example, a function f(r, θ, φ) can be integrated over every point in R³ by the triple integral

$${\displaystyle \underset{0}{\overset{2\pi }{\int }}\underset{0}{\overset{\pi }{\int }}\underset{0}{\overset{\mathrm{\infty }}{\int }}f(r,\theta ,\phi )r^2\sin \theta \mathrm{d}r\mathrm{d}\theta \mathrm{d}\phi .}$$

The del operator in this system leads to the following expressions for the gradient and Laplacian for scalar fields,

$${\displaystyle \begin{array}{rl}\mathrm{\nabla }f& =\frac{\mathrm{\partial }f}{\mathrm{\partial }r}\hat{\mathbf{r}}+\frac{1}{r}\frac{\mathrm{\partial }f}{\mathrm{\partial }\theta }\hat{\mathit{\theta }}+\frac{1}{r\sin \theta }\frac{\mathrm{\partial }f}{\mathrm{\partial }\phi }\hat{\mathit{\phi }},\\ {\mathrm{\nabla }}^{2}f& =\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2\frac{\mathrm{\partial }f}{\mathrm{\partial }r}\right)+\frac{1}{r^2\sin \theta }\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left(\sin \theta \frac{\mathrm{\partial }f}{\mathrm{\partial }\theta }\right)+\frac{1}{r^2{\sin }^2\theta }\frac{{\mathrm{\partial }}^{2}f}{\mathrm{\partial }{\phi }^2}\\ & =\left(\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{r}^2}+\frac{2}{r}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\right)f+\frac{1}{r^2\sin \theta }\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left(\sin \theta \frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\right)f+\frac{1}{r^2{\sin }^2\theta }\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{\phi }^2}f,\end{array}}$$

And it leads to the following expressions for the divergence and curl of vector fields,

$${\displaystyle \mathrm{\nabla }\cdot \mathbf{A}=\frac{1}{r^2}\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r^2{A}_r\right)+\frac{1}{r\sin \theta }\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left(\sin \theta A_{\theta }\right)+\frac{1}{r\sin \theta }\frac{\mathrm{\partial }{A}_{\phi }}{\mathrm{\partial }\phi },}$$

$${\displaystyle \begin{array}{rl}\mathrm{\nabla }\times \mathbf{A}=& \frac{1}{r\sin \theta }\left[\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\left(A_{\phi }\sin \theta \right)-\frac{\mathrm{\partial }{A}_{\theta }}{\mathrm{\partial }\phi }\right]\hat{\mathbf{r}}\\ & +\frac{1}{r}\left[\frac{1}{\sin \theta }\frac{\mathrm{\partial }{A}_r}{\mathrm{\partial }\phi }-\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r{A}_{\phi }\right)\right]\hat{\mathit{\theta }}\\ & +\frac{1}{r}\left[\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(r{A}_{\theta }\right)-\frac{\mathrm{\partial }{A}_r}{\mathrm{\partial }\theta }\right]\hat{\mathit{\phi }},\end{array}}$$

Further, the inverse Jacobian in Cartesian coordinates is

$${\displaystyle J^{-1}=\left(\begin{array}{ccc}{\displaystyle \frac{x}{r}}& {\displaystyle \frac{y}{r}}& {\displaystyle \frac{z}{r}}\\ \\ {\displaystyle \frac{xz}{r^2\sqrt{x^2+y^2}}}& {\displaystyle \frac{yz}{r^2\sqrt{x^2+y^2}}}& {\displaystyle \frac{-\left(x^2+y^2\right)}{r^2\sqrt{x^2+y^2}}}\\ \\ {\displaystyle \frac{-y}{x^2+y^2}}& {\displaystyle \frac{x}{x^2+y^2}}& 0\end{array}\right).}$$

The metric tensor in the spherical coordinate system is ${\displaystyle g=J^{T}J}$.

Distance in spherical coordinates

In spherical coordinates, given two points with φ being the azimuthal coordinate

$${\displaystyle \begin{array}{rl}\mathbf{r}& =(r,\theta ,\phi ),\\ {\mathbf{r}}^{\prime }& =(r^{\prime },{\theta }^{\prime },{\phi }^{\prime })\end{array}}$$

The distance between the two points can be expressed as

$${\displaystyle \begin{array}{rl}\mathbf{D}& =\sqrt{r^2+r^{\prime 2}-2r{r}^{\prime }(\sin \theta \sin {\theta }^{\prime }\cos (\phi -{\phi }^{\prime })+\cos \theta \cos {\theta }^{\prime })}\end{array}}$$

Kinematics

In spherical coordinates, the position of a point or particle (although better written as a triple${\displaystyle (r,\theta ,\phi )}$) can be written as

$${\displaystyle \mathbf{r}=r\hat{\mathbf{r}}.}$$

Its velocity is then

$${\displaystyle \mathbf{v}=\frac{\mathrm{d}\mathbf{r}}{\mathrm{d}t}=\dot{r}\hat{\mathbf{r}}+r\dot{\theta }\hat{\mathit{\theta }}+r\dot{\phi }\sin \theta \hat{\mathit{\phi }}}$$

and its acceleration is

$${\displaystyle \begin{array}{rl}\mathbf{a}=& \frac{\mathrm{d}\mathbf{v}}{\mathrm{d}t}\\ =& \phantom{}\left(\ddot{r}-r{\dot{\theta }}^2-r{\dot{\phi }}^2{\sin }^2\theta \right)\hat{\mathbf{r}}\\ & +\left(r\ddot{\theta }+2\dot{r}\dot{\theta }-r{\dot{\phi }}^2\sin \theta \cos \theta \right)\hat{\mathit{\theta }}\\ & +\left(r\ddot{\phi }\sin \theta +2\dot{r}\dot{\phi }\sin \theta +2r\dot{\theta }\dot{\phi }\cos \theta \right)\hat{\mathit{\phi }}\end{array}}$$

The angular momentum is

$${\displaystyle \mathbf{L}=\mathbf{r}\times \mathbf{p}=\mathbf{r}\times m\mathbf{v}=m{r}^2\left(-\dot{\phi }\sin \theta \hat{\mathit{\theta }}+\dot{\theta }\hat{\mathit{\phi }}\right)}$$

Where ${\displaystyle m}$ is mass. In the case of a constant φ or else θ = π/2, this reduces to vector calculus in polar coordinates.

The corresponding angular momentum operator then follows from the phase-space reformulation of the above,

$${\displaystyle \mathbf{L}=-i\hslash \mathbf{r}\times \mathrm{\nabla }=i\hslash \left(\frac{\hat{\mathit{\theta }}}{\sin (\theta )}\frac{\mathrm{\partial }}{\mathrm{\partial }\varphi }-\hat{\mathit{\varphi }}\frac{\mathrm{\partial }}{\mathrm{\partial }\theta }\right).}$$

The torque is given as

$${\displaystyle \tau =\frac{\mathrm{d}\mathbf{L}}{\mathrm{d}t}=\mathbf{r}\times \mathbf{F}=-m\left(2r\dot{r}\dot{\phi }\sin \theta +r^2\ddot{\phi }\sin \theta +2r^2\dot{\theta }\dot{\phi }\cos \theta \right)\hat{\mathit{\theta }}+m\left(r^2\ddot{\theta }+2r\dot{r}\dot{\theta }-r^2{\dot{\phi }}^2\sin \theta \cos \theta \right)\hat{\mathit{\phi }}}$$

The kinetic energy is given as

$${\displaystyle E_k=\frac{1}{2}m\left[{\left(\dot{r}\right)}^2+{\left(r\dot{\theta }\right)}^2+{\left(r\dot{\phi }\sin \theta \right)}^2\right]}$$

Computational electromagnetics

From Wikipedia, the free encyclopedia
https://en.wikipedia.org/wiki/Computational_electromagnetics

Computational electromagnetics (CEM), computational electrodynamics or electromagnetic modeling is the process of modeling the interaction of electromagnetic fields with physical objects and the environment using computers.

It typically involves using computer programs to compute approximate solutions to Maxwell's equations to calculate antenna performance, electromagnetic compatibility, radar cross section and electromagnetic wave propagation when not in free space. A large subfield is antenna modeling computer programs, which calculate the radiation pattern and electrical properties of radio antennas, and are widely used to design antennas for specific applications.

Background

Several real-world electromagnetic problems like electromagnetic scattering, electromagnetic radiation, modeling of waveguides etc., are not analytically calculable, for the multitude of irregular geometries found in actual devices. Computational numerical techniques can overcome the inability to derive closed form solutions of Maxwell's equations under various constitutive relations of media, and boundary conditions. This makes computational electromagnetics (CEM) important to the design, and modeling of antenna, radar, satellite and other communication systems, nanophotonic devices and high speed silicon electronics, medical imaging, cell-phone antenna design, among other applications.

CEM typically solves the problem of computing the E (electric) and H (magnetic) fields across the problem domain (e.g., to calculate antenna radiation pattern for an arbitrarily shaped antenna structure). Also calculating power flow direction (Poynting vector), a waveguide's normal modes, media-generated wave dispersion, and scattering can be computed from the E and H fields. CEM models may or may not assume symmetry, simplifying real world structures to idealized cylinders, spheres, and other regular geometrical objects. CEM models extensively make use of symmetry, and solve for reduced dimensionality from 3 spatial dimensions to 2D and even 1D.

An eigenvalue problem formulation of CEM allows us to calculate steady state normal modes in a structure. Transient response and impulse field effects are more accurately modeled by CEM in time domain, by FDTD. Curved geometrical objects are treated more accurately as finite elements FEM, or non-orthogonal grids. Beam propagation method (BPM) can solve for the power flow in waveguides. CEM is application specific, even if different techniques converge to the same field and power distributions in the modeled domain.

Overview of methods

The most common numerical approach is to discretize ("mesh") the problem space in terms of grids or regular shapes ("cells"), and solve Maxwell's equations simultaneously across all cells. Discretization consumes computer memory, and solving the relevant equations takes significant time. Large-scale CEM problems face memory and CPU limitations, and combating these limitations is an active area of research. High performance clustering, vector processing, and/or parallelism is often required to make the computation practical. Some typical methods involve: time-stepping through the equations over the whole domain for each time instant; banded matrix inversion to calculate the weights of basis functions (when modeled by finite element methods); matrix products (when using transfer matrix methods); calculating numerical integrals (when using the method of moments); using fast Fourier transforms; and time iterations (when calculating by the split-step method or by BPM).

Choice of methods

Choosing the right technique for solving a problem is important, as choosing the wrong one can either result in incorrect results, or results which take excessively long to compute. However, the name of a technique does not always tell one how it is implemented, especially for commercial tools, which will often have more than one solver.

Davidson gives two tables comparing the FEM, MoM and FDTD techniques in the way they are normally implemented. One table is for both open region (radiation and scattering problems) and another table is for guided wave problems.

Maxwell's equations in hyperbolic PDE form

Maxwell's equations can be formulated as a hyperbolic system of partial differential equations. This gives access to powerful techniques for numerical solutions.

It is assumed that the waves propagate in the (x,y)-plane and restrict the direction of the magnetic field to be parallel to the z-axis and thus the electric field to be parallel to the (x,y) plane. The wave is called a transverse magnetic (TM) wave. In 2D and no polarization terms present, Maxwell's equations can then be formulated as:

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}\overline{u}+A\frac{\mathrm{\partial }}{\mathrm{\partial }x}\overline{u}+B\frac{\mathrm{\partial }}{\mathrm{\partial }y}\overline{u}+C\overline{u}=\overline{g}}$$

where u, A, B, and C are defined as

$${\displaystyle \begin{array}{rl}\overline{u}& =\left(\begin{array}{c}{E}_x\\ E_y\\ H_z\end{array}\right),\\ A& =\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& \frac{1}{ϵ}\\ 0& \frac{1}{\mu }& 0\end{array}\right),\\ B& =\left(\begin{array}{ccc}0& 0& \frac{-1}{ϵ}\\ 0& 0& 0\\ \frac{-1}{\mu }& 0& 0\end{array}\right),\\ C& =\left(\begin{array}{ccc}\frac{\sigma }{ϵ}& 0& 0\\ 0& \frac{\sigma }{ϵ}& 0\\ 0& 0& 0\end{array}\right).\end{array}}$$

In this representation, ${\displaystyle \overline{g}}$ is the forcing function, and is in the same space as ${\displaystyle \overline{u}}$. It can be used to express an externally applied field or to describe an optimization constraint. As formulated above:

$${\displaystyle \overline{g}=\left(\begin{array}{c}{E}_{x,\text{constraint}}\\ E_{y,\text{constraint}}\\ H_{z,\text{constraint}}\end{array}\right).}$$

${\displaystyle \overline{g}}$ may also be explicitly defined equal to zero to simplify certain problems, or to find a characteristic solution, which is often the first step in a method to find the particular inhomogeneous solution.

Integral equation solvers

The discrete dipole approximation

The discrete dipole approximation is a flexible technique for computing scattering and absorption by targets of arbitrary geometry. The formulation is based on integral form of Maxwell equations. The DDA is an approximation of the continuum target by a finite array of polarizable points. The points acquire dipole moments in response to the local electric field. The dipoles of course interact with one another via their electric fields, so the DDA is also sometimes referred to as the coupled dipole approximation. The resulting linear system of equations is commonly solved using conjugate gradient iterations. The discretization matrix has symmetries (the integral form of Maxwell equations has form of convolution) enabling fast Fourier transform to multiply matrix times vector during conjugate gradient iterations.

Method of moments and boundary element method

The method of moments (MoM) or boundary element method (BEM) is a numerical computational method of solving linear partial differential equations which have been formulated as integral equations (i.e. in boundary integral form). It can be applied in many areas of engineering and science including fluid mechanics, acoustics, electromagnetics, fracture mechanics, and plasticity.

MoM has become more popular since the 1980s. Because it requires calculating only boundary values, rather than values throughout the space, it is significantly more efficient in terms of computational resources for problems with a small surface/volume ratio. Conceptually, it works by constructing a "mesh" over the modeled surface. However, for many problems, MoM are significantly computationally less efficient than volume-discretization methods (finite element method, finite difference method, finite volume method). Boundary element formulations typically give rise to fully populated matrices. This means that the storage requirements and computational time will tend to grow according to the square of the problem size. By contrast, finite element matrices are typically banded (elements are only locally connected) and the storage requirements for the system matrices typically grow linearly with the problem size. Compression techniques (e.g. multipole expansions or adaptive cross approximation/hierarchical matrices) can be used to ameliorate these problems, though at the cost of added complexity and with a success-rate that depends heavily on the nature and geometry of the problem.

MoM is applicable to problems for which Green's functions can be calculated. These usually involve fields in linear homogeneous media. This places considerable restrictions on the range and generality of problems suitable for boundary elements. Nonlinearities can be included in the formulation, although they generally introduce volume integrals which require the volume to be discretized before solution, removing an oft-cited advantage of MoM.

Fast multipole method

The fast multipole method (FMM) is an alternative to MoM or Ewald summation. It is an accurate simulation technique and requires less memory and processor power than MoM. The FMM was first introduced by Greengard and Rokhlin and is based on the multipole expansion technique. The first application of the FMM in computational electromagnetics was by Engheta et al.(1992). The FMM has also applications in computational bioelectromagnetics in the Charge based boundary element fast multipole method. FMM can also be used to accelerate MoM.

Plane wave time-domain

While the fast multipole method is useful for accelerating MoM solutions of integral equations with static or frequency-domain oscillatory kernels, the plane wave time-domain (PWTD) algorithm employs similar ideas to accelerate the MoM solution of time-domain integral equations involving the retarded potential. The PWTD algorithm was introduced in 1998 by Ergin, Shanker, and Michielssen.

Partial element equivalent circuit method

The partial element equivalent circuit (PEEC) is a 3D full-wave modeling method suitable for combined electromagnetic and circuit analysis. Unlike MoM, PEEC is a full spectrum method valid from dc to the maximum frequency determined by the meshing. In the PEEC method, the integral equation is interpreted as Kirchhoff's voltage law applied to a basic PEEC cell which results in a complete circuit solution for 3D geometries. The equivalent circuit formulation allows for additional SPICE type circuit elements to be easily included. Further, the models and the analysis apply to both the time and the frequency domains. The circuit equations resulting from the PEEC model are easily constructed using a modified loop analysis (MLA) or modified nodal analysis (MNA) formulation. Besides providing a direct current solution, it has several other advantages over a MoM analysis for this class of problems since any type of circuit element can be included in a straightforward way with appropriate matrix stamps. The PEEC method has recently been extended to include nonorthogonal geometries. This model extension, which is consistent with the classical orthogonal formulation, includes the Manhattan representation of the geometries in addition to the more general quadrilateral and hexahedral elements. This helps in keeping the number of unknowns at a minimum and thus reduces computational time for nonorthogonal geometries.

Cagniard-deHoop method of moments

The Cagniard-deHoop method of moments (CdH-MoM) is a 3-D full-wave time-domain integral-equation technique that is formulated via the Lorentz reciprocity theorem. Since the CdH-MoM heavily relies on the Cagniard-deHoop method, a joint-transform approach originally developed for the analytical analysis of seismic wave propagation in the crustal model of the Earth, this approach is well suited for the TD EM analysis of planarly-layered structures. The CdH-MoM has been originally applied to time-domain performance studies of cylindrical and planar antennas and, more recently, to the TD EM scattering analysis of transmission lines in the presence of thin sheets and electromagnetic metasurfaces, for example.

Differential equation solvers

Finite-Difference Frequency-Domain

Finite-difference frequency-domain (FDFD) provides a rigorous solution to Maxwell’s equations in the frequency-domain using the finite-difference method. FDFD is arguably the simplest numerical method that still provides a rigorous solution. It is incredibly versatile and able to solve virtually any problem in electromagnetics. The primary drawback of FDFD is poor efficiency compared to other methods. On modern computers, however, a huge array of problems are easily handled such as calculated guided modes in waveguides, calculating scattering from an object, calculating transmission and reflection from photonic crystals, calculate photonic band diagrams, simulating metamaterials, and much more.

FDFD may be the best "first" method to learn in computational electromagnetics (CEM). It involves almost all the concepts encountered with other methods, but in a much simpler framework. Concepts include boundary conditions, linear algebra, injecting sources, representing devices numerically, and post-processing field data to calculate meaningful things. This will help a person learn other techniques as well as provide a way to test and benchmark those other techniques.

FDFD is very similar to finite-difference time-domain (FDTD). Both methods represent space as an array of points and enforces Maxwell’s equations at each point. FDFD puts this large set of equations into a matrix and solves all the equations simultaneously using linear algebra techniques. In contrast, FDTD continually iterates over these equations to evolve a solution over time. Numerically, FDFD and FDTD are very similar, but their implementations are very different.

Finite-difference time-domain

Finite-difference time-domain (FDTD) is a popular CEM technique. It is easy to understand. It has an exceptionally simple implementation for a full wave solver. It is at least an order of magnitude less work to implement a basic FDTD solver than either an FEM or MoM solver. FDTD is the only technique where one person can realistically implement oneself in a reasonable time frame, but even then, this will be for a quite specific problem. Since it is a time-domain method, solutions can cover a wide frequency range with a single simulation run, provided the time step is small enough to satisfy the Nyquist–Shannon sampling theorem for the desired highest frequency.

FDTD belongs in the general class of grid-based differential time-domain numerical modeling methods. Maxwell's equations (in partial differential form) are modified to central-difference equations, discretized, and implemented in software. The equations are solved in a cyclic manner: the electric field is solved at a given instant in time, then the magnetic field is solved at the next instant in time, and the process is repeated over and over again.

The basic FDTD algorithm traces back to a seminal 1966 paper by Kane Yee in IEEE Transactions on Antennas and Propagation. Allen Taflove originated the descriptor "Finite-difference time-domain" and its corresponding "FDTD" acronym in a 1980 paper in IEEE Trans. Electromagn. Compat. Since about 1990, FDTD techniques have emerged as the primary means to model many scientific and engineering problems addressing electromagnetic wave interactions with material structures. An effective technique based on a time-domain finite-volume discretization procedure was introduced by Mohammadian et al. in 1991. Current FDTD modeling applications range from near-DC (ultralow-frequency geophysics involving the entire Earth-ionosphere waveguide) through microwaves (radar signature technology, antennas, wireless communications devices, digital interconnects, biomedical imaging/treatment) to visible light (photonic crystals, nanoplasmonics, solitons, and biophotonics). Approximately 30 commercial and university-developed software suites are available.

Discontinuous time-domain method

Among many time domain methods, discontinuous Galerkin time domain (DGTD) method has become popular recently since it integrates advantages of both the finite volume time domain (FVTD) method and the finite element time domain (FETD) method. Like FVTD, the numerical flux is used to exchange information between neighboring elements, thus all operations of DGTD are local and easily parallelizable. Similar to FETD, DGTD employs unstructured mesh and is capable of high-order accuracy if the high-order hierarchical basis function is adopted. With the above merits, DGTD method is widely implemented for the transient analysis of multiscale problems involving large number of unknowns.

Multiresolution time-domain

MRTD is an adaptive alternative to the finite difference time domain method (FDTD) based on wavelet analysis.

Finite element method

The finite element method (FEM) is used to find approximate solution of partial differential equations (PDE) and integral equations. The solution approach is based either on eliminating the time derivatives completely (steady state problems), or rendering the PDE into an equivalent ordinary differential equation, which is then solved using standard techniques such as finite differences, etc.

In solving partial differential equations, the primary challenge is to create an equation which approximates the equation to be studied, but which is numerically stable, meaning that errors in the input data and intermediate calculations do not accumulate and destroy the meaning of the resulting output. There are many ways of doing this, with various advantages and disadvantages. The finite element method is a good choice for solving partial differential equations over complex domains or when the desired precision varies over the entire domain.

Finite integration technique

The finite integration technique (FIT) is a spatial discretization scheme to numerically solve electromagnetic field problems in time and frequency domain. It preserves basic topological properties of the continuous equations such as conservation of charge and energy. FIT was proposed in 1977 by Thomas Weiland and has been enhanced continually over the years. This method covers the full range of electromagnetics (from static up to high frequency) and optic applications and is the basis for commercial simulation tools: CST Studio Suite developed by Computer Simulation Technology (CST AG) and Electromagnetic Simulation solutions developed by Nimbic.

The basic idea of this approach is to apply the Maxwell equations in integral form to a set of staggered grids. This method stands out due to high flexibility in geometric modeling and boundary handling as well as incorporation of arbitrary material distributions and material properties such as anisotropy, non-linearity and dispersion. Furthermore, the use of a consistent dual orthogonal grid (e.g. Cartesian grid) in conjunction with an explicit time integration scheme (e.g. leap-frog-scheme) leads to compute and memory-efficient algorithms, which are especially adapted for transient field analysis in radio frequency (RF) applications.

Pseudo-spectral time domain

This class of marching-in-time computational techniques for Maxwell's equations uses either discrete Fourier or discrete Chebyshev transforms to calculate the spatial derivatives of the electric and magnetic field vector components that are arranged in either a 2-D grid or 3-D lattice of unit cells. PSTD causes negligible numerical phase velocity anisotropy errors relative to FDTD, and therefore allows problems of much greater electrical size to be modeled.

Pseudo-spectral spatial domain

PSSD solves Maxwell's equations by propagating them forward in a chosen spatial direction. The fields are therefore held as a function of time, and (possibly) any transverse spatial dimensions. The method is pseudo-spectral because temporal derivatives are calculated in the frequency domain with the aid of FFTs. Because the fields are held as functions of time, this enables arbitrary dispersion in the propagation medium to be rapidly and accurately modelled with minimal effort. However, the choice to propagate forward in space (rather than in time) brings with it some subtleties, particularly if reflections are important.

Transmission line matrix

Transmission line matrix (TLM) can be formulated in several means as a direct set of lumped elements solvable directly by a circuit solver (ala SPICE, HSPICE, et al.), as a custom network of elements or via a scattering matrix approach. TLM is a very flexible analysis strategy akin to FDTD in capabilities, though more codes tend to be available with FDTD engines.

Locally one-dimensional

This is an implicit method. In this method, in two-dimensional case, Maxwell equations are computed in two steps, whereas in three-dimensional case Maxwell equations are divided into three spatial coordinate directions. Stability and dispersion analysis of the three-dimensional LOD-FDTD method have been discussed in detail.

Other methods

Eigenmode expansion

Eigenmode expansion (EME) is a rigorous bi-directional technique to simulate electromagnetic propagation which relies on the decomposition of the electromagnetic fields into a basis set of local eigenmodes. The eigenmodes are found by solving Maxwell's equations in each local cross-section. Eigenmode expansion can solve Maxwell's equations in 2D and 3D and can provide a fully vectorial solution provided that the mode solvers are vectorial. It offers very strong benefits compared with the FDTD method for the modelling of optical waveguides, and it is a popular tool for the modelling of fiber optics and silicon photonics devices.

Physical optics

Physical optics (PO) is the name of a high frequency approximation (short-wavelength approximation) commonly used in optics, electrical engineering and applied physics. It is an intermediate method between geometric optics, which ignores wave effects, and full wave electromagnetism, which is a precise theory. The word "physical" means that it is more physical than geometrical optics and not that it is an exact physical theory.

The approximation consists of using ray optics to estimate the field on a surface and then integrating that field over the surface to calculate the transmitted or scattered field. This resembles the Born approximation, in that the details of the problem are treated as a perturbation.

Uniform theory of diffraction

The uniform theory of diffraction (UTD) is a high frequency method for solving electromagnetic scattering problems from electrically small discontinuities or discontinuities in more than one dimension at the same point.

The uniform theory of diffraction approximates near field electromagnetic fields as quasi optical and uses ray diffraction to determine diffraction coefficients for each diffracting object-source combination. These coefficients are then used to calculate the field strength and phase for each direction away from the diffracting point. These fields are then added to the incident fields and reflected fields to obtain a total solution.

Validation

Validation is one of the key issues facing electromagnetic simulation users. The user must understand and master the validity domain of its simulation. The measure is, "how far from the reality are the results?"

Answering this question involves three steps: comparison between simulation results and analytical formulation, cross-comparison between codes, and comparison of simulation results with measurement.

Comparison between simulation results and analytical formulation

For example, assessing the value of the radar cross section of a plate with the analytical formula:

$${\displaystyle {\text{RCS}}_{\text{Plate}}=\frac{4\pi A^2}{{\lambda }^2},}$$

where A is the surface of the plate and ${\displaystyle \lambda }$ is the wavelength. The next curve presenting the RCS of a plate computed at 35 GHz can be used as reference example.

Cross-comparison between codes

One example is the cross comparison of results from method of moments and asymptotic methods in their validity domains.

Comparison of simulation results with measurement

The final validation step is made by comparison between measurements and simulation. For example, the RCS calculation and the measurement of a complex metallic object at 35 GHz. The computation implements GO, PO and PTD for the edges.

Validation processes can clearly reveal that some differences can be explained by the differences between the experimental setup and its reproduction in the simulation environment.

Light scattering codes

There are now many efficient codes for solving electromagnetic scattering problems. They are listed as:

• discrete dipole approximation codes,
• codes for electromagnetic scattering by cylinders,
• codes for electromagnetic scattering by spheres.

Solutions which are analytical, such as Mie solution for scattering by spheres or cylinders, can be used to validate more involved techniques.