Transforming Mean and Osculating Elements Using Numerical Methods

Mean element propagation of perturbed two body orbits has as its mathematical basis the averaging theory of nonlinear dynamical systems. Mean elements define an orbit’s long-term evolution characteristics consisting of both secular and long-period effects. Using averaging theory, a near-identity transformation can be found that transforms between the mean elements and their osculating counterparts that augment the mean elements with short period effects. The ability to perform the conversion is necessary so that orbit design conducted in either mean elements or osculating can be effectively converted between each element type. In the present work, the near-identity transformation is found using the Fast Fourier Transform. An efficient method is found that is capable of recovering the mean or osculating elements to first-order.


Introduction
Mean element theories have proven useful for efficient determination of an orbit's long-term characteristics; however, a mean orbit propagation often requires transforming a set of osculating initial conditions to mean to properly initialize the propagation. Conversely, once a mean orbit is obtained converting the results back to osculating elements requires the inverse transformation. Fortunately the mean element theory provides an available mechanism for finding a transformation that is consistent with the applied averaging technique. Indeed, the process to obtain the mean elements using averaging produces the transformation equations directly. This is a well-known result (cf. Sanders et al. [1], Chapter 7 on Averaging Over Angles) that will be reviewed later in this paper, and applied to the present context of perturbed two-body orbits. There is an extensive body of work applying averaging theory to satellite astrodynamics, which has traditionally focused on deriving analytic or semi-analytic theories to find satellite mean orbital elements and the associated mean-to-osculating or osculating-to-mean transformations. Two key contributions that have found wide use are Draper's Semi-analytic Satellite Theory (DSST), documented extensively by Danielson et al. [2], McClain [3], McClain et al. [4], and the work by Guinn [5]. These prior works obtain the transformations by using analytical techniques to compute explicit formulae for the Fourier series that correspond to the short period terms of the orbit being examined. In the present study, numerical techniques based on the Fast Fourier Transform (FFT) are utilized for deriving the Fourier series coefficients. This has the advantage of being able to find the short period terms due to any perturbing acceleration without the need for an explicit derivation of the associated Fourier coefficients. Nor are there any approximations present due to truncating eccentricity expansions. Prior research by Lutzky and Uphoff [6] utilized a numerical integration technique to find the Fourier Series coefficients; however, their work focused on directly integrating the integrals defining the Fourier coefficients. The present study takes advantage of progress in FFTs to more efficiently find the coefficients. Furthermore, unlike Lutzky and Uphoff's work, specific issues with regard to mean motion resonances are examined.
In a prior paper written by the author (Ely [7]) the subject of mean element propagations using numerical averaging is presented. In that work, averaging the osculating element equations of motion to obtain their mean element counterparts is given without derivation. Presently, a near-identity transformation is used to derive the formal conversion between osculating and mean elements, and, as a consequence, the mean element equations of motion are also found; thus, providing a formal justification for the averaging method used in Ely [7]. The near-identity transformation method is then applied to specific orbital perturbations that, in some cases, require special modifications depending on the perturbation's form.

