Electromagnetic monitoring and control of a plurality of nanosatellites

Abstract

A method for monitoring position of and controlling a second nanosatellite (NS) relative to a position of a first NS. Each of the first and second NSs has a rectangular or cubical configuration of independently activatable, current-carrying solenoids, each solenoid having an independent magnetic dipole moment vector, μ1 and μ2. A vector force F and a vector torque are expressed as linear or bilinear combinations of the first set and second set of magnetic moments, and a distance vector extending between the first and second NSs is estimated. Control equations are applied to estimate vectors, μ1 and μ2, required to move the NSs toward a desired NS configuration. This extends to control of N nanosatellites.

Claims

1. A method of monitoring and controlling a second nanosatellite relative to a position of a first nanosatellite, the method comprising: (i) initializing a time interval; (ii) for a present time interval, based on at least measured magnetometer data, estimating a distance vector that extends between a first location associated with a first nanosatellite and a second location associated with a second nanosatellite, the nanosatellites comprising respective first and second plurality of solenoids; (iii) estimating force vectors and torque vectors to move at least the second plurality of solenoids relative to the first nanosatellite; (iv) for at least the second nanosatellite, estimating magnetic dipole moment vectors associated with both force vectors and torque vectors; (v) causing at least the second plurality of solenoids to produce the magnetic dipole moment vectors; (vi) incrementing the present time interval, and returning to step (ii), wherein the first plurality of solenoids is arranged in a rectangular geometry.

2. The method according to claim 1, wherein the first location is on a first solenoid within the first plurality of solenoids, and the second location is on a second solenoid within the second plurality of solenoids.

3. The method according to claim 2, wherein the measured magnetometer data comprises a magnetic field vector value associated with the second solenoid.

4. The method according to claim 3, further comprising using the magnetic field vector value to estimate a magnitude of the distance vector.

5. The method according to claim 1, wherein each of the first plurality of solenoids is a current-carrying solenoid that can be independently activated.

6. The method according to claim 1, further comprising using an RF link to communicate commands from the first nanosatellite to the second nanosatellite to control the second plurality of solenoids.

7. A method of monitoring and controlling a second nanosatellite relative to a position of a first nanosatellite, the method comprising: (i) initializing a time interval; (ii) for a present time interval, based on at least measured magnetometer data, estimating a distance vector that extends between a first location associated with a first nanosatellite and a second location associated with a second nanosatellite, the nanosatellites comprising respective first and second plurality of solenoids, the first location being on a first solenoid within the first plurality of solenoids, the second location being on a second solenoid within the second plurality of solenoids, the measured magnetometer data comprising a magnetic field vector value associated with the second solenoid; (iii) estimating force vectors and torque vectors to move at least the second plurality of solenoids relative to the first nanosatellite; (iv) for at least the second nanosatellite, estimating magnetic dipole moment vectors associated with both force vectors and torque vectors; (v) causing at least the second plurality of solenoids to produce the magnetic dipole moment vectors; (vi) incrementing the present time interval, and returning to step (ii), and (vii) using the magnetic field vector value to estimate a normalized vector value associated with the distance vector.

8. The method according to claim 7, wherein each of the first plurality of solenoids is a current-carrying solenoid that can be independently activated.

9. The method according to claim 7, further comprising using an RF link to communicate commands from the first nanosatellite to the second nanosatellite to control the second plurality of solenoids.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 schematically illustrates use of a single known normalized magnetic dipole moment vector μ^ and corresponding magnetic field strength vector B to estimate a distance vector d extending between first and second nanosatellites, relative to a coordinate system determined by μ and B, and use of this information to estimate electromagnetic forces and torques on a given nanosatellite, according to a first embodiment.

(2) FIGS. 2, 4 and 6 are flow charts of different procedures for estimating distance vectors d and angular orientations of non-reference nanosatellites relative to a reference nanosatellite.

(3) FIG. 3 graphically illustrates variation of a length squared L(ψ−ψ0).sup.2=|μ^−(μ^.Math.α)α|.sup.2) as a function of an angle ψ−ψ0 between the normalized distance vector α and the normalized magnetic moment vector μ^.

(4) FIG. 5 illustrates orientation of α for one embodiment.

(5) FIG. 7 illustrates two solenoid cubes (i=1 and k=2) confined to a plane P(x,y).

(6) FIG. 8 is a general flow chart for practicing the invention.

DESCRIPTION OF THE INVENTION

(7) One challenge in monitoring and controlling one or more nanosatellites, using electromagnetic forces and torques, is determination of all relevant components of the three dimensional distance vectors that naturally arise in description of the physical parameters and resulting forces in such an environment. Beginning with the expression in Eq. (2) for a magnetic field flux that is generated by presence of a known magnetic moment μ, one is first confronted with estimation of the distance vector d (with source at a center of the magnetic moment source or associated solenoid(s)) at the point of measurement of the magnetic field B(d). Equation (2) is analyzed to estimate the vector components of d, using the known magnetic moment vector μ and the measurable field vector B(d). Equation (2) is rewritten as
B′(d)=(4π/μ0)B(d)={μ−(μ.Math.d)d/d.sup.2}/d.sup.3={μ−(μ.Math.α)α}/d.sup.3.  (3)
d=dα,  (4)
where α is a normalized vector (|α|=1) whose direction is parallel to the (unknown) direction of d. The normalized distance vector α is re-expressed as a sum of components,
α=αn1^+bn2^+cn3^,  (5)
where n1^, n2^ and n3^ are mutually perpendicular vectors, each of unit length, and a, b and c are to be estimated. Note that both the scalar value d and the vector α are initially unknown. One suitable set of choices for the vectors n1^, n2^ and n3^ is
n1^=μ^=μ/|μ,  (6A)
n2^={B′−(B′.Math.μ^)μ^}/|B′−(B′.Math.μ^)μ^|,  (6B)
n3^=n1^×n2,  (6C)
(n1^).sup.2=(n2^).sup.2=(n3^).sup.2=1.  (6D)

(8) Using the decompositions set forth in Eqs. (4), (5) and (6), Eq. (3) can be rewritten

