Degenerate Conic

Julian Date

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

1. "Converting Between Julian Dates and Gregorian Calendar Dates", United States Naval Observatory.
2. D. Steel, "Marking Time: The Epic Quest to Invent the Perfect Calendar", John Wiles & Sons, 2000.

Merry Christmas

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  

Fortran Build Tools

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

• Fortran Build Tools - at the Fortran Wiki

Speeding up SPICE

The SPICE Toolkit software is an excellent package of very well-written and well-documented routines for a variety of astrodynamics applications. It is produced by NASA's Navigation and Ancillary Information Facility (NAIF). Versions are available for Fortran 77, C, IDL, and Matlab.

To speed up the execution of SPICE-based programs, there are a few things you can do:

• For the Fortran SPICELIB, recompile it with optimization enabled (say, -O2). The default library released by NAIF is not compiled with optimization.
• Turn off the SPICE traceback system (call TRCOFF()). If you are using a compiler that has a built-in stack trace routine (for example TRACEBACKQQ in the Intel compiler), just include a call to it in the SPICE BYEBYE.F routine. That will give you a stack trace for any fatal errors.
• Converting a non-native binary PCK to native form will also speed up data access somewhat.
• Calls to ephemeris routines where the target and observer bodies are input as strings will be slightly slower than the ones where the inputs are the NAIF ID codes (which are integers). For example, for the geometric state (position and velocity) of one body relative to another, use SPKGEO instead of SPKEZR. Also, if all you require is position and not velocity, use SPKGPS instead of SPKGEO, since it will be a bit faster.

Also note that for applications only requiring the ephemerides of the solar system major bodies, there is an older code from JPL (PLEPH) which is simpler and faster than SPICE, but uses a different format for the ephemeris files.

Rocket Equation

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.

Porkchop Plots

A porkchop plot is a data visualization tool used in interplanetary mission design which displays contours of various quantities as a function of departure and arrival date. Example pork chop plots for 2016 Earth-Mars transfers are shown here. The x-axis is the Earth departure date, and the y-axis is the Mars arrival date. The contours show the Earth departure C3, the Mars arrival C3, and the sum of the two. C3 is the "characteristic energy" of the departure or arrival, and is equal to the \(v_{\infty}^2\) of the hyperbolic trajectory. Each point on the curve represents a ballistic Earth-Mars transfer (computed using a Lambert solver). The two distinct black and blue contours on each plot represent the "short way" (\(<\pi\)) and "long way" (\(>\pi\)) transfers. The contours look vaguely like pork chops, the centers being the lowest C3 value, and thus optimal transfers. The January - October 2016 opportunity will be used for the ExoMars Trace Gas Orbiter, and the March - September 2016 opportunity will be used for InSight.

References

• A. B. Sergeyevsky, G. C. Snyder, R. A. Cunniff, "Interplanetary Mission Design Handbook, Volume I, Part 2: Earth to Mars Ballistic Mission Opportunities, 1990-2005", NASA JPL, September 1983.
• L. M. Burke, R. D. Falck, and M. L. McGuire, "Interplanetary Mission Design Handbook: Earth-to-Mars Mission Opportunities 2026 to 2045", NASA/TM—2010-216764, October 2010.
• Chauncey Uphoff, The History of the Term C3, 2001.

Complex Step Differentiation

Complex step differentiation is a method for estimating the derivative of a function \(f(x)\) using the equation:

• \(f'(x) \approx \frac{ Im[f(x + ih)] }{h}\)

Unlike finite-difference formulas, the complex-step formula does not suffer from roundoff errors due to subtraction, so a very small value of the step size \(h\) can be used, producing a much more accurate derivative. The example shown below is for the function: \(f(x) = e^x + \sin(x)\). The derivative error is compared to three finite-difference formulas:

• \(f'(x) \approx \frac{f(x+h) - f(x)}{h}\) (Two-point forward difference)
• \(f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}\) (Two-point central difference)
• \(f'(x) \approx \frac{f(x-2h) -8 f(x-h) + 8 f(x+h) - f(x+2h) }{12 h}\) (Four-point central difference)

For large values of \(h\), truncation error dominates. The finite-difference error decreases as \(h\) is decreased, until roundoff error begins to dominate (around \(10^{-8}\) for the forward difference, \(10^{-5}\) for the 2-point central difference, and \(10^{-3}\) for the 4-point central difference). The complex-step estimate, however, can produce a derivative estimate to machine precision for any value of \(h < 10^{-8}\).

See also

1. J.R.R.A. Martins, I.M. Kroo, J.J. Alonso, "An Automated Method for Sensitivity Analysis using Complex Variables", AIAA-2000-0689.
2. complex_step.f90 in the Fortran Astrodynamics Toolkit.

Lambert's Problem

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

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. D. Izzo, "Revisiting Lambert’s problem", Celestial Mechanics and Dynamical Astronomy, Oct. 2014.

Fortran & C Interoperability

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

1. C interop to popen, comp.lang.fortran, 12/2/2009.

Midpoint Circle Algorithm

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

1. Jack E. Bresenham, "Algorithms for Computer Control of a Digital Plotter", IBM System Journal, 1965.