Basic Theory of the Near-identity Transformation
The trajectory of a satellite in two-body perturbed motion can be modeled using a variation of parameter formulation of six 1 st -order differential equations of motion (EOMs). In a basic form with only one fast angle, perhaps some slowly varying angles, and no resonances (the case of mean motion resonances will be examined later) the EOMs can be represented as where the osculating elements have been selected to be the equinoctial set {α, λ} with α = [a, h, k, p, q], and the mean longitude λ, which has been singled out because it is the short periodic ('fast') angle in this problem. As seen in Eq. (1), the mean longitude varies with a rate that is of zeroth-order O (1) = O ε 0 because of its dependence on the Keplerian mean motion n (a) = μ a 3 . The small parameter ε defines first-order O (ε) and will be used to order terms in subsequent asymptotic expansions. For orbit applications, ε is dependent on the specific perturbation; two key examples include the zonal harmonic J 2 or the relative distance of a perturbing body. In the formal analysis that follows the presence of ε is an aid in tracking the order of individual terms, but might not be explicitly present in the actual application. The variables θ in Eq. 1 are either long-period or secular and slowly vary with rates that are themselves of first-orderθ = O (ε) or higher relative to the mean motion. The inclusion of θ in Eq. 1 generalizes the presentation given in Ely [7]. Typical examples of this include: the orbital rates of other bodies in their apparent motion around the central body, Venus' rotation rate, slowly varying atmospheres, etc. These slow variables are assumed to be known functions of time, thus, if present in a problem, make Eq. 1 non-autonomous. The EOMs in Eq. 1 are in a standard variation of parameters formulation, and, as proven in Lemma 7.8.3 in Sanders et al. [1], there exists a near-identity transformation that eliminates the mean longitude's short periodic contribution to the EOMs to produce an averaged system of equations that approximates the original system to first-order O (ε) on time scales of order O 1 ε . The set of averaged equations of motion take the form ᾱ where it has been assumed that the osculating variables {α, λ} can be decomposed into mean and short periodic components The mean elements ᾱ,λ , represented with the overbar, exhibit long-period and/or secular behavior, and the short periodic variations around the mean elements { α s , λ s } have been eliminated from Eq. 2 via the near-identity transformation. A few observations about (2) are in order: 1. The mean mean longitudeλ exhibits a combination of secular growth from a leading dependence on the mean motionn, an O (1)-rate, and long-period variations from both g λ and an evolving mean semimajor axisā vian (ā).
2. By assumption, the angles θ only vary with long-period motion or secularly with O (ε)-rates, thus their mean and osculating values are equivalent. That is,θ = θ .
The perturbing functions g α and g λ in Eq. 2 are currently unknown and will now be determined along with the solution for the near-identity transformation. The functional form of the first-order near-identity transformation from mean-toosculating can be represented functionally as where, in what follows, it is assumed that the mean elements at the current epoch are known. Later this assumption will be switched, and the inverse problem (finding the mean elements from the osculating) will be solved. The procedure for finding the solution to the first-order transformation function u 1 is now outlined; a similar procedure for finding v 1 also applies. Substitute (4) into (1) and asymptotically expandα Differentiating u 1 produces Usingα = εg α and collecting the first-order terms leads to A convenient choice for the unknown function g α is the average of f α defined as follows where the function f α is evaluated overλ ranging from λ ,λ + 2π with the mean elements ᾱ,λ and θ held constant at their current values. The underlying orbit used to numerically compute f α conforms to the unperturbed Keplerian orbit defined by ᾱ,λ at the current epoch. Using this definition for g α leads to the following partial differential equation for u 1 Doing the same analysis for the v 1 equation yields where the term u 1 1 represents the first component of the vector u 1 and corresponds to the expression associated with the semimajor axis. Equations 9 and 10 represent the homological partial differential equations for finding the near-identity transformations of u 1 and v 1 , respectively. Solutions to these equations can be obtained by assuming that the functions f α , f λ can be expanded in a complex Fourier series as functions ofλ (11) where i ≡ √ −1 , and the coefficients defined as It follows from Eq. 12 that the zero-frequency terms (j = 0) are simply the average values f α (ᾱ, θ ) ,f λ (ᾱ, θ ) . Note that for real valued perturbation functions the Fourier coefficients with negative indices are complex conjugates of their positive index counterparts, f −j (α, θ ) = f j (α, θ ) * , with the * representing complex conjugation. Substituting (11) into Eq. 9 and integrating yields the following formal solution for u 1 (13) where the integration constant is zero because the average of the transformation has been defined to be zero, that isū 1 ᾱ, θ ,λ ≡ 0. Likewise, substituting (11) into Eq. 10 leads to The near identity transformation can now be found by substituting (13) and (14) into Eq. 4 to obtain The transformation from osculating-to-mean elements is found by transposing (4), substituting the transposition into the right hand side, and performing an asymptotic expansion accurate to first-order as follows where the terms in Eq. 16 are now evaluated with osculating elements.
In the actual applications that follow the small parameter ε will be included as part of the expression for the perturbing functions f α , f λ , in these cases the ordering parameter ε can be set to a value of one in the near identity transformations without loss of generality because its effect will be included in the computation of Fourier series coefficients. It should be emphasized that in Eq. 15 the perturbing functions f α , f λ are evaluated using the current known mean element values ᾱ,λ to obtain the Fourier coefficients f α j , f λ j , while in Eq. 16 the problem is inverted and the current osculating elements are known {α, λ} and the perturbing functions are evaluated using these values to obtain their associated Fourier coefficients. Often with initial value problems the osculating elements are known at the initial epoch, (16) enables an analyst to find the corresponding mean element initial conditions for conducting a mean element study that is consistent with the selected osculating initial conditions. Likewise, in a mean element trajectory propagation the analyst may desire the associated osculating trajectory; Eq. 15 provides a method to obtain this using known mean elements.