(9) $\begin{array}{cc}\begin{array}{c}{B}^{\prime }\left(d\right)/\mu =\left\{B^{\prime }-\left(B^{\prime }.Math.{\mu }^{.Math.}\right){\mu }^{.Math.}+\left(B^{\prime }.Math.{\mu }^{.Math.}\right){\mu }^{.Math.}\right\}/\mu \\ =\left\{{\mu }^{.Math.}-\left({\mu }^{.Math.}.Math.\alpha \right)\alpha \right\}{d}^3,\\ =\{{\mu }^{.Math.}-a\left\{a{\mu }^{.Math.}+bn{2}^{.Math.}+cn{3}^{.Math.}\right\}/d^3,\\ =\{\left(1-{\left(a\right)}^2\right){\mu }^{.Math.}-\left(\mathrm{ab}\right)\left\{B^{\prime }-\left({\mu }^{.Math.}.Math.B^{\prime }\right){\mu }^{.Math.}\right\}/\\ .Math.B^{\prime }-\left({\mu }^{.Math.}.Math.B^{\prime }\right){\mu }^{.Math.}.Math.\}\left(1/d^3\right),\end{array}& \left(7\right)\end{array}$
where the term involving the vector n3^ is dropped here because this vector only appears in one term in Eq. (7) (equivalent to c=0). Equation (7) can be rewritten by collecting all coefficients involving the vector μ^ and all coefficients involving the perpendicular vector B′−(B′.Math.μ^)μ^ in Eq. (7),
(1/μ){B′−((B′.Math.μ^)μ^+(B′.Math.μ^)μ^}={(1−a.sup.2)μ^−ab{B′−(B′.Math.μ^)μ^}/|(B′−B′.Math.μ^)μ^|}(1/d.sup.3).  (8)
Separately collecting all coefficients of each of the two vectors μ^ and B′−(B′.Math.μ^)μ^ yields the relations
(1−a.sup.2)/d.sup.3=(B′.Math.μ^)/μ,  (9)
b=(−d.sup.3/aμ)|B′−((B′.Math.μ^)μ^|,  (10)
1−a.sup.2−b.sup.2=(B′.Math.μ^)d.sup.3/μ,−(1/a.sup.2)|(B′−B′.Math.μ^)μ^|.sup.2d.sup.6/μ.sup.2,   (11)
a.sup.2={B′.sup.2−(B′.Math.μ^).sup.2μ^}d.sup.3//μ(B′.Math.μ^)  (12)
d.sup.3(B′).sup.2/(B′.Math.μ)=μ,  (13)

(10) Equations (7) and (9)-(13) determine the quantities d and a, and the distance vector d. The vector quantities B′(d) and B(d) are measurable for any selected distance vector d.

(11) One suitable procedure is the following, illustrated in a flow chart in FIG. 2. In step 21, a first magnetic moment μ1, associated with the reference and controllable by the associated solenoid current, is selected. In step 22, a magnetic field, B′(d) or B(d), whose source is the magnetic moment μ1, is measured at a selected distance vector d=d12 (not yet estimated), which extends from a center of the reference to a center of a second (non-reference) nanosatellite. In step 23 (optional), the distance vector d12=dα is estimated for a location of the center of the second nanosatellite or associated solenoid relative to the center of the reference. Steps 22 and 23 can be repeated to estimate a distance vector d for location of the center of each non-reference nanosatellite solenoid.

(12) Electromagnetic forces, F12 and F21=−F12 and the corresponding torques, τ12 and τ21=τ12 are computed in this two-moment system, the center locations for each non-reference nanosatellite are incremented for a sequence of time increments Δt.

(13) A second method for estimation of the direction, in a plane P, of the normalized distance vector d uses a combined activation of two distinct solenoids to provide a first hybrid magnetic moment, μ1′=j1μ1+k1μ2, and to provide a second hybrid magnetic moment, μ2′=j2μ1+k2μ2, where μ1 and μ2 are magnetic moments of equal magnitude μ produced by activation of only the first solenoid and by activation of only the second solenoid, respectively, and j1, k1, j2 and k2 are selected coefficients, with j1/k1≠±j2/k2 and |μ1|=|μ2|=μ. Equation (3) is re-expressed for normalized magnetic moments, μ1^ and μ2^, as
B1′(d)=(4π/(μμ0)B1(d)={μ1′^−(μ1′^.Math.α)α}/d.sup.3,  (14)
B2′(d)=(4π/(μμ0)B2(d)={μ2′^−(μ2′^.Math.α)α}/d.sup.3,  (14)
|μ1^|=μ2^|=1.  (16)

(14) The vector lengths |B1′(d)| and |B2′(d)| will generally have different values, because the scalar products (μ1′^.Math.α) and (μ2′^.Math.α) will have different values. In a graphical illustration in FIG. 3, it is assumed that μ1′^.Math.α=cos ψ1 and μ2^.Math.α=cos ψ2 and that an angle value, ψ=ψ0, that minimizes the quantity
L(ψ).sup.2=|μ′^−(μ′^.Math.α)α={1−(μ′^.Math.α).sup.2}=1−cos.sup.2(ψ−ψ0)  (17)
is not yet known. The quantities L(ψ1).sup.2 and L(ψ2).sup.2 are measurable, using Eqs. (14)-(17), and can be re-expressed as

(15) $\begin{array}{cc}\left[\begin{array}{cc}\cos \left(2\psi 1\right)& \sin \left(\psi 1\right)\\ \cos \left(2\psi 2\right)& \sin \left(2\psi 2\right)\end{array}\right]\left[\begin{array}{c}\cos \left(2\psi 0\right)=1-2{L\left(\psi 1\right)}^2\\ \sin \left(2\psi 0\right)=1-2{L\left(\psi 2\right)}^2\end{array}\right]& \left(18\right)\end{array}$
where the angle ψ0 is to be estimated and other quantities are determinable. The 2×2 matrix in Eqs. (18) has an associated determinant, with a value, Det=sin 2(ψ2−ψ1), whose magnitude can be made as large as possible by one of the choices ψ2−ψ1=±π/4, ±3π/4. The solutions of Eqs. (18) are

(16) $\begin{array}{cc}\left[\begin{array}{c}\cos \left(2\psi 0\right)\\ \sin \left(2\psi 0\right)\end{array}\right]\begin{array}{c}=\\ =\end{array}\left[\begin{array}{cc}\sin \left(2\psi 2\right)& -\sin \left(2\psi 1\right)\\ -\cos \left(2\psi 2\right)& \cos \left(2\psi 1\right)\end{array}\right]\left[\begin{array}{c}\left\{1-2{L\left(\psi 1\right)}^2\right\}/\mathrm{Det}\\ \left\{1-2{L\left(\psi 2\right)}^2\right\}/\mathrm{Det}\end{array}\right]& \left(19\right)\end{array}$
The solution angle, ψ=ψ0, provides the angle at which the normalized distance vector α is oriented relative to the direction of the first magnetic moment μ1 in the plane P=P(x,y), as illustrated in FIGS. 1 and 3. This alternative procedure provides a second procedure for estimating the normalized distance vector α. After the distance vector d=dα. has been estimated, Eq. (3) may be used to estimate the scalar distance value d.

(17) FIG. 4 is a flow chart of a second embodiment of a procedure for estimating the normalized distance vector α and the scalar distance value d. In step 41, first and second magnetic moments, μ1 and μ2, and hybrid magnetic moments, μ1′ and μ2′, as linear combinations of μ1 and μ12 (with associated angles ψ1 and ψ2), are provided for the reference satellite. In step 42, the length proportions squared, L(ψ1−ψ0).sup.2=|μ1′^−(μ1′^.Math.α)α|.sup.2 and L(ψ2−ψ0).sup.2=|μ2′^−(μ2′^.Math.α)α|.sup.2, are measured for known and distinct angles ψ=ψ1 and ψ=ψ2 for the hybrid magnetic moments, μ1′ and μ2′. Here, ψ=ψ0 is an angle value, unknown as yet, that minimizes the length squared value L(ψ−ψ0).sup.2.

(18) In step 43, solutions, cos ψ0 and sin ψ0, are estimated for the equations
cos(2ψ1)cos(2ψ0))+sin(2ψ1)sin(2ψ0)=1−2L(ψ1).sup.2,
cos(2ψ2)cos(2ψ0))+sin(2ψ2)sin(2ψ0)=1−2L(ψ2).sup.2.
In step 44, the solution angle, ψ=ψ0, is interpreted as the angle that distance vector d makes with a non-hybrid, magnetic moment vector μ1 or μ2. This indicates a direction in the plane P(x,y) for the normalized distance vector α. In step 45 (optional), the scalar distance value d=|d| is estimated, using, for example, Eq. (3).

(19) The developments of Eqs. (3)-(15), or (16)-(22), assume that the magnetic moment vector μ and the corresponding measured magnetic field vector B(d) are known. This requires knowledge of B(d) and μ for a single (initially unknown) distance vector d at which the field is measured. Alternatively, both the distance vector d and the direction of the moment μ may be unknown initially, although the magnitude, |μ|=μ may be known, through control of the current supplied to the solenoid. In this situation, it is assumed that two solenoids, each with the same magnetic moment magnitude, |μ1|=|μ2|=μ, are oriented at a known, non-zero angle (preferably π/2) relative to each other, that the two solenoids are associated with each other at a common end, and that the electromagnetic forces and motions are confined to a plane P(x,y). The moments are expressed as analogous to Eqs. (3) and (4):
B1(d)=(μ0/4πd.sup.3){μ1−(μ1.Math.d)d/d.sup.2)=(μ0/4πd.sup.3){μ1−(μ1.Math.α)α},  (20)
B2(d)=(μ0/4πd.sup.3){μ2−(μ2.Math.α.Math.d)d/d.sup.2)=(μ0/4πd.sup.3){μ2−(μ2.Math.α)α}  (21)
d=dα,  (22)
where the two solenoids have (approximately) a common distance vector d, whose angular orientation relative to the directions of the first and second solenoid sin the plane P(x,y) is initially unknown. Note that the field vectors B1 and B2 are each perpendicular to the vector d or a:
B1.Math.α=0,  (23A)
B2.Math.α=0,  (23B)
although generally |B1|≠|B2|, because the moments μ1 and μ2 are generally oriented differently relative to the vector α.

(20) Introduce new unit length vectors, which are measurable,
n1^=B1/|B1|,  (24)
n2^=B2/|B2|,  (25)
n2′^={n2^−(n1^.Math.n2)n1^}/|n2^−(n1^.Math.n2)n1^|,  (26)
n3^=n1^×n2′^,  (27)
where n1^, n2′^ and n3′ are mutually perpendicular and each has unit length. The normal forms of a plane P1 with normal vector B1 and of a plane P2 with normal vector B2 are
n1^.Math.s−|B1|=0,  (28)
n2′^.Math.s−|B2|=0,  (29)
where s=(x,y,z) is a coordinate vector in the coordinate system adopted here and |B1| and |B2| are the lengths of a first normal vector (n1^) and a second normal vector (n2′^), extending from the coordinate origin to the first plane P1 and from the coordinate origin to the second plane P2, respectively. Each of the planes, P1 and P2, contains the distance vector d, and the two normal vectors, n1^ and n2′^, (or n1^ and n2^) are not parallel. The intersection of the planes P1 and P2, therefore, defines a common line CL, which is parallel to the distance vector d.

(21) Express the common line CL, lying in the intersection P1ΩP2, as a vector
s=s1n1^+s2n2′^+s3n3^.  (30)
The coefficients s1, s2 and s3 are determined using the normal forms in Eqs. (28) and (29):
n1′^.Math.s−|B1|=s1(n1^.Math.n1^)+s2(n1^.Math.n2′^)+s3(n1^.Math.n3^)−|B1|=0,  (31)
n2′^.Math.s−|B2|=s1(n2′^.Math.n1^)+s2(n2′^.Math.n2′^)+s3(n2′^.Math.n3^)−|B2|=0,  (32)
Taking account of the mutual orthogonality of n1^, n2′^ and n3^, a solution of Eqs. (31) and (32) is
s1=|B1|,  (33A)
s2=|B2|,  (33B)
s3 is not determined.  (33C)
The intersection line CL is expressible as
s=|B1|n1^+|B2|n2′^+s3n3^,  (34)
α=s/|s|,  (35)
which defines a direction of the normalized distance vector α=d/|d|. FIG. 5 illustrates formation of the normalized distance vector α. This requires that B1 and B2 not be parallel, but does not require that B1 and B2 be perpendicular.

(22) The scalar value d of the distance vector d is estimated using the relation
|B1|.sup.2+|B2|.sup.2=(μ0/4πd.sup.3).sup.2{2μ.sup.2−μ.sup.2(cos.sup.2 ψ+sin.sup.2 ψ)}=(μ0/4πd.sup.3).sup.2μ.sup.2,  (36)
where ψ is an angle between the vector α in the plane P(x,y) and the vector direction μ1^. The angle ψ itself, illustrated in FIG. 5, is estimated using the relation

(23) $\begin{array}{cc}\begin{array}{c}B1.Math.B2={\left(\mu 0/4\pi d^3\right)}^2\left\{\mu 1.Math.\mu 2-{\mu }^2\cos \psi \sin \psi \right\}\\ =-{\left(\mu 0/4\pi {\mu }^3\right)}^2\left({\mu }^2/2\right)\sin 2\psi \},\end{array}& \left(37\right)\end{array}$
where it is assumed here that the scalar value d is already estimated and that the magnetic moments μ1 and μ2 are orthogonal to each other so that μ1, μ2=0. These procedures, using Eqs. (34)-(35) or Eqs. (36)-(37), provide third and fourth alternative approaches to estimation of the scalar value d and the normalized distance vector α, which do not require ab initio knowledge of the separate vector directions of the magnetic moments μ1 and μ2, only the requirement that μ1μ2=0.

(24) A suitable procedure for these alternatives is illustrated in a flow chart in FIG. 6. In step 61, the reference causes each non-reference nanosatellite, to generate two or more distinct magnetic moments, μ1 and μ2, by activating two or more solenoids on the non-reference nanosatellite at separate times. In step 62, a magnetic field strength vector, B1 and B2, for the respective magnetic moments is generated and measured at the center of the reference solenoids. In step 63, the distance vector d, extending from the center of the non-reference nanosatellite solenoid(s) to the center of the reference solenoid(s), is estimated, using the common line CL of intersection between the planes, P1 and P2, generated by the plane normals, B1 and B2. In step 64, the directions for the magnetic moments, μ1 and μ2, associated with the non-reference nanosatellite, are estimated. The angles (ψ, π/2−ψ) in the plane P between the distance vector d and the respective magnetic moments, μ1 and μ2, are estimated, in step 65. This determines the location coordinates of the center of the non-reference nanosatellite solenoid(s) and the angular orientation of the non-reference nanosatellite, relative to the center and coordinate axes orientation of the reference. In this embodiment, solenoids that produce the magnetic moments, μ1 and μ2, are activated on a non-reference nanosatellite, the corresponding magnetic field vectors, B1 and B2, thus produced are measured at the reference nanosatellite, and the corresponding angular orientations of the solenoids (indicated by μ1 and μ2) on the non-reference nanosatellite are estimated, using the angle pair (ψ, π/2−ψ) estimated at the reference. In step 66 (optional), the scalar distance value d is estimated, using Eq. (3). The roles of the reference and non-reference nanosatellites can be exchanged.

(25) Application of Distance Vector Information to Satellite Location and Control.

(26) In the development and application of the control law to follow, equations of motion of each of first and second solenoid rectangles (referred to herein as “cubes), relative to each other are developed. The first and second solenoid cubes lie in a common plane P(x,y), illustrated in FIG. 7, and their subsequent motions are expressed in terms of motions only in this common plane. This model of two interacting cubes, the first a reference nanosatellite (i=1) and the second a controlled nanosatellite (k=2) is then extended to one reference satellite and to N−1 controlled nanosatellites (k=2, . . . , N), all lying in the common plane P(x,y) and each with its own activatable solenoids.

(27) The first or reference cube (i=1) is defined by four line segments, indexed as j=1, 2, 3, 4, and by a cube center c1. Each of the four line segments for the first cube has a pair of indices, (ij)=(1,1), (1,2), (1,3) and (1,4), as illustrated in FIG. 7. In a similar manner, the second (controlled) nanosatellite (indexed here as k=2 or as k=2, 3, . . . , N) is defined by each of four line segments (m=1, 2, 3, 4) and by a center c2, with the pair of indices (k,m)=(2,1), (2,2), (2,3) and (2,4). Motions of the second nanosatellite (translation of the center, rotation about the center) are preferable stated with reference to location and angular orientation of the first or reference nanosatellite. Each solenoid (line segment) of each of the first and second nanosatellites is independently activated, and the magnitude and direction of the solenoid electrical current is independently selected.

(28) Each of the four solenoid line segments that defines a nanosatellite solenoid (i=1 and k=2) has a line segment center. A distance vector, d=d(ij;k,m) (also written as d.sub.ijkm in the following) extending between a line segment (i,j) and a line segment (k,m) connects these two line centers in the development to follow. Each pair of solenoid line segments, (i,j) and (k,m), on the first and second cubes, respectively, has an associated pair of magnetic dipole moment vectors, μ.sub.i,j and μ.sub.k,m, independently chosen, and an associated distance vector d.sub.ijkm. Each of the quantities, μ.sub.i,j, μ.sub.km, and d.sub.ijkm will vary with time. A magnetic dipole moment, μ.sub.i,j and/or μ.sub.km may be 0 if that particular solenoid is not presently activated. More generally, each magnetic dipole moment vector, μ.sub.i,j and/or μ.sub.km may have a vector value 0 during one or more time intervals and may have a non-zero vector value during one or more other time intervals. The magnetic dipole moment vectors, μ.sub.i,j and μ.sub.km, will be chosen during each time interval to move the two cubes (more generally, to move the N cubes) toward a desired cube configuration (relative location and relative angular orientation) over the course of time, and the desired cube configuration may change with passage of time. The preceding formalism for two nanosatellites is extended to a first (reference) nanosatellite and N−1 controlled or non-reference nanosatellites, all lying in a common plane P(x,y), where N is a reasonable number (e.g., N=2-20).

(29) The length Δt of all time intervals may be uniform (e.g., Δt=0.1 sec), or the time interval length may itself vary with time (e.g., Δt=0.01-0.5 sec) according to the rate at which the present cube configuration is to approach a desired cube configuration. Movement of the different cubes toward a desired cube configuration is facilitated by temporally varying values of a force vector F.sub.ijkm and/or by a temporally varying electrical torque vector τ.sub.ijkm. The force vectors F.sub.ijkm and/or the torque vectors τ.sub.ijkm vary in response to the temporally changing (non-zero) magnetic dipole moments, μ.sub.i,j and/or μ.sub.km, and the temporally changing solenoid distance vectors d.sub.ijkm.

(30) FIG. 8 is a general flow chart illustrating broadly a procedure for a time interval sequence to estimate distance vectors d, magnetic dipole moment vectors μ, force vectors F, torque vectors τ, and solenoid commands that are presently estimated to be required to move a present cube configuration toward a desired cube configuration. In step 81, an initial time interval, t.sub.n≦t<t.sub.n+1 with n=0, is identified. In step 82, a present distance vector, d=d.sub.n is estimated for each nanosatellite, relative to a first (reference) nanosatellite, using magnetometer measurements and magnetic dipole moment vectors μ.sub.n, if known, for each nanosatellite, for the present time interval (t.sub.n≦t<t.sub.n+1). Step 82 may be implemented using any of the several embodiments discussed in connection with Eqs. (1)-(37), or any other suitable procedure for estimation of a distance vector, d=d.sub.n, for each relevant or active solenoid on each nanosatellite.

(31) In step 83, a force vector F and a torque vector T are estimated for each nanosatellite that are required to move each nanosatellite toward a desired cube configuration, for the present time interval (t.sub.n≦t<t.sub.n+1). In step 84, a magnetic dipole moment vector μ is estimated for each solenoid on each nanosatellite that is required to provide the force vectors F and torque vectors τ, for the present time interval (t.sub.n≦t<t.sub.n+1). In step 85, the required solenoid commands (estimated in step 84) are executed, for the present time interval (t.sub.n≦t<t.sub.n+1). In step 86, the present time interval, t.sub.n≦t<t.sub.n+1, is replaced by the next time interval, t.sub.n+1≦t<t.sub.n+2 (equivalent to incrementing the time interval index n to n+1), and the system returns to step 82.

(32) Implementation of Nanosatellite Control Commands.

(33) Returning to FIG. 7, the force on a given cube is equal to the sum of all the force interactions between the EMs on one cube and the EMs on the other cube. The force on the i.sup.th cube generated by the k.sup.th cube is given by

(34) $\begin{array}{cc}{\mathrm{Fc}}_i=\underset{j=1}{\overset{4}{.Math.}}\underset{k=1,k\ne i}{\overset{\#\mathrm{Cubes}}{.Math.}}\underset{m=1}{\overset{4}{.Math.}}{F}_{\mathrm{ijkm}},& \left(41\right)\end{array}$
where
F.sub.ijkm=M.sub.F(μ.sub.km,d.sub.ijkm)μ.sub.ij  (42)
Using the far field model, we have, from Appendix I,

(35) $\begin{array}{cc}{M}_F\left({\mu }_{\mathrm{km}},d_{\mathrm{ijkm}}\right)=-\frac{3\mu }{4{\pi d}_{\mathrm{ijkm}}^5}\left(\left(d_{\mathrm{ijkm}}.Math.{\mu }_{\mathrm{km}}\right)+\left({\mu }_{\mathrm{km}}.Math.d_{\mathrm{ijkm}}\right)+\left({\mu }_{\mathrm{km}}^T{\mathrm{•d}}_{\mathrm{ijkm}}\right)\left(I-5\frac{d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijkm}}}{d_{\mathrm{ijkm}}^T{\mathrm{•d}}_{\mathrm{ijkm}}}\right)\right)& \left(43\right)\end{array}$
The symbol is the outer product operator.

(36) Similarly, the torque on a given cube is the sum of all the torque interactions between the EMs on one cube and the EMs on the other cube, plus a cross product of the moment arm from the center of the cube to the center of the EM with the force on that EM. Thus, the torque on the i.sup.th cube by the k.sup.th cube is given by

(37) $\begin{array}{cc}\tau c_i=\underset{j=1}{\overset{4}{.Math.}}\underset{k=1,k\ne i}{\overset{\#\mathrm{Cubes}}{.Math.}}\underset{m=1}{\overset{4}{.Math.}}{\tau }_{\mathrm{ijkm}},& \left(44\right)\end{array}$
where
τ.sub.ijkm=M.sub.τ(μ.sub.km,d.sub.ijkm)μ.sub.ij,  (45)

(38) $\begin{array}{cc}{M}_{\tau }\left({\mu }_{\mathrm{km}},d_{\mathrm{ijkm}}\right)=\frac{{\mu }_0}{4\pi }{\left(\frac{3\left({\mu }_{\mathrm{km}}.Math.d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijkm}}^5}{d}_{\mathrm{ijkm}}-\frac{{\mu }_{\mathrm{km}}}{d_{\mathrm{ijkm}}^3}\right)}_x^T+{\mathrm{rc}}_{\mathrm{ij}}\times M_F\left({\mu }_{\mathrm{km}},d_{\mathrm{ijkm}}\right)\mathrm{and}& \left(46\right)\end{array}$
rc.sub.ij the moment arm on from the center of the i.sup.th cube to the center of the j.sup.th EM The symbol ( ).sup.T, in Eq. (46) is the skew matrix.

(39) To be able to control the nanosatellites, Eqs. (44)-(46) must be solved for the control input μ for each nanosatellite. To accomplish this, the μs will be determined in the reference frame that is attached to the body of the reference nanosatellite, referred to as the body frame.

(40) To accomplish this, a transformation from the inertial frame to the body frame is needed. These μs will go through a rotation θ from the inertial frame to the body frame. Also note that the μs are note free to point in any arbitrary direction, the μs are fixed to the body frame, and thus the x-y components of each μ will be isolated by projection matrices, I.sub.x and I.sub.y, associated with the x-axis and the y-axis. This can easily be extended to the z-axis.

(41) To demonstrate this, we will work through the equations for two nanosatellites, which is easily extended to any reasonable number N of nanosatellites. The force acting on nanosatellite 1 (reference) can be written as
Fc.sub.1=M.sub.F.sub.11μ.sub.11+M.sub.F.sub.12μ.sub.12+M.sub.F.sub.13μ.sub.13+M.sub.F.sub.14μ.sub.14  (47)

(42) ${\mathrm{Fc}}_1=\left[\begin{array}{c}{\mathrm{Fc}}_{1x}\\ {\mathrm{Fc}}_{1y}\end{array}\right]$
is the components of the force acting on nanosatellite 1 in the inertial frame,
M.sub.F.sub.11=M.sub.F(μ.sub.21,d.sub.1121)+M.sub.F(μ.sub.22,d.sub.1122)+M.sub.F(μ.sub.23,d.sub.1123)+M.sub.F(μ.sub.24,d.sub.1124),
M.sub.F.sub.12=M.sub.F(μ.sub.21,d.sub.1221)+M.sub.F(μ.sub.22,d.sub.1222)+M.sub.F(μ.sub.23,d.sub.1223)+M.sub.F(μ.sub.24,d.sub.1224),
M.sub.F.sub.13=M.sub.F(μ.sub.21,d.sub.1321)+M.sub.F(μ.sub.22,d.sub.1322)+M.sub.F(μ.sub.23,d.sub.1323)+M.sub.F(μ.sub.24,d.sub.1324),
M.sub.F.sub.12=M.sub.F(μ.sub.21,d.sub.1421)+M.sub.F(μ.sub.22,d.sub.1422)+M.sub.F(μ.sub.23,d.sub.1423)+M.sub.F(μ.sub.24,d.sub.1424)  (48)
Resolving and isolating the components of μ in the body frame Eq. (41) is

(43) 0$\begin{array}{cc}{\mathrm{Fc}}_1=[\begin{array}{cc}{M}_{F_{11}}R\left(\theta \right)I_x+M_{F_{12}}R\left(\theta \right)I_y& M_{F_{13}}R\left(\theta \right)I_x+M_{F_{14}}R\end{array}\left(\theta \right)I_y]\hspace{1em}\left[\begin{array}{c}{\mu }_{11x}\\ {\mu }_{12y}\\ {\mu }_{13x}\\ {\mu }_{14y}\end{array}\right]& \left(49\right)\end{array}$
where μ.sub.ijx and μ.sub.iky represents the μs in the body frame,

(44) $\begin{array}{cc}{I}_x=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],I_y=\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right],\mathrm{and}R\left({\theta }_1\right)=\left[\begin{array}{ccc}\cos {\theta }_1& -\sin {\theta }_1& 0\\ \sin {\theta }_1& \cos {\theta }_1& 0\\ 0& 0& 1\end{array}\right]& \left(50\right)\end{array}$
Applying the same process to the torque equation, we have

(45) $\begin{array}{cc}\tau c_{1z}=[\begin{array}{cc}{M}_{\tau 11}R\left({\theta }_1\right)I_x+M_{\tau 12}R\left({\theta }_1\right)I_y& M_{\tau 13}R\left({\theta }_1\right)I_x+M_{\tau 14}R\left({\theta }_1\right)I_y\end{array}]\hspace{1em}\left[\begin{array}{c}{\mu }_{11x}\\ {\mu }_{12y}\\ {\mu }_{13x}\\ {\mu }_{14y}\end{array}\right]& \left(51\right)\end{array}$
where
M.sub.τij=M.sub.τ(μ.sub.21,d.sub.ij21)+M.sub.τ(μ.sub.22,d.sub.ij22)+M.sub.τ(μ.sub.23,d.sub.ij23)+M.sub.τ(μ.sub.24,d.sub.ij24),

(46) Writing Eqs. (50), (51 and (52) in matrix form, we have

(47) $\begin{array}{cc}[\begin{array}{c}{\mathrm{Fc}}_{1x}\\ {\mathrm{Fc}}_{1y}\\ {\mathrm{Fc}}_{1z}\end{array}]=\hspace{1em}[\begin{array}{cc}{M}_{F_{11}}R\left({\theta }_1\right)I_x+M_{F_{12}}R\left({\theta }_1\right)I_y& M_{F_{13}}R\left({\theta }_1\right)I_x+M_{F_{14}}R\left({\theta }_1\right)I_y\\ M_{F\tau 11}R\left({\theta }_1\right)I_x+M_{F\tau 12}\left({\theta }_1\right)I_y& M_{F\tau 13}R\left({\theta }_1\right)I_x+M_{F\tau 14}R\left({\theta }_1\right)I_y\end{array}]\hspace{1em}[\begin{array}{c}{\mu }_{11x}\\ {\mu }_{12y}\\ {\mu }_{13x}\\ {\mu }_{14y}\end{array}]& \left(53\right)\end{array}$
Using a pseudo inverse operator, we can solve for the μs. Note, we are using four variables to solve for three values. This is an extra degree of freedom in the control input that can be used to optimize another constraint, if desired.
Control Law Signal.

(48) In the previous section we showed that, given a desired force, torque applied to c.sub.1, the distance from each EM on c.sub.2 to each EM on c.sub.1, and the μs of c.sub.2 we can solve for the μs of c.sub.1. Here we will explore a very simple feedback control law to demonstrate in simulation that given a desired relative position and pointing angle, c.sub.1 can be driven to the desired position and angle by control of c.sub.1's μs.

(49) Because we want to exercise relative position control, we will define RelativePosition as the relative position from c.sub.2 to c.sub.1, expressed as a vector. We formulate an error signal between the relative position vector and a desired relative position we call DesiredDistance.
Error=RelativePosition−DesiredDistance.  (54)

(50) For a linear time invariant plant, the proportional part of a PD controller is the Error multiplied by a gain. For an EM plant, this signal would be a force command Because force between two magnetic dipole moments is inversely proportional to the relative distance raised to the 4.sup.th power, we compensate for this non-linearity by multiplying the Error signal by |RealtivePosition|.sup.4. To include damping, we multiply the difference between the velocities by a gain. The control law has the form
pdSignal=Error kp|relativePosition|.sup.4+Error kd.  (55)

(51) For the orientation of the nanosatellite we formulate an angle error signal. This will be a torque command. The torque between two magnetic dipole moments is inversely proportional to the relative distance raised to the 3rd power. Similar to above, we compensate for this by multiplying the error signal by |RelativePosition|.sup.3. To include damping, we multiply the angular velocity by a gain. The control law has the form of
pdSignalw=angleError kw|relativePosition|.sup.3+w kdw  (56)

(52) To solve for the μs, we let Fx=pdSignalx, Fy=pdSignaly and τc.sub.1z=pdSignalw, the control law can be written as

(53) $\begin{array}{cc}\left[\begin{array}{c}{\mu }_{11x}\\ {\mu }_{12y}\\ {\mu }_{13x}\\ {\mu }_{14y}\end{array}\right]={[\begin{array}{cc}{M}_{F_{11}}R\left({\theta }_1\right)I_x+M_{F_{12}}R\left({\theta }_1\right)I_y& M_{F_{13}}R\left({\theta }_1\right)I_x+M_{F_{14}}R\left({\theta }_1\right)I_y\\ M_{F\tau 11}R\left({\theta }_1\right)I_x+M_{F\tau 12}\left({\theta }_1\right)I_y& M_{F\tau 13}R\left({\theta }_1\right)I_x+M_{F\tau 14}R\left({\theta }_1\right)I_y\end{array}]}^{-1}\hspace{1em}\left[\begin{array}{c}{\mathrm{pdSignal}}_x\\ {\mathrm{pdSignal}}_y\\ \mathrm{pdSignalw}\end{array}\right]& \left(57\right)\end{array}$

APPENDIX. CONTROL SIGNAL FACTORED OUT OF FORCE EQUATION

(54) In this Appendix we introduce a matrix form of the force equation. Then we show how we write the force equation in this form. We then conclude with the general case of the force and torque acting on a nanosatellite from N−1 other nanosatellites, where N is the total number of nanosatellites.

(55) A typical form of the force acting on an electromagnet (EM) due to another electromagnet using the far field model is

(56) $\begin{array}{cc}{F}_1=-\frac{3{\mu }_0}{4\pi d^5}\left(\left({\mu }_1{\mathrm{•\mu }}_2\right)d+\left({\mu }_1•d\right){\mu }_2+\left({\mu }_2•d\right){\mu }_1-5\frac{\left({\mu }_1•d\right)\left({\mu }_2•d\right)}{d^2}d\right)& \left(61\right)\end{array}$
where
μ.sub.0 is the permeability of free space
μ.sub.1 is the magnetic dipole of EM 1
μ.sub.2 is the magnetic dipole of EM 2
d is the distance between the dipoles and
F.sub.1 is the force acting on EM 1

(57) The identity (x.Math.y)z=(zy).Math.x is used to write Eq. (61) as

(58) $\begin{array}{cc}{F}_1=-\frac{3{\mu }_0}{4\pi d^5}\left(\left(d.Math.{\mu }_2\right)+\left({\mu }_2.Math.d\right)+\left({\mu }_2^T•d\right)\left(I-5\frac{d.Math.d}{d^T•d}\right)\right).Math.{\mu }_1,& \left(62\right)\end{array}$
where is an outer product. Equation (62) can be written as a matrix times a vector, such as,
F.sub.1=M.sub.F(μ.sub.2,d).Math.μ.sub.1,  (63)
where

(59) $\begin{array}{cc}{M}_F\left({\mu }_2,d\right)=-\frac{3{\mu }_0}{4\pi d^5}\left(\left(d.Math.{\mu }_2\right)+\left({\mu }_2.Math.d\right)+\left({\mu }_2^T•d\right)\left(I-5\frac{d.Math.d}{d^T•d}\right)\right)& \left(64\right)\end{array}$
The value of writing the force equation in this form is that for a give force it is a simple process to solve for the control input μ.sub.1.
Re-Expression of the Force Equation.

(60) First, we will verify the identity (x.Math.y) z=(zy).Math.x, and then will show the substitution that gives the matrix form of the force equation. We first verify that
(x.Math.y)z=(zy).Math.x,  (65)
where

(61) $\begin{array}{cc}x=\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right],y=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right],\mathrm{and}z=\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right]& \left(66\right)\end{array}$
Multiplying out the left hand side of Eq. (65), we obtain

