Nonsingular Geopotential Models
The gravitational potential of a non-homogeneous celestial body at radius \(r\), latitude \(\phi\) , and longitude \(\lambda\) can be represented in spherical coordinates as:
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)
else
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
References
- S. Pines. "Uniform Representation of the Gravitational Potential and its Derivatives", AIAA Journal, Vol. 11, No. 11, (1973), pp. 1508-1511.
- 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.
- 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 86-FM-15, April 1986.
- W. M. Lear, "The Programs TRAJ1 and TRAJ2", JSC Internal Note 87-FM-4, 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", JSC-CN-23097, January 24, 2011.