Computing the Fourier Coefficients
In a departure from other techniques that find the Fourier coefficients analytically (cf. Danielson et al. [2] or Guinn [5]) the coefficients in the present work are found numerically using the Discrete Fourier Transform (DFT). As developed in Sundararajan [8] and adapted presently for the perturbing functions discussed previously, the DFTs F α j ,F λ j at frequency multiple j for continuous vector functions where N is the number of samples taken of f α , f λ . Again, the overbars on F α j ,F λ j are used to represent DFTs computed using mean element evaluations of the osculating perturbing functions. The DFTs F α j ,F λ j can be related to the Fourier series coefficients f α j , f λ j by approximating their integral definitions in Eq. 12 using a rectangular integration rule that has divided the 2π interval into N rectangles each with a width of 2πk N, that is and similarly Comparing the results in Eq. 18 with the definition for the DFTs in Eq. 17 it is found that where the indices have been limited to j ≤ N 2 because, as Sundararajan [8] in Chapter 12 shows, larger values of j would introduce aliasing effects and the approximation will fail as proven by the Nyquist sampling theorem. The range in Eq. 20 has been extended to the negative indices using the property of real valued signals that the Fourier series coefficients and the DFT terms with negative indices are complex conjugates of their positive index counterparts [8]. Indeed, in proving this result and comparing to the DFT definition it can be shown which is useful for relating the results returned by the numerical DFT routines with the desired Fourier coefficients.
Summarizing, N samples recover N/2 frequency multiples of jλ. If a particular problem's associated Fourier series possesses a large number of short periodic terms with significant amplitudes, then a large N will be required to recover them satisfactorily. In the examples that follow some heuristics will be discussed for choosing an appropriate N.
Substituting (20) into Eqs. 13 and 14 yields approximate relationships for the mean-to-osculating transformations using the DFTs as follows The computation of the inverse transformations from osculating-to-mean given in Eq. 16 follows the same procedure as outlined above, and is presented now for completenessᾱ where the DFTs are evaluated with the osculating elements, hence, no overbar is present. The DFTs are found efficiently using a Fast Fourier Transform (FFT) algorithm that is usually implemented in most mathematical processing packages. In the present case, the FFTs are computed using the GNU Scientific Library package [9]. Each call to the FFT routine computes a set of the DFTs. For example, in a mean-to-osculating transformation the following DFTs are found when the FFT routine is provided a set of numerical samples of the functions f α , f λ evaluated with the current mean orbital elementsᾱ and the mean longitudeλ k at N times between λ ,λ + 2π using a sampling interval of 2π/N.
Example Consider an orbiter at Venus with conditions similar to the Magellan orbiter as given in Table 1. This problem includes perturbations from the zonal and tesseral harmonics and the Sun. Since Venus' rotational rate is approximately 243 days there is no possibility of a significant mean motion resonance and hence the sidereal angle can be treated as a secular slow variable. Likewise, since Venus' orbital period is 0.6 years, the apparent mean longitude of the Sun's perturbation can also be treated as  [7] is used to find the mean elements as a function of time, which are then converted to osculating elements using Eq. 22. The one day propagation is done in two ways, 1) standard propagation of the osculating equations of motion; 2) the transformation technique just described to convert a mean orbit to osculating.
In a first case the Sun's perturbations are turned off, hence the only active perturbations are from the zonal and tesseral harmonics. The leading harmonic perturbation is due to J 2 with a magnitude of ∼ 4.4E − 6; however, many of the zonals and tesserals included in the 10 x 10 field have magnitudes approaching J 2 and provide significant contributions to the transformations. In the preceding theory all of the zonal and tesseral harmonics are treated as first-order perturbations, and identification of an explicit coefficient for ε is not necessary because evaluating the DFTs includes the full effect of these included terms, thus ε has been set to one. An illustration of the recovery of the semi-major axis is shown in Fig. 1 where the numerically averaged and integrated mean elements are shown in red ('mean'), the 1 st -order numerical conversion of the mean elements to osculating is shown in green ('mean2osc'), and the direct numerical integration of the osculating equations of motion are shown in blue ('osc'). The order of the FFT used is N = 128 (this selection will be discussed shortly). The recovery of the osculating elements from the mean is quite good, with no perceptible differences seen between 'mean2osc' and 'osc' (i.e., the 'green' line overlies the 'blue').
To see this comparison more clearly and to examine the effect that the selected FFT order has on the results, differences of semi-major axis for orders N = 32, 64, 128, 256 are shown in Fig. 2. The results illustrate the importance of selecting a proper order to ensure a sufficient spectrum is recovered with the FFT. In this case, N = 32 yields semimajor axis differences on the order of 80 m or less while N = 128 produces differences that are significantly less at 6 cm or less. Increasing the order to N = 256 yields no substantial improvement; thus N = 128 yields a sufficient order in this example. Currently, the process for selecting a proper order is done by trial and error; however, a starting point is to examine the highest order term in the included gravity field and, via the sampling theorem, a minimum FFT order would need to be at least twice that order to properly capture the highest order gravity harmonic signature, in the present case that is 20. Other factors that also influence this selection are the rate at which the magnitudes of the gravity fields Fourier coefficients decrease in value, as well as the orbits initial eccentricity where the presence of significant harmonic terms is directly correlated to eccentricity [12]. Furthermore, FFT evaluations are most efficient for orders that are powers of 2. These considerations have led to the selection of 128 for the current problem. Fortunately, once an appropriate order has been found for a particular problem it remains relatively constant when analyzing nearby orbits for parametric studies.
In the next case the Sun's perturbations are activated in addition to the zonal/tesseral harmonics. The inclusion of these terms yields nearly identical results as with the prior case, and is illustrated in Fig. 3. Indeed, there is no visually perceptible difference in the results relative to the prior case. An explicit computation of the difference in the two errors (without the Sun as shown in the lower left plot of Fig. 2 and with the Sun as shown in the right plot of Fig. 3) at one day yields a value of ∼ 2 mm, small enough that the difference between these errors is insignificant. This case illustrates that even with multiple and independent O (ε) slow angles (in this case Venus' sidereal rotation and its orbit around the Sun) the theory yields satisfactory results without special treatment. As will be shown later, if additional angles in the problem (beyond the mean longitude) vary with rates that approach O (1), then other methods are required.

