Degenerate Conic

The New Features of Fortran 202x

The glacially slow pace of Fortran language development continues! From the smoke-filled room of the Fortran Standards Committee, the next Fortran standard is coming along. It has not yet been given an official name but is currently referred to as "Fortran 202x" (so hopefully we will see it before the end of this decade). Recall that Fortran 90 was originally referred to as "Fortran 198x".

A document (from John Reid) describing the new features is available here and is summarized below:

• The limit on line length has been increased to ten thousand characters, and the limit on total statement length (with all continuation lines) has been increased to a million characters.
• Allocatable-length character variables were added in Fortran 2003, but were ignored by all intrinsic procedures, until now! Better late than never, I suppose. In Fortran 202x, intrinsic routines that return strings can now be passed an allocatable string and it will be auto-allocated to the correct size. This is actually a pretty great addition, since it fixes an annoyance that we've had to live with for 20 years. Of course, what we really need is a real intrinsic string class. Maybe in 20 more years we can have that.
• typeof and classof. No idea what these are for. Probably a half-baked generics feature? What we really need is a real generics feature. Since that is in work for Fortran 202y I don't know why these have been added now. Probably we will discover in 10 years that these aren't good for much, but will be stuck with them until the end of time (see parameterized derived types).
• Conditional expressions. This adds the dreadful C-style (a ? b : c) conditional expression to Fortran. Somebody thought this was a good idea? This abomination can also be passed to an optional procedure argument, with a new .nil. construct to signify not present.
• Improvements to the half-baked binary, octal, and hexadecimal constants.
• Two string parsing routines (split and tokenize) were added. These are basically Fortran 90 style string routines that I guess are OK, but why bother? Or why limit to just these two? See above comment about adding a real string class.
• Trig functions that accept input in degrees have been added (e.g. sind, cosd, etc). These are already commonly available as extensions in many compilers, so now they will just be standardized.
• Half revolution trig function have also been added (acospi, asinpi, etc.). I have no idea what these are good for.
• A selected_logical_kind function was added. This was a missing function that goes along with the selected_real_kind and selected_int_kind functions that were added in Fortran 2003. Logical type parameters (e.g., logical8, logical16, etc. were also added to the intrinsic iso_fortran_env module).
• Some updates to the system_clock intrinsic.
• Some new functions in the ieee_arithmetic module for conformance with the 2020 IEEE standard.
• Update to the c_f_pointer intrinsic procedure to allow for specifying the lower bound.
• New procedures were added to iso_c_binding for conversion between Fortran and C strings. This is nice, since everyone who ever used this feature had to write these themselves. So, now they will be built in.
• A new AT format code, which is like A, but trims the trailing space from a string. Another nice little change that will save a lot of trim() functions.
• Opening a file now includes a new option for setting the leading zeros in real numbers (print them or suppress them).
• Public namelists can now include private variables. Sure, why not?
• Some updates for coarrays (coarrays are the built-in MPI like component of Fortran).
• A new simple procedure specification. A simple procedure is sort of a super pure procedure that doesn't reference any variables outside of the procedure. Presumably, this would allow some additional optimizations from the compiler.
• A new "multiple subscript" array feature. The syntax for this is just dreadful (using the @ character in Fortran for the first time).
• Integer arrays can now be used to specify the rank and bounds of an array. I tried to read the description of this, but don't get it. Maybe it's useful for something.
• Can now specify the rank of a variable using an integer constant. This is nice, and makes the rank parts of the language a little more consistent with the other parts. Again, this is something that is cleaning up a half-baked feature from a previous standard.
• Reduction specifier for do concurrent, which is putting more OpenMP-like control into the actual language (but just enough so people will still complain that is isn't enough).
• enumeration types have been added for some reason. In addition, the half-baked Fortran 2003 enum has been extended so it can now be used to create a new type. I'm not entirely sure why we now have two different ways to do this.

So, there you have it. There are some nice additions that current users of Fortran will appreciate. It will likely be years before compilers implement all of this though. Fortran 202x will be the 3rd minor update of Fortran since Fortran 2003 was released 20 years ago. Unfortunately, there really isn't much here that is going to convince anybody new to start using Fortran (or to switch from Python or Julia). What we need are major additions to bring Fortran up to speed where it has fallen behind (e.g., in such areas as generics, exception handling, and string manipulation). Or how about some radical upgrades like automatic differentiation or built-in units for variables?

References

• J. Reid, The new features of Fortran 202x, ISO/IEC JTC1/SC22/WG5 N2194, March 21, 2022
• JTC1/SC22/WG5
• INCITS PL22.3, the US Fortran Programming Language Standards Technical Committee
• The new features of Fortran 2015 (degenerateconic.com)

Packs

There are a lot of classic Fortran libraries from the 1970s to the early 1990s whose names end in "pack" (or, more specifically, "PACK"). I presume this stands for "package", and it seems to have been a common nomenclature for software libraries back in the day. Here is a list:

Note that some of these dates are my best guesses for the original release date. I can update the table if I get better information.

LAPACK is the only one of these "packs" that has seen continuous development since it began. But the majority of these codes were abandoned by their original developers and organizations decades ago, and they have been frozen ever since, a symptom of the overall decline in the Fortran ecosystem. Though in a state of benign neglect, they yet remain useful to the present day. For example, SciPy includes wrappers to MINPACK and QUADPACK (and lots of other FORTRAN 77 code). However, the only "upstream" for many of these packages is Netlib, a sourcecode graveyard from another time that has no collaborative component. So, mostly any updates anybody makes are never distributed to anyone else. Frequently, these codes are converted to other languages (e.g, by f2c) and fixes or updates are made in the C code that never get back into the Fortran code. I've mentioned before how Fortran users end up calling Fortran wrappers to C code that was originally converted from FORTRAN 77 code that was never modernized.

The problem is that these classic codes, while great, are not perfect. They are written in the obsolete FORTRAN 77 fixed-form style, which nobody wants anything to do with nowadays, but that continues to poison the Fortran ecosystem. They are littered with now-unnecessary and confusing spaghetti code constructs such as GOTOs and arithmetic IF statements. They are not easy to incorporate into modern codes (there is no package manager for Netlib). Development could have continued up to the present day, and each of these libraries could have state of the art, modern Fortran implementations of both the classic and the latest algorithms. Well it's not too late. We have the internet now, and ways to collaborate on code (e.g, GitHub). We can restart development of some of these libraries:

• Convert them to free-form source
• Clean up the formatting
• Remove obsolete and now-unnecessary clutter
• Clean up the doc strings so we can use FORD to auto-generate online documentation
• Put them in a module so there is no chance of calling them with the wrong input types
• Update them for multiple real kinds
• Add automated unit testing
• Add C APIs so that they can be called from C (and, more importantly, any language that is interoperable with C, such as Python)
• Start adding newer algorithms developed in recent decades

Consider QUADPACK. This library has been used for decades, and translated into numerous other programming languages. The FORTRAN 77 version is still used in SciPy. Note that recently some minor bugs where found in the decades-old code. The edits actually made it back to Netlib through some opaque process, but the comment was made "I do not expect a lots of edits on the package. It feels like a bug here and there." Unfortunately, this has been the perception of these classic libraries: the old Fortran code is used because it works but nobody really wants to maintain it. Well, I've started a complete refactor of QUADPACK, and have already added new features as well as new methods not present in the original library. Also consider FFTPACK and MINPACK, which were recently brought under the fortran-lang group on GitHub and development restarted as community projects. Our (ambitious) goal is to try to get SciPy to use these modernized libraries instead of the old ones. Other libraries can also move in this direction. Note: we don't have forever. SciPy recently replaced FFTPACK with a new library rewritten in C++. It's inevitable that all the old FORTRAN 77 code is going to be replaced. No one wants to work with it, and for good reason.

It's time for the Fortran community to show examples of how modern Fortran upgrades to these libraries are a better option than rewriting them from scratch in another language. With a little spit and polish, these classic libraries can have decades more life left in them. There's no need for our entire scientific programming heritage to be rewritten every few years in yet another new programming language (C, C++, Matlab, Octave, R, Python, Julia, Rust, Go, Dart, Chapel, Zig, etc, etc, etc). These libraries can be used as-is in other languages such as Python, but modernizing also makes them easier to use in actual Fortran applications. We already have the new versions available for use via the Fortran Package Manager (FPM), so adding them to your Fortran project is now a simple addition to your fpm.toml file:

[dependencies]
minpack = { git="https://github.com/fortran-lang/minpack.git" }
quadpack = { git="https://github.com/jacobwilliams/quadpack.git" }

Now isn't that better than manually downloading unversioned source files from an ftp server set up before many of the people reading this were born?

Update: See also this thread at the fortran-lang Discourse.

References

• Jack Dongarra, Gene Golub, Eric Grosse, Cleve Moler, Keith Moore, "Netlib and NA-Net: building a scientific computing community", IEEE Annals of the History of Computing, 30 (2): 30-41, 2008
• Cleve Moler, A Brief History of MATLAB, MathWorks.
• J. J. Dongarra, C. B. Moler, EISPACK: A package for solving matrix eigenvalue problems.
• J. J. Moré, B. S. Garbow, and K. E. Hillstrom, User Guide for MINPACK-1, Argonne National Laboratory Report ANL-80-74, Argonne, Ill., 1980.
• Robert Piessens, Elise de Doncker-Kapenga, Christoph W. Überhuber, David Kahaner, "QUADPACK: A subroutine package for automatic integration", Springer-Verlag, ISBN 978-3-540-12553-2, 1983

Modified Equinoctial Elements

Last time, I mentioned an alternative orbital element set different from the classical elements. Modified equinoctial elements \((p,f,g,h,k,L)\) are another alternative that are applicable to all orbits and have nonsingular equations of motion (except for a singularity at \(i = \pi\)). They are defined as [1]:

• \(p = a (1 - e^2)\)
• \(f = e \cos (\omega + \Omega)\)
• \(g = e \sin (\omega + \Omega)\)
• \(h = \tan (i/2) \cos \Omega\)
• \(k = \tan (i/2) \sin \Omega\)
• \(L = \Omega + \omega + \nu\)

Where \(a\) is the semimajor axis, \(p\) is the semi-latus rectum, \(i\) is the inclination, \(e\) is the eccentricity, \(\omega\) is the argument of periapsis, \(\Omega\) is the right ascension of the ascending node, \(L\) is the true longitude, and \(\nu\) is the true anomaly. \((f,g)\) are the x and y components of the eccentricity vector in the orbital frame, and \((h,k)\) are the x and y components of the node vector in the orbital frame.

The conversion from Cartesian state (position and velocity) to modified equinoctial elements is given in this Fortran subroutine:

subroutine cartesian_to_equinoctial(mu,rv,evec)

use iso_fortran_env, only: wp => real64

implicit none

real(wp),intent(in)               :: mu    !! central body gravitational parameter (km^3/s^2)
real(wp),dimension(6),intent(in)  :: rv    !! Cartesian state vector
real(wp),dimension(6),intent(out) :: evec  !! Modified equinoctial element vector: [p,f,g,h,k,L]

real(wp),dimension(3) :: r,v,hvec,hhat,ecc,fhat,ghat,rhat,vhat
real(wp) :: hmag,rmag,p,f,g,h,k,L,kk,hh,s2,tkh,rdv

r = rv(1:3)
v = rv(4:6)

rdv      = dot_product(r,v)
rhat     = unit(r)
rmag     = norm2(r)
hvec     = cross(r,v)
hmag     = norm2(hvec)
hhat     = unit(hvec)
vhat     = (rmag*v - rdv*rhat)/hmag
p        = hmag*hmag / mu
k        = hhat(1)/(1.0_wp + hhat(3))
h        = -hhat(2)/(1.0_wp + hhat(3))
kk       = k*k
hh       = h*h
s2       = 1.0_wp+hh+kk
tkh      = 2.0_wp*k*h
ecc      = cross(v,hvec)/mu - rhat
fhat(1)  = 1.0_wp-kk+hh
fhat(2)  = tkh
fhat(3)  = -2.0_wp*k
ghat(1)  = tkh
ghat(2)  = 1.0_wp+kk-hh
ghat(3)  = 2.0_wp*h
fhat     = fhat/s2
ghat     = ghat/s2
f        = dot_product(ecc,fhat)
g        = dot_product(ecc,ghat)
L        = atan2(rhat(2)-vhat(1),rhat(1)+vhat(2))

evec = [p,f,g,h,k,L]

end subroutine cartesian_to_equinoctial

Where cross is a cross product function and unit is a unit vector function (\(\mathbf{u} = \mathbf{v} / |\mathbf{v}|\)).

The inverse routine (modified equinoctial to Cartesian) is:

subroutine equinoctial_to_cartesian(mu,evec,rv)

use iso_fortran_env, only: wp => real64

implicit none

real(wp),intent(in)               :: mu   !! central body gravitational parameter (km^3/s^2)
real(wp),dimension(6),intent(in)  :: evec !! Modified equinoctial element vector : [p,f,g,h,k,L]
real(wp),dimension(6),intent(out) :: rv   !! Cartesian state vector

real(wp) :: p,f,g,h,k,L,s2,r,w,cL,sL,smp,hh,kk,tkh
real(wp) :: x,y,xdot,ydot
real(wp),dimension(3) :: fhat,ghat

p = evec(1)
f = evec(2)
g = evec(3)
h = evec(4)
k = evec(5)
L = evec(6)

kk      = k*k
hh      = h*h
tkh     = 2.0_wp*k*h
s2      = 1.0_wp + hh + kk
cL      = cos(L)
sL      = sin(L)
w       = 1.0_wp + f*cL + g*sL
r       = p/w
smp     = sqrt(mu/p)
fhat(1) = 1.0_wp-kk+hh
fhat(2) = tkh
fhat(3) = -2.0_wp*k
ghat(1) = tkh
ghat(2) = 1.0_wp+kk-hh
ghat(3) = 2.0_wp*h
fhat    = fhat/s2
ghat    = ghat/s2
x       = r*cL
y       = r*sL
xdot    = -smp * (g + sL)
ydot    =  smp * (f + cL)

rv(1:3) = x*fhat + y*ghat
rv(4:6) = xdot*fhat + ydot*ghat

end subroutine equinoctial_to_cartesian

The equations of motion are:

subroutine modified_equinoctial_derivs(mu,evec,scn,evecd)

use iso_fortran_env, only: wp => real64

implicit none

real(wp),intent(in)               :: mu    !! central body gravitational parameter (km^3/s^2)
real(wp),dimension(6),intent(in)  :: evec  !! modified equinoctial element vector : [p,f,g,h,k,L]
real(wp),dimension(3),intent(in)  :: scn   !! Perturbation (in the RSW frame) (km/sec^2)
real(wp),dimension(6),intent(out) :: evecd !! derivative of `evec`

real(wp) :: p,f,g,h,k,L,c,s,n,sqrtpm,sl,cl,s2no2w,s2,w

p      = evec(1)
f      = evec(2)
g      = evec(3)
h      = evec(4)
k      = evec(5)
L      = evec(6)
s      = scn(1)
c      = scn(2)
n      = scn(3)
sqrtpm = sqrt(p/mu)
sl     = sin(L)
cl     = cos(L)
s2     = 1.0_wp + h*h + k*k
w      = 1.0_wp + f*cl + g*sl
s2no2w = s2*n / (2.0_wp*w)

evecd(1) = (2.0_wp*p*c/w) * sqrtpm
evecd(2) = sqrtpm * ( s*sl + ((w+1.0_wp)*cl+f)*c/w - g*n*(h*sl-k*cl)/w )
evecd(3) = sqrtpm * (-s*cl + ((w+1.0_wp)*sl+g)*c/w + f*n*(h*sl-k*cl)/w )
evecd(4) = sqrtpm * s2no2w * cl
evecd(5) = sqrtpm * s2no2w * sl
evecd(6) = sqrt(mu*p)*(w/p)**2 + sqrtpm * ((h*sl-k*cl)*n)/w

end subroutine modified_equinoctial_derivs

The 3-element scn vector contains the components of the perturbing acceleration (e.g. from a thrusting engine, other celestial bodies, etc.) in the directions perpendicular to the radius vector, in the "along track" direction, and the direction of the angular momentum vector.

Integration using modified equinoctial elements will generally be more efficient than a Cowell type integration of the Cartesian state elements. You only have to worry about the one singularity (note that there is also an alternate retrograde formulation where the singularity is moved to \(i = 0\)).

See also

• Gooding Universal Elements [degenerateconic.com]
• modified_equinoctial_module in Fortran Astrodynamics Toolkit

References

1. M. J. H. Walker, B. Ireland, Joyce Owens, "A Set of Modified Equinoctial Orbit Elements" Celestial Mechanics, Volume 36, Issue 4, p 409-419. (1985)
2. Walker, M. J. H, "Erratum - a Set of Modified Equinoctial Orbit Elements" Celestial Mechanics, Volume 38, Issue 4, p 391-392. (1986)
3. Broucke, R. A. & Cefola, P. J., "On the Equinoctial Orbit Elements" Celestial Mechanics, Volume 5, Issue 3, p 303-310. (1972)
4. D. A. Danielson, C. P. Sagovac, B. Neta, L. W. Early, "Semianalytic Satellite Theory", Naval Postgraduate School Department of Mathematics, 1995
5. P. J. Cefola, "Equinoctial orbit elements - Application to artificial satellite orbits", Astrodynamics Conference, 11-12 September 1972.

Gooding Universal Elements

There are various different 6-element sets that can be used to represent the state (position and velocity) of a spacecraft. The classical Keplerian element set (\(a\), \(e\), \(i\), \(\Omega\), \(\omega\), \(\nu\)) is perhaps the most common, but not the only one. The Gooding universal element set [1] is:

• \(\alpha\) : twice the negative energy per unit mass (\(2 \mu / r - v^2\) or \(\mu/a\))
• \(q\) : periapsis radius magnitude (a.k.a. \(r_p\))
• \(i\) : inclination
• \(\Omega\) : right ascension of the ascending node
• \(\omega\) : argument of periapsis
• \(\tau\) : time since last periapsis passage

The nice thing about this set is that it is always defined, even for parabolas and rectilinear orbits (which other orbital element sets tend to ignore). A table from [1] shows the relationship between \(\alpha\) and \(q\) for the different types of conics:

One advantage of this formulation is for trajectory optimization purposes, since if \(\alpha\) is an optimization variable, the orbit can transition from an ellipse to a hyperbola while iterating without encountering any singularities (which is not true if semimajor axis and eccentricity are being used, since \(a \rightarrow \infty\) for a parabola and then becomes negative for a hyperbola).

Gooding's original paper from 1988 includes Fortran 77 code for the transformations to and from position and velocity, coded with great care to avoid loss of accuracy due to numerical issues. Two other papers [2-3] also include some components of the algorithm. Modernized versions of these routines can be found in the Fortran Astrodynamics Toolkit. For example, a modernized version of ekepl1 (from [2]) for solving Kepler's equation:

pure function ekepl1(em, e)
  !!  Solve kepler's equation, `em = ekepl - e*sin(ekepl)`,
  !!  with legendre-based starter and halley iterator

  implicit none

  real(wp) :: ekepl1
  real(wp),intent(in) :: em
  real(wp),intent(in) :: e

  real(wp) :: c,s,psi,xi,eta,fd,fdd,f

  real(wp),parameter :: testsq = 1.0e-8_wp

  c = e*cos(em)
  s = e*sin(em)
  psi = s/sqrt(1.0_wp - c - c + e*e)

  do
    xi = cos(psi)
    eta = sin(psi)
    fd = (1.0_wp - c*xi) + s*eta
    fdd = c*eta + s*xi
    f = psi - fdd
    psi = psi - f*fd/(fd*fd - 0.5_wp*f*fdd)
    if (f*f < testsq) exit
   end do

   ekepl1 = em + psi

end function ekepl1

The Gooding universal element set can be found in NASA's Copernicus tool, and is being used in the design and optimization of the Artemis I trajectory.

See also

• Lambert's Problem [degenerateconic.com]
• gooding_module in Fortran Astrodynamics Toolkit

References

1. R. H. Gooding, "On universal elements, and conversion procedures to and from position and velocity" Celestial Mechanics 44 (1988), 283-298.
2. A. W. Odell, R. H. Gooding, "Procedures for solving Kepler's equation" Celestial Mechanics 38 (1986), 307-334.
3. R. H. Gooding, A. W. Odell. "The hyperbolic Kepler equation (and the elliptic equation revisited)" Celestial Mechanics 44 (1988), 267-282.
4. J, Williams, J. S. Senent, and D. E. Lee, "Recent improvements to the copernicus trajectory design and optimization system", Advances in the Astronautical Sciences 143, January 2012 (AAS 12-236)
5. G. Hintz, "Survey of Orbit Element Sets", Journal of Guidance Control and Dynamics, 31(3):785-790, May 2008

Fortran Filesystem Blues

One of my major pet peeves about Fortran is that it contains virtually no high-level access to the file system. The file system is one of those things that the Fortran standard pretends doesn't exist. It's one of those idiosyncratic things about the Fortran standard, like how it uses the word "processor" instead of "compiler", and doesn't acknowledge the existence of source code files so it doesn't bother to recommend a file extension for Fortran (gotta preserve that backward compatibility just in case punch cards come back).

There are various hacky workarounds to do various things that poor Fortran programmers have had to use for decades. Basically, all you have is OPEN, CLOSE, READ, WRITE, and INQUIRE. Here are a few examples:

Deleting a file

Here's a standard Fortran routine that can be used to delete a file:

function delete_file(name) result(istat)

  implicit none

  character(len=*),intent(in) :: name
  integer :: istat

  integer :: iunit

  open(name=name,newunit=iunit,status='OLD',iostat=istat)
  if (istat==0) close(iunit,status='DELETE',iostat=istat)

end function delete_file

As you can see, it's just a little trick, using a feature of CLOSE to delete a file after closing it (but first, we had to open it). Of course, note that the error codes are non-standard (the Fortran standard also doesn't consider this information important and lets the different compilers do whatever they want -- so of course, they all do something different). And of course, there is no exception handling in Fortran either (a post for another time).

