Small-angle rotation method for star tracker orientation

Abstract The Wahba problem is the task of constrained optimization seeking the matrix from SO(3), which maximally converges (based on the least squares criterion) two sequences of unit vectors. The solution of this task is vital for satellite attitude determination using star trackers. An iterative method for solving the Wahba problem is proposed. Each iteration of the proposed method is reduced to sequential rotation of the vectors and solving the system of linear algebraic equations. Usage of the method implies that the corresponding vectors of both sequences are located sufficiently close to each other. Two variants of the method are proposed, having linear and quadratic convergence. The Wahba problem solution is interpreted in terms of finding the angular velocity of a system of material points, which have certain angular momentum. Taking into consideration the characteristics of state-of-the-art star trackers, one to two iterations are sufficient for finding the optimal solution using the small-angle rotation method. The primary advantage of the proposed method as compared with classical methods based on calculation of eigenvectors and singular decomposition is the simplicity of its implementation.


Introduction
Let us assume that matrix R belongs to a special rotation group: R ∈ SOð3Þ ⇔ R ∈ R 3×3 : R T R ¼ I; detðRÞ ¼ 1; (1) where I is a unit matrix with dimension 3. Let us assume that there are two sequences of unit vectors v and w in three-dimensional (3-D) space v i ∈ R 3 w i ∈ R 3 : jv i j ¼ 1; jw i j ¼ 1; (2) where i ¼ 1; : : : ; n, and n is the number of vectors in the sequence. We shall denote the function (3) as loss function where k i are weight coefficients, such that P n i¼1 k i ¼ 1. The norm of vector x hereinafter is understood to be the Euclidean norm It is necessary to find the matrix R opt from SO(3) which minimizes the loss function The Wahba problem has important applied significance because its solution is necessary in stellar orientation systems for satellite attitude determination. Stellar orientation systems are intended for determination of satellite attitude in the inertial (geocentric) reference system, and presently they are the most precise among existing 2 satellite orientation systems.
General functional diagram of a star tracker is shown in Fig. 1. Light from stars is passing through the optical path of the device and is projected on a photodetector. Using special algorithms, 3,4 stars are segregated on the background of photodetector noise and optical distortions, and their directional cosines are determined in the star tracker reference system.
If vectors w are interpreted as the directional cosines of stars in the inertial geocentric stellar coordinate system and vectors v are interpreted as the directional cosines of the same stars but in the star tracker reference system, then matrix R opt obtained as the result of solving task (5) will determine the attitude of the satellite with reference to the Earth (if the universal time is known). Therefore, hereinafter we shall name R opt as optimal orientation matrix. Coefficients k i are chosen depending on the accuracy of determination of star coordinates on the photosensitive matrix, which is usually proportional to the star energy.
The values of vectors v are known not accurately because of star orientation instrument errors in the course of determining star positions. Such errors are caused by optical perversions (lens distortion and aberration), noise, and discrete structure of photosensitive matrix. The values of vectors w are also determined with errors due to inaccuracy of existing star catalogues.
Let v true and w true be the true (without introduced error) directional cosines of stars in the inertial reference system and the star tracker reference system, and let R true be the matrix determining the rotation, which aligns vectors v true and w true It is possible to consider that vectors v and w are where ε v i ε w i ∼ Nð0; σÞ. If jε v j and jε w j are small, then it holds that R opt ≈ R true . In practice, the error of star determination by instrument is substantially larger than the star catalogue error, jε v j ≫ jε w j.
There are two approaches to solve Wahba problem (5): direct minimization of loss function and exact analytical solution. Since the first formulation of the Wahba problem, numerous methods for its analytical solution have been proposed. 5 Many of these methods seem quite simple, but actual computations come up against complications, e.g., finding the square root from an ill-conditioned matrix.
On the other hand, not much attention in the literature has been devoted to the possibility of direct minimization of the function LðRÞ. The reason is that the application of minimization methods requires the existence of sufficiently close initial approximation and calculation of function LðRÞ values at every step. The laboriousness of calculating the values of function LðRÞ grows with increase of the number of vectors n.
When solving the Wahba problem in connection with the task of a satellite attitude control, it is necessary to take into consideration the special features of the application domain and the Fig. 1 Functional diagram of a star tracker. 1-stellar sky, 2-optical system, 3-electronic module with photodetector array, 4-computing module, and 5-star tracker. resultant constraints on the character of input data. Usage of these constraints allows the creation of more efficient methods for solving the Wahba problem as compared with the generalized problem formulation. Two such constraints on the values of input data should be noted.
First, the number n of stars for determining orientation usually is not over 50 because of optical characteristics of star trackers and distribution of stars over the celestial sphere. Actually, not all the stars are used for recognition, but only 10 to 15 of the brightest stars, because their coordinates have the highest accuracy.
Second, when considering the solution of the Wahba problem in connection with the task of satellite attitude control, one should take into account the fact that the errors of vectors v and w are small as compared with the minimal distance between two arbitrary vectors of the sequence.
Inequalities in Eq. (8) are fulfilled since star catalogues are compiled providing a condition that a distance between the nearest stars by a few times exceeds permissible error ε v i . If the star position error significantly exceeds the possible one (due to proton influence, pixel defects, etc.), then such stars will not be detected by identification algorithm, since they are rejected according to the criterion of mutual angular distances.
Provided that Eq. (8) is fulfilled, it is possible to calculate easily and with sufficient accuracy the initial approximation for the matrix R 0 ≈ R opt using the triangle attitude determination (TRIAD) method. 6,7 Using the initial approximation R 0 , it is possible to rotate vectors v so that they are sufficiently close to vectors w: The result from this is that Hence, when solving the Wahba problem for the task of satellite orientation, it may be assumed that Eq. (10) applies.

