The gravitational potential of a nonhomogeneous celestial body at radius \(r\), latitude \(\phi\) , and longitude \(\lambda\) can be represented in spherical coordinates as:
\( U = \frac{\mu}{r} [ 1 + \sum\limits_{n=2}^{n_{max}} \sum\limits_{m=0}^{n} \left( \frac{r_e}{r} \right)^n P_{n,m}(\sin \phi) ( C_{n,m} \cos m\lambda + S_{n,m} \sin m\lambda ) ] \)Where \(r_e\) is the radius of the body, \(\mu\) is the gravitational parameter of the body, \(C_{n,m}\) and \(S_{n,m}\) are spherical harmonic coefficients, \(P_{n,m}\) are Associated Legendre Functions, and \(n_{max}\) is the desired degree and order of the approximation. The acceleration of a spacecraft at this location is the gradient of this potential. However, if the conventional representation of the potential given above is used, it will result in singularities at the poles \((\phi = \pm 90^{\circ})\), since the longitude becomes undefined at these points. (Also, evaluation on a computer would be relatively slow due to the number of sine and cosine terms).
A geopotential formulation that eliminates the polar singularity was devised by Samuel Pines in the early 1970s [1]. Over the years, various implementations and refinements of this algorithm and other similar algorithms have been published [27], mostly at the Johnson Space Center. The Fortran 77 code in Mueller [2] and Lear [45] is easily translated into modern Fortran. The Pines algorithm code from Spencer [3] appears to contain a few bugs, so translation is more of a challenge. A modern Fortran version (with bugs fixed) is presented below:
subroutine gravpot(r,nmax,re,mu,c,s,acc) use, intrinsic :: iso_fortran_env, only: wp => real64 implicit none real(wp),dimension(3),intent(in) :: r !position vector integer,intent(in) :: nmax !degree/order real(wp),intent(in) :: re !body radius real(wp),intent(in) :: mu !grav constant real(wp),dimension(nmax,0:nmax),intent(in) :: c !coefficients real(wp),dimension(nmax,0:nmax),intent(in) :: s ! real(wp),dimension(3),intent(out) :: acc !grav acceleration !local variables: real(wp),dimension(nmax+1) :: creal, cimag, rho real(wp),dimension(nmax+1,nmax+1) :: a,d,e,f integer :: nax0,i,j,k,l,n,m real(wp) :: rinv,ess,t,u,r0,rhozero,a1,a2,a3,a4,fac1,fac2,fac3,fac4 real(wp) :: ci1,si1 real(wp),parameter :: zero = 0.0_wp real(wp),parameter :: one = 1.0_wp !JW : not done in original paper, ! but seems to be necessary ! (probably assumed the compiler ! did it automatically) a = zero d = zero e = zero f = zero !get the direction cosines ess, t and u: nax0 = nmax + 1 rinv = one/norm2(r) ess = r(1) * rinv t = r(2) * rinv u = r(3) * rinv !generate the functions creal, cimag, a, d, e, f and rho: r0 = re*rinv rhozero = mu*rinv !JW: typo in original paper rho(1) = r0*rhozero creal(1) = ess cimag(1) = t d(1,1) = zero e(1,1) = zero f(1,1) = zero a(1,1) = one do i=2,nax0 if (i/=nax0) then !JW : to prevent access ci1 = c(i,1) ! to c,s outside bounds si1 = s(i,1) else ci1 = zero si1 = zero end if rho(i) = r0*rho(i1) creal(i) = ess*creal(i1)  t*cimag(i1) cimag(i) = ess*cimag(i1) + t*creal(i1) d(i,1) = ess*ci1 + t*si1 e(i,1) = ci1 f(i,1) = si1 a(i,i) = (2*i1)*a(i1,i1) a(i,i1) = u*a(i,i) do k=2,i if (i/=nax0) then d(i,k) = c(i,k)*creal(k) + s(i,k)*cimag(k) e(i,k) = c(i,k)*creal(k1) + s(i,k)*cimag(k1) f(i,k) = s(i,k)*creal(k1)  c(i,k)*cimag(k1) end if !JW : typo here in original paper ! (should be GOTO 1, not GOTO 10) if (i/=2) then L = i2 do j=1,L a(i,ij1) = (u*a(i,ij)a(i1,ij))/(j+1) end do end if end do end do !compute auxiliary quantities a1, a2, a3, a4 a1 = zero a2 = zero a3 = zero a4 = rhozero*rinv do n=2,nmax fac1 = zero fac2 = zero fac3 = a(n,1) *c(n,0) fac4 = a(n+1,1)*c(n,0) do m=1,n fac1 = fac1 + m*a(n,m) *e(n,m) fac2 = fac2 + m*a(n,m) *f(n,m) fac3 = fac3 + a(n,m+1) *d(n,m) fac4 = fac4 + a(n+1,m+1)*d(n,m) end do a1 = a1 + rinv*rho(n)*fac1 a2 = a2 + rinv*rho(n)*fac2 a3 = a3 + rinv*rho(n)*fac3 a4 = a4 + rinv*rho(n)*fac4 end do !gravitational acceleration vector: acc(1) = a1  ess*a4 acc(2) = a2  t*a4 acc(3) = a3  u*a4 end subroutine gravpot
References
 S. Pines. “Uniform Representation of the Gravitational Potential and its Derivatives“, AIAA Journal, Vol. 11, No. 11, (1973), pp. 15081511.
 A. C. Mueller, “A Fast Recursive Algorithm for Calculating the Forces due to the Geopotential (Program: GEOPOT)”, JSC Internal Note 75FM42, June 9, 1975.
 J. L. Spencer, “Pines Nonsingular Gravitational Potential Derivation, Description and Implementation“, NASA Contractor Report 147478, February 9, 1976.
 W. M. Lear, “The Gravitational Acceleration Equations”, JSC Internal Note 86FM15, April 1986.
 W. M. Lear, “The Programs TRAJ1 and TRAJ2”, JSC Internal Note 87FM4, April 1987.
 R. G. Gottlieb, “Fast Gravity, Gravity Partials, Normalized Gravity, Gravity Gradient Torque and Magnetic Field: Derivation, Code and Data“, NASA Contractor Report 188243, February, 1993.

R. A. Eckman, A. J. Brown, and D. R. Adamo, “Normalization of Gravitational Acceleration Models“, JSCCN23097, January 24, 2011.
Leave a Reply