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]]
degs, rads = 180./pi, pi/180.
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)
answers = []
for X0 in X0s:
answer, info = ODEint(deriv, X0, time, full_output=True)
answers.append(answer)
if 1 == 1:
names = 'x(t)', 'y(t)', 'z(t)', 'r(t)-a0'
answer = answers[1] # just choose one for the plot
posn = answer.T[:3]
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()
for answer in answers:
posn = answer.T[:3]
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.
