Analytical expressions for acceleration due to zonal harmonics of a gravitational field?

Question

Wikipedia's Geopotential_model; The deviations of Earth's gravitational field from that of a homogeneous sphere discusses the expansion of the potential in spherical harmonics. The first few zonal harmonics ($\theta$ dependence only) are seen after the monopole term in

$$u = -\frac{GM}{r} - \sum_{n=2} J^0_n \frac{P^0_n(\sin \theta)}{r^{n+1}}$$

where $P^0_n$ are Legendre polynomials. I want to calculate the first three terms for $J_2, J_3, J_4$ by hand. I have

$$P^0_2(\sin \theta) = \frac{1}{2}(3 \sin^2 \theta - 1)$$

$$P^0_3(\sin \theta) = \frac{1}{2}(5 \sin^3 \theta - 3 \sin \theta)$$

$$P^0_4(\sin \theta) = \frac{1}{8}(35 \sin^4 \theta - 30 \sin^2 \theta + 3)$$

Since these terms are cylindrically symmetric I can write

$$\sin^2(\theta) = \frac{x^2+y^2}{r^2} = \frac{x^2+y^2}{x^2+y^2+z^2} $$

The $J_2$ term in the potential is then:

$$u_{J_2} = -J_2 \frac{1}{2} \frac{1}{r^3} \frac{3x^2 + 3y^2 - r^2}{r^2} = -J_2 \frac{1}{2} \frac{1}{r^5} (2x^2 + 2y^2 - z^2)$$

and the acceleration from this would be the negative gradient $-\nabla u$ or

$$\mathbf{a_{J_2}} = -\nabla u_{J_2}$$

Using this Wolfram Alpha link to make sure I don't make errors taking derivatives, I get (after a slight adjustment)

$$a_x = J_2 \frac{x}{r^7} \left( \frac{9}{2} z^2 - 3(x^2 + y^2) \right)$$

$$a_y = J_2 \frac{y}{r^7} \left( \frac{9}{2} z^2 - 3(x^2 + y^2)\right)$$

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

and these look very similar to but not the same as the results in Wikipedia's Geopotential_model; The deviations of Earth's gravitational field from that of a homogeneous sphere:

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

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

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

I'm close but I can't reproduce Wikipedia's result here. Once I'm confident with the process I can continue for the $J_3$ and $J_4$ terms and start doing numerical integration of orbits.

Answer 1

Let's look at @mmeent's comment suggesting that the spherical coordinates used in the linked Wikipedia article set the polar angle equal to zero at the equator rather than the pole.

where spherical coordinates (r, θ, φ) are used, given here in terms of cartesian (x, y, z) for reference

While that link shows $\theta = 0$ at the "north pole" (how I've usually seen spherical coordinates defined) the equations directly below that line do indeed define $\theta = 0$ to be the equator with $z=0$:

$$x = r \cos \theta \cos \phi$$

$$y = r \cos \theta \sin \phi$$

$$x = r \sin \theta$$

$$\sin^2(\theta) = \frac{z^2}{r^2} = \frac{z^2}{x^2+y^2+z^2} $$

then (noting that in the original question I'd put a minus sign where none existed):

$$u_{J_2} = +J_2 \frac{1}{2} \frac{1}{r^3} \frac{3z^2 - r^2}{r^2} = J_2 \frac{1}{2} \frac{1}{r^5} (2z^2 - (x^2 + y^2))$$

and using $\mathbf{a_{J_2}} = -\nabla u_{J_2}$ and Wolfram Alpha I get:

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

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

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

which agrees with Wikipedia.