## FtyAyH

an implicit method can be used directly, yielding k-1 k-I

Such methods for linear differential equations are known as correct-only methods, for reasons that will be apparent shortly. If we are integrating a system of equations, A is a matrix and Eq. (17-13) requires evaluation of a matrix inverse. (The 1 in the first bracket becomes the identity matrix.)

* An alternative procedure is to utilize the last k values of yK and /„ that have been evaluated, regardless of step size changes, rather than requiring the back values to be at evenly spaced points. This requires inversion of a kxk matrix at each step to find the coefficients, and the stability properties of these methods are less well understood than those of the conventional methods considered here [Yong, 1974]

Implicit methods can be made to have smaller truncation errors and better stability properties than explicit methods [Lambert, 1973], so it is desirable to have some method of finding fn+, to enable an implicit method to be used for nonlinear equations. The usual procedure is to use an explicit method, known as a predictor, to calculate y„+v Then fn+1 is evaluated and an implicit method, known as a corrector in this application, is used to obtain a refined value of/n+1, followed by a second evaluation of fn+i using the new yn+i. Methods in this general class are called predictor-corrector methods. The mode of application described above is the PECE (predict-evaluate-correct-evaluate) mode; it requires two function evaluations per step. It is also possible' to apply the corrector more than once after a single use of the predictor, but analysis indicates that the PECE mode is preferable, and that decreased truncation error is better achieved by decreasing the step size than by multiple applications of the corrector [Lambert, 1973].

It can be shown that no convergent fc-step method can have order greater than k +1 if k is odd or k + 2 if k is even [Henrici, 1962]. However, methods of order k + 2 have poor stability properties, so A: +1 is the optimal order for practical A:-step methods. The most commonly used A:-step methods are the Adams methods, which are defined by choosing <**_,= -! and aj=0 for j k—\. These have good stability properties and reduced computer storage requirements compared with other multistep methods. One explicit Adams method, the Adams-Bashforth, has the form jw.-A+A's'0/.+1+7-* 07-14)

and is of order k. An Adams-Moulton method is an implicit Adams method of order k +1, given by y.+i-y.+ir'k 0/-1+7-* <17-15)

A widely used predictor-corrector pair is a /»-step Adams-Bashforth predictor followed by a (p— l)-step Adams-Moulton corrector; both steps are of order p. One advantage of this pair is that the difference between the predicted and the corrected values of yn+l gives an estimate of truncation error and can be used for step size control. This is in contrast to Runge-Kutta methods, for which step size changes are relatively easy, but estimates of truncation error are difficult to obtain. The fourth-order Adams-Bashforth-Moulton pair is given by Predictor (explicit):

Corrector (implicit):

A-. =yn + ¿(9/„+, + 19/. - 5/„_, +/„_2) (17-16b)

This is only an example; higher order methods are widely used and, unlike higher order Runge-Kutta methods, cost only additional storage space and not additional function evaluations. (See, for example, Hull, et al., [1972] or Enright and Hull [1976].)

If values of y are needed at intermediate times and an Adams integrator is being used, it is convenient to employ an interpolation algorithm based on values of /„ rather than^„, because the former already have to be stored for the integration routine.

In either the Runge-Kutta or the predictor-corrector calculations, some of the function evaluations may be done approximately rather than exactly (by not recalculating torques, for example) to save computational effort. These are called pseudoevaluations, and are represented by E*, so one often reads of a PECE* mode of a predictor-corrector.

In choosing an integration method, the factors of programming complexity, computer storage requirements, execution time, and computational accuracy must all be considered. For a specific application where the characteristic frequencies of the system are known to be nearly constant, a fixed-step method is indicated. If the step size is limited by variations in the driving terms rather than by integration error (noisy input and/or low-accuracy requirements) or if function evaluations are relatively inexpensive, a Runge-Kutta method is preferred. If, on the other hand, the integration step is set by integration error (smooth input, high accuracy), or function evaluations are expensive, a predictor-corrector method is better. Adams methods are favored in this class because they combine good stability properties with relatively low computer storage requirements and programming complexity.

If the characteristic frequencies of the system are not constant, a variable-step method should be used. A complete integration package in this class must include an algorithm for automatic step size variation, based on an estimate of local truncation error. Because predictor-corrector methods provide an automatic estimate of local truncation error, they are the preferred variable-step methods. The best general-purpose integration methods currently available are packages with variable-step and variable-order Adams-Bashforth-Moulton integrators (see Hull, et al„ [1972] and Enright and Hull [1976], which also include comparative tests of integration packages using a wide variety of test cases).

