Nonsingular Geopotential Models

sphereThe gravitational potential of a non-homogeneous 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 [2-7], mostly at the Johnson Space Center. The Fortran 77 code in Mueller [2] and Lear [4-5] 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)
        ci1 = zero
        si1 = zero
    end if

    rho(i)   = r0*rho(i-1)
    creal(i) = ess*creal(i-1) - t*cimag(i-1)
    cimag(i) = ess*cimag(i-1) + t*creal(i-1)
    d(i,1)   = ess*ci1 + t*si1
    e(i,1)   = ci1
    f(i,1)   = si1
    a(i,i)   = (2*i-1)*a(i-1,i-1)
    a(i,i-1) = 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(k-1) + s(i,k)*cimag(k-1)
            f(i,k) = s(i,k)*creal(k-1) - c(i,k)*cimag(k-1)
        end if

        !JW : typo here in original paper
        ! (should be GOTO 1, not GOTO 10)
        if (i/=2) then
            L = i-2
            do j=1,L
                a(i,i-j-1) = (u*a(i,i-j)-a(i-1,i-j))/(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


  1. S. Pines. “Uniform Representation of the Gravitational Potential and its Derivatives“, AIAA Journal, Vol. 11, No. 11, (1973), pp. 1508-1511.
  2. A. C. Mueller, “A Fast Recursive Algorithm for Calculating the Forces due to the Geopotential (Program: GEOPOT)”, JSC Internal Note 75-FM-42, June 9, 1975.
  3. J. L. Spencer, “Pines Nonsingular Gravitational Potential Derivation, Description and Implementation“, NASA Contractor Report 147478, February 9, 1976.
  4. W. M. Lear, “The Gravitational Acceleration Equations”, JSC Internal Note 86-FM-15, April 1986.
  5. W. M. Lear, “The Programs TRAJ1 and TRAJ2”, JSC Internal Note 87-FM-4, April 1987.
  6. 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.
  7. R. A. Eckman, A. J. Brown, and D. R. Adamo, “Normalization of Gravitational Acceleration Models“, JSC-CN-23097, January 24, 2011.

Tagged with: , , , ,

Leave a Reply

Your email address will not be published. Required fields are marked *