Degenerate Conic

Conversion Factors

Conversion Factors

• 1 lbm = 0.45359237 kg
• 1 lbf = 4.4482216152605 N
• 1 ft = 0.3048 m
• 1 mile = 1.609344 km
• 1 nmi = 1.852 km
• 1 slug = 1 \(\mathrm{lbf} ~ \mathrm{s}^2 / \mathrm{ft} \approx\) 14.5939029 kg

References

• NIST Special Publication 811: Guide for the use of the International System of Units.
• Metric system cartoons

Fortran Code

module conversion_factors

use, intrinsic :: iso_fortran_env, only: wp => real64 !double precision

implicit none

real(wp),parameter :: one = 1.0_wp

real(wp),parameter :: lbm2kg = 0.45359237_wp     ! exact
real(wp),parameter :: lbf2N = 4.4482216152605_wp !
real(wp),parameter :: ft2m = 0.3048_wp           !
real(wp),parameter :: mile2km = 1.609344_wp      !
real(wp),parameter :: nmi2km = 1.852_wp          !
real(wp),parameter :: slug2kg = lbf2N/ft2m       ! ~ 14.593902937206362

real(wp),parameter :: kg2lbm = one/lbm2kg        ! ~ 2.2046226218487757
real(wp),parameter :: N2lbf = one/lbf2N          ! ~ 0.2248089430997105
real(wp),parameter :: m2ft = one/ft2m            ! ~ 3.280839895013123
real(wp),parameter :: km2mile = one/mile2km      ! ~ 0.621371192237334
real(wp),parameter :: km2nmi = one/nmi2km        ! ~ 0.5399568034557235
real(wp),parameter :: kg2slug = one/slug2kg      ! ~ 0.06852176585679176

end module conversion_factors

json-fortran 4.0.0

I just tagged the 4.0.0 release of json-fortran. This is the first release with Unicode support (thanks to Izaak Beekman). Who says there are no good open source Fortran projects on the internet?

The Unicode build of the library is optional (and only enabled using the preprocessor directive USE_UCS4). Currently, this only works with Gfortran. It doesn't yet work with the Intel Fortran Compiler, which is lagging behind on Unicode support. The Fortran standard supports Unicode via the selected_char_kind function, which can be used to specify the character set used for a character string, like so:

integer,parameter :: u = selected_char_kind('ISO_10646')
character(kind=u,len=11) :: string = u_'Hello World'

Earth-Mars Free Return

Let's try using the Fortran Astrodynamics Toolkit and Pikaia to solve a real-world orbital mechanics problem. In this case, computing the optimal Earth-Mars free return trajectory in the year 2018. This is a trajectory that departs Earth, and then with no subsequent maneuvers, flies by Mars and returns to Earth. The building blocks we need to do this are:

1. An ephemeris of the Earth and Mars [DE405]
2. A Lambert solver [1]
3. Some equations for hyperbolic orbits [given below]
4. An optimizer [Pikaia]

The \(\mathbf{v}_{\infty}\) vector of a hyperbolic planetary flyby (using patched conic assumptions) is simply the difference between the spacecraft's heliocentric velocity and the planet's velocity:

• \(\mathbf{v}_{\infty} = \mathbf{v}_{heliocentric} - \mathbf{v}_{planet}\)

The hyperbolic turning angle \(\Delta\) is the angle between \(\mathbf{v}_{\infty}^{-}\) (the \(\mathbf{v}_{\infty}\) before the flyby) and the \(\mathbf{v}_{\infty}^{+}\) (the \(\mathbf{v}_{\infty}\) vector after the flyby). The turning angle can be used to compute the flyby periapsis radius \(r_p\) using the equations:

• \(e = \frac{1}{\sin(\Delta/2)}\)
• \(r_p = \frac{\mu}{v_{\infty}^2}(e-1)\)

So, from the Lambert solver, we can obtain Earth to Mars, and Mars to Earth trajectories (and the heliocentric velocities at each end). These are used to compute the \(\mathbf{v}_{\infty}\) vectors at: Earth departure, Mars flyby (both incoming and outgoing vectors), and Earth return.

There are three optimization variables:

• The epoch of Earth departure [modified Julian date]
• The outbound flight time from Earth to Mars [days]
• The return flyby time from Mars to Earth [days]

There are also two constraints on the heliocentric velocity vectors before and after the flyby at Mars:

• The Mars flyby must be unpowered (i.e., no propulsive maneuver is performed). This means that the magnitude of the incoming (\(\mathbf{v}_{\infty}^-\)) and outgoing (\(\mathbf{v}_{\infty}^+\)) vectors must be equal.
• The periapsis flyby radius (\(r_p\)) must be greater than some minimum value (say, 160 km above the Martian surface).

For the objective function (the quantity that is being minimized), we will use the sum of the Earth departure and return \(\mathbf{v}_{\infty}\) vector magnitudes. To solve this problem using Pikaia, we need to create a fitness function, which is given below:

subroutine obj_func(me,x,f)

!Pikaia fitness function for an Earth-Mars free return

use fortran_astrodynamics_toolkit

implicit none

class(pikaia_class),intent(inout) :: me !pikaia class
real(wp),dimension(:),intent(in) :: x !optimization variable vector
real(wp),intent(out) :: f !fitness value

real(wp),dimension(6) :: rv_earth_0, rv_mars, rv_earth_f
real(wp),dimension(3) :: v1,v2,v3,v4,vinf_0,vinf_f,vinfm,vinfp
real(wp) :: jd_earth_0,jd_mars,jd_earth_f,&
dt_earth_out,dt_earth_return,&
rp_penalty,vinf_penalty,e,delta,rp,vinfmag

real(wp),parameter :: mu_sun = 132712440018.0_wp !sun grav param [km3/s2]
real(wp),parameter :: mu_mars = 42828.0_wp !mars grav param [km3/s2]
real(wp),parameter :: rp_mars_min = 3390.0_wp+160.0_wp !min flyby radius at mars [km]
real(wp),parameter :: vinf_weight = 1000.0_wp !weights for the
real(wp),parameter :: rp_weight = 10.0_wp ! penalty functions

!get times:
jd_earth_0 = mjd_to_jd(x(1)) !julian date of earth departure
dt_earth_out = x(2) !outbound flyby time [days]
dt_earth_return = x(3) !return flyby time [days]
jd_mars = jd_earth_0 + dt_earth_out !julian date of mars flyby
jd_earth_f = jd_mars + dt_earth_return !julian date of earth arrival

!get earth/mars locations (wrt Sun):
call get_state(jd_earth_0,3,11,rv_earth_0) !earth at departure
call get_state(jd_mars, 4,11,rv_mars) !mars at flyby
call get_state(jd_earth_f,3,11,rv_earth_f) !earth at arrival

!compute lambert maneuvers:
call lambert(rv_earth_0,rv_mars, dt_earth_out*day2sec, mu_sun,v1,v2) !outbound
call lambert(rv_mars, rv_earth_f,dt_earth_return*day2sec,mu_sun,v3,v4) !return

!compute v-inf vectors:
vinf_0 = v1 - rv_earth_0(4:6) !earth departure
vinfm = v2 - rv_mars(4:6) !mars flyby (incoming)
vinfp = v3 - rv_mars(4:6) !mars flyby (outgoing)
vinf_f = v4 - rv_earth_f(4:6) !earth return

!the turning angle is the angle between vinf- and vinf+
delta = angle_between_vectors(vinfm,vinfp) !turning angle [rad]
vinfmag = norm2(vinfm) !vinf vector mag (incoming) [km/s]
e = one/sin(delta/two) !eccentricity
rp = mu_mars/(vinfmag*vinfmag)*(e-one) !periapsis radius [km]

!the constraints are added to the fitness function as penalties:
if (rp>=rp_mars_min) then
    rp_penalty = zero !not active
else
    rp_penalty = rp_mars_min - rp
end if
vinf_penalty = abs(norm2(vinfm) - norm2(vinfp))

!fitness function (to maximize):
f = - ( norm2(vinf_0) + &
    norm2(vinf_f) + &
    vinf_weight*vinf_penalty + &
    rp_weight*rp_penalty )

end subroutine obj_func

The lambert routine called here is simply a wrapper to solve_lambert_gooding that computes both the "short way" and "long way" transfers, and returns the one with the lowest total \(\Delta v\). The two constraints are added to the objective function as penalties (to be driven to zero). Pikaia is then initialized and called by:

program flyby

use pikaia_module
use fortran_astrodynamics_toolkit

implicit none

integer,parameter :: n = 3

type(pikaia_class) :: p
real(wp),dimension(n) :: x,xl,xu
integer :: status,iseed
real(wp) :: f

!set random number seed:
iseed = 371

