Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

Oct 24, 2021

Atmosphere Models

us-standard-atmosphere

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

$$ d = \frac{1}{2} \rho v^2 A C_D$$

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:

See also

Oct 18, 2021

Great Circle

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:

$$ d = r \arccos\bigl(\sin\phi_1\cdot\sin\phi_2+\cos\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)\bigr) $$

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:

$$ d = 2 r \arcsin \sqrt{\sin^2\left(\frac{\Delta\phi}{2}\right)+\cos{\phi_1}\cdot\cos{\phi_2}\cdot\sin^2\left(\frac{\Delta\lambda}{2}\right)} $$

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

Using the Vincenty equation:

$$ d = r \arctan \frac{\sqrt{\left(\cos\phi_2\cdot\sin(\Delta\lambda)\right)^2+\left(\cos\phi_1\cdot\sin\phi_2-\sin\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)\right)^2}}{\sin\phi_1\cdot\sin\phi_2+\cos\phi_1\cdot\cos\phi_2\cdot\cos(\Delta\lambda)} $$

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:

Algorithm Procedure Result (meters)
Law of Cosines gcdist 6.3784205037462689
Haversine gcdist_haversine 6.3781369999999997
Vincenty gcdist_vincenty 6.3781369999999997

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:

Algorithm Procedure Result (meters)
Law of Cosines gcdist 2272779.3057236285
Haversine gcdist_haversine 2272779.3057236294
Vincenty gcdist_vincenty 2272779.3057236290

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:

Algorithm Procedure Result (meters)
Law of Cosines gcdist 20037508.342789244
Haversine gcdist_haversine 20037508.342789244
Vincenty gcdist_vincenty 20037508.342789244

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:

Algorithm Procedure Result (meters)
Law of Cosines gcdist 20037508.342789244
Haversine gcdist_haversine 20037508.342789244
Vincenty gcdist_vincenty 20037508.252588764

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

Oct 17, 2021

Visual Studio Code + Fortran

vscode 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:

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:

Oct 17, 2021

Pelican

pelican

Since 2014, this website has always been a WordPress site. Until now! I just converted it to a static site, powered by Pelican. Pelican is a static site generator, written in Python. I converted all the posts from the old site to markdown using the wordpress-export-to-markdown tool, which then required a little massaging to get into shape for deployment.

To install Pelican, all I had to do was:

pip install "pelican[markdown]"
pip install pelican-render-math

The site is configured using a pelicanconf.py script. Once that is set up, building the site is as simple as:

pelican content

And that's it. It generates only HTML files, which can then be uploaded to the web server. For that I followed the instructions here to set up a GitHub Action to do this automatically. I also made an .htaccess file to redirect all the old WordPress URLs to the new ones, so old links on other sites to this one will still work.

All in all, it works pretty well. I did lose the comments, but I think it's probably worth it since it makes maintenance so much better. The site is now a set of simple markdown files, which I can easily edit locally, and even view the results in Pelican's simple built-in web server. Everything is backed up to GitHub (currently this is in a private repo). And I don't have to worry about keeping WordPress up to date, databases getting corrupted, or a PHP update breaking everything. Pelican themes and layout options are not as sophisticated as WordPress, but I will continue to tweak things over time.

See also

Oct 02, 2021

The Prehistory of Fortran-Lang

The first Fortran program I ever wrote (Jan. 1997).

The first Fortran program I ever wrote (Jan. 1997).

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).

The FORTRAN book from my first (and only) programming class. I still have it.

The FORTRAN book from my first (and only) programming class. I still have it.

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).

[from archive.org]

Compaq Visual Fortran [from archive.org]

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.

The favicon for this website

The favicon for this website

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.

foss

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.

logos

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?

The Fortran-lang logo.

The Fortran-lang logo.

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

Sep 25, 2021

FortranCon 2021

fortran_logo

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

F

See also

Jun 13, 2021

JSON

json-fortran-logo