Extensions to the Basic Theory
The preceding basic theory applies to perturbations from zonal harmonics, drag, solar radiation pressure, and, in certain cases, with tesseral harmonics and other body perturbations. In addition to the satellite mean longitude λ, the tesseral harmonics are functionally dependent on the central body's sidereal angle φ and other body perturbations are dependent on the other body's mean longitude λ O . These dependencies introduce extensions to the basic theory that will now be outlined.

Tesseral Harmonic Perturbations
The equations of motion with the tesseral harmonics have the following structure ⎡ ⎣α λ where the sidereal rate equationφ ≡ ω is to first order a constant and has increased the dimension of the problem by one. It also follows from the form ofφ thatφ = φ. The current analysis can benefit by examining the multiple periodic Fourier series expansion of the tesseral harmonic contribution to a central body's gravity field, where the superscript T is used to represent the tesseral harmonics, and M is the maximum order of the tesserals in the gravity field being examined (which is equal to or less than the maximum degree N of the field). Note that the inclusion of terms m < 0 is not traditional (cf. Danielson [2], Section 2.7.1 or Kaula [12], Chapter 3); however, the more generic form of the Fourier series expansion presented here includes these terms to aid in identification with the DFT. Without loss of generality ε has been set to 1, as the small parameter is included in the Fourier coefficients f T jm (α, θ ). In an analytical averaging theory, the coefficients f T jm (α, θ ) are often infinite expansions, as well, that need to be carefully truncated. For a numerical approach to averaging and near-identity transformation the coefficients are not needed analytically; however, of importance is the relationship between the competing angles in the problem (the mean longitude and the sidereal angle), which is exposed by the formal expansion in Eq. 25. The presence of the sidereal angle introduces three distinct perturbation problems that will be separately considered: 1. Resonant Case: The mean longitude rateλ (∼ mean motion n) is O(1) and the sidereal rateφ is O(1) and the initial conditions are such that a deep resonance or shallow resonance exists. Deep resonance exists when the mean longitude rate and the sidereal rate are nearly commensurate. That is, there exists a rational ratio of two integers Q and P that satisfẏ λ φ ∼ Q P → Q : P resonance with Q and P integers.
For instance, at Earth a resonance of Q : P = 2 : 1 is a ∼ 12 hr orbit. A shallow (or near) resonance occurs whenλ/φ are close to Q/P in some sense (specific for each resonance). For the resonance to be significant Q and P should be smaller integers. That is a 2:1 resonance that excites gravity harmonics of order 2 will have a more significant effect than a 30:1 resonance that excites the typically much smaller order 30 harmonics. 2. Non-resonant Case: The mean longitude rateλ (∼ mean motion n) is O (1) and the sidereal rateφ is O(1) and the initial conditions are such that no significant resonances exist (i.e., an orbit that is outside of a shallow resonance region). In this case the integers that satisfy (26) are too large to be significant. The conditions that need to exist to make this determination are the subject of Kolmogorov-Arnold-Moser (KAM) theory of nonlinear dynamical systems theory [13]. 3. Adiabatic (or Two Time Scales) Case: The mean longitude rateλ (∼ mean motion n) is O(1) and the sidereal rateφ is O(ε). This case matches the conditions developed in the previous section on the basic theory where the sidereal angle is now a particular instance of θ . The sidereal angle rate is obviously nonresonant (the integers which satisfy (26) are too large to be significant), and is slow enough that it can be considered a constant during an averaging interval. This applies to Venus (as illustrated in the prior example); sometimes the Moon, and often other bodies that are tidally locked in a 1:1 spin resonance.
A detailed discussion of each of these cases follows.

