The 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)
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.
Julian date is a count of the number of days since noon on January 1, 4713 BC in the proleptic Julian calendar. This epoch was chosen by Joseph Scaliger in 1583 as the start of the "Julian Period": a 7,980 year period that is the multiple of the 19-year Metonic cycle, the 28-year solar/dominical cycle, and the 15-year indiction cycle. It is a conveniently-located epoch since it precedes all written history. A simple Fortran function for computing Julian date given the Gregorian calendar year, month, day, and time is:
function julian_date(y,m,d,hour,minute,sec)
use, intrinsic :: iso_fortran_env, only: wp => real64
implicit none
real(wp) :: julian_date
integer,intent(in) :: y,m,d ! Gregorian year, month, day
integer,intent(in) :: hour,minute,sec ! Time of day
integer :: julian_day_num
julian_day_num = d-32075+1461*(y+4800+(m-14)/12)/4+367*&
(m-2-(m-14)/12*12)/12-3*((y+4900+(m-14)/12)/100)/4
julian_date = real(julian_day_num,wp) + &
(hour-12.0_wp)/24.0_wp + &
minute/1440.0_wp + &
sec/86400.0_wp
end function julian_date
References
- "Converting Between Julian Dates and Gregorian Calendar Dates", United States Naval Observatory.
- D. Steel, "Marking Time: The Epic Quest to Invent the Perfect Calendar", John Wiles & Sons, 2000.
program main
implicit none
integer :: i,j,nstar,nspace,ir
character(len=:),allocatable :: stars,spaces
integer,parameter :: n = 10
integer,parameter :: total = 1 + (n-1)*2
do j=1,200
call system('clear')
nstar = -1
do i=1,n
nstar = nstar + 2
nspace = (total - nstar)/2
stars = repeat('*',nstar)
spaces = repeat(' ',nspace)
if (i>1) then
ir = max(1,ceiling(rand(0)*len(stars)))
stars(ir:ir) = ' '
end if
write(*,'(A)') spaces//stars//spaces
end do
spaces = repeat(' ',(total-1)/2)
write(*,'(A)') spaces//'*'//spaces
write(*,'(A)') ''
end do
end program main
Let's face it, make is terrible. It is especially terrible for large modern Fortran projects, which can have complex source file interdependencies due to the use of modules. To use make with modern Fortran, you need an additional tool to generate the proper file dependency. Such tools apparantly exist (for example, makemake, fmkmf, sfmakedepend, and Makedepf90), but I've never used any of them. Any Fortran build solution that involves make is a nonstarter for me.
If you are an Intel Fortran user on Windows, the Visual Studio integration automatically determines the correct compilation order for you, and you never have to think about it (this is the ideal solution). However if you are stuck using gfortran, there are still various decent opensource solutions for building modern Fortran projects that you can use:
- SCons - A Software Construction Tool. I used SCons for a while several years ago, and it mostly worked, but I found it non-trivial to configure, and the Fortran support was flaky. Eventually, I stopped using it. Newer releases may have improved, but I don't know.
- foraytool [Drew McCormack] (formerly called TCBuild) - This one was specifically designed for Fortran, works quite well and is easy to configure. However, it does not appear to be actively maintained (the last release was over four years ago).
- FoBiS [Stefano Zaghi] - Fortran Building System for Poor Men. This is quite new (2014), and was also specifically designed for Fortran. The author refers to it as "a very simple and stupid tool for automatically building modern Fortran projects". It is trivially easy to use, and is also quite powerful. This is probably the best one to try first, especially if you don't want to have to think about anything.
With FoBiS, if your source is in ./src, and you want to build the application at ./bin/myapp, all you have to do is this: FoBiS.py build -s ./src -compiler gnu -o ./bin/myapp. There are various other command line flags for more complicated builds, and a configuration file can also be used.
See also
The rocket equation describes the basic principles of a rocket in the absence of external forces. Various forms of this equation relate the following fundamental parameters:
- the engine thrust (\(T\)),
- the engine specific impulse (\(I_{sp}\)),
- the engine exhaust velocity (\(c\)),
- the initial mass (\(m_i\)),
- the final mass (\(m_f\)),
- the burn duration (\(\Delta t\)), and
- the effective delta-v (\(\Delta v\)) of the maneuver.
The rocket equation is:
- \(\Delta v = c \ln \left( \frac{m_i}{m_f} \right)\)
The engine specific impulse is related to the engine exhaust velocity by: \(c = I_{sp} g_0\), where \(g_0\) is the standard Earth gravitational acceleration at sea level (defined to be exactly 9.80665 \(m/s^2\)). The thrust is related to the mass flow rate (\(\dot{m}\)) of the engine by: \(T = c \dot{m}\). This can be used to compute the duration of the maneuver as a function of the \(\Delta v\):
- \(\Delta t = \frac{c m_i}{T} \left( 1 - e^{-\frac{\Delta v}{c}} \right)\)
A Fortran module implementing these equations is given below.
module rocket_equation
use,intrinsic :: iso_fortran_env, only: wp => real64
implicit none
!two ways to compute effective delta-v
interface effective_delta_v
procedure :: effective_dv_1,effective_dv_2
end interface effective_delta_v
real(wp),parameter,public :: g0 = 9.80665_wp ![m/s^2]
contains
pure function effective_dv_1(c,mi,mf) result(dv)
real(wp) :: dv ! effective delta v [m/s]
real(wp),intent(in) :: c ! exhaust velocity [m/s]
real(wp),intent(in) :: mi ! initial mass [kg]
real(wp),intent(in) :: mf ! final mass [kg]
dv = c*log(mi/mf)
end function effective_dv_1
pure function effective_dv_2(c,mi,T,dt) result(dv)
real(wp) :: dv ! delta-v [m/s]
real(wp),intent(in) :: c ! exhaust velocity [m/s]
real(wp),intent(in) :: mi ! initial mass [kg]
real(wp),intent(in) :: T ! thrust [N]
real(wp),intent(in) :: dt ! duration of burn [sec]
dv = -c*log(1.0_wp-(T*dt)/(c*mi))
end function effective_dv_2
pure function burn_duration(c,mi,T,dv) result(dt)
real(wp) :: dt ! burn duration [sec]
real(wp),intent(in) :: c ! exhaust velocity [m/s]
real(wp),intent(in) :: mi ! initial mass [kg]
real(wp),intent(in) :: T ! engine thrust [N]
real(wp),intent(in) :: dv ! delta-v [m/s]
dt = (c*mi)/T*(1.0_wp-exp(-dv/c))
end function burn_duration
pure function final_mass(c,mi,dv) result(mf)
real(wp) :: mf ! final mass [kg]
real(wp),intent(in) :: c ! exhaust velocity [m/s]
real(wp),intent(in) :: mi ! initial mass [kg]
real(wp),intent(in) :: dv ! delta-v [m/s]
mf = mi/exp(dv/c)
end function final_mass
end module rocket_equation
References
- J.E. Prussing, B.A. Conway, "Orbital Mechanics", Oxford University Press, 1993.
Lambert's problem is to solve for the orbit transfer that connects two position vectors in a given time of flight. It is one of the basic problems of orbital mechanics, and was solved by Swiss mathematician Johann Heinrich Lambert.
A standard Fortran 77 implementation of Lambert's problem was presented by Gooding [1]. A modern update to this implementation is included in the Fortran Astrodynamics Toolkit, which also includes a newer algorithm from Izzo [2].
There can be multiple solutions to Lambert's problem, which are classified by:
- Direction of travel (i.e., a "short" way or "long" way transfer).
- Number of complete revolutions (N=0,1,...). For longer flight times, solutions exist for larger values of N.
- Solution number (S=1 or S=2). The multi-rev cases can have two solutions each.
In the example shown here, there are 6 possible solutions:
References
- 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.
- D. Izzo, "Revisiting Lambert’s problem", Celestial Mechanics and Dynamical Astronomy, Oct. 2014.
The ISO_C_BINDING intrinsic module and the BIND attribute introduced in Fortran 2003 are very handy for producing standard and portable Fortran code that interacts with C code. The example given here shows how to use the popen, fgets, and pclose routines from the C standard library to pipe the result of a shell command into a Fortran allocatable string.
module pipes_module
use,intrinsic :: iso_c_binding
implicit none
private
interface
function popen(command, mode) bind(C,name='popen')
import :: c_char, c_ptr
character(kind=c_char),dimension(*) :: command
character(kind=c_char),dimension(*) :: mode
type(c_ptr) :: popen
end function popen
function fgets(s, siz, stream) bind(C,name='fgets')
import :: c_char, c_ptr, c_int
type (c_ptr) :: fgets
character(kind=c_char),dimension(*) :: s
integer(kind=c_int),value :: siz
type(c_ptr),value :: stream
end function fgets
function pclose(stream) bind(C,name='pclose')
import :: c_ptr, c_int
integer(c_int) :: pclose
type(c_ptr),value :: stream
end function pclose
end interface
public :: c2f_string, get_command_as_string
contains
!**********************************************
! convert a C string to a Fortran string
!**********************************************
function c2f_string(c) result(f)
implicit none
character(len=*),intent(in) :: c
character(len=:),allocatable :: f
integer :: i
i = index(c,c_null_char)
if (i<=0) then
f = c
else if (i==1) then
f = ''
else if (i>1) then
f = c(1:i-1)
end if
end function c2f_string
!**********************************************
! return the result of the command as a string
!**********************************************
function get_command_as_string(command) result(str)
implicit none
character(len=*),intent(in) :: command
character(len=:),allocatable :: str
integer,parameter :: buffer_length = 1000
type(c_ptr) :: h
integer(c_int) :: istat
character(kind=c_char,len=buffer_length) :: line
str = ''
h = c_null_ptr
h = popen(command//c_null_char,'r'//c_null_char)
if (c_associated(h)) then
do while (c_associated(fgets(line,buffer_length,h)))
str = str//c2f_string(line)
end do
istat = pclose(h)
end if
end function get_command_as_string
end module pipes_module
A example use of this module is:
program test
use pipes_module
implicit none
character(len=:),allocatable :: res
res = get_command_as_string('uname')
write(*,'(A)') res
res = get_command_as_string('ls -l')
write(*,'(A)') res
end program test
References
- C interop to popen, comp.lang.fortran, 12/2/2009.
The Midpoint circle algorithm is a clever and efficient way of drawing a circle using only addition, subtraction, and bit shifts. It is based on the Bresenham line algorithm developed by Jack Bresenham in 1962 at IBM.
The algorithm was also independently discovered by Apple programmer Bill Atkinson in 1981 when developing QuickDraw for the original Macintosh.
A Fortran implementation is given below (which was used to draw the circle shown here, which has a radius of 7 pixels):
subroutine draw_circle(x0, y0, radius, color)
implicit none
integer,intent(in) :: x0, y0, radius, color
integer :: x,y,err
x = radius
y = 0
err = 1-x
do while (x >= y)
call color_pixel( x + x0, y + y0, color)
call color_pixel( y + x0, x + y0, color)
call color_pixel( -x + x0, y + y0, color)
call color_pixel( -y + x0, x + y0, color)
call color_pixel( -x + x0, -y + y0, color)
call color_pixel( -y + x0, -x + y0, color)
call color_pixel( x + x0, -y + y0, color)
call color_pixel( y + x0, -x + y0, color)
y = y + 1
if (err<0) then
err = err + 2 * y + 1
else
x = x - 1
err = err + 2 * (y - x + 1)
end if
end do
end subroutine draw_circle
References
- Jack E. Bresenham, "Algorithms for Computer Control of a Digital Plotter", IBM System Journal, 1965.
Fortran90.org has a nice page showing side-by-side comparisons of the same code in both NumPy and Fortran. Aside from a few weird quirks (see the array slice examples), Python code should be easily understandable by a Fortran programmer.
Fortran 2003 was a significant update to Fortran 95 (probably something on the order of the update from C to C++). It brought Fortran into the modern computing world with the addition of object oriented programming features such as type extension and inheritance, polymorphism, dynamic type allocation, and type-bound procedures.
The following example shows an implementation of an object oriented class to perform numerical integration of a user-defined function using a Runge-Kutta method (in this case RK4). First, an abstract rk_class is defined, which is then extended to create the RK4 class by associating the deferred step method to rk4. Other RK classes could also be created in a similar manner. In order to use the class, it is only necessary to set n (the number of variables in the state vector) and f (the user-supplied derivative function), and then call the integrate method.
module rk_module
use, intrinsic :: iso_fortran_env, wp => real64 !double precision reals
implicit none
real(wp),parameter :: zero = 0.0_wp
!main integration class:
type,abstract,public :: rk_class
private
!user specified number of variables:
integer :: n = 0
!user-specified derivative function:
procedure(deriv_func),pointer :: f => null()
contains
procedure,non_overridable,public :: integrate !main integration routine
procedure(step_func),deferred :: step !the step routine for the selected method
end type rk_class
!extend the abstract class to create an RK4 method:
! [all we need to do is set the step function]
type,extends(rk_class),public :: rk4_class
contains
procedure :: step => rk4
end type rk4_class
interface
subroutine deriv_func(me,t,x,xdot) !derivative function
import :: rk_class,wp
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),dimension(me%n),intent(out) :: xdot
end subroutine deriv_func
subroutine step_func(me,t,x,h,xf) !rk step function
import :: rk_class,wp
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),intent(in) :: h
real(wp),dimension(me%n),intent(out) :: xf
end subroutine step_func
end interface
contains
subroutine integrate(me,t0,x0,h,tf,xf)
! main integration routine
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t0
real(wp),dimension(me%n),intent(in) :: x0
real(wp),intent(in) :: h
real(wp),intent(in) :: tf
real(wp),dimension(me%n),intent(out) :: xf
real(wp) :: t,dt,t2
real(wp),dimension(me%n) :: x
logical :: last
if (h==zero) then
xf = x0
else
t = t0
x = x0
dt = h
do
t2 = t + dt
last = ((dt>=zero .and. t2>=tf) .or. & !adjust last time step
(dt<zero .and. t2<=tf)) !
if (last) dt = tf-t !
call me%step(t,x,dt,xf)
if (last) exit
x = xf
t = t2
end do
end if
end subroutine integrate
subroutine rk4(me,t,x,h,xf)
! Take one Runge Kutta 4 integration step: t -> t+h (x -> xf)
implicit none
class(rk4_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),intent(in) :: h
real(wp),dimension(me%n),intent(out) :: xf
!local variables:
real(wp),dimension(me%n) :: f1,f2,f3,f4
real(wp) :: h2
!parameters:
real(wp),parameter :: half = 0.5_wp
real(wp),parameter :: six = 6.0_wp
h2 = half*h
call me%f(t,x,f1)
call me%f(t+h2,x+h2*f1,f2)
call me%f(t+h2,x+h2*f2,f3)
call me%f(t+h,x+h*f3,f4)
xf = x + h*(f1+f2+f2+f3+f3+f4)/six
end subroutine rk4
end module rk_module
An example use of this module follows. Here, we are integrating the equations of motion of a spacecraft in Earth orbit, using two-body assumptions: \(\ddot{\mathbf{r}} = -\frac{\mu}{r^3}\mathbf{r}\)
The rk4_class is extended to create a new spacecraft class, which contains the data required by the derivative routine (in this case, simply the central body gravitational parameter mu). A more complicated version could also include spherical harmonic gravity, solar radiation pressure, atmospheric drag, engine thrust, etc. Note that the select type statement is necessary in the derivative routine in order to access the data from the extended type.
program main
use rk_module
!spacecraft propagation type:
! extend the rk class to include data used in the deriv routine
type,extends(rk4_class) :: spacecraft
real(wp) :: mu = zero ! central body gravitational
! parameter (km^3/s^2)
end type spacecraft
integer,parameter :: n=6 !number of state variables
type(spacecraft) :: s
real(wp) :: t0,tf,x0(n),dt,xf(n)
!constructor:
s = spacecraft(n=n,f=twobody,mu=398600.436233_wp) !main body is Earth
!initial conditions:
x0 = [10000.0_wp,10000.0_wp,10000.0_wp,& !initial state [x,y,z] (km)
1.0_wp,2.0_wp,3.0_wp] ! [vx,vy,yz] (km/s)
t0 = zero !initial time (sec)
dt = 10.0_wp !time step (sec)
tf = 1000.0_wp !final time (sec)
call s%integrate(t0,x0,dt,tf,xf)
write(*,'(A/,*(F15.6/))') 'Final state:',xf
contains
subroutine twobody(me,t,x,xdot)
! derivative routine for two-body orbit propagation
implicit none
class(rk_class),intent(inout) :: me
real(wp),intent(in) :: t
real(wp),dimension(me%n),intent(in) :: x
real(wp),dimension(me%n),intent(out) :: xdot
real(wp),dimension(3) :: r,v,a_grav
real(wp) :: rmag
select type (me)
class is (spacecraft)
r = x(1:3)
v = x(4:6)
rmag = norm2(r)
a_grav = -me%mu/rmag**3 * r !accel. due to gravity
xdot(1:3) = v
xdot(4:6) = a_grav
end select
end subroutine twobody
end program main
The result is:
Final state: 10667.963305 11658.055962 12648.148619 0.377639 1.350074 2.322509
See also
- Object-Oriented Programming in Fortran 2003 Part 1: Code Reusability, PGI Insider, Feb. 2011.
- Object-Oriented Programming in Fortran 2003 Part 2: Data Polymorphism, PGI Insider, June 2011.
- Fortran Astrodynamics Toolkit (Github)
