# Brouwer-Lyddane mean semi major axis bias

I'm using the software GMAT (General Mission Analysis Tool) (website, YouTube) for propagating spacecraft and calculating their Brouwer-Lyddane (BL) mean semi-major axis (SMA). I found out something interesting that I don't understand.

I propagate a spacecraft, with the initial orbital element:

• BL mean: [a,e,i,RAAN,AOP,MA] = [7000 km,0,45 deg,0,0,0]
• osculating elements: [a,e,i,RAAN,AOP,MA] = [7004.718 km , 0.000898, 45.019 deg, 0, 0, 0]
• and the epoch: 01 Jan 2000 12:00:00.000

Where the aforementioned osculating elements and the BL mean elements correspond to each other, by using the Brouwer-Lyddane Short transformation.

The propagator uses a Runge-Kutta 8th order integrator, with a force model of the Earth’s gravitational field JGM-2 with zonal and tesseral harmonics up to degree 21 and order 21. I don't use drag, sun radiation pressure, or third body gravity.

The produced BL mean SMA is:

Where the y axis is the BL SMA in [km], and the x axis is the ElapsedDays

When I propagate the same spacecraft with the same initial orbital elements, same propagator, but with a different epoch: 01 Jan 2000 00:00:00.000, I get the following BL SMA:

Notice that while the first epoch produces a significant bias value of approximately 130 meters, the second epoch does not.

Does any one has any explanation for this bias? Why is it epoch-depended? Why is it constant in time?

• In each case, are you specifying the BL mean elements or the osculating? – Chris Jul 1 '17 at 13:06
• Both. The osculating elements [7000 km,0,45 deg,0,0,0] correspond to the mean orbital elements [7004.718 km , 0.000898, 45.019 deg, 0, 0, 0] via the Brouwer-Lyddane Short transformation. Both cases are initialized with the same orbital elements. – Orbital Fun Jul 1 '17 at 14:05
• – uhoh Jul 1 '17 at 17:25
• I'm not familliar with the software, but if you start the satellite at the same everything except time (12 hours difference) then the satellite will be starting in a different gravitational potentials for the two cases. I am guessing that if you ran it 100 times spacing your starting Epochs evenly over 24 hours, you'd see their "altitude biases" would be spread above and below, depending on if you started in a gravity peak or valley. There should be no expectation that the mean should be equal to the instantaneous starting value, you're starting randomly. – uhoh Jul 1 '17 at 17:39
• That is what I also thought happened. It still does not explain why is this bias constant (and different for each epoch). I assumed that for long periods of time, the BL-SMA should be sometimes above and sometimes below the true value. But, as you can see, it is not. Depending on the epoch I choose, the BL-SMA is constantly over-biased or under-biased. – Orbital Fun Jul 1 '17 at 18:38

I'll expand on my comment:

I'm not familliar with the software, but if you start the satellite at the same everything except time (12 hours difference) then the satellite will be starting in a different gravitational potentials for the two cases. I am guessing that if you ran it 100 times spacing your starting Epochs evenly over 24 hours, you'd see their "altitude biases" would be spread above and below, depending on if you started in a gravity peak or valley. There should be no expectation that the mean should be equal to the instantaneous starting value, you're starting randomly.

To keep it simple, I'll use only the $J_2$ multipole component of the Earth's gravitational potential in addition to the monopole term $GM_E$. I've written a short python script using numbers and equations from Wikipedia (links shown within the script).

I've constructed a 90 degree inclination orbit, at an altitude of about 901 kilometers and propagated for 14 orbits or about 1 day. (I chose this altitude so I could set the inclination to at 99 degrees and confirm it was sun-synchronous by propagating for 91 days and watching the ascending node move from the x-axis to the y-axis, thereby checking that the gradient of the $J_2$ potential term is OK.)

I computed a semi-major axis and orbital speed $a_0$ and $v_0$ for a circular orbit in the monopole field $GM_E/r$ and used those for starting conditions.

The first plot shows $x(t), y(t), z(t)$ and $r(t)-a_0$ for 24 hours (14 orbits) starting on the equator with phase = 0. The calculation was repeated for phase = 30, 60, 90 degrees, such that the last one started above the pole, also at the same altitude of 901 kilometers.

In the second plot, I show the four plots of $r(t)-a_0$ which start at the same altitude, but different locations over the Earth. Right away you can see that the "bias" you speak of depends on the starting location and the particular gravitational potential that exists there.

Starting over the equator with phase = 0 degrees, where one is sitting lower in the Earth's gravity well, the speed is a little low for a circular orbit, and so the satellite is actually at apoapsis, dropping to a slightly lower altitude a half-period later.

However, starting over the pole with phase = 90 degrees, higher in Earth's gravity well, the speed is a little high for circular orbit, and so this is now periapsis, and the altitude rises a half-period later.

summary: By choosing an easier to understand gravitational field, the behavior of the "bias" is more systematic here. In your question, you're using a very complex, bumpy-looking gravity field which is constantly rotating. By starting at the same place in space, but different epochs, you were actually rotating the earth, moving different gravitational potentials to the starting location of your satellite.