Creating a directory

There is absolutely no standard ability in Fortran to create or delete a directory. Again, Fortran basically pretends that directories don't exist. Note that the Intel Fortran Compiler provides a super useful non-standard extension of the INQUIRE statement to allow it to be used to get information about directories. It's a pretty ridiculous state of affairs (Fortran added a bazillion IEEE routines in Fortran 2018 for some reason that nobody needed, but it still doesn't have something like this that everybody has needed for decades). Intel's IFPORT portability module provides many useful (non-standard) routines for accessing the file system (for example DELFILESQQ for deleting files and DELDIRQQ for deleting a directory). I use these all the time and the fact that they are non-standard and not present in other compilers is a major source of annoyance.

For some file or directory operations, you can always resort to using system calls, and thus have to provide different methods for different platforms (of course Fortran provides no standard way to get any information about the platform, so you have to resort to non-standard, possibly compiler-specific, preprocessor directives to do that). For example, to create a directory you could use:

subroutine make_directory(name)

  implicit none

  character(len=*),intent(in) :: name

  call execute_command_line ('mkdir '//name)

end subroutine make_directory

That one was easy because "mkdir" works on macOS, Window, and Linux. See previous post about getting command line calls into strings, which can also be useful for some things like this.

Getting a list of files

How about getting a list of all the files that match a specific pattern? Again, this can be done if you are using the Intel compiler using the IFPORT module GETFILEINFOQQ routine:

function get_filenames(pattern, maxlen) result(list)

  use ifport

  implicit none

  character(len=*),intent(in) :: pattern
  integer,intent(in) :: maxlen
  character(len=maxlen),dimension(:),allocatable :: list

  integer(kind=int_ptr_kind()) :: handle
  type(file$infoi8) :: info
  integer :: length

  allocate(list(0))
  handle = file$first
  do
    length = GETFILEINFOQQ(trim(pattern), info, handle)
    if ((handle==file$last) .or. (handle==file$error)) exit
    list = [list, info%name]
  end do

end function get_filenames

The file$infoi8 type is some super-weird non-standard Intel thing. Note that this example only returns the file names, not the full paths. To get the full path is left as an exercise to the reader. Also note that Fortran doesn't provide a good string class, so we just set a maximum string length as an argument to this function (but note that Intel is only returning a 255-character string anyway).

Final thoughts

There are various libraries out there that can do some of this. For example, M_system (a Fortran module interface for calling POSIX C system routines). Hopefully, a comprehensive set of filsystem routines will eventually make it into the Fortran Standard Library.

See also

• Portability Library Routines [Intel]
• getfileinfoqq -- Intel Portability Function that returns information about the specified file. [Intel]
• Interface to POSIX I/O API [Fortran stdlib issue]

Atmosphere Models

The drag force magnitude on a spacecraft can be computed using the following equation:

Which can be easily implemented in the following Fortran function:

pure function drag(density, v, area, cd)

    use iso_fortran_env, only: wp => real64

    implicit none

    real(wp) :: drag !! drag force magnitude (N)
    real(wp), intent(in) :: density !! atmospheric density (kg/m^3)
    real(wp), intent(in) :: v       !! relative velocity magnitude (m/s)
    real(wp), intent(in) :: area    !! spacecraft cross-sectional area (m^2)
    real(wp), intent(in) :: cd      !! spacecraft drag coefficient

    drag = 0.5_wp * density * v**2 * area * cd

end function drag

The density \(\rho\) is computed using an atmosphere model. Here is a listing of some atmosphere model codes available in Fortran:

• 1976 Standard Atmosphere
• INPE-atmosphere-models -- My archive of models originally from this page (which now seems to be gone). Includes five models:
  • Jacchia 70
  • Jacchia-Roberts
  • Jacchia 77
  • LaFontaine-Hughes
  • Jacchia-Lineberry
• MSIS86
• MSISE90
• NRLMSISE00
• NRLMSIS 2.0
• Earth Global Reference Atmospheric Model 2010 (Earth GRAM 2010)

See also

• Fox & Mcdonald, "Introduction To Fluid Mechanics", Wiley
• Drag Coefficient [Wikipedia]
• U.S. Standard Atmosphere, 1976, U.S. Government Printing Office, Washington, D.C., 1976.
• ModelWeb Catalogue and Archive [NASA GSFC]
• W. W. Vaughan, "Guide to Reference and Standard Atmospheric Models", AIAA G-003C-2010, January 2010

Great Circle

If we assume that the Earth is a sphere of radius \(r\), the shortest distance between two sites (\([\phi_1,\lambda_1]\) and \([\phi_2,\lambda_2]\)) on the surface is the great-circle distance, which can be computed in several ways:

Law of Cosines equation:

This equation, while mathematically correct, is not usually recommended for use, since it is ill-conditioned for sites that are very close to each other.

Using the haversine equation:

The haversine one is better, but can be ill-conditioned for near antipodal points.

Using the Vincenty equation:

This one is accurate for any two points on a sphere.

Code

Fortran implementations of these algorithms are given here:

Law of Cosines:

pure function gcdist(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: clat1,clat2,slat1,slat2,dlon,cdlon

clat1 = cos(lat1)
clat2 = cos(lat2)
slat1 = sin(lat1)
slat2 = sin(lat2)
dlon = long1-long2
cdlon = cos(dlon)

d = r * acos(slat1*slat2+clat1*clat2*cdlon)

end function gcdist

Haversine:

pure function gcdist_haversine(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: dlat,clat1,clat2,dlon,a

dlat  = lat1-lat2
clat1 = cos(lat1)
clat2 = cos(lat2)
dlon  = long1-long2

a = sin(dlat/2.0_wp)**2 + clat1*clat2*sin(dlon/2.0_wp)**2
d = r * 2.0_wp * asin(min(1.0_wp,sqrt(a)))

end function gcdist_haversine

Vincenty:

pure function gcdist_vincenty(r,lat1,long1,lat2,long2) &
              result(d)

implicit none

real(wp) :: d !! great circle distance from 1 to 2 [km]
real(wp),intent(in) :: r     !! radius of the body [km]
real(wp),intent(in) :: long1 !! longitude of first site [rad]
real(wp),intent(in) :: lat1  !! latitude of the first site [rad]
real(wp),intent(in) :: long2 !! longitude of the second site [rad]
real(wp),intent(in) :: lat2  !! latitude of the second site [rad]

real(wp) :: c1,s1,c2,s2,dlon,clon,slon

c1   = cos(lat1)
s1   = sin(lat1)
c2   = cos(lat2)
s2   = sin(lat2)
dlon = long1-long2
clon = cos(dlon)
slon = sin(dlon)

d = r*atan2(sqrt((c2*slon)**2+(c1*s2-s1*c2*clon)**2), &
    (s1*s2+c1*c2*clon))

end function gcdist_vincenty

Test cases

Using a radius of the Earth of 6,378,137 meters, here are a few test cases (using double precision (real64) arithmetic):

Close points

The distance between the following very close points \(([\phi_1,\lambda_1], [\phi_2,\lambda_2]) =([0, 0.000001], [0.0, 0.0])\) (radians) is:

Houston to New York

The distance between the following two points \(([\phi_1,\lambda_1], [\phi_2,\lambda_2]) =([29.97, -95.35], [40.77, -73.98])\) (deg) is:

Antipodal points

The distance between the following two antipodal points \(([\phi_1,\lambda_1], [\phi_2,\lambda_2]) = ([0, 0], [0.0, \pi/2])\) (radians) is:

Nearly antipodal points

The distance between the following two nearly antipodal points \(([\phi_1,\lambda_1], [\phi_2,\lambda_2]) = ([10^{-8}, 10^{-8}], [0.0, \pi/2])\) (radians) is:

Other methods

For a more accurate result, using an oblate-spheroid model of the Earth, the Vincenty inverse geodetic algorithm can be used (the equation above is just a special case of the more general algorithm when the ellipsoidal major and minor axes are equal). For the Houston to New York case, the difference between the spherical and oblate-spheroid equations is 282 meters.

See also

• GCDIST, a Fortran 77 version of the Law of Cosines formula, 1996.
• J. Williams, Direct and Inverse, January 31, 2016 [degenerateconic.com]
• GIS FAQ Q5.1: Great circle distance between 2 points [Movable Type Scripts]

Visual Studio Code + Fortran

With the right plugins, Microsoft Visual Studio Code can be turned into a very good modern Fortran IDE. Here are some of my favorites:

• FORTRAN IntelliSense (hansec.fortran-ls). To use this, you must pip install fortran-language-server.
• vscode-modern-fortran-formatter (yukiuuh.vscode-modern-fortran-formatter). To use this, you must pip install fprettify.
• Modern Fortran (krvajalm.linter-gfortran)
• Fortran Breakpoint Support (ekibun.fortranbreaker)

Historically, good full-featured IDEs for modern Fortran have been hard to find. Intel Fortran integrates with the full MS Visual Studio, and I used that for many years. But, except for syntax highlighting and automatic determination of compile order, it is basically terrible. I always have to disable the "IntelliSense" features since it brings the interface to a crawl (it's been like this for years and they don't seem to care enough to fix it). Unfortunately, Intel's excellent Fortran debugger also only works from within Visual Studio, so if you need to use it, you are stuck in that environment.

I've pretty much converted to use VSCode now for all Fortran coding (on both Mac and Windows). In the past, I have also used TextWrangler, Sublime Text, and Atom.

See also:

• vscode-icons -- A nice set of icons for various file formats, although the Fortran icon is terrible. I submitted a pull request to use the new one from fortran-lang.org. As of Oct 2021, it has yet to be merged.
• "Parse File…" with lagging visual studio [Intel Fortran Forum]. Note: I believe one of my comments on this post was censored by the admins. I distinctly remember asking if we were a bunch of suckers for continuing to give money to a multi-billion dollar corporation (Intel) when they can't be bothered to fix their defective product.
• Fortran IDEs at fortran-lang.org.
• Is Visual Studio Code Really The Best Code Editor?, February 28, 2021. [tabnine.com]

The Prehistory of Fortran-Lang

At the recent FortranCon2021, I was amused to see in the presentation "The State of Fortran 2021" that the first bullet point on the slide "Formation of Fortran-Lang" was:

August 2019 Conversations on Twitter between Ondřej Čertík, Milan Curcic, and Jacob Williams bring out common perceived shortcomings in the Fortran ecosystem

Laurence Kedward, et. al., "The State of Fortran 2021", Fortran2021, Sept. 2021

So, now, over two years later, I think it's finally time to tell my side of the story. 😀

I was first introduced to Fortran 25 years ago, in 1997 (my only previous programming experience was with BASIC on Apple ]['s). We used the WATFOR FORTRAN 77 compiler on DOS machines in an undergraduate programming class for engineers. Little did I know (or little did the university know) that Fortran 90 had been published years earlier (and indeed Fortran 95 was published at the end of that very semester).

Over the course of my academic career, I encountered and used Fortran a few times, mostly FORTRAN 77, and mostly to solve homework problems or do class projects. I came to use IMSL, Numerical Recipes in FORTRAN 77, Harwell Subroutine Library, and had some minimal exposure to Fortran 90 one semester when asked to modernize some old code by one of the professors for some project that I've forgotten about (something to do with gravity anomalies). But mostly I used a lot of Matlab, which is very popular among engineering students (due to cheap student licenses). I also used a little MathCad, some Mathematica, REALBasic (since renamed Xojo), and Visual Basic (classic). There was then no exposure in the aerospace engineering curriculum of two major US universities to any sort of formal software development practices (for all I know this is still the case). I never once encountered version control, unit testing, continuous integration, etc. Maybe those weren't as common then? (Subversion didn't come out until 2000, and Git until 2005).

When I entered professional life in late 2006, it wasn't long before I became the lead developer of a Fortran tool at NASA called Copernicus. Originally created by one of my former professors, Copernicus is a tool for spacecraft trajectory optimization and design. While working on Copernicus, I was introduced to Fortran 90/95 in earnest. Copernicus was originally developed on Windows with Compaq Visual Fortran (which, by the time I came along, had recently become defunct) and one of my first projects was to convert the code over to use the Intel Fortran compiler (which we have used ever since). During this period I first became aware of modern-ish (Fortran 95) concepts such as user derived types, operator overloading, function overloading, etc. Eventually Fortran 2003 came along, which was quite a revelation, and as soon as non-buggy implementations of the new features appeared in the Intel compiler (it took a while) I was using it like crazy. It's hard to imagine programming Fortran now without modern features such as type bound procedures. FORTRAN 77 (and indeed, even Fortran 95), to me, now seems like something out of the dark ages. And so I came to really appreciate the Modern Fortran programming language, and wondered why more people didn't use it (or seem to even know about it).

The fact was, all was not well in the larger Fortran community. Really, there was very little community. Most of the code I encountered on the internet was still FORTRAN 77, like the last 25 years never happened. There were also few if any Fortran blogs (there were a few of exceptions: Fortran Dev [now defunct], Fortran in a C World [now also defunct], and Steve Lionel's Doctor Fortran column at Intel [still updated infrequently]). One of the bright spots was the Intel Fortran Compiler forum, which was always full of helpful folks (Steve, of course, the most helpful of all). There were very few websites to speak of that had much useful Modern Fortran content. Fortran.com is a mostly-useless site, and fortran.org has been held hostage by a domain name squatter for years. Ondřej Čertík's site fortran90.org is good (but of course it is folly to put the standard name in the URL since that will eventually be out of date). Fortranwiki.org, administrated by Jason Blevins is also a good site.

On June 26, 2014, I started this blog, with the intent to focus on modern Fortran and astrodynamics algorithms. Earlier that year, I had also discovered GitHub and started to work on some projects there as well. At work, we were still using Subversion for version control at this time, but I wanted to learn more about git. The blog and my presence on GitHub was an experiment to attempt to publicly develop an ecosystem of Modern Fortran libraries (either by creating new ones, or refactoring old codes that hadn't been updated in decades that nobody else seemed to want to update), as well as to advocate for Fortran and show what it could do. On GitHub I encountered others who were also working on Fortran tools, namely Stefano Zaghi (author of lots of great libraries and tools, including FoBiS, the Fortran Building System for poor people, which is very dear to my heart), Zaak Beekman (who helped me greatly with JSON-Fortran and introduced me to continuous integration), Chris MacMackin (author of the ford automatic documentation generator for Fortran, which I contributed a very very small amount to), and a few others. Zaak also started a GitHub group called the Fortran F/OSS Programmers Group, which I also participated in and view as sort of a precursor to the fortran-lang group. Zaak and I coauthored a paper in 2018 on Modern Fortran applications in my field of spacecraft trajectory optimization, which I consider to be my magnum opus.

In April 2017, I was approached by Manning Publications, who wanted to discuss the possibility of writing a Fortran book. I didn't really have time to do that so I declined. But luckily, they also approached Milan Curcic, who did end up writing it (Modern Fortran, 2020) and it turned out excellent.

I currently have 67 repositories on GitHub, the vast majority of them being Fortran libraries, and they are used by lots of people all over the world in various fields of science and engineering. But one thing I didn't really have the aptitude or know-how to do was to form an online community. Years ago I considered buying the fortran.io domain when I noticed it was still up for grabs, but it was eventually bought by Nick Doiron around 2016 (see FORTRAN.io -- finally, a Fortran Web Framework) which made the rounds somewhat on the tech internet, basically just as a novelty. I guess the .io fad for tech domain names is over now anyway?

There just wasn't a website that really served as a focal point for the entire Fortran community. It was a huge image problem. This was one of the things I lamented in my (infamous?) August 2019 post. Luckily, Milan had purchased fortran-lang.org sometime in 2018, and that is what he activated for the new site in April 2020. Milan and Ondřej have already written about their amazing efforts to rejuvenate the Fortran community and ecosystem (I first became aware of Ondřej in April 2019 when he announced LFortran, which totally blew my mind). The fortran-lang website, the Fortran Standard Library, the Fortran Package Manager, and everything else that has come out these activities has been spectacular. I have actually contributed very little to all of this, although I did help with the logo (see this previous post for the story on that). I also continue to toil away on my libraries like before (recently getting a good number of them set up to compile with the new Fortran Package Manager).

So, now you know my side of the story. The future of Fortran is looking bright, and I'm pleased to have contributed in some small way. Now get off my lawn.

See also

• J. Williams, R. D. Falck, and I. B. Beekman, "Application of Modern Fortran to Spacecraft Trajectory Design and Optimization", 2018 Space Flight Mechanics Meeting, 8–12 January 2018, AIAA 2018-145.
• J. Williams, End of Life [degenerateconic.com] August 3, 2019
• J. Williams, New Blood [degenerateconic.com] April 26, 2020
• M. Curcic, First year of Fortran-lang [medium.com] Dec 21, 2020
• O. Čertík, Resurrecting Fortran [ondrejcertik.com] March 12, 2021

FortranCon 2021

FortranCon 2021 (September 23-24, 2021), the second international conference dedicated to the Fortran programming language, has concluded. Once again, it was an amazing conference, with lots of great talks, including a fortran-lang minisymposium. This year there were about 160 registered attendees and about 80 people on average on Zoom and Slack at any given time.

Lists of talks

• Alfio Lazzaro Welcome
• Arjen Markus Handling dimensions and units of measure in a program
• Mohd Furquan Some adventures with compile time evaluation
• Alexey Cherkaev Fortran in Hydrometallurgy: Practice and Research
• Leandro Martinez Invited Presentation: A computational chemist perspective on Fortran vs. Julia programming
• Marshall Ward Bit-reproducible methods of dimensional and rotational verification testing
• Jeff Hammond Standard Fortran on GPUs and its utility in quantum chemistry codes
• Damian Rouson Keynote: Fortran at the Intersection
• Ondřej Čertík LFortran: Interactive LLVM-based Fortran Compiler for Modern Architectures
• Alessandro Fanfarillo HIPFort: Present and Future Directions for Portable GPU Programming in Fortran
• Sarah Tipler SPACE: Three Decades of Spacecraft Power Systems Analysis with Fortran
• Sebastian Costamagna SEACAL: Cross Platform Parallel Hybrid Code for Calculations of Ship Motions
• Rohit Goswami f2py: Two Decades Later
• Nicolas Abreu Netto Test automation on a Fortran Legacy Software: the AnaTest framework
• Jeff Hammond The Case for Asynchronous Task Parallelism in Fortran
• Brad Richardson Your Requirements Specification as an Executable Test Suite
• Laurence Kedward The State of Fortran
• Thirumalai Shaktivel Finish AST generation in LFortran
• Gagandeep Singh Supporting Arrays and Allocatables in LFortran
• Rohit Goswami Implementing Fortran Standardese within LFortran
• Nathaniel Shaffer What's new in the Fortran Standard Library?
• Aman Godara Improving Strings Support in Fortran
• Chetan Karwa et al. Linked lists for stdlib
• Sebastian Ehlert Fortran package manager
• Jakub Jelínek Handling Compiler Flags in fpm
• Tiziano Müller Wrap-up

See also

• FortranCon 2021 Website
• FortranCon 2020 [previous post]
• Countdown to Fortran 2023 Release (Is this real? There's no way to know...)