(62) $\begin{array}{cc}\left(x.Math.y\right)z=\left(x_1{y}_1+x_2{y}_2+x_3{y}_3\right)\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\end{array}\right]=\left[\begin{array}{c}{z}_1{y}_1{x}_1+z_1{y}_2{x}_2+z_1{y}_3{x}_3\\ z_2{y}_1{x}_1+z_2{y}_2{x}_2+z_2{y}_3{x}_3\\ z_3{y}_1{x}_1+z_3{y}_2{x}_2+z_3{y}_3{x}_3\end{array}\right]& \left(67\right)\end{array}$

(63) Multiplying out the right hand side we get the same vector. We have shown that
(x.Math.y)z=(zy).Math.x  (68)
Using the identity (68) to bring μ.sub.1 to the right side of each dot product, we have

(64) 0$\begin{array}{cc}{F}_1=-\frac{3{\mu }_0}{4\pi d^5}\left(\left(d.Math.{\mu }_2\right){\mu }_1+\left({\mu }_2.Math.d\right){\mu }_1+\left({\mu }_2•d\right){\mu }_1-5\frac{\left({\mu }_2•d\right)}{d2}\left(d.Math.d\right){\mu }_1\right)& \left(69\right)\end{array}$
Now factoring out μ.sub.1,

