How could one determine the acceleration of a satellite given its ephemeris data, taking into account the earth's oblateness?

Question

I am looking to determine the acceleration vector of a satellite orbiting the earth at a certain point in time given its ephemeris data for that same point in time. I am taking into account the effects of an oblate earth, but not of third-body forces, aerodynamic forces, or solar radiation pressure.

To find the ephemeris data for a satellite, I use JPL's HORIZONS Web-Interface (https://ssd.jpl.nasa.gov/horizons.cgi), and I use the VECTORS ephemeris type, which gives me the following information for a given point in time (key included as well): (I apologize for including this information for anyone that is already familiar with this tool, as I only wanted to make the data format I am working with more accessible)

Information given for one point in time:
JDTDB
X     Y     Z
VX    VY    VZ
LT    RG    RR

Key (units can be changed):
JDTDB    Julian Day Number, Barycentric Dynamical Time
  X      X-component of position vector (au)
  Y      Y-component of position vector (au)
  Z      Z-component of position vector (au)
  VX     X-component of velocity vector (au/day)                           
  VY     Y-component of velocity vector (au/day)                           
  VZ     Z-component of velocity vector (au/day)                           
  LT     One-way down-leg Newtonian light-time (day)
  RG     Range; distance from coordinate center (au)
  RR     Range-rate; radial velocity wrt coord. center (au/day)

How would I solve for the acceleration vector given this information?

I know that assuming a spherical earth, the magnitude of the acceleration is given by

|Ā| = GM / R^2 (where R is the distance between the satellite and the center of the earth),

but would that work for an oblate earth? In addition how can I determine the direction of the vector? Assuming a spherical earth, the direction of the acceleration vector would be opposite to that of the position vector, but I do not think that would work either.

I appreciate the help! Also, if my problem is beyond the scope of a StackExchange question, I would appreciate references to source material that I could look into.

Answer 1

Units

I aways use kilograms, meters and seconds (MKS) to avoid getting tripped up by mixed units. Horizons offers both AU & AU/day and km & km/sec options and the first thing I do when importing ephemeris data is convert the km & km/sec to m & m/sec!

This is from this JPL Horizons tutorial

Central force ("monopole term")

The main term in the spherical harmonic expansion of Earth's gravity that's associated with oblateness is called $J_2$. It's the second term, the first term is the central force "the main gravity" that I also call the monopole term.

The acceleration of a body in the gravitation field of another body of standard gravitational parameter $GM$ can be written:

$$\mathbf{a_{Central}} = -GM \frac{\mathbf{r}}{|r|^3},$$

where $r$ is the vector from the body $M$ to the body who's acceleration is being calculated. Remember that the acceleration of each body depends only on the mass of the other body, even though the force depends on both masses, because the first mass cancels out by $a=F/m$.

Written out as components that would be

$$\mathbf{a_{Central}} = -GM \frac{x}{|r|^3}\mathbf{\hat{x}} -GM \frac{y}{|r|^3}\mathbf{\hat{y}} -GM \frac{z}{|r|^3}\mathbf{\hat{z}}$$

In Python I write:

a = -GM * X * np.sqrt((X**2).sum())**-1.5

where X is a position vector.

This is from this answer to How to calculate the planets and moons beyond Newtons's gravitational force? as is the quote below

Oblateness ($J_2$ only):

I'm just using the math from Wikipedia's article on the Geopotential Model with a very important-to-remember approximation; I am assuming that the oblateness is in the plane of the ecliptic — that the rotational axis of the oblate body is in the $\mathbf{\hat{z}}$ direction, perpendicular to the ecliptic. Don't forget that this is an approximation! [...]

$$\mathbf{r} = x \mathbf{\hat{x}} + y \mathbf{\hat{y}} + z \mathbf{\hat{z}} $$

$$a_x = J_2 \frac{x}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$

$$a_y = J_2 \frac{y}{|r|^7} (6z^2 - 1.5(x^2+y^2)) $$

$$a_z = J_2 \frac{z}{|r|^7} (3z^2 - 4.5(x^2+y^2)) $$

The following should be added to the Newtonian monopole term:

$$\mathbf{a_{J2}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}} $$

Unless you are doing very precise work, it's probably safe to assume that the oblateness is in the plane of the ecliptic; that Earth's axis is still pointing in the same direction that it was in the year 2000. It does slowly precess but how to deal with that mathematically is beyond the scope of this answer.