Resonant Case
For deep resonance or shallow resonance, the commensurability of the mean motion with the planet rotation rate introduces a new slow variable that can be revealed by introducing the resonant node ψ. It is defined as This definition is consistent with the equinoctial elements in that it is non-singular for orbits that are circular or have inclinations below 180 • . However, for missions where repeating groundtracks are relevant, a more useful angle is one that accounts for the precession of the orbital plane and the line of apsides when orbiting an oblate central body. This form of the node, called the stroboscopic node ψ s by Gedeon [14], is defined using where it should be noted that this becomes undefined for equatorial orbits (in this event use of the resonant node, given in Eq. 27, should be made).
Using the resonant node ψ definition given in Eq. 27, the transformation of Eq. 24 into the new variable yields ⎡ ⎣α ψ where, in resonance or near resonance, Q and P have been selected to satisfy the following propertyψ Note that the change of variables is generic to the system with all perturbations active, no restriction to considering only tesseral harmonics in Eq. 29 has been made. Hence if other perturbations (i.e., zonal harmonics, other bodies, drag, and solar pressure) are active they will participate in the new system of equations as indicated in Eq. 29. This change of variables has eliminated the O (1) dependence on the isolated mean motion, and now all of the components of the right hand side in Eq. 29 are O (ε) or higher. The near-identity transformation for this system takes the form Equation 29 has the same structural form as Eq. 1, thus the same procedure for finding the homological equations can be used with the following result for the α elementsn and for the mean longitudē The expansion for the node reveals that Hence, only the mean longitude homological equation v 1 needs to be solved and then scaled by P to get w 1 for the node. In Ely [7] it was shown that the formal Fourier series expansion for the perturbing rate equations due to the tesseral harmonics with the change of variable to the resonant node takes the form where for convenience of notation exp () has been used for the exponential function rather than e () , f T jm (ᾱ, θ ) ≡ f α,T jm (ᾱ, θ ) , f λ,T jm (ᾱ, θ ) , andψ (evaluated using the current mean values ᾱ,λ ) is itself considered a slow quantity that is held fixed in later FFT evaluations. Additionally, the variables θ might represent slowly varying orientation terms, or some other long periodic gravity effect. The surviving harmonics from the averaging procedure include Subtracting (36) from Eq. 35 produces the Fourier series expression with the average With these observations (37) can be recast into a generic one-dimensional form as follows where k ≡ jQ − mP , k = 0 follows from jQ = mP , and f T k (ᾱ, θ ,ψ) identifies the formal Fourier coefficients that are functionally dependent on only the long-period variables as a result of the change in variables. This conforms more closely to the previous results and the DFT can be readily applied. Using Eq. (41) the transformation function u 1 can be recast with the mean-to-osculating numerical transformation on α taking the form The DFT is computed by evaluating f α,T ᾱ,θ ,λ k , φ k = f α,T ᾱ,θ ,λ k , Pλ k −ψ /Q where the sidereal angle has been constrained using φ = Pλ −ψ /Q and the mean mean longitudeλ k takes N values between λ ,λ + 2πQ using a sampling interval of 2πQ/N. The span 2πQ is present becauseλ must vary this amount for the independent angleλ/Q to vary a full 2π period. A few notes on these evaluations will follow, but first performing a similar set of operations for v 1 yields the mean-to-osculating numerical transformation on λ to be To perform the operations indicated in Eqs. (42) and (43) it should be noted that: 1. The evaluation period of 2πQ is, in fact, a full sidereal period and is sufficient to recover the m-'daily' oscillatory terms. Because the FFT needs to capture frequencies that include multiples of the mean longitude and multiples of the m-'dailies', the order of the FFT needed to capture all of these increases as compared to the zonal-only case. This will be illustrated later in an example. 2. The subtraction of the average rate terms of the Fourier series in Eq. (37) is not explicitly done, when the FFT returns the coefficients setting the components that correspond to the average term to zero yields the equivalent result. That is, computation of the mean rate is never required for the osculating/mean transformation -only the osculating rates are needed as functions of the mean or osculating elements -depending on the direction of the transformation.
The process for finding the osculating-to-mean transformation is the same as before; transpose Eqs. (42) and (43), and replace the mean elements with their osculating counterparts in the function evaluations.
In the preceding analysis if the stroboscopic node ψ s had been selected for analysis rather than the resonant node ψ the results in Eqs. (42) and (43) would still apply; however, the sidereal angle would now be constrained to be φ = P λ − −ψ s Q + . The choice of using ψ s or ψ depends on the application being considered. In the example to follow ψ s has been selected since the orbit is near polar, if the orbit being investigated were near equatorial ψ would have been a preferred choice.
Example Consider a Mars orbiter that is in a Q : P = 13 : 1 resonance with Mars' rotation rate with initial conditions and other parameters as indicated below in Table 2. These elements are similar to those for the Mars Reconnaissance Orbiter; however, the semimajor axis has been adjusted to put the orbit in a deep resonance. Two cases will be examined. In both cases, perturbations from the Sun, drag, and solar radiation pressure are included. In the first case only the zonal harmonics are active, and in the second the tesserals are active along with the zonals. The inclination recovery for both of these cases is shown in Fig. 4, where the highest FFT order N necessary to sufficiently recover the short periodic terms is shown. In the case with only the zonal harmonics present, a FFT of order 128 proves sufficient; however, the inclusion of the tesserals increases the necessary order to 1024. The reason for this can be seen by the presence of a significant m-'daily' signature when the tesseral harmonics are active. That is, there is a daily period for the short periodic inclination that is not present in the zonal harmonic case. To adequately capture periods on the order of a Martian day or less requires the FFT to carry more terms than is needed when sampling over only an orbital period, in this case 1024. It turns out that the number of terms needed is a function of the particular resonance with higher resonances typically requiring more terms than lower order resonances. For instance, a simulation of a 2:1 resonance case (not shown) indicates that a 64 th -order FFT is sufficient.  . 4 Comparison of inclination recovery of the mean-to-osculating transformation with zonal harmonics (Case 1) and zonal plus tesseral harmonics (Case 2) and the highest order of N necessary for the FFT to sufficiently recover the short periodic terms.