(65) $\begin{array}{cc}{F}_1=-\frac{3{\mu }_0}{4\pi d^5}\left(\left(d.Math.{\mu }_2\right)+\left({\mu }_2.Math.d\right)+\left({\mu }_2•d\right)-5\frac{\left({\mu }_2•d\right)}{d2}\left(d.Math.d\right)\right){\mu }_1& \left(70\right)\end{array}$
Then factoring (μ.sub.2□d), we have

(66) $\begin{array}{cc}{F}_1=-\frac{3{\mu }_0}{4\pi d^5}\left(\left(d.Math.{\mu }_2\right)+\left({\mu }_2.Math.d\right)+\left({\mu }_2^T•d\right)\left(I-5\frac{\left(d.Math.d\right)}{d^T•d}\right)\right).Math.{\mu }_1.& \left(71\right)\end{array}$
Expressions for M.sub.F and M.sub.T.

(67) The notation that will be used here is the following.

(68) μ.sub.ij is the magnetic dipole moment of the j.sup.th EM on the cube.

(69) d.sub.ijkm is distance vector from the center of the m.sup.th EM on the k.sup.th cube to the center of the j.sup.th EM on the cube.

(70) F.sub.ijkm is the force on the i.sup.th cube, j.sup.th EM due to the k.sup.th cube m.sup.th EM.

