Degenerate Conic | Converting Legacy C++ Code To Modern Fortran

Degenerate Conic

Algorithms • Modern Fortran Programming • Orbital Mechanics

May 25, 2026

Converting Legacy C++ Code To Modern Fortran

"I am bound Upon a wheel of fire, that mine own tears Do scald like molten lead." - Shakespeare, "King Lear", Act IV, Scene 7

Typically, we see examples of poor saps who think that legacy Fortran codes must be converted to "modern" languages such as C++.<br>Well, let's try to do the reverse! Let's try converting some "legacy" C++ code to modern Fortran. The example I will show is the Jacchia-Roberts atmosphere model from GMAT (General Mission Analysis Tool). This model is widely used for satellite orbit simulations to model satellite drag. There is no modern Fortran implementation publicly available as far as I can tell, so let's create one.<br>And just to see what happens, I'll do this with assistance from AI, particularly Claude Sonnet (initially 4.5, but also trying 4.6) in GitHub Copilot. What could go wrong?

All About the Vibes

First, we start with the C++ files we want to convert. These files are at least 20 years old, with heritage going back even further. The main ones are:

JacchiaRobertsAtmosphere.cpp - the main implementation of the Jacchia-Roberts density model.

SolarFluxReader.cpp - the algorithm to read in a space weather file. There are a couple different file formats supported here, but I'm only going to worry about the CSSI one.

A1Mjd.cpp - Some time conversion routines. This will be important later when trying to diagnose some validation issues.

Using AI is a very interesting and also frustrating experience.<br>It's like working with someone who is a genius, but is also an idiot and sometimes a compulsive liar.<br>The initial conversion got me 80% of the way there but was subtly wrong in various ways.<br>I started with a basic prompt ("Convert this C++ file into modern Fortran") and went from there.<br>There was a lot of back and forth and manual interventions to fix little issues (syntax errors, etc.)<br>In one instance, I tried to get it to identify the source of a bug, and it "thought" for 15 minutes until I finally just gave up and looked at the code myself and fixed it instantly.<br>Sometimes, it would get fixated on one approach that wasn't correct, and I'd have to keep telling it to ignore that approach and try something else (e.g., I had to figure out the correct way to do something, and once I gave it more directions, it was fine and could do it).<br>One thing that was annoying was that it kept trying to generate code from first principles rather than performing a straight up conversion of the C++ code to Fortran. I had to keep reminding it that what I wanted was a faithful translation.<br>At one point, I switched to Sonnet 4.6 and it also found a couple of code typos that Sonnet 4.5 had generated.

The Result

See jacchia-roberts-fortran for the end result. The library provides a class that can be used to initialize the model with a given space weather file (or constant values can be specified). Then a method is provided for computing the density.<br>An example use case is:

program example

use jacchia_roberts_module, only: jacchia_roberts_type<br>use jacchia_roberts_kinds, only: ip, dp

implicit none

type(jacchia_roberts_type) :: jr<br>integer(ip) :: status<br>real(dp) :: density

! example inputs:<br>real(dp),parameter :: rad_earth = 6356.766_dp ! Earth polar radius (km)<br>character(len=*),parameter :: sw_file = 'data/SpaceWeather-All-v1.2.txt' ! space weather file to load<br>real(dp),parameter :: utc_mjd = 59215.5_dp ! MJD: Jan 1, 2021 12:00 UTC<br>real(dp),parameter :: alt_km = 200.0_dp ! geodetic altitude (km)<br>real(dp),parameter :: lat = 20.0_dp ! geodetic latitude (deg)<br>real(dp),dimension(3),parameter :: r = [7000.0_dp, 0.0_dp, 0.0_dp] ! spacecraft position vector (km)<br>real(dp),dimension(3),parameter :: sun = [1.0_dp, 0.0_dp, 0.0_dp] ! sun direction unit vector

! initialize the model:<br>call jr%initialize(rad_earth, filename=sw_file, status=status)

! compute the density:<br>density = jr%density(alt_km, r, sun_vector, lat, utc_mjd)

end program example

The idea is, first you would initialize the class, and then during your simulation, call the model to get the atmospheric density using the current epoch and state of the spacecraft. A flow chart of the process would look a little like this:

The resultant Fortran code is very clean compared to the C++ code. Consider this example (the rho_cor function). The original C++ code:

Real JacchiaRobertsAtmosphere::rho_cor(Real height, Real a1_time, Real geo_lat,<br>GEOPARMS *geo)<br>Real geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha;<br>Real sin_lat, eta_lat;

// Compute geomagnetic activity correction<br>if (height 200.0)<br>geo_cor = 0.012 * geo->tkp + 0.000012 * exp(geo->tkp);<br>else<br>geo_cor = 0.0;

// Compute semiannual variation correction<br>f = (5.876e-7 * pow(height, 2.331) + 0.06328) *<br>exp(-0.002868 * height);<br>day_58 = (a1_time - 6204.5)/365.2422; // @todo - should this use...