Non-resonant Case
In this case, the appropriate averaging operator is a double average over the two angles with a result that the non-resonant tesseral harmonics produce no long-period or secular effects on the orbit [7]. That is, the tesserals introduce only short period terms that are multiples of the mean longitude and the sidereal angle. In this case, the near-identity transformation must be capable of recovering both of these angles. To facilitate the analysis, the functional form of the EOMs with non-resonant tesserals active can be represented as where the mean longitude and the sidereal angle are both considered 'fast' variables. The procedure for finding the near-identity transformation follows the same argument as before. However, the Fourier series expansions no longer have a single independent variable; they are multiple periodic Fourier expansions. It is a simple exercise to show that the homological equations governing the system in Eq. (44) take the form The mean rates f α (ᾱ, θ ) ,f λ (ᾱ, θ ) would now be computed using double averages over both (λ, φ) and, when considering only the non-resonant tesseral harmonics with no zonals present, it can be shown thatf α (ᾱ, θ ) = 0 andf λ (ᾱ, θ ) = 0 [7]. Regarding the other perturbations, there is no dependence on the sidereal angle, so, with no loss in generality, the subsequent analysis will be restricted to the tesseral harmonics only. The multiple periodic Fourier series expansion of the tesseral harmonics in Eq. (25) applies with no transformation. A direct substitution of Eq. (25) into Eq. (45) and integrating leads to the following solution for u 1 Comparing Eq. (46) to Eq. (13) an immediate realization can be made, the multidimensional result in Eq. (45) can become singular (hence divergent) for terms with jn−mω ≈ 0, while the scalar case in Eq. (13) is never singular. This result illustrates the famous 'small divisor' problem in perturbation theory and is the raison d'être for the celebrated Kolmogorov-Arnol'd-Moser (KAM) theorem of nonlinear dynamical systems theory [13]. The restriction in this case of being non-resonant has a more formal definition in the KAM theorem, which states for the selected initial conditions if there exists constants γ > 0 and τ > 0 that satisfy the inequality then the series in Eq. (46) is convergent and non-singular. Satellite orbits that are sufficiently far from resonance should satisfy the condition in Eq. (47) with determining sufficiency typically found via analysis of selected cases. If the analyst is in doubt whether the condition is satisfied for his/her problem, then the transformation to the stroboscopic node in Case 1 will always work; however, the efficiency of a mean element propagation may suffer (via the numerical integrator having to compensate for stroboscopic node that has too large a frequency). Continuing under the assumption that Eq. (47) is satisfied, the transformation function for the mean longitude v 1 is found by substituting the result in Eq. (46) and Eq. (25) into Eq. (45) and integrating, which yields (48) Now the near-identity transformation from mean-to-osculating elements can be found by evaluating the transformation at the current mean elementsᾱ, mean mean longitudeλ, and sidereal angle φ. The multiple periodic Fourier coefficients in Eq. (25) with independent mean longitudeλ and sidereal angle φ can be found using the GNU Scientific Library's multiple periodic FFT capabilities via sampling the functions f α,T (ᾱ, θ , λ, φ) , f λ,T (ᾱ, θ , λ, φ) from λ ,λ + 2π with steps 2π/N and [φ, φ + 2π ) with steps 2π/M * (where typically M * should be selected to be ≥ M). This results in the following transformations and the relationship between the DFTs and the Fourier coefficients conform to To first-order, the transformation from osculating to mean elements can be found by transposing Eq. (49) and evaluating the right hand side with osculating elements.
Example Here the semimajor axis from the prior example has been increased by 50 km to place the orbit well outside of the 13:1 resonance, but still far from other major resonances (such as the 12:1 resonance). The double angle non-resonant tesseral harmonic conversion just described applies. The FFT evaluations sample the mean longitude over an orbital period using the same 128 th order as was done for the zonal-only case, and the sidereal angle is sampled over its period at 64 th order. The mean-to-osculating and osculating inclination histories and errors for a one-day propagation are shown in Fig. 5. The transformation error (right) is on par (and slightly better) than the resonant case shown previously. However, this comes at a larger computational expense. Each two angle FFT requires 8192 (= 128 x 64) function calls whereas the prior resonance case needed only 1024 calls to achieve similar error performance.

