How do SDP4's "Deep space" corrections to SGP4 account for the Sun's and Moon's gravity?

Question

The Simplified Perturbations model SGP4 is used to calculate Earth satellite state vectors (position and velocity) using standard ephemeris data encoded as TLEs (Two Line Elements). According to Wikipedia:

Current code libraries have merged SGP4 and SDP4 algorithms into a single codebase handling the range of orbital periods which are usually referred to generically as SGP4.

where SDP4 is the deep-space partner to the original SGP4, using only the simplest drag model but now also accounting for other perturbative effects, including the gravitational perturbations of the Moon and the Sun (as well as resonant effects near 1 and 2 orbits per day).

Published TLEs are calculated specifically to work with the appropriate SGP predictor. According to the original 1980/1988 version of Spacetrack Report No. 3, Models for Propagation of NORAD Element Sets:

All space objects are classified by NORAD as near-Earth (period less than 225 minutes) or deep-space (period greater than or equal 225 minutes). Depending on the period, the NORAD element sets are automatically generated with the near-Earth or deep-space model. The user can then calculate the satellite period and know which prediction model to use.

In SGP4 the initialization uses the TLE's mean motion to set a flag that determines which propagation method is used later in the execution. For example something along the lines of:

if ((2*pi / satrec.no) >= 225.0)
  {
    satrec.method = 'd';
    satrec.isimp  = 1;

or

IF((TWOPI/XNODP/XMNPDA) .GE. .15625) IDEEP=1

where 0.15625 is exactly 225/(24*60).

QUESTION: Can someone explain how SDP4 mathematically estimates the gravitational perturbations from the Sun and Moon? Does it contain a "mini-ephemeris" for the relative positions of the Sun, Earth, Moon system as a function of epoch, or at least their average periods, and propagate the satellite's motion including these forces, or does it use some average perturbation model?

note: I'm not looking for a general answer like "it uses perturbation theory", I'd like to know roughly how SGP4 actually does it.

Just for one particular example, in January the Sun will pull in one direction, but in July it will pull in the opposite direction. If the orbit is highly elliptical, does this matter for the perturbation calculation? Does it matte if the Sun pulls in the direction of periapsis, apoapsis or to the side?

SGP4 is also discussed in the 2006 report Revisiting Spacetrack Report #3: Rev 2.

Answer 1

I found Hujsak 1979 -- its title is "A Restricted Four Body Solution for Resonating Satellites Without Drag", but the linked pdf shows it is also Spacetrack Report #1. Four body means Earth, Moon, Sun, and satellite. There is a lot of integration going on: some things are averaged three or four times, over different periods with respect to different masses simultaneously. It's a typescript, not the final published version, so it can be challenging to read.

Hoots 1981, on the other hand, is the final version published in Celestial Mechanics, so I would start reading with this; click "print" to see more than one page at a time.

Hoots 1980, unlike the other links, is paywalled (only the first page is free). However, once I eventually got access, I discovered that this particular article states explicitly, "This paper is declared a work of the U.S. Government and therefore is in the public domain", so perhaps a complaint to the publisher (AIAA) is in order. A detailed description of its contents is now available in my answer to Differences between SGP8 and the standard SGP4? Is it ever used in practice?

Hujsak and Hoots together I haven't found, but since it is listed as "Aerospace Defense Command Space Computational Center Program Documentation", I didn't expect to.

There are also two master's theses from 1993 at the Naval Postgraduate School which describe parallelizing the SGP family of propagators: Ostrom and Brewer both revisit the same material as Spacetrack #3, but with more explanatory comments on which equation does what.

Answer 2

If you want a complete answer about how SGP4 (or SDP4) does it, you should read the original source: The Space Track Report #3. For the actual theory behind it, you should read the references by Hujsak and Hoots, for which I could not have access.

During the description of SDP4 routines, it says:

"At this point SDP4 calls the initialization section of DEEP which calculates all initialized quantities needed for the deep-space perturbations (see Section Ten).

(...)

Here SDP4 calls the periodics section of DEEP which adds the deep-space lunar and solar periodics to the orbital elements (see Section Ten). From this point on, it will be assumed that n, e, I, ω, Ω, and M are the mean motion, eccentricity, inclination, argument of perigee, longitude of ascending node, and mean anomaly after lunar-solar periodics have been added."

Section 10, however, provides the code but no equations for the deep-space routines. The Fortran code is full of constants and (apparently at least) dummy variables, which pollute the code. But there is also a well-known implementation in Matlab available on Celestrak's website, which I find more readable.

From the Matlab code, I can infer that there are no mini-ephemeris for either the sun or moon, as there is no computation involving julian date or universal time going on anywhere. In Earth resonance effects, the sidereal time is used, though. One comment in dspace routine claims:

"This procedure provides deep space contributions to mean elements for perturbing third body. These effects have been averaged over one revolution of the Sun and Moon. For Earth resonance effects, the effects have been averaged over no [sic] revolutions of the satellite. (Mean motion)."

Checking the dsinit function, you'll see that the mean elements are modified by adding their averaged time perturbation multiplied by the elapsed time since the TLE epoch, such as in:

em = em + dedt * t

Where the mean eccentricity is augmented by the perturbing effect. The factor dedt depends on several computations, but there is no trigonometric series involved (as in the past, computing trigonometric functions was done with parsimony).

So, you could check the references for a more accurate answer, but I'll leave with the conclusion that they add perturbations by some (very well scaled) average perturbing rate on mean elements, which afterwards are converted to osculating position and velocity with usual methods.