Degenerate Conic

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++. 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. 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. It's like working with someone who is a genius, but is also an idiot and sometimes a compulsive liar. The initial conversion got me 80% of the way there but was subtle wrong in various ways. I started with a basic prompt ("Convert this C++ file into modern Fortran") and went from there. There was a lot of back and forth and manual interventions to fix little issues (syntax errors, etc.) 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. 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). 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. 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. An example use case is:

program example

use jacchia_roberts_module, only: jacchia_roberts_type
use jacchia_roberts_kinds,  only: ip, dp

implicit none

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

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

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

! compute the density:
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,
                                       GEOPARMS *geo)
{
   Real geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha;
   Real sin_lat, eta_lat;

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

   // Compute semiannual variation correction
   f = (5.876e-7 * pow(height, 2.331) + 0.06328) *
      exp(-0.002868 * height);
   day_58 = (a1_time - 6204.5)/365.2422; // @todo - should this use GmatConstants value?
   tausa = day_58 + 0.09544 * (
         pow( 0.5*(1.0 + sin(2*GmatMathConstants::PI*day_58 +6.035)), 1.65 )  - 0.5);
   alpha = sin(4.0*GmatMathConstants::PI*tausa + 4.259);
   g = 0.02835 + (0.3817 + 0.17829 * sin(2*GmatMathConstants::PI*tausa + 4.137)) *
         alpha;
   semian_cor = f * g;

   // Compute seasonal latitudinal variation
   sin_lat = sin(geo_lat);
   eta_lat = sin(2.0*GmatMathConstants::PI*day_58 + 1.72) * sin_lat * fabs(sin_lat);
   slat_cor = 0.014 * (height - 90.0) * eta_lat *
      exp(-0.0013 * (height - 90.0) * (height - 90.0));

   return pow(10.0, geo_cor + semian_cor + slat_cor);
}

And the Fortran equivalent:

pure function rho_cor(height, utc_mjd, geo_lat, geo) result(correction)
    real(dp), intent(in) :: height !! Spacecraft altitude (km)
    real(dp), intent(in) :: utc_mjd !! UTC Modified Julian Date
    real(dp), intent(in) :: geo_lat !! Geodetic latitude of spacecraft (radians)
    type(geoparms_type), intent(in) :: geo !! Geomagnetic parameters
    real(dp) :: correction !! Correction factor (multiplicative)

    real(dp) :: geo_cor, semian_cor, slat_cor, f, g, day_58, tausa, alpha, sin_lat, eta_lat

    ! Compute geomagnetic activity correction
    if (height < 200.0_dp) then
        geo_cor = 0.012_dp * geo%tkp + 0.000012_dp * exp(geo%tkp)
    else
        geo_cor = 0.0_dp
    end if

    ! Compute semiannual variation correction
    f      = (5.876e-7_dp * height**2.331_dp + 0.06328_dp) * exp(-0.002868_dp * height)
    day_58 = (utc_mjd - 36204.0_dp) / 365.2422_dp ! years since Jan 1, 1958 midnight
    tausa  = day_58 + 0.09544_dp * &
            ((0.5_dp * (1.0_dp + sin(2.0_dp * PI * day_58 + 6.035_dp)))**1.65_dp - 0.5_dp)
    alpha  = sin(4.0_dp * PI * tausa + 4.259_dp)
    g      = 0.02835_dp + (0.3817_dp + 0.17829_dp * sin(2.0_dp * PI * tausa + 4.137_dp)) * alpha
    semian_cor = f * g

    ! Compute seasonal latitudinal variation
    sin_lat  = sin(geo_lat)
    eta_lat  = sin(2.0_dp * PI * day_58 + 1.72_dp) * sin_lat * abs(sin_lat)
    slat_cor = 0.014_dp * (height - 90.0_dp) * eta_lat * &
                exp(-0.0013_dp * (height - 90.0_dp)**2)

    correction = 10.0_dp**(geo_cor + semian_cor + slat_cor)

end function rho_cor

The main differences in the Fortran are:

• The real kinds use the dp kind variable, allowing us to easily compile the code with a specified real kind (e.g., single, double, or quad precision). This is one of the neat features of Fortran, and how I write all my codes.
• The crazy GSFC MJD (a1_time) is replaced with real MJD, so we had to modify that one constant (the epoch of Jan 1, 1958).
• Normal syntactical differences: abs instead of fabs, no semicolons, ** operator instead of the pow function, then/endif instead of {} symbols, etc. The C++ code using (height - 90.0) * (height - 90.0) instead of (height - 90.0_dp)**2 is also kind of funny.

I find the Fortran version easier to read without the extra C++ clutter, but maybe it's just because i'm more used to it? This aspect of Fortran (it being easy to read) isn't always appreciated or discussed in these "language war" type discussions, but it's very real in my opinion, especially for code is that mostly math.

Here's an example output density plot (generated using my pyplot-fortran library) for three epochs near the max/min/average parts of the solar cycle:

Validation

One of the most interesting parts of the process was using AI to help validate the result. I had no unit tests to start with, but I had two sources of "truth":

• GMAT itself, which I could run to generate density data for different altitudes, epochs, etc. and compare it to the Fortran results.
• An older Fortran implementation of Jacchia-Roberts. This was a Fortran 77 implementation by Valdemir Carrara (Instituto Nacional de Pesquisas Espaciais, INPE) from the 80s, which I obtained from his website, which is no longer there (I archived the codes on GitHub in 2019 thinking I would do something with them one day, and now that day has arrived). This was also another AI opportunity, since the space weather component of this code was not preserved, so I had Copilot generate a wrapper so that it could call the modern Fortran module we had produced instead.