Review of Methods for Solution of the Wahba Problem
A short but recent review of methods for solving the Wahba problem is given in Ref. 8. An earlier but detailed analysis of methods for solving the Wahba problem was made by Markley and Mortari in Ref. 5. The basic results related to this problem were obtained in 1980s and 1990s; 9-13 however, due to rapid progress in space instruments development, in particular, in star sensors, this topic is still actively discussed during the last decade. [5][6][7][14][15][16][17][18][19][20][21] The majority of classic methods for solution of the Wahba problem are based on the usage of matrix B (sometimes called the attitude profile matrix 6 ), which is formed on the basis of vectors v and w: One of the first analytical solutions of the Wahba problem has been proposed by Farrell and Stuelpnagel, 22 who have shown that R opt coincides with orthogonal matrix W of polar decomposition B ¼ WS, where S is a symmetrical matrix.
According to Farrell, matrix W may be expressed using matrix B as where kAk F ¼ trðA T AÞ is the Frobenius norm.
According to Ref. 22, the W matrix may be expressed through the B matrix as Though Eq. (13) appears quite simple, its practical application is connected with considerable difficulties. During solution of the task of satellite orientation, the vectors of sequence v corresponding to the directional vectors of the stars are located sufficiently close to each other (usually within the cone of 10 angular degrees or smaller), due to the size of the star tracker's field of view (FOV). For this reason, matrix B T B is often an ill-conditioned matrix, and the process of computing its square root, for instance, with the Denman-Beaver iteration 23 is either diverging or has large error.
Matrix W of polar decomposition may be expressed in terms of singular value decomposition (SVD) of matrix as Equation (14) is the basis of the method proposed by Markley. 11 The fast optimal attitude matrix (FOAM) and slower optimal matrix algorithm (SOMA) methods by Markley 10 are also based on the idea of SVD of matrix B, which in the mentioned methods is reduced to solving the system of nonlinear algebraic equations. According to the computational experiments performed by Markley, 10 determination of the optimal orientation matrix using the SVD method has higher accuracy than the FOAM and the SOMA methods; however, the latter two methods are faster.
The most popular group of methods for solving Wahba problem are the methods (usually called q-method or Davenport method) calculating the quaternion of the respective orientation matrix R opt , as the eigenvector corresponding to maximal eigenvalue λ max of matrix K: where s ¼ TrðBÞ; Vector b may also be expressed using components of matrix B Matrix K is a symmetrical matrix; therefore, efficient algorithms exist for calculation of eigenvectors of the matrix, e.g., Jacobi method or the methods based on QR decomposition with preliminary reduction to tridiagonal form. 24 However, implementation and debugging of such algorithms may be quite labor consuming, considering that coding often has to be performed using low-level programming languages or erasable programmable logic device (EPLD). It is necessary to keep in mind that the specified task has to be solved in real time. Besides, the capacity of spaceborne processors is significantly constrained. Therefore, despite small dimension of matrix K, calculation of its eigenvectors may be laborious.
To avoid straightforward calculation of eigenvectors of matrix K, Shuster and Oh 9 have proposed the quaternion estimator (QUEST) algorithm, which gives an explicit formula for finding the eigenvector corresponding to the maximal eigenvalue λ max , provided that this value is known. The maximal eigenvalue is searched for as the root of a nonlinear algebraic equation, which the authors propose to solve using Newton-Raphson method. However, such method of finding eigenvalues is the least stable method, therefore, direct search of eigenvectors of matrix K gives more accurate results.
Other important results obtained by Schuster are as follows: he was the first to consider the problem from the statistical point of view and has shown that the loss function as a random value is distributed by chi-square distribution law with 2n − 3 degrees of freedom: Schuster also has shown that the covariance matrix P of rotational angle error equals to: Attitude angle errors are the values ϕ ¼ ðϕ 1 ; ϕ 2 ; ϕ 3 Þ, such that where ½ϕ× is a skew-symmetric matrix determined by elements of vector ϕ. Markley 11 has shown that in case of a large number of observations, the following approximate estimate of the covariance matrix will be valid The results obtained by Schuster and Markley show that the presented estimates of the Wahba problem are consistent. However, the issue of unbiasedness of the proposed estimates remains open. In order to determine what is the unbiased estimate for SO(3), at first it is necessary to clarify the concepts of mathematical expectation and the average value for the elements of this manifold. Two approaches for determination of the concept of average value for the rotation group-Euclidean and Riemannian-are given, for example, in Refs. 25 and 26. The approaches to solve the Wahba problem that are closest to the approach proposed in this work were presented by Park and Kim, 27 Mortari et al., 28 and Mortari. 13 Authors of Ref. 27 were solving the task of constrained minimization similar to Wahba problem for position fixing using global positioning system. The methodology proposed by them is using the exponential relation between Lie group SO(3) and Lie algebra soð3Þ for determining the gradient of the loss function and for sequential search of its minimum using Newton's method and steepest descent method.
The "energy approach algorithm" 13 proposed by Mortari has a common feature with the proposed, in this article, small-angle rotation (SAR) method: based on the closeness of vectors v and w, it substitutes the loss function by an approximate function-an energy function. Despite the similarity of the basic approach, the energy approach algorithm is principally different from the SAR because it is based on a computation of eigenvalues or SVD, while the SAR algorithm is based on solving the system of algebraic equations.

