Abstract:
A method and system for locating the three dimensional coordinates of a moving or stationary object in real time. The three dimensional coordinates of an object in half space or full space are determined based upon the time of arrival or phase of the wave front measured by a plurality of receiver elements and an established vector magnitudes proportional to the measured time of arrival or phase at each receiver element. The coordinates of the object are calculated by solving a matrix equation or a set of closed form algebraic equations.

Description:
CONTRACTUAL ORIGIN OF THE INVENTION 
     The United States Government has rights in this invention pursuant to Contract No. DE-AC07-76ID01570 between the U.S. Department of Energy and EG&amp;G Idaho, Inc. 
    
    
     This is a continuation-in-part of application Ser. No. 08/069,647 filed Jun. 1, 1993 now abandoned, which was a continuation-in-part of application Ser. No. 07/896,867 filed Jun. 12, 1992, now abandoned. 
    
    
     BACKGROUND OF THE INVENTION 
     This invention relates to a method for the remote tracking in real-time of the spatial coordinates of a stationary or moving object in three dimensional space. 
     Interferometry methods for the location of an emitter source are known which require the knowledge of several parameters such as, waveform modulation, emitter signal angle of arrival or time of arrival, or measurements of angular rates between the emitter source and the measurements site. It is known in the art to use interferometry and tracking systems which utilized amplitude and phase comparison relationships to determine either the direction or spatial coordinates of a remote object relative to that of the receiver or an arbitrary reference. Conventional direction and range tracking systems typically utilize planar propagating electromagnetic waves. Conventional planar wave propagating systems become less accurate as the distance between the object under observation and the receiver decreases. For close range applications (which are generally considered to be where the location of the object relative to the sensory array is less than one hundred times the largest distance between sensor elements) the margin of error associated with conventional techniques is often unacceptable. Also, additional expense is required if spatial coordinates must be determined in addition to direction finding capability. 
     In order to obtain the actual coordinates of the emitter location using conventional plane wave propagating techniques, either triangulation methods, or range and direction detection methods are commonly used. Often, unacceptable system costs and limitations in tracking may result from additional measurement time requirements in such systems. If three dimensional rather than two dimensional tracking is required, time and cost requirements can be amplified even further. Due to the relatively large physical dimensions of interferometer sensors and triangulation baseline, serious limitations in performance and application often occur with the use of the planar wave interferometry. When the distance separating the receiver and the object decreases to a point where the accuracy of the tracking measurements is degraded to an unacceptable state, a different technique must be considered. 
     Optical and sub-optical frequency electromagnetic (EM) wave energies as well as acoustical wave energies, have been most commonly used. Systems such as optical interferometers utilized for direction tracking require mechanical servo mechanisms and are not fault tolerant, in that realignment and calculation of the target coordinates after power failure is required. 
     U.S. Pat. No. 4,788,548 of Hammerquist discloses a passive radar ranging system capable of determining the two dimensional range of an emitter source using four co-linear antenna receiver elements. However, the technique of Hammerquist is limited to two dimensional range determination. 
     It is an object of this invention to provide a method for the fault tolerant, continuous, real-time, high resolution determination of three dimensional spatial coordinates of objects. 
     It is another object of this invention to provide a method for determining in real-time close range three dimensional coordinates of objects using a broad spectrum of frequencies. 
     It is still another object of this invention to provide a method for determining the three dimensional full or half space coordinates of an object when the waveform is unmodulated, randomly modulated, or pulsed. 
     It is still a further object of this invention to provide a method for determining spatial coordinates or tracking of an object at low frequencies to penetrate poor climate and atmospheric conditions. 
     Additional objects, advantages and novel features of the invention will become apparent to those skilled in the art upon examination of the following and by practice of the invention. 
     SUMMARY OF THE INVENTION 
     To achieve the foregoing and other objects, this invention comprises a method and system for locating the three dimensional coordinates of a moving or stationary object in real-time, based upon the known geometric configuration of a propagating wave front. In one embodiment, the three dimensional coordinates of an object in half space are established based upon the phase of the wave front measured by a minimum of four receiver elements positioned in a planar antenna array. In another embodiment, the three dimensional coordinates in full space are calculated from the phase of the wave front measured by a minimum of five receiver elements in an configuration having no more that three sensors in a plane. The system can operate using optical and sub-optical frequency electromagnetic wave energies as well as acoustical wave energies. Additionally, the system provides exceptionally accurate tracking of objects at close range, exceeding that of conventional planar wave techniques. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The present invention is illustrated in the accompanying drawings where: 
     FIG. 1 is a graphical illustration of the position of the emitter and the antenna elements for the three dimensional half space coordinate determination or tracking. 
     FIG. 2 is a schematic block diagram of the half space coordinate determination or tracking system. 
     FIG. 3 is a graphical illustration showing the position of the emitter and the antenna elements for the three dimensional full space coordinate determination or tracking. 
    
    
     DETAILED DESCRIPTION OF THE INVENTION 
     The present invention utilizes a propagating wave emanating from a source. In a preferred embodiment of the present invention for two dimensional tracking, the wave front is propagating radially outward from the object under observation. For three dimensional tracking, it is assumed that the wave front propagates spherically outward from the object. In a preferred embodiment, the present invention utilizes, but is not limited to signal frequencies for acoustic waves (longitudinal) of between 200 Hz to 50 kHz, and electromagnetic (traverse) radio frequencies of between 27 MHz to 23 GHz. The present invention is also capable of utilizing pulsed or single event emission sources. 
     Referring now to FIG. 1, the configuration of the present invention is shown for the coordinate determination or tracking of an object in three dimensional half space. The small circles in FIG. 1 represent sensors elements of either longitudinal or transverse waves. The small filled circles identify the location where a ray cast from the source through the sensor would pierce a sphere representing a feature in the propagating wave front. To track the position or coordinates of an object in three dimensional half space, a planar antenna array consisting of a minimum of four elements, depicted in FIG. 1 as elements A, B, C and D, are used for the receiver. In the simplest configuration, all of the receiver elements are located in the XY plane and are spatially offset from one another. An emitter, represented as E and located at coordinates Sx, Sy, Sz, is shown in FIG. 1 as emitting a propagating spherical equiphase wave front. The emitter (E) and receiver elements (A, B, C, and D) are positioned with respect to each other in the Cartesian coordinate system with receiver element D located at the origin. FIG. 1 shows a graphical representation of a spherical wave front. For one instant in time, the spherical wave front is shown passing through element D and is represented by the single spherical surface intersecting receiver element D. Unequal signal propagation times between the emitter and each of the receiver elements results in signal phase differences, which are referenced to the propagating signal wave front arriving at receiver element D. With this antenna array configuration, all three Cartesian coordinates Sx, Sy, and Sz of the emitter location may be determined. Three phase difference measurements are used to determine the emitter coordinates. In FIG. 1, a, b, and c are electrical lengths corresponding to the phase differences between receiver element D and receiver element locations A, B, and C, respectively. These radii dimensions dynamically change as the coordinates of the emitter change. 
     The relationship between the signal phase and lengths a, b, and c are expressed as ##EQU1## where λ is the free space wavelength, ψ a , ψ b , and ψ c  are the three phase differences measured in radians, Δt a , Δt b , and Δt c  are the measured differences in times of arrival for the wave front at each sensor element with respect to that at element D, and υ is the velocity of the propagating wave. Equation (1) illustrates the relationship between the time and phase differences of the wave (i.e., either one may be calculated from the measured value of the other). 
     Four vector quantities, EA, EB, EC, and ED are used to describe the vectors originating at the emitter and extending to each of the array elements. Vector loop relationships can be derived where the vector path begins and ends at the same point. A vector extending from the emitter to an array element, from the array element to the wave front tangency point and from the tangency point back to the emitter, is an example of a vector loop. The expressions representing the configuration in FIG. 1 are: 
     
         EA=(Ax-Sx)i+(Ay-Sy)j-(Az-Sz) k=(R+a)ea                     (2) 
    
     
         EB=(Bx-Sx)i+(By-Sy)j-(Bz-Sz) k=(R+b)eb                     (3) 
    
     
         EC=(Cx-Sx)i+(Cy-Sy)j-(Cz-Sz) k=(R+c)e                      (4) 
    
     and 
     
         ED=-Sxi-Syj-Szk=(R)ed                                      (5) 
    
     The terms Ax, Bx, Cx, Ay, By, Cy, Az, Bz, and Cz are the X, Y and Z spatial displacements of elements A, B, and C relative to the reference element D. The terms ea, eb and ec are the direction unit vectors extending from the emitter to receiver elements A, B, and C, respectively. To solve for the four unknowns, Sx, Sy, Sz and R, only the vector magnitudes need to be considered and the equations may be expressed as follows: 
     
         |EA|.sup.2 =(Ax-Sx).sup.2 +(Ay-Sy).sup.2 +(Az-Sz) .sup.2 =(R+a).sup.2                                              (6) 
    
     
         |EB|.sup.2 =(Bx-Sx).sup.2 +(By-Sy).sup.2 +(Bz-Sz) .sup.2 =(R+b).sup.2                                              (7) 
    
     
         |EC|.sup.2 =(Cx-Sx).sup.2 +(Cy-Sy).sup.2 +(Cz-Sz) .sup.2 =(R+c).sup.2                                              (8) 
    
     and 
     
         |ED|.sup.2 =Sx.sup.2 +Sy.sup.2 +Sz.sup.2 =(R).sup.2(9) 
    
     For planar array configurations, equations (6) through (9) are used to derive the following set of linear equations which describe the three unknowns, R, Sx, and Sy: ##EQU2## 
     A solution for Sz in terms of Sx, Sy, and R is now obtained by using equation 9 and is expressed as ##EQU3## and only the positive real root is acceptable. The position of E is assumed to always reside in the positive half space, therefore, Sz is always positive. However, if the sensing array is not planar then both positive and negative half-spaces are resolved and Sz may be either positive of negative. From the foregoing, it can be seen that the solution to the coordinates of the emitter object in three dimensional Cartesian space is given by Equations (10) and (11). Note that from Equations (10) and (11) it is not necessary to calculate the range of the emitter from the reference in order to solve for the x and y coordinates of the emitter source, thereby providing an instantaneous method for calculating the two dimensional coordinates of the emitter source. 
     For rapid calculation and efficient implementation of the solution, Equations (10) and (11) can be solved in the closed form as represented by the following Equations (12) through (28). ##EQU4## 
     The implementation of the three dimensional half space embodiment of the present invention is shown in FIG. 2. This embodiment consists of four antenna elements designated 12, 14, 16 and 18 located at points A, B, C and D of FIG. 1 respectively. Antenna element D is established as a reference location. The antenna elements are aligned in a planar array and detect either a longitudinal or transverse wave front. Additionally, the wave front can be a unmodulated, randomly modulated, or pulsed energy signal. 
     The energy collected at antenna elements 12, 14, 16 and 18 is fed to a phase converter system 20 to be translated into time or distance. Energy collected at antenna elements 12, 14, and 16 is transferred to a corresponding plurality of phase detectors 22, 24, and 26 while the energy collected at the reference antenna element 18 is also input into phase detectors 22, 24, and 26. The output of the phase detectors is transmitted to a plurality of phase/distance converters 32, 34, and 36 respectively. The phase/distance converters convert the phase to either a distance parameter, or alternately to a time parameter, as shown in Equation (1). 
     Output from the phase/distance converters is then transmitted to the position processor 40 for processing the information as described in the closed form solutions set forth in Equations (12) through (28). Alternately, using a digital computer, the position processor could be used to solve either Equations (10) and (11) or the closed form Equations (12) through (28). However, this process is slightly more time consuming. For rapid, real time calculations, the closed form solution is preferred for hardware implementation. As shown in FIG. 2, the output from each phase/distance converter is transmitted to multipliers 42, 44, 46 and 48. The sensor elements y-coordinates, Ay, By, and Cy, as well as constants K 1 , K 2  and K 3  are also provided to multiplier 42. Constants K 4 , K 5 , and K 6  are provided to multiplier 44. The sensor elements x-coordinates Ax, Bx and Cx, and constants K 7 , K 8 , and K 9  are provided to multiplier 46. Multiplier 48 inputted information consists of constants K 10 , K 11  and K 12 . The multipliers 42, 44, 46, and 48 perform the multiplication function of Equations (12) through (28) and the output of the multipliers is fed into a corresponding number of adders, denoted in FIG. 2 as 52, 54, 56 and 58 respectively. Adders 52, 54, 56, and 58 perform the addition differencing functions of Equations (12) through (28). 
     Adder 54 provides the denominator for the Sx, Sy and R calculations and the denominator information is transmitted to dividers 60, 62 and 64. The numerator data for dividers 60, 62 and 64 is provided from the output of adder 52, 56 and 58 respectively. Dividers 60, 62, and 64 perform the division function of Equations (12) through (14) and yield the solutions for Sx, Sy and R respectively. The solutions for Sx, Sy and R is then transmitted to multiplier 66 which performs the squaring functions of Equation (15) while the subsequent addition is performed by adder 68. Finally the Sz coordinate is calculated by taking the square root of the sum (R 2  -Sx 2  -Sy 2 ). Sx, Sy and Sz represent the x, y and z coordinates of the object. 
     Referring now to FIG. 3, the configuration of the present invention is shown for the coordinate determination or tracking of an object in three dimensional full space. To track the position or coordinates of an object in three dimensional full space, a non-planar array consisting of a minimum of five elements, depicted in FIG. 3 as elements A, B, C, D, and E, are used for the antenna receiver. 
     An emitter, represented as S and located at coordinates X s , Y s , and Z s , is shown in FIG. 3 as emitting a propagating spherical wave front. Using the five element sensor array, a closed form solution for finding the location of the emanating point energy source emitting a transverse or longitudinal wave front can be derived. For purposes of the derivation it is assumed that element D is a reference sensor to which all other sensor elements may be compared. The small circles in FIG. 3 represent sensors of either longitudinal or transverse waves. The small filled circles identify the location where a ray cast from the source through the sensor would pierce a sphere representing a feature in the propagating wave front. In some cases the ray pierces the sphere before reaching the sensor (e.g., elements A, B, and C), and in others it pierces the sphere after passing the sensor (e.g., element E). The following sensor position vectors can be written for each sensor element: 
     
         A=Axi+Ayj+Azk                                              (29) 
    
     
         B=Bxi+Byj+Bzk                                              (30) 
    
     
         C=Cxi+Cyj+Czk                                              (31) 
    
     
         D=Dxi+Dyj+Dzk                                              (32) 
    
     
         E=Exi+Eyj+Ezk                                              (33) 
    
     The source position vector is defined as: 
     
         S=Sxi+Syj+Szk                                              (34) 
    
     From these position vectors, the following relative positions vectors for the sensors with respect to the energy source can be written: 
     
         SA=(Ax-Sx)i+(Ay-Sy)j+(Az-Sz)k                              (35) 
    
     
         SB=(Bx-Sx)i+(By-Sy)j+(Bz-Sz)k                              (36) 
    
     
         SC=(Cx-Sx)i+(Cy-Sy)j+(Cz-Sz)k                              (37) 
    
     
         SD=(Dx-Sx)i+(Dy-Sy)j+(Dz-Sz)k                              (38) 
    
     
         SE=(Ex-Sx)i+(Ey-Sy)j+(Ez-Sz)k                              (39) 
    
     The distance from the propagating energy front or wave to each of the sensors as the front passes the reference sensor D is denoted by a lower case letter which is the same as the upper designation for each sensor (i.e., a goes with A). These lengths are measured values that carry a sign with them. Equations 40 through 43 show how the lengths are calculated from phase differences or differences in time of arrival with respect to the reference sensor. ##EQU5## where λ is the free space wavelength, ψa, ψb, ψc and ψe are the three phase differences measured in radians, αt a , αt b , αt c  and αt e  are the measured time of arrival for the wave front at each sensor element, and υ is the velocity of the propagating wave. Equations (40) through (43) illustrate the relationship between the time and phase differences of the wave (i.e., either one may be calculated from the measured value of the other). 
     From Equations (40) through (43) and the vectors shown in FIG. 3, a set of vector loop equations can be written. The loops are degenerate in the sense that the vectors are collinear or have the same unit vector. The assumptions is made that positiveness is in the direction of energy propagation except for the measured quantities which are signed values. The signed values are summed in the following vector loop equations: 
     
         |SA|sa+asa-|SD|sa=0    (44) 
    
     
         |SB|sb+bsb-|SD|sb=0    (45) 
    
     
         |SC|sc+csc-|SD|sc=0    (46) 
    
     
         |SE|se+ese-|SD|se=0    (47) 
    
     From these vector loop expressions, the magnitude expressions are as follows: 
     
         |SA|+a-|SD|=0          (48) 
    
     
         |SB|+b-|SD|=0          (49) 
    
     
         |SC|+c-|SD|=0          (50) 
    
     
         |SE|+e-|SD|=0          (51) 
    
     
         |Sd=R                                             (52) 
    
     By expanding the magnitude equations and simplifying, the following results can be derived: 
     
         [(Ax-Sx).sup.2 +(Ay-Sy).sup.2 +(Az-Sz).sup.2 ].sup.1/2 +a-[(Dx-Sx).sup.2 +(Dy-Sy).sup.2 +(Dz-Sz).sup.2 ].sup.1/2 =0                (53) 
    
     
         (Ax-Sx).sup.2 +(Ay-Sy).sup.2 +(Az-Sz).sup.2 =(Dx-Sx).sup.2 +(Dy-Sy).sup.2 +(Dz-Sz).sup.2 -2a[(Dx-Sx).sup.2 +(Dy-Sy).sup.2 +(Dz-Sz).sup.2 ].sup.1/2 +a.sup.2                                                  (54) 
    
     In equation (53), the radical term can be replaced with R which is the unknown range or distance from the reference sensor to the source. By combining terms and reorganizing, the following equations with unknowns S x , S y , S z  and R can be written: 
     
         Ax.sup.2 -2AxSx+Sx.sup.2 +Ay.sup.2 -2AySy +Sy.sup.2 +Az.sup.2 -2AzSz+Sz.sup.2 =Dx.sup.2 2DxSx+Sx.sup.2 +Dy.sup.2 -2DySy+Sy.sup.2 +Dz.sup.2 -2DzSz+Sz.sup.2 -2aR+a.sup.2                    (55) 
    
     
         Ax.sup.2 -2AxSx+Ay.sup.2 -2AySy+Az.sup.2 -2AzSz=Dx.sup.2 -2DxSx+Dy.sup.2 -2DySy+Dz.sup.2 -2DzSz-2aR+a.sup.2                        (56) 
    
     
         2(Dx-Ax)Sx+2(Dy-Ay)Sy+2(Dz-Az)Sz+2aR=Dx.sup.2 +Dy.sup.2 +Dz.sup.2 -Ax.sup.2 -Ay.sup.2 -Az.sup.2 +a.sup.2                              (57) 
    
     Expanding Equations (55) through (57) similarly yields the following equations respectively: 
     
         2(Dx-Bx)Sx+2(Dy-By)Sy+2(Dz-Bz)Sz+2bR=Dx.sup.2 +Dy.sup.2 +Dz.sup.2 -Bx.sup.2 -By.sup.2 -Bz.sup.2 +b.sup.2                              (58) 
    
     
         2(Dx-Cx)Sx+2(Dy-Cy)Sy+2(Dz-Cz)Sz+2cR=Dx.sup.2 +Dy.sup.2 +Dz.sup.2 -Cx.sup.2 -Cy.sup.2 -Cz.sup.2 +c.sup.2                              (59) 
    
     
         2(Dx-Ex)Sx+2(Dy-Ey)Sy+2(Dz-Ez)Sz+2eR=Dx.sup.2 +Dy.sup.2 +Dz.sup.2 -Ex.sup.2 -Ey.sup.2 -Ez.sup.2 +e.sup.2                              (60) 
    
     From Equations (58) through (60), a matrix representation of the closed form solutions may be written and computer coded. Equation (61) is the matrix representation for five sensors in an arbitrary configuration. ##EQU6## 
     Note that in Equation (61) it is not necessary to calculate the range of the emitter from the reference in order to solve for the x, y and z coordinates of the emitter source. It is also important to note that a closed form solution for Equation (61) can be used to solve for the three dimensional coordinates of an object in full space. As in FIG. 2, the closed form solution for the full space three dimension coordinates may be implemented. Appendix A provides an example of the closed form solution to Equation (61). 
     As in conventional interferometry, the phase difference measurement for either the two or three dimensional tracking, is unambiguous only if the phase angle ψ does not exceed 2π. Several known methods may be used to solve the phase ambiguity problems, such as using multiple continuous wave signals differing slightly in frequency or signal modulation techniques. 
     The foregoing description of a preferred embodiment of the invention has been presented for purposes of illustration and description. It is not intended to be exhaustive or to limit the invention to the precise form disclosed, and obviously many modifications and variations are possible in light of the above teaching. The embodiments described explain the principles of the invention and practical application and enable others skilled in the art to utilize the invention in various embodiments and with various modifications as are suited to the particular use contemplated. It is intended that the scope of the invention be defined by the claims appended hereto. 
     
                                           APPENDIX A__________________________________________________________________________Sx = - 1/2{(((Dz - Cz) Ey + (Cy - Dy) Ez + Dy Cz - Cy Dz) b +((-Dz + Bz) Ey + (Dy - By) Ez + By Dz - Dy Bz)c + %l e)a.sup.2 +(((Dz - Az) Ey + (-Dy + Ay) Ez + Dy Az - Ay Dz) c + %8 e) b.sup.2 + {((Cz - Dz) Ey + (Dy - Cy) Ez - Dy Cz + Cy Dz) b.sup.2 +((Dz - Bz) Ey + (By - Dy) Ez + Dy Bz - By Dz) c.sup.2 + %2 e.sup.2 + %,Ex.sup.2 +%1 Ey.sup.2 + {((Dz - Cz) Bx.sup.2 + (Dz - Cz) By.sup.2 + (Dz - Cz)Bz.sup.2 + %5 Bz - Dz Cx.sup.2 -Dz Cy.sup.2 - Dz Cz.sup.2 + %4 Cz}Ey + %1 Ez.sup.2 + {(Cy - Dy)Bx.sup.2 +(Cy - Dy)By.sup.2 +%6 By + (Cy - Dy) Bz.sup.2 + Cx.sup.2 Dy + Dy Cy.sup.2 + %3 Cy + DyCz.sup.2 }Ez +(Dy Cz - Cy Dz) Bx.sup.2 + (Dy Cz - Cy Dz) By.sup.2 +{Dz Cx.sup.2 + Dz Cy.sup.2 + Dz Cz.sup.2 + %3 Cz}By + (Dy Cz - Cy Dz)Bz.sup.2 +{-Cx.sup.2 Dy - Dy Cy.sup.2 + %4 Cy - Dy Cz.sup.2 }Bz) a + {((-Dz + Az) Ey + (Dy - Ay) Ez - Dy Az + Ay Dz) c .sup.2 +((Cz - Dz) Ay + (Dy -  Cy) Az - Dy Cz + Cy Dz) e.sup.2 + %8 Ex.sup.2 + %8Ey.sup.2 + {Dz Cz.sup.2 + %3 Cz) Ey + %8 Ez.sup.2 + ((Dy - Cy) Ax.sup.2 + (Dy - Cy)Ay.sup.2 + %5 Ay +(Dy - Cy) Az.sup.2 - Cx.sup.2 Dy - Dy Cy.sup.2 + %4 Cy - Dy Cz.sup.2) Ez(-Dy Cz + Cy Dz) Ax.sup.2 + (-Dy Cz + Cy Dz) Ay.sup.2 +(-Dz Cx.sup.2 - Dz Cy.sup.2 - Dz Cz.sup.2 + %4 Cz) Ay + (-Dy Cz + Cy Dz)Az.sup.2 +(Cx.sup.2 Dy Cy.sup.2 + %3 Cy + Dy Cz.sup.2 )Az)b + %7 e c.sup.2 + (((Dz - Bz)Ay + (By - Dy)Az +  Dy Bz - By Dz)e.sup.2 + %7 Ex.sup.2 + %7Ey.sup.2 + ((Dz - Bz) Ax.sup.2 + (Dz - Bz) Ay.sup.2 + (Dz - Bz) Az.sup.2 +(By.sup.2 + Bz.sup.2 - Dy.sup.2 - Dz.sup.2 - Dxhu 2 + Bx.sup.2) Az -Bx.sup.2 Dz - By.sup.2 Dz - Bz.sup.2 Dz - Bz.sup.2 Dz +%4 Bz) Ey + %7 Ez.sup.2 + ((By - Dy) Ax.sup.2 + (By - Dy) Ay.sup.2 +(Dy.sup.2 - By.sup.2 - Bz.sup.2 - Bx.sup.2 + Dz.sup.2 + Dx.sup.2) Ay +(By - Dy) Az.sup.2 + Dy Bx.sup.2 + Dy By.sup.2 +%3 By + Dy Bz.sup.2) Ez + (Dy Bz -  By Dz) Ax.sup.2 + (Dy Bz - By Dz)Ay.sup.2 +(Bx.sup.2 Dz + By.sup.2 Dz + Bz.sup.2 Dz + %3 Bz) Ay + (Dy Bz - By Dz)Az.sup.2 +(-Dy Bx.sup.2 - Dy By.sup.2 + %4 By - Dy Bz.sup.2) Az) c + (%2 Ax.sup.2 +%2 Ay.sup.2 +(Cz - Dz) Bx.sup.2 + (Cz - Dz)By.sup.2 + (Cz- Dz)Bz.sup.2 + %6 Bz + DzCx.sup.2 + Dz Cy.sup.2 +Dz Cz.sup.2 + %3 Cz) Ay + %2 Az.sup.2 + ((Dy - Cy) Bx.sup.2 + (Dy - Cy)Byhu 2 + %5 By +(Dy - Cy) Bz.sup.2 - Cx.sup.2 Dy - Dy Cy.sup.2 + %4 Cy - Dy Cz.sup.2) Az+(-Dy Cz + Cy Dz) Bx.sup.2 + (-Dy Cz + Cy Dz) Byhu 2 +(-Dz Cx.sup.2 - Dz Cy.sup.2 - Dz Cz.sup.2 + %4 Cz) By + (-Dy Cz + Cy Dz)Bz.sup.2 +(Cx.sup.2 Dy + Dy Cy.sup.2 + %3 Cy + Dy Cz.sup.2)Bz) e)/((%2 Ex +((Cz - Dz) Bx + (Dx - Cx) Bz - Dx Cz + Cx Dz)Ey +((Dy - Cy) Bx + (-Dx + Cx)By + Dx Cy - Cx Dy) Ez + (-Dy Cz + Cy Dz) Bx +(-Cx Dz + Dx Cz) By + (Cx Dy - Dx Cy) Bz) a . . . +(((Cz - Dz) Ay + (Dy - Cy) Az - Dy Cz + Cy Dz) Ex +((Dz - Cz) Ax + (-Dx + Cx) Az -  Cx Dz + Dx Cz) Ey +((Cy - Dy) Ax + (Dx - Cx) Ay + Cx Dy - Dx Cy) Ez + (Dy Cz - Cy Dz) Ax +(-Dx Cz + Cx Dz) Ay + (Dx Cy - Cx Dy) Az) b +(((Dz - Bz) Ay + (By - Dy) Az + Dy Bz - By Dz) Ex +((-Dz + Bz) Ax + (Dx - Bx) Az - Dx Bz + Bx Dz) Ey +((Dy -By) Ax + (Bx - Dx) Ay + Dx By - Bx Dy) Ez + (By Dz - Dy Bz) Ax +(Dx Bz - Bx Dz) Ay + (-Dx By + Bx Dy) Az) c + (%l Ax +((Dz - Cz) Bx + (-Dx + Cx) Bz - Cx Dz + Dx Cz) Ay +((Cy - Dy) Bx + (Dx - Cx) By + Cx Dy - Dx Cy) Az + (Dy Cz - Cy Dz) Bx +(-Dx Cz + Cx Dz) By + (Dx Cy - Cx Dy) Bz) e)%I: = (Cz - Dz)By + (Dy - )Bz - Dy Cz + Cy Dz%2: = (Dz - Cz) By + (Cy - Dy) Bz + Dy Cz - Cy Dz%3: = -Dx.sup.2 - Dz.sup.2 - Dy.sup.2%4: = Dy.sup.2 + Dz.sup.2 + Dx.sup.2%5: = -Dx.sup.2 - Dz.sup.2 + Cy.sup.2 - Dy.sup.2 + Cz.sup.2 + Cx.sup.2%6: = Dy.sup.2 + Dx.sup.2 - Cz.sup.2 - Cx.sup.2 + Dz.sup.2 - Cy.sup.2%7: = (-Dz + Bz) Ay + (Dy - By) Az + By Dz - Dy Bz%8: =  (Dz - Cz) Ay + (Cy - Dy) Az + Dy Cz - Cy DzSy =  1/2((((Dz - Cz) Ex + (-Dx + Cx) Ez - Cx Dz + Dx Cz) b +((-Dz + Bz)Ex + (Dx - Bx)Ez - Dx Bz + Bx Dz) c + %4 e) a.sup.2 +(((Dz - Az) Ex + (-Dx + Ax) Ez + Az Dx - Dz Ax) c + %3 e) b.sup.2 +(((Cz - Dz) Ex + (Dx - Cx) Ez - Dx Cz + Cx Dz) b.sup.2 +((Dz - Bz) Ex + (Bx - Dx) Ez + Dx Bz - Bx Dz) c.sup.2 + %l e.sup.2 + %4Ex.sup.2 +((Dz - Cz) Bx.sup.2 + (Dz - Cz) By.sup.2 + (Dz - Cz) Bz.sup.2 +  %7 Bz -Dz CX.sup.2 - Dz cy.sup.2 -Dz Cz.sup.2 + %6 Cz} Ex + %4 Ey.sup.2 + %4 Ez.sup.2 + {(-Dx + Cx)Bx.sup.2 + %8 Bx +(-Dx + Cx) By.sup.2 + (-Dx + Cx) Bz.sup.2 + Cx.sup.2 Dx + %5 Cx + DxCz.sup.2 + Dx Cy.sup.2 } Ez +(-Cx Dz + Dx Cz) Bx.sup.2 + {Dz Cx.sup.2 + Dz Cy.sup.2 + Dz Cz.sup.2 + %5Cz} Bx +(-Cx Dz + Dx Cz) By.sup.2 + (-Cx Dz + Dx Cz) Bz.sup.2 +{-Cx.sup.2 Dx + %6 Cx - Dx Cy.sup.2 - Dx Cz.sup.2 } Bz} a +{((-Dz + Az) Ex + (Dx -  Ax) Ez + Dz Ax - Az Dx) c.sup.2 +((Cz - Dz) Ax + (Dx - Cx) Az - Dx Cz + Cx Dz) e.sup.2 + %3 Ex.sup.2 +{(Cz - Dz) Ax.sup.2 + (Cz - Dz) Ay.sup.2 + (Cz - Dz) Az.sup.2 + %8 Az +Dz Cx.sup.2 + Dz CY.sup.2 +Dz Cz.sup.2 + %5 Cz} Ex + %3 Ey.sup.2 + %3 Ez.sup.2 + {(Dx - Cx) Ax.sup.2+ %7 Ax +(Dx - Cx) Ay.sup.2 + (Dx - Cx) Az.sup.2 - Cx.sup.2 Dx + %6 Cx - DxCy.sup.2 - Dx Cz.sup.2 } Ez +. . .(-Dx Cz + Cx Dz) Ax.sup.2 + {-Dz CX.sup.2 - Dz Cy.sup.2 - Dz Cz.sup. 2 +%6 Cz} Ax +(-Dx Cz + Cx Dz) Ay.sup.2 + (-Dx Cz + Cx Dz) Az.sup.2 +{Cx.sup.2 Dx + %5 Cx + Dx Cz.sup.2 + Dx Cy.sup.2 } Az} b + %2 e C.sup.2+{((Dz - Bz) Ax + (Bx - Dx) Az + Dx Bz - Bx Dz) e.sup.2 + %2 Ex.sup.2 +{(Dz - Bz) Ax.sup.2 +(Dz - Bz) Ay.sup.2 + (Dz - Bz) Az.sup.2 + {By.sup.2 + Bz.sup.2 - Dy.sup.2- Dz.sup.2 - Dx.sup.2 + Bx.sup.2 } Az -Bx.sup.2 Dz - By.sup.2 Dz - Bz.sup.2 Dz + %6 Bz} Ex + %2 Ey.sup.2 + %2Ez.sup.2 + {(Bx -  Dx) Ax.sup.2 +{Dy.sup.2 - By.sup.2 - Bz.sup.2 - Bx.sup.2 + Dz.sup.2 } Ax + (Bx - Dx)Ay.sup.2 + (Bx - Dx) Az.sup.2 +Dx Bx.sup.2 + %5 Bx + Dx By.sup.2 + Dx Bz.sup.2 } Ez + (Dx Bz - Bx Dz)Ax.sup.2 +{Bx.sup.2 Dz + By.sup.2 Dz + Bz.sup.2 Dz + %5 Bz} Ax + (Dx Bz - Bx Dz)Ay.sup.2 +(Dx Bz - Bx Dz) Az.sup.2 + {-Dx Bx.sup.2 + %6 Bx - Dx By.sup.2 - DxBz.sup.2 } Az} c +{%l Ax.sup.2 + {(Cz - Dz) Bx.sup.2 + (Cz - Dz) By.sup.2 + (Cz - Dz)Bz.sup.2 + %8 Bz + Dz CX.sup.2 +Dz Cy.sup.2 + Dz Cz.sup.2 + %5 Cz} Ax + %1, Ay.sup.2 + %l Az.sup.2 + {(Dx- Cx) Bx.sup.2 + %7 Bx +(Dx - Cx) By.sup.2 + (Dx - Cx) Bz.sup.2 - Cx.sup.2 Dx + %6 Cx - DxCy.sup.2 - Dx Cz.sup.2 } Az +(-Dx Cz + Cx Dz) Bx.sup.2 + {-Dz Cx.sup.2 - Dz Cy.sup.2 - Dz Cz.sup.2 +%6 Cz} Bx +(-Dx Cz + Cx Dz) By.sup.2 + (-Dx Cz + Cx Dz) Bz.sup.2 +{CX.sup.2 Dx + %5 Cx + Dx Cz.sup.2 + Dx Cy.sup.2 }Bz} e}/((Dz - Cz) By + (Cy - Dy) Bz + Dy Cz - Cy Dz) Ex + %4 Ey +((Dy - Cy) Bx + (- Dx + Cx) By + Dx Cy - Cx Dy) Ez + (-Dy Cz + Cy Dz) Bx+(-Cx Dz + Dx Cz) By + (Cx Dy - Dx Cy) Bz) a +(((Cz - Dz) Ay + (Dy - Cy) Az - Dy Cz + Cy Dz) Ex + %3 Ey +((Cy - Dy) Ax + (Dx - Cx) Ay + Cx Dy - Dx Cy) Ez + (Dy Cz - Cy Dz) Ax +(- Dx Cz + Cx Dz) Ay + (Dx Cy - Cx Dy) Az) b +((((Dz - Bz) Ay + (By - Dy) Az + Dy Bz - By Dz) Ex + %2 Ey +((Dy - By) Ax + (Bx - Dx) Ay + Dx By - Bx Dy) Ez + (By Dz - Dy Bz) Ax +(Dx Bz - Bx Dz) Ay + (-Dx By + Bx Dy) Az) c +(((Cz - Dz) By + (Dy - Cy) Bz - Dy Cz + Cy Dz) Ax + %1 Ay +((Cy - Dy) Bx + (Dx - Cx) By + Cx Dy - Dx Cy) Az + (Dy Cz - Cy Dz) Bx +(- Dx Cz + Cx Dz) By + (Dx Cy - Cx Dy) Bz) e)%l: = (Dz - Cz) Rx + (-Dx + Cx) Bz - Cx Dz + Dx Cz%2: = (-Dz + Bz) Ax + (Dx - Bx) Az - Dx Bz + Bx Dz%3: = (Dz - Cz) Ax + (-Dx + Cx) Az - Cx Dz + Dx Cz%4: = (Cz - Dz) Bx + (Dx - Cx) Bz - Dx Cz + Cx Dz%5: = -Dx.sup. 2 - Dz.sup.2 - Dy.sup.2%6: = Dy.sup.2 + Dz.sub.2 + Dx.sup.2%7: = -Dx.sup.2 - Dz.sup.2 + Cy.sup.2 - Dy.sup.2 + Cz.sup.2 + Cx.sup.2%8: = Dy.sup.2 + Dx.sup.2 - Cz.sup.2 - Cx.sup.2 + Dz.sup.2 - Cy.sup.2Sz= 1/2{(((-Dy + Cy) Ex + (Dx - Cx) Ey - Dx Cy + Cx Dy) b +((Dy - By) Ex + (Bx - Dx) Ey + Dx By - Bx Dy) c + %1 e) a.sup.2 +(((-Dy + Ay) Ex + (Dx - Ax) Ey + Ax Dy - Dx Ay) c + %8 e) b.sup.2 +{((Dy - Cy) Ex + (-Dx + Cx) Ey - Cx Dy +  Dx Cy) b.sup.2 +((By - Dy) Ex + (Dx - Bx) Ey - Dx By + Bx Dy) c.sup.2 + %2 e.sup.2 + %1Ex.sup.2 +{(-Dy + Cy) Bx.sup.2 + (-Dy + Cy) By.sup.2 + %5 By + (-Dy + Cy) Bz.sup.2+ Dy Cx.sup.2 + Dy Cy.sup.2 +%4 Cy + Dy Cz.sup.2 } Ex + %1 Ey.sup.2 + {(Dx - Cx) Bx.sup.2 + %6 Bx +(Dx - Cx) By.sup.2 +(Dx - Cx) Bz.sup.2 - Dx Cx.sup.2 + %3 Cx - Dx Cz.sup.2 - Dx Cy.sup.2 } Ey+ %1 Ez.sup.2 +(-Dx Cy + Cx Dy) Bx.sup.2 + {-Dy Cx.sup.2 - Dy Cy.sup.2 + %3 Cy - DyCz.sup.2 } Bx +(-Dx Cy + Cx Dy) Bz.sup.2 } a +{((Dy - Ay) Ex + (-Dx + Ax) Ey - Ax Dy + Dx Ay) c.sup.2 +((Dy - Cy) Ax + (-Dx + Cx) Ay - Cx Dy + Dx Cy) e.sup.2 + %8 Ex.sup.2 +{(Dy - Cy) Ax.sup.2 + (Dy - Cy) Ay.sup.2 + %6 Ay + (Dy - Cy) Az.sup.2 -Dy Cx.sup.2 - Dy Cy.sup.2 +%3 Cy - Dy Cz.sup.2 } Ex + %8 Ey.sup.2 + {(-Dx + Cx) Ax.sup.2 + %5 Ax +(-Dx + Cx) Ay.sup.2 +(-Dx + Cx) Az.sup.2 + Dx Cx.sup.2 + %4 Cx + Dx Cz.sup.2 + Dx Cy.sup.2 }Ey + %8 Ez.sup.2 +(-Cx Dy + Dx Cy) Ax.sup.2 + {Dy Cx.sup.2 + Dy Cy.sup.2 + %4 Cy + DyCz.sup.2 } Ax +(-Cx Dy + Dx Cy) Ay.sup.2 + {-Dx Cx.sup.2 + %3 Cx - Dx Cz.sup.2 - DxCy.sup.2 } Ay +(-Cx Dy + Dx Cy) Az.sup.2 } b + %7 e c.sup.2 +{((By - Dy) Ax + (Dx - Bx) Ay + Bx Dy - Dx By) e.sup.2 + %7 Ex.sup.2 +{(By - Dy) Ax.sup.2 +(By - Dy) Ay.sup.2 + {-Bx.sup.2 + Dy.sup.2 + Dx.sup.2 - By.sup.2 +Dz.sup.2 - Bz.sup.2 } Ay + (By - Dy) Az.sup.2 +Dy Bx.sup.2 + Dy By.sup.2 + %4 By + Dy Bz.sup.2 } Ex + %7 Ey.sup.2 + {(Dx- Bx) Ax.sup.2 +{By.sup.2 30  Bz.sup.2 - Dy.sup.2 + Bx.sup.2 - Dz.sup.2 - Dx.sup.2 } Ax +(Dx - Bx) Ay.sup.2 + (Dx - Bx) Az.sup.2 -Dx Bx.sup.2 + %3 Bx - Dx By.sup.2 - Dx Bz.sup.2 } Ey + %7 Ez.sup.2 + (BxDy - Dx By) Ax.sup.2 +{-Dy Bx.sup.2 - Dy By.sup.2 + %3 By - Dy Bz.sup.2 } Ax + (Bx Dy - Dx By)Ay.sup.2 +{Dx Bx.sup.2 %4 Bx + Dx Bz.sup.2 + Dx By.sup.2 } Ay + (Bx Dy - Dx By)Az.sup.2 }  c + {%2 Ax.sup.2 +{(Dy - Cy) Bx.sup.2 + (Dy - Cy) By.sup.2 + %6 By + (Dy - Cy) Bz.sup.2 -Dy Cx.sup.2 - Dy Cy.sup.2 +%3 Cy - Dy Cz.sup.2 } Ax + %2 Ay.sup.2 + {(-Dx + Cx) Bx.sup.2 + %5 Bx +(-Dx + Cx) By.sup.2 +(-Dx + Cx) Bz.sup.2 + Dx Cx.sup.2 + %4 Cx + Dx Cz.sup.2 + Dx Cy.sup.2 }Ay + %2 Az.sup.2 +(-Cx Dy + Dx Cy) Bx.sup.2 + {Dy Cx.sup.2 + Dy Cy.sup.2 + %4 Cy + DyCz.sup.2 } Bx +(-Cx Dy + Dx Cy) By.sup.2 + {-Dx Cx.sup.2 + %3 Cx - Dx Cz.sup.2 - DxCy.sup.2 } By +(-Cx Dy + Dx Cy) Bz.sup.2 } e}/((((-Dz + Cz) By + (Dy - Cy) Bz + Cy Dz - Dy Cz) Ex +((Dz - Cz) Bx + (-Dx + Cx) Bz - Cx Dz + Dx Cz) Ey + %2 Ez +(-Cy Dz + Dy Cz) Bx + (Cx Dz - Dx Cz) By + (-Cx Dy + Dx Cy) Bz) a + . ..(((Dz - Cz) Ay + (-Dy + Cy) Az - Cy Dz + Dy Cz) Ex +((-Dz + Cz) Ax + (Dx - Cx) Az + Cx Dz - Dx Cz) Ey +((Dy - Cy) Ax + (-Dx + Cx) Ay - Cx Dy + Dx Cy) Ez + (Cy Dz - Dy Cz) Ax +(-Cx Dz + Dx Cz) Ay + (-Dx Cy + Cx Dy) Az) b +(((Bz - Dz) Ay + (Dy - By) Az + By Dz - Dy Bz) Ex +((Dz - Bz) Ax + (Bx - Dx) Az - Bx Dz + Dx Bz) Ey +((By - Dy) Ax + (Dx - Bx) Ay + Bx Dy - Dx By) Ez + (-By Dz + Dy Bz) Ax +(-Dx Bz + Bx Dz) Ay + (Dx By - Bx Dy) Az) c +(((Dz - Cz) By + (-Dy + Cy) Bz - Cy Dz + Dy Cz) Ax +((-Dz + Cz) Bx + (Dx - Cx) Bz + Cx Dz - Dx Cz) Ay + %1 Az +(Cy Dz - Dy Cz) Bx + (-Cx Dz + Dx Cz) By + (-Dx Cy + Cx Dy) Bz) e) +%1: = (Dy -  Cy) Bx + (-Dx + Cx) By - Cx Dy + Dx Cy%2: = (-Dy + Cy) Bx + (Dx - Cx) By - Dx Cy + Cx Dy%3: = Dx.sup.2 + Dz.sup.2 + Dy.sup.2%4: = -Dy.sup.2 - Dz.sup.2 - Dx.sup.2%5: = Dx.sup.2 - Cx.sup.2 + Dz.sup.2 - Cy.sup.2 + Dy.sup.2 - Cz.sup.2%6: = -Dy.sup.2 - Dx.sup.2 + Cx.sup.2 + Cz.sup.2 + Cy.sup.2 - Dz.sup.2%7: = (Dy - By) Ax + (Bx - Dx) Ay - Bx Dy + Dx By%8: = (-Dy + Cy) Ax + (Dx - Cx) Ay - Dx Cy + Cx Dy__________________________________________________________________________