AI seems to get "unit test happy" sometimes, spewing out numerous long unit tests that do something simple and then prints up elaborate tables with unicode icons and paragraphs of text that tells you what it just tested. That seemed a bit excessive, and I ended up deleting many of the tests it generated, but kept a good number of them. Overall it was quite useful for helping to generate and evaluate tests. When tests failed (e.g., the new results didn't match the other two data sets) I would tell it to investigate to locate the source of the discrepancy, which it could sometimes do successfully, and sometimes not. The AI went on a flight of fancy for a while trying to figure out a timing discrepancy. It did help me debug it, but it couldn't quite figure it out. At one point I asked if it could go to GitHub and look at the full GMAT code that was there, it told me it could and provided an unconvincing reply suspiciously quickly that it had done so. I guess it was lying to me (Copilot I thought we were friends). I ended up downloading the extra file manually, putting it in the workspace, and directing the AI's attention to it, which seemed to help.

During the validation process, the AI did help me identify three possible issues with the original C++ code. I link the three bug reports I submitted here:

1. PrepareKpData uses wrong row index for F10.7a - This is one that the AI detected that I would never have noticed in a million years. It looks like a copy/paste type bug but I'm not 100% sure if it's a real issue or not.
2. Legacy UTC calculation used by JacchiaRobertsAtmosphere - This is one that I eventually noticed while trying to debug timing discrepancies during the validation process comparing the results from GMAT to the new version. GMAT (and other GSFC tools) use this crazy alternate definition of modified Julian date (MJD). So while trying to figure out how to convert what it was outputting in a real MJD, I realized they were using an approximation of UTC (it turns out the correct UTC value output from GMAT is slightly different from what is being used internally by the atmosphere model). It looks like something that was a holdover from the earlier tools, but I'm not sure if there's a good reason to keep using it. According to the git history when this file was first committed, the file was originally created for the GSS project in 1995 (I presume this is referring to the Generalized Support Software project at GSFC).
3. Low-altitude Jacchia-Roberts bugs - This one is very interesting. There appear to be bugs in the GMAT code for low (<125 km) altitudes. Note that this part of the code is there but not used, since an exception is raised for this case, even though the Jacchia model should support the low altitude ranges (in Roberts' paper, he says his results "are identical" to Jacchia's original equations in the 90 - 125 km altitude regime). According to the comments and git history on SourceForge, this file was committed in 2004 and the source was the Windows Swingby code (Swingby was another legacy NASA tool that was one of the progenitors of STK/Astrogator). The roots function comments says it was converted from existing Fortran code (likely from GTDS, which included the Jacchia-Roberts model). However, the polynomial deflation technique was written in 1993 based on the "Numerical Recipes in C" book (note that Numerical Recipes talks about the need to "polish" a root found by this method, something not implemented in this code. Maybe that's where they went wrong?) Is it possible that Swingby had this same bug all those years ago? Or maybe my Fortran conversion is the buggy one? In any event, I replaced this algorithm in the Fortran version with an alternate root finding approach (based on the one in the INPE code) which seems to work fine. I also tested routines from polyroots-fortran and those gave the same results (I might end up just replacing this custom routine with the external dependency).

Legacy Schmegacy

The earlier Fortran tool that some of this code was based on was the Goddard Trajectory Determination System (GTDS). This tool, like many Fortran tools, had a long heritage, but was thrown in the trash in the end, never modernized but replaced with C++ rewrites. Incidentally, GTDS and Swingby were developed by Computer Sciences Corporation (CSC), a company that existed from 1959 to 2017, which was cofounded by Roy Nutt, one of the members of the original team that created Fortran. So, now we have come full circle on this wheel of fire with a new modern Fortran port.

Conclusions

The end result of our AI-assisted C++ to Fortran conversion produced something that I think is likely quite useable. Overall I believe it did save me time compared to how long it would have taken me to do manually. But, getting this over the finish line definitely required manual intervention and domain knowledge. The process found a few potential bugs in the original code that need to be looked at by actual experts, since I don't trust the AI not a gaslight me with plausible-sounding reasons for them. Having the AI write the unit tests was also a huge time saver and help with debugging the implementation. In the end, for a task of this magnitude, I think AI can get you something like 95% of the way there, but then that last 5% needs to be completed and checked very carefully by somebody who more or less knows what they are doing. Blindly "vibe coding" is dangerous for technical applications such as this, since it's very easy for an AI to provide plausible looking procedures that are totally (or even worse, subtly) wrong.

Note: No AI was used to write the text of this article.

References

1. Jacchia, L. G., "New Static Models of the Thermosphere and Exosphere with Empirical Temperature Profiles", Smithsonian Astrophysical Observatory Special Report No. 313, 1970.
2. Roberts, C. E., Jr., "An Analytic Model for Upper Atmosphere Densities Based Upon Jacchia's 1970 Models", Celestial Mechanics, Vol. 4, pp. 368-377, 1971.
3. GMAT: General Mission Analysis Tool
4. Archive of legacy INPE atmosphere models github/@jacobwilliams
5. Atmosphere Models [degenerateconic.com] 2021-10-24
6. "Does the Programming Language Still Matter When AI Writes Your Code? I Tried Python and Came Back to C#" medium.com/@binnmti, Mar 16, 2026.
7. J. Carrico & E. Fletcher, "Software Architecture and Use of Satellite Tool Kit’s Astrogator Module for Libration Point Orbit Missions", Libration Point Orbits and Applications, pp. 471-487 (2003)
8. A. C. Long, et al., "Goddard Trajectory Determination System (GTDS) Mathematical Theory Revision 1", FDD/552-89/001, Computer Sciences Corporation, July 1989.