Basic Definitions
Let soð3Þ be the set of skew-symmetric matrices The exponent of a square matrix A is the square matrix equaling the sum of infinite series The properties of the matrix exponent are verified easily expðAÞ expð−AÞ ¼ I; det½expðAÞ > 0: (29) In that case, if matrix M ∈ soð3Þ, then according to Eqs. (27) and (28) exp ðMÞ T expðMÞ ¼ expð−MÞ expðMÞ ¼ I: Taking into consideration Eq. (29), this means that expðMÞ ∈ SOð3Þ. That is, the exponent of a skew-symmetric matrix is a rotation matrix ∀ M ∈ soð3Þ ⇒ expðMÞ ∈ SOð3Þ: (31) This remarkable fact has deep connection with the theory of Lie groups and algebra. 29 For vector x ¼ ðx 1 ; x 2 ; x 3 Þ, we shall denote, hereinafter, as ½x× the skew-symmetric matrix: Let ½a; b be the operation of cross product, then the following relations hold true

SAR Algorithm of the First Order (Angular Momentum Method)
Let us assume that a system of vectors v and w exists, satisfying condition (10) specified in Sec. 1. As already mentioned in Sec. 1, after determining the approximate attitude using TRIAD method, 6,7 using transformations (9), it is possible to achieve the fulfillment of condition (10).
Let ω ¼ ðω 1 ; ω 2 ; ω 3 Þ be the rotation vector of the coordinate system. Then by reason of expression (31), it is possible to replace LðRÞ ¼ L½expðωÞ ¼ LðωÞ, and in this case, the function L may be represented in the form As mentioned in Sec. 1, the rotation angle ω opt is small. Therefore, it is possible to limit the expansion into series (26) by the first several terms.
in this case Taking into account Eq. (33) where const means a constant independent from ω. Thus, the task of constrained minimization of the loss function LðωÞ may be reduced to the task of unconstrained minimization of the quadratic form f 1 ðωÞ.
The necessary condition for achieving the extreme point by a function is equality of its gradient to zero. Due to symmetric property of matrix ½v× 2 , gradient f 1 ðωÞ equals to The search for ω satisfying the equality ∇f 1 ðωÞ ¼ 0 is reduced to solve the system of algebraic equations and b is determined in Eq. (19).
It is easy to note that matrix A coincides with the inertia tensor of material points determined by vectors v i and having masses k i . Hence, the solution of Eq. (40) has a vivid physical interpretation-it is necessary to find such angular velocity ω for which the inertia momentum of a rigid system of material points defined by vectors v and having masses k i will equal to the inertia momentum of the same system, if unit velocity w i is applied to each material point. At that, f 1 ðωÞ function may be regarded as the value of kinetic energy of the system of material points.
Considering the mentioned mechanical interpretation, the SAR method of the first order may be also called the "angular momentum method." It is worth mentioning that Mortari in Ref. 13 also approximated the loss function by a function that he regarded as the kinetic energy function-energy of n compressed springs with coefficients of rigidity k i . In all the rest, Mortari's approach has nothing common with the solution proposed here.

