Runge-Kutta Nyström methods can be used to integrate second-order ordinary differential equations of the form:
- (a) \(\ddot{\mathbf{x}}=f(t,\mathbf{x})\)
There are also versions for the more general form:
- (b) \(\ddot{\mathbf{x}}=f(t,\mathbf{x},\dot{\mathbf{x}})\)
For second-order systems, Nyström methods are generally more efficient than the more familiar Runge-Kutta methods. The following algorithm is a 4th order Nyström method from Reference [1], encapsulated into a simple Fortran module. It integrates equations of type (b). Several others are also given in the same paper. A high-accuracy variable step size Nyström method for systems of type (a) can also be found in Reference [2].
module nystrom
use, intrinsic :: iso_fortran_env, wp=>real64
implicit none
private
! derivative function prototype
abstract interface
function deriv(n,t,x,xd) result(xdd)
import :: wp
integer,intent(in) :: n
real(wp),intent(in) :: t
real(wp),dimension(n),intent(in) :: x
real(wp),dimension(n),intent(in) :: xd
real(wp),dimension(n) :: xdd
end function deriv
end interface
public :: nystrom_lear_44
contains
subroutine nystrom_lear_44(n,df,t0,x0,xd0,dt,xf,xdf)
! 4th order Nystrom integration routine
! Take a step from t0 to t0+dt (x0,xd0 -> xf,xdf)
implicit none
integer,intent(in) :: n !number of equations
procedure(deriv) :: df !derivative routine
real(wp),intent(in) :: t0 !initial time
real(wp),dimension(n),intent(in) :: x0,xd0 !initial state
real(wp),intent(in) :: dt !time step
real(wp),dimension(n),intent(out) :: xf,xdf !final state
real(wp) :: t
real(wp),dimension(n) :: x,xd,k1,k2,k3,k4
real(wp),parameter :: s5 = sqrt(5.0_wp)
real(wp),parameter :: delta2 = (5.0_wp-s5)/10.0_wp
real(wp),parameter :: delta3 = (5.0_wp+s5)/10.0_wp
real(wp),parameter :: delta4 = 1.0_wp
real(wp),parameter :: a1 = (3.0_wp-s5)/20.0_wp
real(wp),parameter :: b1 = 0.0_wp
real(wp),parameter :: b2 = (3.0_wp+s5)/20.0_wp
real(wp),parameter :: c1 = (-1.0_wp+s5)/4.0_wp
real(wp),parameter :: c2 = 0.0_wp
real(wp),parameter :: c3 = (3.0_wp-s5)/4.0_wp
real(wp),parameter :: alpha1 = 1.0_wp/12.0_wp
real(wp),parameter :: alpha2 = (5.0_wp+s5)/24.0_wp
real(wp),parameter :: alpha3 = (5.0_wp-s5)/24.0_wp
real(wp),parameter :: alpha4 = 0.0_wp
real(wp),parameter :: ad1 = (5.0_wp-s5)/10.0_wp
real(wp),parameter :: bd1 = -(5.0_wp+3.0_wp*s5)/20.0_wp
real(wp),parameter :: bd2 = (3.0_wp+s5)/4.0_wp
real(wp),parameter :: cd1 = (-1.0_wp+5.0_wp*s5)/4.0_wp
real(wp),parameter :: cd2 = -(5.0_wp+3.0_wp*s5)/4.0_wp
real(wp),parameter :: cd3 = (5.0_wp-s5)/2.0_wp
real(wp),parameter :: beta1 = 1.0_wp/12.0_wp
real(wp),parameter :: beta2 = 5.0_wp/12.0_wp
real(wp),parameter :: beta3 = 5.0_wp/12.0_wp
real(wp),parameter :: beta4 = 1.0_wp/12.0_wp
if (dt==0.0_wp) then
xf = x0
xdf = xd0
else
t = t0
x = x0
xd = xd0
k1 = dt*df(n,t,x,xd)
t = t0+delta2*dt
x = x0+delta2*dt*xd0+a1*dt*k1
xd = xd0+ad1*k1
k2 = dt*df(n,t,x,xd)
t = t0+delta3*dt
x = x0+delta3*dt*xd0+dt*(b1*k1+b2*k2)
xd = xd0+bd1*k1+bd2*k2
k3 = dt*df(n,t,x,xd)
t = t0+dt
x = x0+dt*xd0+dt*(c1*k1+c3*k3)
xd = xd0+cd1*k1+cd2*k2+cd3*k3
k4 = dt*df(n,t,x,xd)
xf = x0+dt*xd0+dt*(alpha1*k1+alpha2*k2+alpha3*k3)
xdf = xd0+beta1*k1+beta2*k2+beta3*k3+beta4*k4
end if
end subroutine nystrom_lear_44
end module nystrom
References
- W. M. Lear, "Some New Nyström Integrators", 78-FM-11, JSC-13863, Mission Planning and Analysis Division, NASA Johnson Space Center, Feb. 1978.
- R. W. Brankin, I. Gladwell, J. R. Dormand, P. J. Prince, and W. L. Seward, "Algorithm 670: a Runge-Kutta-Nyström code", ACM Transactions on Mathematical Software, Volume 15 Issue 1, March 1989. [code]
For many years, I have used the closed-form solution from Heikkinen [1] for converting Cartesian coordinates to geodetic latitude and altitude. I've never actually seen the original reference (which is in German). I coded up the algorithm from the table given in [2], and have used it ever since, in school and also at work. At one point I compared it with various other closed-form algorithms, and never found one that was faster and still had the same level of accuracy. However, I recently happened to stumble upon Olson's method [3], which seems to be better.
A Fortran implementation of Olson's algorithm is:
pure subroutine olson(rvec, a, b, h, lon, lat)
use, intrinsic :: iso_fortran_env, wp=>real64 !double precision
implicit none
real(wp),dimension(3),intent(in) :: rvec !position vector [km]
real(wp),intent(in) :: a ! geoid semimajor axis [km]
real(wp),intent(in) :: b ! geoid semiminor axis [km]
real(wp),intent(out) :: h ! geodetic altitude [km]
real(wp),intent(out) :: lon ! longitude [rad]
real(wp),intent(out) :: lat ! geodetic latitude [rad]
real(wp) :: f,x,y,z,e2,a1,a2,a3,a4,a5,a6,w,zp,&
w2,r2,r,s2,c2,u,v,s,ss,c,g,rg,rf,m,p,z2
x = rvec(1)
y = rvec(2)
z = rvec(3)
f = (a-b)/a
e2 = f * (2.0_wp - f)
a1 = a * e2
a2 = a1 * a1
a3 = a1 * e2 / 2.0_wp
a4 = 2.5_wp * a2
a5 = a1 + a3
a6 = 1.0_wp - e2
zp = abs(z)
w2 = x*x + y*y
w = sqrt(w2)
z2 = z * z
r2 = z2 + w2
r = sqrt(r2)
if (r < 100.0_wp) then
lat = 0.0_wp
lon = 0.0_wp
h = -1.0e7_wp
else
s2 = z2 / r2
c2 = w2 / r2
u = a2 / r
v = a3 - a4 / r
if (c2 > 0.3_wp) then
s = (zp / r) * (1.0_wp + c2 * (a1 + u + s2 * v) / r)
lat = asin(s)
ss = s * s
c = sqrt(1.0_wp - ss)
else
c = (w / r) * (1.0_wp - s2 * (a5 - u - c2 * v) / r)
lat = acos(c)
ss = 1.0_wp - c * c
s = sqrt(ss)
end if
g = 1.0_wp - e2 * ss
rg = a / sqrt(g)
rf = a6 * rg
u = w - rg * c
v = zp - rf * s
f = c * u + s * v
m = c * v - s * u
p = m / (rf / g + f)
lat = lat + p
if (z < 0.0_wp) lat = -lat
h = f + m * p / 2.0_wp
lon = atan2( y, x )
end if
end subroutine olson
For comparison, Heikkinen's algorithm is:
pure subroutine heikkinen(rvec, a, b, h, lon, lat)
use, intrinsic :: iso_fortran_env, wp=>real64 !double precision
implicit none
real(wp),dimension(3),intent(in) :: rvec !position vector [km]
real(wp),intent(in) :: a ! geoid semimajor axis [km]
real(wp),intent(in) :: b ! geoid semiminor axis [km]
real(wp),intent(out) :: h ! geodetic altitude [km]
real(wp),intent(out) :: lon ! longitude [rad]
real(wp),intent(out) :: lat ! geodetic latitude [rad]
real(wp) :: f,e_2,ep,r,e2,ff,g,c,s,pp,q,r0,u,v,z0,x,y,z,z2,r2,tmp,a2,b2
x = rvec(1)
y = rvec(2)
z = rvec(3)
a2 = a*a
b2 = b*b
f = (a-b)/a
e_2 = (2.0_wp*f-f*f)
ep = sqrt(a2/b2 - 1.0_wp)
z2 = z*z
r = sqrt(x**2 + y**2)
r2 = r*r
e2 = a2 - b2
ff = 54.0_wp * b2 * z2
g = r2 + (1.0_wp - e_2)*z2 - e_2*e2
c = e_2**2 * ff * r2 / g**3
s = (1.0_wp + c + sqrt(c**2 + 2.0_wp*c))**(1.0_wp/3.0_wp)
pp = ff / ( 3.0_wp*(s + 1.0_wp/s + 1.0_wp)**2 * g**2 )
q = sqrt( 1.0_wp + 2.0_wp*e_2**2 * pp )
r0 = -pp*e_2*r/(1.0_wp+q) + &
sqrt( max(0.0_wp, 1.0_wp/2.0_wp * a2 * (1.0_wp + 1.0_wp/q) - &
( pp*(1.0_wp-e_2)*z2 )/(q*(1.0_wp+q)) - &
1.0_wp/2.0_wp * pp * r2) )
u = sqrt( (r - e_2*r0)**2 + z2 )
v = sqrt( (r - e_2*r0)**2 + (1.0_wp - e_2)*z2 )
z0 = b**2 * z / (a*v)
h = u*(1.0_wp - b2/(a*v) )
lat = atan2( (z + ep**2*z0), r )
lon = atan2( y, x )
end subroutine heikkinen
Olson's is about 1.5 times faster on my 2.53 GHz i5 laptop, when both are compiled using gfortran with -O2 optimization (Olson: 3,743,212 cases/sec, Heikkinen: 2,499,670 cases/sec). The level of accuracy is about the same for each method.
References
- M. Heikkinen, "Geschlossene formeln zur berechnung raumlicher geodatischer koordinaten aus rechtwinkligen Koordinaten". Z. Ermess., 107 (1982), 207-211 (in German).
- E. D. Kaplan, "Understanding GPS: Principles and Applications", Artech House, 1996.
- D. K. Olson, "Converting Earth-Centered, Earth-Fixed Coordinates to Geodetic Coordinates," IEEE Transactions on Aerospace and Electronic Systems, 32 (1996) 473-476.
Introducing json-fortran, an easy-to-use JSON API written in modern Fortran. As far as I know, it's currently the only publicly-available production-ready JSON API for Fortran that does not involve an interface to C libraries. The code is hosted at GitHub, and is released under a BSD-style license.
Fortran users may find JSON quite useful as a configuration file format. Typically, large and complicated Fortran codes (say, in the fields of climate modeling or trajectory optimization) require data to be read from configuration files. JSON has many advantages over a roll-your-own file format, or a Fortran namelist, namely:
- It's a standard.
- It's human-readable and human-editable.
- API's exist for many other programming languages.
Consider the example of a program to propagate a spacecraft state vector. The required information is the initial time t0, step size dt, final time tf, central body gravitational parameter mu, and the 6-element state vector x0. The JSON configuration file might look something like this:
{
"t0": 0.0,
"dt": 1.0,
"tf": 86400.0,
"mu": 398600.4418,
"x0": [ 10000.0,
10000.0,
10000.0,
1.0,
2.0,
3.0
]
}
The code would look something like this:
program propagate
use json_module
implicit none
real(wp) :: t0, dt, tf, mu
real(wp),dimension(:),allocatable :: x0
type(json_file) :: config
logical :: found
!load the file:
call config%load_file('config.json')
!read in the data:
call config%get('t0',t0,found)
call config%get('dt',dt,found)
call config%get('tf',tf,found)
call config%get('mu',mu,found)
call config%get('x0',x0,found)
!propagate:
! ...
end program propagate
Of course, real production code would have more graceful error checking. For example, to make sure the file was properly parsed and that each variable was really present. These features are included in the API, including helpful error messages if something goes wrong.