(71) τ.sub.ijkm is the torque on the i.sup.th cube, j.sup.th EM due to the k.sup.th cube m.sup.th EM.

(72) rc.sub.ij is the moment arm from the center of the i.sup.th cube to the center of the j.sup.th EM.

(73) Using the above notation, the force and torque equations become

(74) $\begin{array}{cc}{\mathrm{Fc}}_i=\underset{j=1}{\overset{4}{.Math.}}\underset{k=1,k\ne i}{\overset{\mathrm{Cubes}}{.Math.}}\underset{m=1}{\overset{4}{.Math.}}{F}_{\mathrm{ijkm}}& \left(72\right)\end{array}$

(75) $\begin{array}{cc}\tau c_i=\underset{j=1}{\overset{4}{.Math.}}\underset{k=1,k\ne i}{\overset{\mathrm{Cubes}}{.Math.}}\underset{m=1}{\overset{4}{.Math.}}{\tau }_{\mathrm{ijkm}}& \left(73\right)\end{array}$
where

(76) $\begin{array}{cc}{F}_{\mathrm{ijkm}}=-\frac{3{\mu }_0}{4\pi d_{\mathrm{ijkm}}^5}\left(\left({\mu }_{\mathrm{ij}}•{\mu }_{km}\right)d_{\mathrm{ijkm}}+\left({\mu }_{\mathrm{ij}}•d_{\mathrm{ijkm}}\right){\mu }_{km}+\left({\mu }_{km}•d_{\mathrm{ijkm}}\right){\mu }_{\mathrm{ij}}-5\frac{\left({\mu }_{\mathrm{ij}}•d_{\mathrm{ijkm}}\right)\left({\mu }_{km}•d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijkm}}^2}{d}_{\mathrm{ijkm}}\right)& \left(74\right)\\ {\tau }_{\mathrm{ijkm}}=\frac{{\mu }_0}{4\pi }\left({\mu }_{\mathrm{ij}}\times \left(\frac{3\left({\mu }_{km}.Math.d_{\mathrm{ijmk}}\right)}{d_{\mathrm{ijkm}}^5}{d}_{\mathrm{ijkm}}-\frac{{\mu }_{km}}{d_{\mathrm{ijkm}}^3}\right)\right)+{\mathrm{rc}}_{\mathrm{ij}}\times F_{\mathrm{ijmk}}& \left(75\right)\end{array}$
The above equations will be re-expressed to be linear in μ.sub.ij.