SAR Algorithm of the Second Order
It is logical to expect from a solution of the Wahba problem that in case of permutation of the sequences of vectors v and w, the sought-for rotation matrix A would have to be inversed. For the SAR, it means that the SAR vector ω changes its sign. However, for the SAR of the first order, in case of permutation of v and w, the SAR vector, in general case, changes not only its sign but also the absolute value as well. The reason is that matrix A (41) depends only upon vectors v but does not depend upon vectors w.
A modification of the SAR method invariant to permutation of vectors and giving the more accurate estimate of ω than Eq. (40) is presented in this section.
Let us transform the loss function L into the gain function, as was done in Ref. 6: Thus, the task of minimizing the loss function was transformed into the task of maximizing the function gðωÞ. Using uncomplicated transformations, it is possible to demonstrate that Replacing the exponent by the first three terms of relevant series, it is possible to obtain the approximation of the function gðωÞ ≈ f 2 ðωÞ In that case Finding ω satisfying the equality ∇f 2 ðωÞ ¼ 0 is reduced to solving the system of algebraic equations then matrix A looks as follows: Matrix A is connected with attitude profile matrix B by the following equation: where scalar s and matrix S are expressed through attitude profile matrix B, as shown in Eqs. (17) and (18).

Matrix Corresponds to the Small-Rotation Vector
When the SAR vector ω is found, the sought-for approximation of matrix R opt is obtained as There are many methods for calculating matrix exponent, 30 but for the case of a skewsymmetric matrix, the simplest way is to use the Rodriguez formula: 29 where Θ ¼ jωj, μ ¼ ω∕Θ.
Instead of using the matrix form for setting the rotation in R 3 , it is preferable to use the quaternion presentation, which is more convenient as it requires less memory space and allows reduction of the number of operations necessary for vector rotation. Calculation of the orientation quaternion corresponding to rotation ω is easier than calculation of the corresponding matrix: Though usage of quaternions is more convenient from the computational point of view than usage of matrices, we shall continue using the matrix notations in this work for the sake of convenience of theoretical calculations. If required, the matrix notations may be easily substituted by quaternion analogs.