!bounds:
xl = [jd_to_mjd(julian_date(2018,1,1,0,0,0)), 100.0_wp,100.0_wp]
xu = [jd_to_mjd(julian_date(2018,12,31,0,0,0)),400.0_wp,400.0_wp]

call p%init(n,xl,xu,obj_func,status,&
            ngen = 1000, &
            np = 200, &
            nd = 9, &
            ivrb = 0, &
            convergence_tol = 1.0e-6_wp, &
            convergence_window = 200, &
            initial_guess_frac = 0.0_wp, &
            iseed = iseed )

call p%solve(x,f,status)

end program flyby

The three optimization variables are bounded: the Earth departure epoch must occur sometime in 2018, and the outbound and return flight times must each be between 100 and 400 days.

Running this program, we can get different solutions (depending on the value set for the random seed, and also the Pikaia settings). The best one I've managed to get with minimal tweaking is this:

• Earth departure epoch: MJD 58120 (Jan 2, 2018)
• Outbound flight time: 229 days
• Return flight time: 271 days

Which is similar to the solution shown in [2].

References

1. R.H. Gooding, "A procedure for the solution of Lambert's orbital boundary-value problem", Celestial Mechanics and Dynamical Astronomy, Vol. 48, No. 2, 1990.
2. Dennis Tito's 2021 Human Mars Flyby Mission Explained. SPACE.com, February 27, 2014.
3. J.E. Prussing, B.A. Conway, "Orbital Mechanics", Oxford University Press, 1993.

Pikaia

I've started a new project on Github for a new modern object-oriented version of the Pikaia genetic algorithm code. This is a refactoring and upgrade of the Pikaia code from the High Altitude Observatory. The original code is public domain and was written by Paul Charbonneau & Barry Knapp (version 1.2 was released in 2003). The new code differs from the old code in the following respects:

• The original fixed-form source (FORTRAN 77) was converted to free-form source.
• The code is now object-oriented Fortran 2003/2008. All user interaction is now through a class, which is the only public entity in the module.
• All real variables are now 8 bytes (a.k.a. "double precision"). The precision can also be changed by the user by modifying the wp parameter.
• The random number generator was replaced with calls to the intrinsic Fortran random_number function.
• There are various new options (e.g., a convergence window with a tolerance can be specified as a stopping condition, and the user can specify a subroutine for reporting the best solution at each generation).
• Mapping the variables to be between 0 and 1 now occurs internally, rather than requiring the user to do it.
• Can now include an initial guess in the initial population.

This is another example of how older Fortran codes can be refactored to bring them up to date with modern computing concepts without doing that much work. In my opinion, the object-oriented Fortran 2003 interface is a lot better than the old FORTRAN 77 version. For one, it's easier to follow. Also, data can now easily be passed into the user fitness function by extending the pikaia class and putting the data there. No common blocks necessary!

References

1. Original Pikaia code and documentation [High Altitude Observatory]

Fortran Stream I/O