Adiabatic (or Two-Time Scales) Case
In the case whereφ λ there is no possibility of any significant resonances; hence, Eq. (25) can be examined directly without transformation. Since the sidereal angle's rate is slow it leads to long-period terms and the angle is an instance of θ that is found in Eq. (1). Because the functional form of this case matches the original form outlined in the Basic Theory section the results in Eq. (22) apply for the mean-toosculating transformation, and the results in Eq. (23) apply for the osculating-tomean transformation.

Other Body Perturbations
As with the tesseral harmonics, the Fourier expansion of the other body perturbations results in a multiple periodic expansion that formally can be represented as An example of this expansion is found in McClain [3] in which the explicit dependence on both the satellite's mean longitude λ and the other body's mean longitude λ O (of its apparent orbit around the central body) are revealed. Compared to Eq. (25), this result is very similar to the tesseral case with a few subtle differences that will now be explored. For most practical situations there are no significant resonances between the two mean longitude rates. Examination of Eq. (51) shows that there exist some combinations of {j, k} that result in a resonance. But, as j or k increase, the coefficients f O jk (α, θ ) diminish in magnitude for a convergent Fourier series. These higher resonances are not significant in most situations on time scales of order O (1/ε) (again these statements can be made more precise using the KAM theorem). Hence, resonances with other bodies will not be considered further. With this elimination, there are two distinct cases that need examination.