Evaluation of Solution Accuracy
It is necessary to determine the criteria for evaluation of the accuracy of the obtained solution. Different approaches to this issue are possible.
In technical applications, the error of star trackers is measured most often by three Euler angles by which matrix R 1 must be additionally rotated in order to be aligned with matrix R 2 . Euler angles set the sequential rotation around axes OX, OY, OZ and are commonly used in mechanics and computer graphics and are called roll, pitch, and yaw.
For star sensors, axis OZ traditionally coincides with the optical axis of the device; therefore, the largest orientation errors are achieved for the third angle yaw (rotation around line of sight), whereas the roll and pitch errors are usually equaling to each other.
Euler angles are convenient for practical usage. Usually, a satellite is always facing the Earth with one of its sides (where scanning equipment or radio transmitter is located). The location of the star tracker within the satellite is known; therefore, the orientation error measured in Euler angles determines the error of directing the satellite on the Earth.
However, Euler angles are extremely inconvenient for theoretical evaluations, as they are not invariant to the change of coordinate system. Besides, it is desirable to express the combined error by one number rather than by three numbers.
For evaluation of a difference of rotation matrices, usually two metrics are used-Riemannian and Euclidean. In the first case, the set of rotation matrices SO(3) is regarded as Riemannian metric manifold with metric 25 The Euclidean metric is based on the Frobenius norm The value d R varies between 0 and π, whereas d E varies between 0 and 2. The set of all rotation matrices SOð3Þ is homeomorphic to a 3-D sphere in four-dimensional space, and each matrix may be interpreted as a point on this sphere. In such an interpretation, d R may be regarded as angle, whereas d E may be regarded as chord between respective points on the sphere.
The Riemannian distance dðR 1 ; R 2 Þ defines the minimal absolute value of the angle by which the coordinate system R 1 must be rotated around an arbitrary axis in order to align it with the coordinate system R 2 .
If the rotation is performed around one and the same axis for the angles φ 1 and φ 2 , then for the matrices R 1 ðφ 1 Þ and R 2 ðφ 2 Þ, the value of Riemannian metric d R ðR 1 ; R 2 Þ equals to jφ 1 − φ 2 j, naturally complying with our perception of orientation error.
Considering that accuracy in optics and in astronomy is usually measured by angular minutes (or angular seconds), the value d determining accuracy will be expressed hereinafter in angular minutes.
It is worth mentioning that in case of close matrices, Euclidean metric and Riemannian metric are close to each other Based on geometrical interpretation of metrics, relation (56) means equality of arc length and chord length for small angles. As star trackers are high-precision instruments with accuracy of the order of angular seconds, it may be assumed that relation (56) holds for them. In order to highlight this fact, hereinafter, we shall denote the metric d R simply as d. 5 Step-by-Step Form of the Algorithm Input data for the algorithm are: R 0 -initial approximation of orientation matrix; w, v-sequences of unit 3-D vectors of dimensionality n; ε-required accuracy of solution (in radians).

Output data:
R-optimal orientation matrix.
All matrices and vectors have the dimension equal to 3. The general scheme of the SAR is as follows. Gain function gðωÞ is approximated by quadratic form f 2 ðωÞ. By solving the system of linear Eq. (42), the minimum of function f 2 ðωÞ is found. When the respective rotation angles ω are found which minimize the function f 2 ðωÞ, the rotation matrix R is calculated on its basis. Vectors v are rotated by matrix R and approach to w.
These operations are repeated until rotation angle ω becomes negligible in the framework of the task being solved.
The steps of iterative algorithm for solving Wahba problem using the SAR method: 1. Determining the initial approximation R 0 using TRIAD method.
: : : ; n. 4. Determining, on the basis of formula (42), the parameters of the system of linear algebraic equations A and b using data w and v. 5. Solving the system of linear algebraic equations Aω ¼ b. 6. According to formula (52), obtaining rotation matrix R ω from vector ω. 7. R ¼ R ω R. 8. If the condition kR ω k 2 F < ε 2 is satisfied, then go to step 9, else to step 3. 9. End of algorithm.
If the rotation is set by quaternion, the algorithm will be similar. All the steps of the SAR algorithm are sufficiently simple for coding including low-level programming languages. The most source-consuming among these steps are extraction of square root and solving system of linear equations. Binary arithmetic has fast and effective methods for calculating square root. Implementation of Gauss method is also rather simple even using lowlevel programming languages. Moreover, solution of 3-D system may be easily implemented directly using Cramer formula (for such dimensions, it has a negligible influence on execution time and accuracy). At the same time, finding eigen numbers, even by means of the simplest Jacobi method, requires nontrivial logic which is not simple for implementation using assembler or EPLD without saying about more complicated methods. 24 Methods of singular matrix decomposition are also complicated for implementation.
The first-order SAR algorithm was implemented in software of 329K star tracker, 31 which successfully operates on board the data relay satellite Luch-5A since December 2011 and on board the satellite Luch-5B since November 2012. Implementation of the algorithm using digital signal processor (DSP) has shown that it is simple for coding and program debugging. 6 Simulation Results

