Trouble deriving rectangular components of acceleration of satellite in orbit around Earth with J2 consideration

Question

I am trying to derive the rectangular components of acceleration for a satellite in orbit, with Earth oblateness in consideration, in order to use the RK4 method to find the updated position and velocity of the satellite. I am aware that there are equations already out there, but I want to know how they were derived. However, my solutions don't exactly look like those equations - more importantly, the power of r in those sets of equations and in the wikipedia force equations are not the same. These equations have $\frac{1}{r^7}$, while mine have $\frac{1}{r^6}$

I will show how I derived the equations, please let me know if I've made an error!

I am using the following equation that I grabbed from the NASA Flight Test Cases document and Fonte 1993 (Implementation of 50x50 Gravity Field Model):

$$U = -\frac{\mu }{r}\left [ 1+\sum_{n=2}^{\infty}\sum_{m=0}^{n}\left(\frac{ R_{e}}{r}\right)^{n}P_{n,m}(sin\phi )\left(C_{n,m}cos(m\lambda) +S_{n,m}sin(m\lambda ) \right)\right]$$

I knew that for zonal harmonics $m=0$, so the equation looks like $$U = -\frac{\mu }{r}\left[1+\sum_{n=2}^{\infty}\left(\frac{R_{e}}{r}\right)^{n}P_{n,0}(sin\phi )(C_{n,0}cos(0) )\right]$$

I will then solve for $n=2$ knowing that $C_{n,0}=J_{2}$

There is a section below equation 10 under the deviations of Earth's gravitational field from that of a homogenous sphere section of the Geopotential Wiki page that says

$$u=\frac{J_{2}P_{2}^{0}(sin\theta )}{r^3}=J_{2}\frac{1}{r^3}\frac{1}{2}(3sin^{2}\theta-1)=J_{2}\frac{1}{r^5}\frac{1}{2}(3z^{2}-r^{2})$$

So I deduced that

$$P_{2,0}=\frac{1}{r^2}\frac{1}{2}(3z^{2}-r^{2})$$

So the equation looks like

$$U = -\frac{\mu }{r}\left[1+\left(\frac{R_{e}}{r}\right)^{2}(J_{2} )\left(\frac{1}{r^2}\frac{1}{2}\left(3z^{2}-r^{2}\right)\right)\right]$$

I then multiplied everything out and put the equation in this form

$$U = -\mu\left[\frac{1 }{r}+\frac{3}{2}J_{2}R_{E}^{2}z^{2}\frac{1 }{r^5}-\frac{1}{2}J_{2}R_{E}^{2}\frac{1 }{r^3}\right]$$

Knowing that $F_{x} = -\frac{\partial U}{\partial x}$, and using $\frac{\partial }{\partial x}r^{-5}=\frac{\partial }{\partial x}(x^{2}+y^{2}+z^{2})^{-5}=-\frac{10x}{r^6}$ etc., and then rearranging, I got the following for the x (and y) component:

$$F_{x} = -\frac{\partial U}{\partial x} = -\frac{\mu x}{r^2}\left[2+\frac{3J_{2}R_{E}^{2}}{r^2}\left(\frac{5z^{2}}{r^2}-1\right)\right]$$

Which looks somewhat familiar to this, but is not it exactly. I've tried doing that math three times now and get the same answer.

Sidenote: I asked OP of that thread, and he said that he got those equations from Section 7A of this competition

Can anyone tell me what I'm doing wrong?

Answer 1

Start with the OP's expression for the gravitational potential

$$U = -\mu\left[\frac{1 }{r}+\frac{3}{2}J_{2}R_{E}^{2}z^{2}\frac{1 }{r^5}-\frac{1}{2}J_{2}R_{E}^{2}\frac{1 }{r^3}\right]$$

can be rewritten as

$$U = -\mu\left[\frac{1 }{r} + \frac{1}{2}J_{2}R_{E}^{2} \left(3\frac{z^2}{r^5}-\frac{1}{r^3}\right)\right]$$

where $U$ is the reduced potential (the potential of a body of mass $m$ is then $mU$) and use

$$\mathbf{a}=\nabla U = \frac{\partial U}{\partial x} \mathbf{\hat{x}} + \frac{\partial U}{\partial y} \mathbf{\hat{y}} \frac{\partial U}{\partial z} \mathbf{\hat{z}} = a_x \mathbf{\hat{x}} + a_y \mathbf{\hat{y}} + a_z \mathbf{\hat{z}}$$

to get the acceleration. I'll do the x component with a little help from Wolfram alpha:

$$\frac{\partial}{\partial x} \frac{1}{r} = -\frac{x}{r^3},$$

$$\frac{\partial}{\partial x} \frac{z^2}{r^5} = -\frac{5xz^2}{r^7},$$

$$\frac{\partial}{\partial x} \frac{1}{r^3} = -\frac{3x}{r^5}.$$

Put those back into $\partial U/ \partial x$ and you get:

$$-\mu\left[ -\frac{x}{r^3} + \frac{x}{2}J_{2}R_{E}^{2} \left(-3 \frac{5z^2}{r^7} +\frac{3}{r^5}\right)\right]$$

$$-\mu\left[ -\frac{x}{r^3} + \frac{x}{2}J_{2}R_{E}^{2} \left(-3 \frac{5z^2}{r^7} +\frac{3r^2}{r^7}\right)\right]$$

$$-\mu\left[ -\frac{x}{r^3} + \frac{x}{2}J_{2}R_{E}^{2} \left(-3 \frac{5z^2}{r^7} +\frac{3(x^2 + y^2 + z^2)}{r^7}\right)\right]$$

At this point I'll leave it as an exercise for the reader to combine terms and arrive at

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

using $r^2=x^2+y^2+z^2$, the math and links in this answer plus the relationship between the dimensional and dimensionless (normalized, unitless) forms of $J_2$ described in For the mathematical relationship between J2 (km^5/s^2) and dimensionless J2 - which one is derived from the other? noting that $\mu$ is (at least in this case) just another name for the product $GM$.