(77) Using the relationship (x.Math.y) z=(zy).Math.x=(zx).Math.y to factor out μ.sub.ij from the dot products, we obtain

(78) $\begin{array}{cc}{F}_{\mathrm{ijkm}}=-\frac{3{\mu }_0}{4\pi d_{\mathrm{ijkm}}^5}\left(\left(d_{\mathrm{ijkm}}.Math.{\mu }_{km}\right)+\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)+\left({\mu }_{km}•d_{\mathrm{ijkm}}\right)-5\frac{\left({\mu }_{km}•d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijkm}}^2}\left(d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijkm}}\right)\right){\mu }_{\mathrm{ij}}& \left(76\right)\end{array}$
Factoring out (μ.sub.km□d.sub.ijkm) we obtain

(79) $\begin{array}{cc}{F}_{\mathrm{ijkm}}=-\frac{3{\mu }_0}{4\pi d_{\mathrm{ijkm}}^5}\left(\left(d_{\mathrm{ijkm}}.Math.{\mu }_{km}\right)+\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)+\left({\mu }_{km}•d_{\mathrm{ijkm}}\right)\left(I-5\frac{\left(d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijkm}}\right)}{\left(d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijmk}}\right)}\right)\right){\mu }_{\mathrm{ij}}& \left(77\right)\end{array}$
We define the matrix