Convergence of the Method
The SAR of the first order and the second order is obtained using approximation of matrix exponent by its decomposition into exponential series up to the first and the second power, respectively [Eqs. (36) and (44)]. Thus, approximation by means of loss functions f 1 ðωÞ and f 2 ðωÞ has, respectively, an error of orders ω and ω 2 . Since minimum of functions f 1 ðωÞ and f 2 ðωÞ is seeking exactly at each step according to formulae (40) and (47), so the final error of function LðωÞ minimum will also have order of ω and ω 2 .
In order to test the accuracy and convergence of the proposed method, its modeling in MATLAB system has been performed. Data of Tables 1-4 are obtained by means of averaging not less than 100,000 experiments.
Simulation was performed as follows. Rotation matrix R true was set in random manner. n unit vectors w i were uniformly placed in a random manner in the FOV (spherical square 2W) of the star sensor. Some errors were added in coordinates of each vector to x and y components of each vector.
where ε ∼ Nð0; σÞ-normally distributed random values. In this case, it was not allowed that star angular distances to be less than 10 SD σ of the added error. Such restriction is quite reasonable since the adjacent stars are not included usually at the guide star catalogues.
Using the SAR and singular decomposition method, estimates of the matrix R opt were obtained for the sequence w and v: R opt1 R opt2 . According to formula (49), the errors Δ 1 ¼ dðR opt1 ; R true Þ and Δ 2 ¼ dðR opt2 ; R true Þ were calculated. Table 1 Average error of estimate R opt using the singular decomposition method and the smallangle rotation (SAR) method of the first order in angular minutes.

Method
Step   The mentioned sequence of operations was performed for two SAR methods: the first and second order.
The values Δ 1 and Δ 2 and their difference Δ ¼ Δ 2 − Δ 1 for the SAR methods of the first and second order are given in Tables 1 and 3. The values of square root from the loss function L are given in Tables 2 and 4.
The following parameters were used during modeling: the side of the square FOV 2W ¼ 20 angular degrees, number of stars n ¼ 15, SD σ ¼ 10.0 arc min. It should be noted that the SD value used during the simulation is much larger than actual values. This was done in order to increase the number of steps executed by the algorithm before reaching the optimal point. Solution by TRIAD method was used as the initial approximation.
As may be seen from Tables 1-4, the singular decomposition method provides more accurate estimate than the SAR though the difference becomes negligibly small after several iterations. Since the exact value R opt is not known, the estimate R opt may be approximated (with the accuracy of singular decomposition by method 32 ) by the estimate R opt of the SVD method. In this case, the third row of the table may be interpreted as the error of the SAR method.
In order to better understand the convergence behavior of the SAR, Figs. 2 and 3 show the dependence of the difference between accuracies of the methods from the number of iterations (in logarithmic scale). Every line in Figs. 2 and 3 corresponds to a single experiment (differing from Tables 1-4 whose data are obtained by means of averaging).
The simulation results confirm the conclusion about the method's convergence rate. As may be seen from Fig. 2 and Table 1, SAR of the first order has linear convergence. Based on Fig. 3 and Table 3, the conclusion may be made that SAR of the second order has super-linear convergence. Graphs in Fig. 3 have form which differs from parabola because the difference between the SAR and the SVD method is stabilized starting from the third step. This is caused by the error of number representation with the type "double" and the error of the singular decomposition algorithm implemented by LAPACK software package. 32 Table 4 Value of loss function LðRÞ for the estimate using the singular decomposition method and the SAR method of the second order in angular minutes.

Method
Step  At any rate, experimental data are indicative of the fact that the second-order method is by far more effective than that of the first order; therefore, the former is more expedient to be used for calculations. Given this, its computational complexity is just very little higher than that of the first order method. A matrix calculation for the second-order method will require more multiplication operations than it will for the first-order method. Otherwise, computational complexity of the algorithms is identical.
Simulation shows that if the error of initial approximation of orientation is of the order of 1 angular minute, then after one iteration of the SAR, the achieved resulting accuracy is of the order of 10 −5 angular minute, abundantly meeting the requirements to the errors of state-of-theart star trackers. In case of necessity, two iterations of the method may be performed. Performing more than three iterations of the SAR makes no sense because in this case, the accuracy of the method exceeds the accuracy of data presentation using the type double.