There are a lot a JSON libraries available for Fortran nowadays. The complete list, as far as I know is:

  • fson Fortran 95 JSON Parser. JSON-Fortran is a fork of this library. (started 2012).
  • YAJL-Fort A Modern Fortran Interface to YAJL. This one is an interface to a C library, and not pure Fortran. See also petaca. (started 2013).
  • JSON-Fortran A Fortran 2008 JSON API. This is my library. As far as I know, this was the first production-ready JSON parser written in modern Fortran. (started 2014).
  • fortjson JSON library written in Fortran 2003. Designed with portability across HPC architectures in mind. (started 2018).
  • jsonff JSON for Fortran. (started 2019).

Example

Here's an example of using JSON-Fortran to read in some data from a JSON string:

program test

 use json_module, wp => json_RK, ip => json_IK

 implicit none

 type(json_file) :: json
 integer(ip) :: t
 real(wp),dimension(:),allocatable :: x

 call json%deserialize('{"t": 1, "x": [2.0, 3.0, 4.0]}')

 call json%get('t', t); write(*,*) t
 call json%get('x', x); write(*,*) x

end program test

This prints:

           1
   2.00000000000000        3.00000000000000        4.00000000000000

Note that JSON-Fortran has all sorts of options for:

  • Using different real (real32, real64, and real128) and integer kinds (int8, int16, int32, and int64).
  • Controlling the JSON format for printing, including number of spaces for indenting, printing vectors on one line only, or full minification with no extra whitespace or line breaks.
  • Support for comments in a JSON file.
  • Multiple ways to get data from a JSON structure, such as RFC 6901 "JSON Pointer" paths or JSONPath) "bracket-notation".
  • Graceful handing of unusual inputs such as NaN or Infinity.
  • Thread safe error handling.

See also

Dec 18, 2020

oneAPI

intel-logo-vector

Intel recently released its new oneAPI Toolkit, which had been in beta for a while. According to their documentation, oneAPI includes various products:

The big news here is that this is all now entirely free to use. They are all available for Windows, Linux, and macOS. The Fortran compiler also now includes all of the Fortran 2018 standard. Apparently, Intel will still sell you a license, but it seems the only reason to buy one would be to get the Premier Support.

See also

Oct 30, 2020

Arbitrary Precision

numbers2

MPFUN2015 is a Fortran arbitrary precision library by David H. Bailey. An archive of the source code can be found on his website here (note that I also have a mirror of the source at GitHub if you want to peruse it). This library seems to be one of the few freely-available production-ready arbitrary precision Fortran libraries available that I can find, although arbitrary precision Fortran code has existed for quite a while (e.g., see Reference 1 from 1971). MPFUN2015 is pretty easy to use, but I think could use a little bit of modernization.

See also

  1. L. C. Maximon, Fortran Program for Arbitrary Precision Arithmetic, National Bureau of Standards, NBS 10-563, April 1, 1971.
  2. D. H. Bailey, A Thread-Safe Arbitrary Precision Computation Package (Full Documentation), August 23, 2019.

Jul 03, 2020

FortranCon 2020

fortrancon2020

FortranCon 2020 was the first international conference entirely focused on the Fortran programming language. It took place in Zurich, Switzerland, but because of the current pandemic, it was almost entirely virtual. It was a terrific two-day conference with a lot of amazing presentations from people all around the world.

Here is the list of talks:

I presented a talk on Copernicus, the spacecraft trajectory design and optimization software that I develop for NASA. We recently completed the 5.0 release of Copernicus, which included some major upgrades and a new GUI. The core of Copernicus, of course, is written in Fortran. Copernicus is being used to design the trajectories that will return humans to the Moon.

A mind numbing chart from my talk on Copernicus at FortranCon 2020.

A mind numbing chart from my talk on Copernicus at FortranCon 2020.

Many thanks to the organizers of FortranCon 2020. I hope this is just the first of many FortranCon's to come.

See also

Update (July 26, 2020): added the missing links

← Previous Next → Page 2 of 14