Adiabatic Single Average Case
The other body mean longitude rate is sufficiently slowλ O λ such thatλ = O (ε). As noted in Ely [7], this case typically involves a single average over the satellite's mean longitude λ, and results in a functional form that matches the basic theory previously outlined with λ O ∈ θ . Hence, the results in Eq. (22) apply for the meanto-osculating transformation, and the results in Eq. (23) apply for the osculating-tomean transformation.

Non-Adiabatic Double Average Case
Here, the other body mean longitude rate is still slowλ O < λ; however, the order of the rate is not distinctly between zeroth and first; that is, O(1)<λ O < O(ε). In this scenario a second average over the other body's mean longitude λ O may be warranted. This is similar in form to the case with non-resonant tesseral harmonics, and the present transformation follows the same argument as was presented for the nonresonant tesseral harmonics. The resulting transformation for the other body double averaged perturbations takes the form To first-order, the transformation from osculating to mean elements can be found by transposing Eq. (53) and evaluating the right hand side with osculating elements.

Conclusions
The preceding work has shown the utility of a numerical FFT-based approach for converting between mean element orbits and their osculating counterparts. The method is efficient and does not require the explicit computation of an analytic Fourier series; rather the coefficients are obtained via direct calculations of the osculating element rate functions. This generality allows the method to be applied to a large set of perturbing accelerations, which for this study includes the zonal and tesseral harmonics, other bodies, drag, and solar radiation pressure. The following cases have been examined,