Method Stability against Single Errors
An important issue in star tracker attitude estimation is the presence of inaccurate data. Radiation, pixel failures, etc., can cause one or more of the star vectors to have noise values that are significantly larger than others. Methods such as the SVD method and Davenport's q-method are very robust to this, while faster methods such as QUEST have more problems dealing with this.
Simulation was performed to verify stability of the SAR against errors of such kind. The errors of 60 arc sec, i.e., significantly larger than usual errors for stars, were added to directing vectors of two stars from 15. Since SAR accuracy significantly depends on accuracy of initial approximation, the most interesting for simulation is the case when two vectors having maximum errors are chosen to estimate initial approximation. The simulation data are presented in Table 5. The data of Table 5 were obtained by means of averaging not less than 50,000 experiments. Fig. 3 Dependence of the difference between accuracies of the SAR method of the second order and the SVD method Δ ¼ jΔ 1 − Δ 2 j. OX axis, number of iterations of SAR method. OY axis, difference between accuracies of the two methods in angular minutes (logarithmic scale). Taking into account the characteristics of modern star trackers (with the FOV of angular radius W ≈ 10 deg ), SD for star position error is not able to exceed 1 arc min since else it will not be identified by the star identification algorithm according to criterion of mutual angular distances. Proceeding from this estimation, errors of 60 arc sec were added to the first two stars' positions at the simulation. The SD σ of errors for other stars is presented in Table 5 along with the simulation results.
The given Table 5 shows that the SAR algorithm is stable against random single errors and possesses, concerning this criterion, properties close to SVD algorithm. The difference in accuracy for SVD and SAR methods is comparable with the data presentation error.

Time of Algorithm's Execution
The purpose of this section is to compare the computational complications of the SAR algorithm and the classical algorithms of solving the Wahba problem based either on SVD or on calculation of eigenvectors.
The components requiring the most computational efforts at each iteration of the SAR are rotation of n vectors with subsequent computation of matrix A and evaluation of rotation vector ω . From the computational point of view, these tasks are reduced to sequential multiplication of n vectors by square matrix or quaternions and subsequent solution of the system of linear algebraic equations.
Although finding the eigenvectors is similar to solving the system of linear equations and is formally characterized by asymptotic 33 complexity Oðm 3 Þ, the actual number of operations required for finding the eigenvectors is much higher. The reason is that the Jacobi method, similar to the methods of matrix reduction to standard forms and QR decomposition, 24  It is obvious that computational complexity of all the algorithms mentioned in the review is OðnÞ depending on the number of stars, since these algorithms are related to calculating attitude profile matrix B [Eq. (11)]. Moreover, besides the calculation of the matrix B rotation of n vectors is required for SAR algorithm. Due to additional computational consumption for the vectors' rotation, SAR algorithm has linear coefficient, which is larger than SVD and q-method have.
It is also obvious that time of the algorithm operation linearly depends on quantity of iterations.
To compare algorithmic complexity of the proposed algorithms with state-of-the-art algorithms, evaluation of their execution time in PC was performed. SVD, q-method, QUEST, and SAR of the second order for one and two iterations were chosen for the comparison. The algorithms were implemented using MATLAB 7/0 and a computer with the processor Intel Core Duo 2.00 GHz. Time of the algorithms' execution versus number of stars and the number of iterations are presented in Figs. 4 and 5.
Interpretive program MATLAB is known for its slow operation, so the absolute values of execution time should be regarded very accurately. Real time of algorithm's execution depends significantly on the method of its implementation using one or another processor. For comparison, time of SAR algorithm's execution at a single iteration for 10 to 15 stars using 40-MHz DSP processor is equal to 10 μs approximately. For modern space class DSP, this time will be tens of milliseconds.
It is seen from Fig. 4 that the methods SVD, q-method, and QUEST have the same linear coefficient of increase depending on the number of stars. The SAR method has larger linear coefficient comparing with SVD, the q-method, and QUEST and therefore its execution time will rise a greater extent at increasing number of stars. However, as was mentioned above, not all stars are used for recognition but usually only 10 to 15 of the brightest stars because their coordinates have the highest accuracy.
The fact, that SAR algorithms' execution time rises a greater extent at increasing number of stars comparing with SVD, q-method, and QUEST is a fundamental property of the algorithms. At the same time, the constant constituent of the algorithm's execution time (the graph value for two stars) may significantly vary depending on implementation of the algorithm using one or another processor.
The SAR algorithm's execution time versus quantity of the method iterations is presented in Fig. 5. The number of stars in FOV is 10. Execution time for the algorithms SVD, the q-method, and QUEST (shown by horizontal lines) is given for comparison. As should be expected, the  dependence has linear form. For quantity of iterations more than two, time of algorithms' execution for SAR is essentially larger than for SVD, q-method, and QUEST. However, as was mentioned in the preceding section, to meet current accuracy requirements, the method's one or two iterations are quite sufficient.