Fortran 2003 introduced many new features to the language that make things a lot easier than they were in the bad old days. Deferred-length character strings and stream I/O are two examples. Deferred-length strings were a huge addition, since they allow the length of character variables to be resized on-the-fly, something never before possible in Fortran. (Note: if you happen to come across a Fortran 90 monstrosity called ISO_VARYING_STRING, avoid it like the plague. That is not what I'm talking about.) Stream I/O is also quite handy, allowing you to access files as a stream of bytes or characters. The following example uses both features to create a subroutine that reads in the contents of a text file into an allocatable string. Note that it uses the size argument of the intrinsic inquire function to get the file size before allocating the string.

subroutine read_file(filename,str)  
!  
! Reads the contents of the file into the allocatable string str.  
! If there are any problems, str will be returned unallocated.  
!  
implicit none

character(len=*),intent(in) :: filename  
character(len=:),allocatable,intent(out) :: str

integer :: iunit,istat,filesize

open(newunit=iunit,file=filename,status='OLD',&  
form='UNFORMATTED',access='STREAM',iostat=istat)

if (istat==0) then  
    inquire(file=filename, size=filesize)  
    if (filesize>0) then  
        allocate( character(len=filesize) :: str )  
        read(iunit,pos=1,iostat=istat) str  
        if (istat/=0) deallocate(str)  
        close(iunit, iostat=istat)  
    end if  
end if

end subroutine read_file  

References

1. Stream Input/Output in Fortran [Fortran Wiki]
2. Text file to allocatable string [Intel Fortran Forum] 03/03/2015

Multidimensional B-Spline Interpolation

I just started a new modern Fortran software library called bspline-fortran, which is for multidimensional (multivariate) b-spline interpolation of data defined on a regular grid. It is available on GitHub, and released under a permissive BSD-style license.

It seems impossible to find code for higher than 3D spline interpolation on the internet. Einspline only has 1D-3D, as do the NIST Core Math Library DBSPLIN and DTENSBS routines. Python's SciPy stops at 2D (Bivariate splines). At first glance, the SciPy routine map_coordinates seems to do what I want, but not really since there doesn't seem to be any way to specify the abscissa vectors (x, y, z, etc.). Note: maybe there is a way, but looking at the less-than-great documentation, the minimally-commented C code, and then reading several confusing Stack Overflow posts didn't help me much. Apparently, I'm not the only one confused by this routine.

After much searching, I finally just took the 2D and 3D NIST routines (written in 1982), refactored them significantly so I could better understand what they were doing, and then extended the algorithm into higher dimensions (which is actually quite easy to do for b-splines). The original FORTRAN 77 code was somewhat hard to follow, since it was stuffing matrices into a giant WORK vector, which required a lot of obscure bookkeeping of indices. For example, I replaced this bit of code:

    IZM1 = IZ - 1  
    KCOLY = LEFTY - KY + 1  
    DO 60 K=1,KZ  
      I = (K-1)*KY + 1  
      J = IZM1 + K  
      WORK(J) = DBVALU(TY(KCOLY),WORK(I),KY,KY,IDY,YVAL,INBV1,WORK(IW))  
60  CONTINUE  

with this:

kcoly = lefty - ky + 1  
do k=1,kz  
    temp2(k) = dbvalu(ty(kcoly:),temp1(:,k),ky,ky,idy,yval,inbv1,work)  
end do  

which makes it a lot easier to deal with and add the extra dimensions as necessary.

The new library contains routines for 2D, 3D, 4D, 5D, and 6D interpolation. It seems like someone else may find them useful, I don't know why they don't seem to be available anywhere else. Eventually, I'll add object-oriented wrappers for the core routines, to make them easier to use.

Comments

I just came across this article from 2005 written by Jef Raskin (probably most famous for having initiated the Macintosh project at Apple) on source code documentation [1]. I agree with much of what he recommends. I have never believed that code can be "self-documenting" without comments, and generally will start writing the comments before I start the actual code. I also like some aspects of Knuth's "literate programming" concept [2].

I do find automatic documentation generators useful, however. There aren't a lot of solutions for modern Fortran projects though. I've not tried to use Doxygen, but am told that support for modern Fortran features is somewhat lacking [3]. ROBODoc is pretty good, although it does not perform any code analysis, and only extracts specific comment blocks that you indicate in the source. A brand new tool called ford is very promising, since it was created with modern Fortran in mind.

I also have an idea that the best code (including comments) should be understandable to a reasonably intelligent person who is not familiar with the programming language. Some languages (Fortran, Python) make this easier than others (C++, Haskell). There is a good article here, discussing the "beauty" of some C++ code (although I disagree with his views on comments). I was struck by the discussion on how to write functions so that more information is provided to the reader of the code [4]. The example given is this:

int idSurface::Split( const idPlane &plane, const float epsilon,  
                      idSurface **front, idSurface **back,  
                      int *frontOnPlaneEdges, int *backOnPlaneEdges ) const;  

Of course, if you are not familiar with the C++ *, **, and & characters and the const declaration, there is no information you can derive from this. It is complete gibberish. Modern Fortran, on the other hand, provides for the declaration of an intent attribute to function and subroutine arguments. The equivalent Fortran version would look something like this:

integer function split(me,plane,eps,front,back,&  
                       frontOnPlaneEdges,backOnPlaneEdges)

class(idSurface),intent(in) :: me  
type(idPlane), intent(in) :: plane  
real, intent(in) :: eps  
type(idSurface), intent(out) :: front  
type(idSurface), intent(out) :: back  
integer, intent(out) :: frontOnPlaneEdges  
integer, intent(out) :: backOnPlaneEdges  

Here, it should be more obvious what is happening (for example, that plane is an input and front is an output). The use of words rather than symbols definitely improves readability for the non-expert (although if taken to the extreme, you would end up with COBOL).

References

1. J. Raskin, "Comments are More Important than Code", Queue, Volume 3 Issue 2, March 2005, Pages 64-65.
2. D. E. Knuth, "Literate Programming", The Computer Journal (1984) 27 (2): 97-111.
3. F03/08 supporting Documentation tools, comp.lang.fortran.
4. S. McGrath, "The Exceptional Beauty of Doom 3's Source Code", kotaku.com, Jan. 14, 2013.

IAU Rotation Models

The IAU Working Group on Cartographic Coordinates and Rotational Elements (WGCCRE) is the keeper of official models that describe the cartographic coordinates and rotational elements of planetary bodies (such as the Earth, satellites, minor planets, and comets). Periodically, they release a report containing the coefficients to compute body orientations, based on the latest data available. These coefficients allow one to compute the rotation matrix from the ICRF frame to a body-fixed frame (for example, IAU Earth) by giving the direction of the pole vector and the prime meridian location as functions of time. An example Fortran module illustrating this for the IAU Earth frame is given below. The coefficients are taken from the 2009 IAU report [Reference 1]. Note that the IAU models are also available in SPICE (as "IAU_EARTH", "IAU_MOON", etc.). For Earth, the IAU model is not suitable for use in applications that require the highest possible accuracy (for that a more complex model would be necessary), but is quite acceptable for many applications.

module iau_orientation_module

use, intrinsic :: iso_fortran_env, only: wp => real64

implicit none

private

!constants:
real(wp),parameter :: zero = 0.0_wp
real(wp),parameter :: one = 1.0_wp
real(wp),parameter :: pi = acos(-one)
real(wp),parameter :: pi2 = pi/2.0_wp
real(wp),parameter :: deg2rad = pi/180.0_wp
real(wp),parameter :: sec2day = one/86400.0_wp
real(wp),parameter :: sec2century = one/3155760000.0_wp

public :: icrf_to_iau_earth

contains

!Rotation matrix from ICRF to IAU_EARTH
pure function icrf_to_iau_earth(et) result(rotmat)

implicit none

real(wp),intent(in) :: et !ephemeris time [sec from J2000]
real(wp),dimension(3,3) :: rotmat !rotation matrix

real(wp) :: w,dec,ra,d,t
real(wp),dimension(3,3) :: tmp

d = et * sec2day !days from J2000
t = et * sec2century !Julian centuries from J2000

ra = ( - 0.641_wp * t ) * deg2rad
dec = ( 90.0_wp - 0.557_wp * t ) * deg2rad
w = ( 190.147_wp + 360.9856235_wp * d ) * deg2rad

!it is a 3-1-3 rotation:
tmp = matmul( rotmat_x(pi2-dec), rotmat_z(pi2+ra) )
rotmat = matmul( rotmat_z(w), tmp )

end function icrf_to_iau_earth

!The 3x3 rotation matrix for a rotation about the x-axis.
pure function rotmat_x(angle) result(rotmat)

implicit none

real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]