I now predict that if you turn off the multipole field in JGM-2 and use a spherically symmetric potential, your calculations will show a flat altitude, no variations and no biases.

closing thoughts: Realistic orbits are not perfect conics, and so they and their Keplerian elements do not represent realistic orbits. They are only approximations to reality, and so are not right even though they are close.

Keplerian elements were used when people were writing with feathers using light from burning animal fat (if they were not busy being burned at the stake themselves). They are a mixed blessing in the 21st century when everything has so many more digits.

def deriv(X, t):

rr, vel = X.reshape(2, -1)

acc0, acc2 = accs(rr)

return np.hstack([vel, acc0 + acc2])

def accs(rr):

x,   y,   z   = rr
xsq, ysq, zsq = rr**2

rsq = (rr**2).sum()
rm3 = rsq**-1.5
rm7 = rsq**-3.5

acc0  = -GM_earth * rr * rm3

# https://en.wikipedia.org/wiki/Geopotential_model#The_deviations_of_Earth.27s_gravitational_field_from_that_of_a_homogeneous_sphere
acc2x = x * rm7 * (6*zsq - 1.5*(xsq + ysq))
acc2y = y * rm7 * (6*zsq - 1.5*(xsq + ysq))
acc2z = z * rm7 * (3*zsq - 4.5*(xsq + ysq))

acc2  = J2_earth * np.hstack((acc2x, acc2y, acc2z))

return acc0, acc2

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

halfpi, pi, twopi = [f*np.pi for f in [0.5, 1, 2]]

GM_earth = 3.986004418E+14 #  m^3/s^2  https://en.wikipedia.org/wiki/Standard_gravitational_parameter
J2_earth = 1.7555E+25      #  m^5/s^2  https://en.wikipedia.org/wiki/Geopotential_model

a0  = (901 + 6378) * 1E+03 # meters (roughly) for a n=14 sun-synchronous orbit

v0  = np.sqrt(GM_earth/a0)  # m/s  (vis-viva)
T   = twopi * a0 / v0       # sec

print T/60, ' minutes'

phases = [rads*d for d in [0, 30, 60, 90]]

X0s = []
for phase in phases:
sp, cp = np.sin(phase), np.cos(phase)

X0 = np.hstack(( [ a0*cp, 0, a0*sp],   # initial positions
[-v0*sp, 0, v0*cp] )) # initial velocitys
X0s.append(X0)

n, days = 14, 1.  # orbits/day, days

time = np.linspace(0, n*days*T, int(n*days*100)+1)

for X0 in X0s:

answer, info = ODEint(deriv, X0, time, full_output=True)

if 1 == 1:

names  = 'x(t)', 'y(t)', 'z(t)', 'r(t)-a0'
rdiff  = np.sqrt((posn**2).sum(axis=0)) - a0
things = [thing for thing in posn] + [rdiff]

plt.figure()
for i, (thing, name) in enumerate(zip(things, names)):
plt.subplot(4, 1, i+1)
plt.plot(time/(24.*3600), thing)   # x-axis in days, not seconds
plt.title(name, fontsize=16)
plt.xlim(0, days+0.01)  # a little more than n days
plt.show()

if 1 == 1:

plt.figure()

rdif   = np.sqrt((posn**2).sum(axis=0)) - a0

plt.plot(time/(24.*3600), rdif)   # x-axis in days, not seconds

plt.xlim(0, days+0.1)  # a little more than n days
plt.text(1.01, -6000, u' 0\u00B0', fontsize=15 )
plt.text(1.01,  2000, u'30\u00B0', fontsize=15 )
plt.text(1.01, 16000, u'60\u00B0', fontsize=15 )
plt.text(1.01, 22000, u'90\u00B0', fontsize=15 )
plt.title('r(t)-a0', fontsize=16)
plt.show()


below: a quick check running for 91 days with an inclination of 99 degrees to confirm the orbit is roughly sun-synchronous in order to make sure there are at least no gross errors with the way I typed the equations for the gradient of the $J_2$ component, or the numerical values. The script above does not allow for inclination, but you'd include it by mixing the $v_y$ and $v_z$ values when X0 is calculated.

• Thank you for the informative response. After I read your response, and also consulted with other space enthusiasts, I think I understand now. But still, let me rephrase it in my own words: In complex gravitational problems, the BL mean elements, as well as the osculating elements, are functions, and not constants. So, by selecting different epochs or different initial true anomalies, we change the gravity or the initial phase with respect to the beginning of the period. So, the initial BL SMA does not necessarily corresponds to the true mean value. but rather to a some value in the period. – Orbital Fun Jul 3 '17 at 18:42
• @OrbitalFun it's a really good question and has (obviously) given me a lot to think about! Thanks also for your answers (1, 2) to the other question... you've given more to think about and I need to dig in there next. – uhoh Jul 4 '17 at 4:37