Conditionality of the Equation System Solution
Unlike the mentioned methods of finding the eigenvectors, solving a system of linear equations using the Gauss method does not require the calculation of square root, and that is the advantage of the SAR. However, employing other method of solution than the Gauss method may necessitate calculation of square root. For example, since matrix A is a symmetric matrix, the search for system may be solved using the Cholesky method, which, in general case, is more stable than the Gauss method.
To compare stability of the Gauss method and the Cholesky method, a computational experiment was performed. Both methods were programmed in C++ language using Visual Studio 2010. Numbers were represented with double precision. For device FOV angular radius ranging from 1 to 10 deg, matrices A were generated, and SAR vectors ω were calculated by the Gauss method and the Cholesky method. The computational experiment showed that the difference between solutions obtained by the two methods is within 10 −12 arc min, which is substantially smaller then the accuracy requirements raised to the present-day star trackers. Therefore, application of the Gauss method instead of the Cholesky method in the context of the task being solved is quite reasonable.
The largest difference in the accuracy of the Gauss method and the Cholesky method was observed for the smallest FOV. The reason is that the closer the stars (and the corresponding vectors) are located to each other, the more degenerated the task of finding the minimum of LðωÞ function. The following fact holds that the smaller the FOVof the star tracker, the more elongated function L is.
The problem of degeneracy of the search task is illustrated in Fig. 6 showing the dependence of the loss function L from parameters ω 2 and ω 3 . In order to enable the presentation of the dependence in 3-D space, it is supposed that ω 1 value is known precisely. Figure 6 distinctly shows that diagram L is elongated in shape, and its level lines are similar to ellipses. In optimization theory, 34 such functions are named as "ravine functions." If L function is approximated by a quadratic form, as was done in Eqs. (36) and (44), then the ratio of maximum and minimum singular values of matrix A of the quadratic form may serve as the measure of elongation of L. If matrix A is positively defined, then the square root of the above mentioned ratio is the ratio of the lengths of maximum and minimum semiaxes of Fig. 6 Dependence of loss function L from parameters ω 2 and ω 3 . the ellipsoid determined by the corresponding quadratic form. The ratio of maximum and minimum singular values is the matrix condition number condðAÞ for Euclidean norm. Therefore, the size of the device FOV is connected with the elongation of L function, which, in its turn, is connected with the condition number of matrix A and, respectively, with accuracy of finding ω.
Based on the mathematic simulation, a following observation was made, in case of disposition of stars in a circular FOVof angular radius W, the condition number of matrix A is proportional to squared cotangent of the FOV angular radius W of the star tracker.
The author did not find the proof of this fact, but its truth is confirmed by simulation results for the several particular cases. During simulation, the stars were arranged in the FOV of the device as shown in Fig. 7. In total, 3 possible relative positions of stars were considered, for the number of stars n ¼ 2; 3; 4. Simulation results are shown in Fig. 8. The indicated dependence is valid also for other quantity and mutual positions of the stars.
It is worth mentioning that the problem of degradation of conditionality when the FOV is decreasing is not only the specific problem of the SAR but is also typical for the Wahba problem in general.

Conclusions
The SAR method proposed in this work is an iterative sequential method for solving the Wahba problem, having linear and quadratic convergence. For the data represented by double precision, the optimal solution is achieved after two iterations, with an error of initial approximation many  times higher than actually possible. For actual data, one or two iterations of the method are sufficient. Each iteration of the method is reduced to sequential multiplication of n vectors by a matrix or to multiplication of quaternions and solving of the system of linear algebraic equations with dimensionality 3. The method is stable against single errors of individual vectors.
The primary advantage of the proposed method as compared with classical methods based on calculation of eigenvectors and singular decomposition is the simplicity of implementation.
The SAR method of the first order was implemented in software of wide-angle star tracker 329K. The trackers of this type successfully operate already more than 1 year on board of the Russian data relay geostationary satellites Luch-5A and Luch-5B.

Future Works
It would be useful to determine the theoretical dependence of conditionality of the task being solved from size of the FOV, which, as experimentally shown in this work, is proportional to the squared cotangent of the FOV's angle of the star tracker.