real(wp) :: c,s

c = cos(angle)
s = sin(angle)

rotmat = reshape([one, zero, zero, &
                  zero, c, -s, &
                  zero, s, c],[3,3])

end function rotmat_x

!The 3x3 rotation matrix for a rotation about the z-axis.
pure function rotmat_z(angle) result(rotmat)

implicit none

real(wp),dimension(3,3) :: rotmat !rotation matrix
real(wp),intent(in) :: angle !angle [rad]

real(wp) :: c,s

c = cos(angle)
s = sin(angle)

rotmat = reshape([ c, -s, zero, &
                   s, c, zero, &
                   zero, zero, one ],[3,3])

end function rotmat_z

end module iau_orientation_module

References

1. B. A. Archinal, et al., "Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2009", Celest Mech Dyn Astr (2011) 109:101-135.
2. J. Williams, Fortran Astrodynamics Toolkit - iau_orientation_module [GitHub]

SOFA

The Standards of Fundamental Astronomy (SOFA) library is a very nice collection of routines that implement various IAU algorithms for fundamental astronomy computations. Versions are available in Fortran and C. Unfortunately, the Fortran version is written to be of maximum use to astronomers who frequently travel back in time to 1977 and compile the code on their UNIVAC (i.e., it is fixed-format Fortran 77). To assist 21st century users, I uploaded a little Python script called SofaMerge to Github. This script can make your life easier by merging all the individual source files into a single Fortran module file (with a few cosmetic updates along the way). It doesn't yet include any fixed to free format changes, but that could be easily added later.

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

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.