(80) $\begin{array}{cc}{M}_F\left({\mu }_{km},d_{\mathrm{ijkm}}\right)=-\frac{3{\mu }_0}{4\pi d_{\mathrm{ijkm}}^5}\left(\left(d_{\mathrm{ijkm}}.Math.{\mu }_{km}\right)+\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)+\left({\mu }_{km}•d_{\mathrm{ijkm}}\right)\left(I-5\frac{\left(d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijkm}}\right)}{\left(d_{\mathrm{ijkm}}.Math.d_{\mathrm{ijkm}}\right)}\right)\right)& \left(78\right)\end{array}$
The force equation is now written as
F.sub.ijkm=M.sub.F(μ.sub.km,d.sub.ijkm)μ.sub.ij  (79)

(81) Now factor out μ.sub.ij from the torque equation. The cross product can be written in matrix form as

(82) $\begin{array}{cc}a\times b={\left[a\right]}_{x}b=\left[\begin{array}{ccc}0& -a_3& a_2\\ a_3& 0& -a_1\\ -a_2& a_1& 0\end{array}\right]\left[\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right]={\left[b\right]}_x^{T}a=\left[\begin{array}{ccc}0& b_3& -b_2\\ -b_3& 0& b_1\\ b_2& -b_1& 0\end{array}\right]\left[\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right]& \left(80\right)\end{array}$
We obtain

