Facebook Linkedin Twitter
Posted Sun Dec 5, 2021 •  Reading time: 7 minutes

Calculating Satellite Orbits with SGP4 and Go

Every Christmas, I catch myself looking up at the night sky, trying to spot a bearded man pulled by a crew of flying deer. Most years, I’m disappointed and resort to soft drink advertising for my Christmas miracles. What I do see, however, are little lights streaking across the dark background of the endless void that is our universe. A firefly? No, its path is too straightforward! A plane? Too far away! A radioactive bird? I certainly hope not!

What I see — and what you, too, may spot when you look up at night — are satellites. According to the Union of Concerned Scientists, there are 4,550 artificial satellites orbiting the Earth at the time of writing (and likely a few more once you read this!). These range from large objects as the well-known International Space Station to tiny little devices as small as a carton of milk.

With so many things out there, how can we make sure we don’t hit anything when we launch new satellites? Hypothetically, how would someone tasked to fly from the North Pole to millions of children around the Earth within 24 hours plan their route without crashing into one of the thousands of objects that orbit our planet?

Thankfully, the North American Aerospace Defense Command (NORAD) keeps track of everything for us in their CelesTrak database. Don’t you love a good database?!

For every object that orbits the Earth1, NORAD provides us with a two-line element set (or TLE) that describes the object’s orbital parameters. For example, this is what the TLE for the International Space Station Zarya module looks like:

1 25544U 98067A   21332.52730935  .00007275  00000-0  14154-3 0  9994
2 25544  51.6436 250.1268 0004444 253.3057 278.3644 15.48699150314068

The two-row, 69-column ASCII format was originally intended for punch cards! It’s not exactly protocol buffers, is it? But it’s also not rocket science2: Wikipedia has a detailed explainer on what each of the elements mean.

There’s a few bits of identifier information, then you have a bit of information about the epoch, which is the timestamp of your data. What we’re interested in is all the orbital information.

You might ask yourself why we have all these different numbers here. In high school physics, many of us are taught about Kepler’s laws of planetary motion, where one body orbits another in an elliptical orbit, which is relatively simple. As always, real life is a bit more complex than that, and we learn that our high school physics doesn’t cover all the details.

For one, the Sun and Moon have some gravitational pull on satellites as well, which we must take into account. Then, Earth’s atmosphere can cause some drag, decreasing satellite speed. We also have effects of radiation! And finally, the Earth is not actually a sphere, but a flat, thin disc.3 Just kidding, it’s a sphere, but not a perfect one. We call the effects of these factors perturbations.

And all of these perturbations make calculating satellite orbits by hand a bit difficult. Thankfully, we have computers! In 1988, we got the SGP4 simplified perturbations model, originally implemented for FORTRAN IV. The model allows taking our orbital parameters and calculating a near-exact4 position of our satellite.

As you may have noticed, FORTRAN IV is no longer the language of choice, so over time some other implementations have appeared. For our purposes, we will turn to go-satellite, a Go port of the official C++ implementation.5

Let’s find out how to use it:

go get github.com/joshuaferrara/go-satellite
touch main.go

First off, we will want to load our TLE data. The library exposes two functions for this, satellite.ParseTLE and satellite.TLEToSat. The first one only parses data and returns a satellite.Satellite type struct, while the second actually runs satellite.sgp4 on that data, initializing a bunch of stuff about that satellite.6 Let’s see what we get with our ISS TLE:

// create a new satellite from our TLE data
s := satellite.TLEToSat(
    "1 25544U 98067A   21332.52730935  .00007275  00000-0  14154-3 0  9994",
    "2 25544  51.6436 250.1268 0004444 253.3057 278.3644 15.48699150314068",

fmt.Printf("%+v\n", s)
    Line1:1 25544U 98067A   21332.52730935  .00007275  00000-0  14154-3 0  9994
    Line2:2 25544  51.6436 250.1268 0004444 253.3057 278.3644 15.48699150314068
    whichconst: {

So much data! You are probably wondering what wgs84 is about: For SGP4, we have to specificy a geodetic system, think of it like a reference frame for Earth. Remember that Earth is not a perfect sphere! Our little blue planet has a lot of quirks and features that we must take into account. WGS 84 is the most recent model, from 1984, and it’s also used by GPS. Alternatives in SGP4 are WGS 72 and WGS 72 Old, but who knows what they are.7

Anyway, in our output (that I have kindly formatted for you), we have tons of information about our satellite. Much of that information is directly copied (e.g., bstar:0.00014154, which is 14154-3 in line 1, columns 55-61 in our TLE and describes the B* Drag Coefficient), while other is calculated based on our input.

With our satellite data in hand, let’s see where it is right now!

// get the current date and time
t := time.Now().UTC()

// propagate the satellite position
pos, vec := satellite.Propagate(s, t.Year(), int(t.Month()), t.Day(), t.Hour(), t.Minute(), t.Second())

fmt.Printf("Position: %+v\n", pos)
fmt.Printf("Velocity Vector: %+v\n", vec)
Position: {X:2928.7746627666966 Y:6134.6035243137985 Z:-111.53131173902264}
Velocity Vector: {X:-4.329668312607247 Y:1.9613689132693668 Z:-6.007449872040121}

Hm, that’s not very helpful! To propagate the satellite position, we have to give our propagator a date and time. We’re using the current date and time, but you could use anything: Maybe you want to know where a satellite is at midnight! Whatever you do, make sure to use time in UTC!

But again, all we get are some numbers that we can’t really put into perspective. Our position is actually not longitude and latitude that we’re used to: We have what are called Earth centered inertial (ECI) coordinates, meaning we have X, Y, and Z coordinates around the planet that don’t rotate with the planet. Thankfully, we have the satellite.ECIToLLA function that allows us to convert ECI to latitude, longitude, and altitude:

// calculate the GST date
gst := satellite.GSTimeFromDate(t.Year(), int(t.Month()), t.Day(), t.Hour(), t.Minute(), t.Second())

// convert Earth Centered Inertial coordinates to Lat/Long
altitude, velocity, ret := satellite.ECIToLLA(pos, gst)

fmt.Printf("Altitude: %+v\n", altitude)
fmt.Printf("Velocity: %+v\n", velocity)
fmt.Printf("Lat/Long: %+v\n", ret)
Altitude: 424.089617385127
Velocity: 7.6549672964590245
Lat/Long: {Latitude:-0.35366624503253075 Longitude:1.309783961510942}

To do the calculation, we have to account for Earth rotation, which ECI does not. And in order to do that, we need to take our time t again and convert it into Greenwich Sidereal Time (GST8). At this point, my understanding of the whole timezone shenanigans ends, but thankfully we have the satellite.GSTimeFromDate function that does all of that for us.

Next, we put our position vector and GST date into satellite.ECIToLLA and receive an altitude of 424km and a latitude/longitude pair. A quick check reveals that the ISS is indeed located at a 424km altitude, so it seems that everything worked out great! It appears to be located somewhere off the West African coast, so no chance to see it at the moment…

Or is it? Well, the whole naming thing can be a bit confusing, but what we got were latitude and longitude in radians, not in degrees as we are used to. To get proper coordinates, we first need to convert our data using satellite.LatLongDeg9:

// convert our Lat and Long into degrees
deg := satellite.LatLongDeg(ret)

fmt.Printf("Lat/Long: %+v\n", deg)
{Latitude:-20.263583196603626 Longitude:75.04509306850244}

There we go, right above the Indian Ocean! But if the ISS is over the ocean, what did I see flying across the night sky…

So, a quick recap of what we have achieved today: We learned how to access publicly available TLE data and what that even means! We saw how we could parse that TLE data with Go and how we to use that to calculate satellite positions! Beyond that, the go-satellite library offers a few other useful features, such as a downloader for CelesTrak or a function to calculate the angle between an observer on Earth and a satellite, so you know exactly in which direction to look.

With your newly acquired knowledge, what will you build next? A quick and efficient web service for satellite locations? A terminal app? A Kubernetes plugin?10

P.S: As a final note before we part, let me share one of my favorite things on the topic of the North American Aerospace Defense Command: Every year around Christmas, NORAD tracks the position of one of our most valuable assets: https://www.noradsanta.org/

  1. Well, at least for every object that isn’t top secret Area 51 stuff! The U indicator in line 1, column 8 actually indicates that the object is Unclassified. But you will probably not come across any classified objects here. ↩︎

  2. Okay, technically… ↩︎

  3. For the record: GitHub Copilot wrote this joke! ↩︎

  4. In this case, there is a ~1km error at epoch, growing about 1-3km per day. That’s why TLE data is often updated! ↩︎

  5. As a side note, can we appreciate that the author Joshua Ferrara claims that they ported the C++ code to Go as one of their first Go projects?! Meanwhile, I still don’t fully understand what a context is supposed to be after two years… ↩︎

  6. You will notice that this library really is a port, and not a rewrite: The API is not exactly written in a style that is familiar to Go devs. But it’s a good start! ↩︎

  7. Again, Wikipedia has a great overview ↩︎

  8. Not to be confused with Gulf Standard Time! ↩︎

  9. Or whatever your favorite radians-to-degrees converter is, you do you! ↩︎

  10. If you come up with a Kubernetes plugin that uses SGP4, please reach out to me! ↩︎