Approximate Closed-Form Solution for the Kinematic Equation. As described in Section 17.1.1, when a set of gyros is part of the attitude determination hardware and the method of gyro modeling is used, only the kinematic equation need be integrated for attitude propagation. Wilcox [1967] and Iwens and Farrenkopf [1971] have presented a method of processing gyro data which yields an approximate closed-form solution to the kinematic equation: If we assume that the gyro data is telemetered or sampled at a fixed rate and that the angular velocity vector in body coordinates is constant over the sampling interval, then a closed-form solution to the kinematic equation (Eq. (17-12)) is

where T is the sampling interval (T= tn+, - /„); i2n is Eq. (17-3) evaluated at time tn; q(tn) is the attitude quaternion at time /„; and q(tn+i) is the propagated attitude quaternion at time /n + 1. The validity of Eq. (17-17) as a solution to the kinematic equation can be established by differentiation.

Equation (17-17) can be rewritten in a more convenient form for numerical computation by evaluating the matrix exponential using the procedures in Appen-

where

and 1 is the identity matrix.

For spacecraft gyros the sampling period is typically 100 to 400 nis. Because of the simplicity of Eq. (17-18), the closed-form solution has been used as, the kinematic integrator for onboard computers [Fish and Chmielewski, 1977]. The gyro system is usually a rate-integrating gyro package which provides an average angular velocity vector over the sampling period (see Section 6.5). The term inside the brackets in Eq. (17-18) is then an orthogonal rotation which retains the normalization of the propagated attitude quaternion.

It remains to illustrate the general validity of the closed-form solution and the computational error. To assess the latter, the quaternion q0„+i) is expanded in a Taylor series about the time tn

By repeated use of the kinematic equation (Eq. (17-2)), the Taylor series can be rewritten as

'-¿Û.B.+ ¿B.Û.J T>q(tn)+±T>nnq(tn)+ • • • (17-21)

The series of terms in the first bracket on the right-hand side of Eq. (17-21) is the Taylor series expansion of expf|i2B7"]. The remaining terms constitute the error introduced for a sampling period T in assuming a constant body rate equal to u„.

In general, the rates are not constant, and information from a rate-integrating gyro package can be used to form the 4x4 matrix

The terms in Eq. (17-21) can be rearranged to yield

The first of the two terms on the right-hand side is the Taylor series expansion of the closed-form expression

which differs from Eq. (17-17) in using time-averaged rate information rather than instantaneous body rates. It can also be written in the form of Eq. (17-18), with fi replacing Q„. Equation (17-23) shows that the error in this clostd-form expression is of order T3 and vanishes if 8n8n = or equivalently, if the vectors w„ and_tb„ are parallel. Thus, the order T3 correction to the closed-form expression using S2 is zero if the axis of rotation is fixed, even though the rates may be time dependent.

17.2 Environmental Torques

As described in Section 17.1, attitude prediction requires a model of the environmental disturbance torques acting on the spacecraft. To numerically integrate Euler's equations, the torque must be modeled as a function of time and the spacecraft's position and attitude. As listed in Table 1-2, the dominant sources of attitude disturbance torques are the Earth's gravitational and magnetic fields, solar radiation pressure, and aerodynamic drag.

Any nonsymmetrical object of finite dimensions in orbit is subject to a gravitational torque because of the variation in the Earth's gravitational force over the object. This gravity-gradient torque results from the inverse square gravitational force field; there would be no gravitational torque in a uniform gravitational field. General expressions for the gravity-gradient torque on a satellite of arbitrary shape have been calculated for both spherical [Nidey, 1960; Roberson, 1961; Hultquist, 1961] and nonspherical [Roberson, 1958b] Earth models. For most applications, it is sufficient to assume a spherical mass distribution for the Earth. If more accuracy is required, this may be obtained from the general potential function for the Earth given in Section 5.2. Alternatively, the effect of the Earth's oblateness can be accounted for in the motion of the orbital plane [Hultquist, 1961; Holland and Sperling, 1969].

In this section, we assume that the spacecraft's moment-of-inertia tensor is known for some arbitrary body reference frame whose origin need not coincide with the spacecraft's center of mass and that the spacecraft is orbiting a spherical Earth. The gravitational force dF, acting on a spacecraft mass element dm, located at a position R, relative to the geocenter is

Ki where /i = GM9 is the Earth's gravitational constant. The torque about the spacecraft's geometric center due to a force, dF„ at a position, r,, relative to the spacecraft's geometric center (see Fig. 17-1) is dN( = r, X dF( = (p + rj) XdF, (17-26)

The vector p is measured from the geometric center to the center of mass and the vector r| is measured from the center of mass to the mass element dm,.. The

Fig. 17-1. Coordinate System for the Calculation of Gravity-Gradient Torque

gravity-gradient torque on the entire spacecraft is obtained by integrating Eq. (17-26) to obtain

The geocentric position vector for the ith mass element can be expressed in terms of the geocentric position vector of the origin of the body reference frame, R,, as

For a practical artificial satellite R, = RJ + p + rJ»p + i<; therefore