(83) 0$\begin{array}{cc}{\tau }_{\mathrm{ijkm}}={\frac{{\mu }_0}{4\pi }\left[\frac{3\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijmk}}^5}{d}_{\mathrm{ijkm}}-\frac{{\mu }_{km}}{d_{\mathrm{ijkm}}^3}\right]}_x^T{\mu }_{\mathrm{ij}}+{\mathrm{rc}}_{\mathrm{ij}}\times F_{\mathrm{ijkm}}& \left(81\right)\end{array}$
Substituting M.sub.F for F.sub.ijkm, we obtain

(84) $\begin{array}{cc}{\tau }_{\mathrm{ijkm}}=\left({\frac{{\mu }_0}{4\pi }\left[\frac{3\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijkm}}^5}{d}_{\mathrm{ijkm}}-\frac{{\mu }_{km}}{d_{\mathrm{ijkm}}^3}\right]}_x^T+{\left[{\mathrm{rc}}_{\mathrm{ij}}\right]}_x{M}_F\left({\mu }_{km},d_{\mathrm{ijkm}}\right)\right){\mu }_{\mathrm{ij}}.& \left(82\right)\end{array}$
We define

(85) $\begin{array}{cc}{M}_{\tau }\left({\mu }_{km},d_{\mathrm{ijkm}}\right)={\frac{{\mu }_0}{4\pi }\left[\frac{3\left({\mu }_{km}.Math.d_{\mathrm{ijkm}}\right)}{d_{\mathrm{ijkm}}^5}{d}_{\mathrm{ijkm}}-\frac{{\mu }_{km}}{d_{\mathrm{ijkm}}^3}\right]}_x^T+{\left[{\mathrm{rc}}_{\mathrm{ij}}\right]}_x{M}_F\left({\mu }_{km},d_{\mathrm{ijkm}}\right),& \left(83\right)\end{array}$
τ.sub.ijkm=M.sub.τ(μ.sub.km,d.sub.ijkm)μ.sub.ij.  (84)

(86) Equations (79) and (84) are referred to as a matrix form of the force and torque equations respectively.