brainwagon

Occasionally I get to talk to hams who are just getting into using amateur satellites, and many of them ask the quite reasonable question “How do I figure out when the next pass occurs?”. For most of them, I suggest that they simply use a program like predict, which is probably what most people expect from a satellite prediction program.

I used just such a program for quite a while, until I first worked NH7WN in Hawaii, myself standing in my front yard, he standing on a beach in Hawaii, and I asked myself the quite normal question: “How often can I work Hawaii from here in California?” The (frankly excellent, don’t get me wrong) predict program didn’t provide a convenient way to answer that question.

Faced with this, I decided to write my own code for satellite prediction, so I could easily adapt it to answering questions like these. Because I like Python, and find it convenient for scripting applications, I chose Python as the implementation language. The algorithm I used was James Miller, G3RUH’s Plan 13 which was originally written in BASIC, and was relatively easy to port. I ended up creating a simple Python library implementation of Plan 13, which I’ve used for all my own satellite prediction needs. The nice thing about Python is that it has a bunch of useful libraries you can use (I used the sqlite3 library to build a database of orbital elements, and urllib to download them automatically from celestrak), and it’s been quite adaptable.

Still there were a few warts with it, and I haven’t had the chance to tidy it up. But recently I discovered that the PyEphem library includes a lot of the same functionality, and with a better version of the orbital model that I implemented. And, it’s pretty easy to use. As an example, here is a small chunk of code that figures out the next three passes of the ISS, using the current (for today) orbital elements. Even if you aren’t an amazing programmer, you can probably figure out how it works.

#!/usr/bin/python

import sys
import math
import ephem

iss = ephem.readtle("ISS (ZARYA)",
	"1 25544U 98067A   09270.78646569  .00012443  00000-0  87997-4 0  6860",
	"2 25544  51.6377 140.0905 0009007 135.9273 312.2213 15.74420558622113")

obs = ephem.Observer()
obs.lat = '38.0'
obs.long = '-122.0'

for p in range(3):
	tr, azr, tt, altt, ts, azs = obs.next_pass(iss)
	while tr < ts :
		obs.date = tr
		iss.compute(obs)
		print "%s %4.1f %5.1f" % (tr, math.degrees(iss.alt), math.degrees(iss.az))
		tr = ephem.Date(tr + 60.0 * ephem.second)
	print
	obs.date = tr + ephem.minute

Running it produces a dump of times, altitude and azimuths for three passes. Here’s a resulting run.

2009/9/28 04:12:18 -0.0 268.5
2009/9/28 04:13:18  1.4 257.6
2009/9/28 04:14:18  2.4 244.9
2009/9/28 04:15:18  2.5 231.3
2009/9/28 04:16:18  1.8 218.3
2009/9/28 04:17:18  0.4 206.7

2009/9/28 17:27:22 -0.1 169.1
2009/9/28 17:28:22  2.1 159.4
2009/9/28 17:29:22  4.0 147.1
2009/9/28 17:30:22  5.2 132.6
2009/9/28 17:31:22  5.1 117.1
2009/9/28 17:32:22  3.9 102.7
2009/9/28 17:33:22  1.9  90.7

2009/9/28 19:00:26  0.0 227.2
2009/9/28 19:01:26  3.8 226.9
2009/9/28 19:02:26  9.2 226.4
2009/9/28 19:03:26 17.8 225.3
2009/9/28 19:04:26 35.8 222.0
2009/9/28 19:05:26 80.4 152.7
2009/9/28 19:06:26 37.7  57.0
2009/9/28 19:07:26 18.5  53.3
2009/9/28 19:08:26  9.6  52.2
2009/9/28 19:09:26  4.1  51.8
2009/9/28 19:10:26  0.2  51.5

I did notice that this simple program may not work properly if the ISS is currently visible from your location, but I am sure that with just a little more work we could figure that out. I’ll probably port my existing software to this framework, because it also includes a bunch of useful functions (like figuring out when the sun rises and sets, and other fun things) which would be a pain for me to add. So here’s my radical suggestion: if you need to compute the location of satellites, write some scripts of your own using pyephem! It’s the homebrewer way. 🙂

Addendum: A few more minutes of work added some more interesting outputs…

#!/usr/bin/python

import sys
import math
import ephem

iss = ephem.readtle("ISS (ZARYA)",
	"1 25544U 98067A   09270.78646569  .00012443  00000-0  87997-4 0  6860",
	"2 25544  51.6377 140.0905 0009007 135.9273 312.2213 15.74420558622113")

obs = ephem.Observer()
obs.lat = '38.0'
obs.long = '-122.0'

for p in range(3):
	tr, azr, tt, altt, ts, azs = obs.next_pass(iss)
	print """Date/Time (UTC)       Alt/Azim	  Lat/Long	Elev"""
	print """====================================================="""
	while tr < ts:
		obs.date = tr
		iss.compute(obs)
		print "%s | %4.1f %5.1f | %4.1f %+6.1f | %5.1f" % \
			(tr, 
			 math.degrees(iss.alt), 
			 math.degrees(iss.az), 
			 math.degrees(iss.sublat), 
			 math.degrees(iss.sublong), 
			 iss.elevation/1000.)
		tr = ephem.Date(tr + 20.0 * ephem.second)
	print
	obs.date = tr + ephem.minute

And some output:

Date/Time (UTC)       Alt/Azim	  Lat/Long	Elev
=====================================================
2009/9/28 04:12:18 | -0.0 268.4 | 34.9 -145.6 | 335.6
2009/9/28 04:12:38 |  0.5 265.0 | 34.0 -144.5 | 335.7
2009/9/28 04:12:58 |  1.0 261.3 | 33.2 -143.4 | 335.9
2009/9/28 04:13:18 |  1.5 257.4 | 32.3 -142.3 | 336.0
2009/9/28 04:13:38 |  1.9 253.4 | 31.4 -141.3 | 336.2
2009/9/28 04:13:58 |  2.2 249.1 | 30.5 -140.2 | 336.3
2009/9/28 04:14:18 |  2.4 244.7 | 29.5 -139.2 | 336.5
2009/9/28 04:14:38 |  2.6 240.2 | 28.6 -138.3 | 336.6
2009/9/28 04:14:58 |  2.6 235.7 | 27.7 -137.3 | 336.8
2009/9/28 04:15:18 |  2.5 231.1 | 26.7 -136.3 | 337.0
2009/9/28 04:15:38 |  2.4 226.7 | 25.8 -135.4 | 337.1
2009/9/28 04:15:58 |  2.1 222.3 | 24.8 -134.5 | 337.3
2009/9/28 04:16:18 |  1.7 218.1 | 23.9 -133.6 | 337.5
2009/9/28 04:16:38 |  1.3 214.0 | 22.9 -132.7 | 337.7
2009/9/28 04:16:58 |  0.9 210.2 | 21.9 -131.8 | 337.9
2009/9/28 04:17:18 |  0.3 206.6 | 20.9 -131.0 | 338.1

Date/Time (UTC)       Alt/Azim	  Lat/Long	Elev
=====================================================
2009/9/28 17:27:24 |  0.0 168.7 | 19.0 -118.4 | 346.7
2009/9/28 17:27:44 |  0.7 165.7 | 20.0 -117.5 | 346.4
2009/9/28 17:28:04 |  1.4 162.5 | 21.0 -116.7 | 346.2
2009/9/28 17:28:24 |  2.1 158.9 | 22.0 -115.8 | 345.9
2009/9/28 17:28:44 |  2.8 155.1 | 22.9 -115.0 | 345.6
2009/9/28 17:29:04 |  3.5 151.0 | 23.9 -114.1 | 345.3
2009/9/28 17:29:24 |  4.1 146.6 | 24.9 -113.2 | 345.0
2009/9/28 17:29:44 |  4.6 141.9 | 25.8 -112.3 | 344.7
2009/9/28 17:30:04 |  4.9 137.0 | 26.8 -111.3 | 344.4
2009/9/28 17:30:24 |  5.2 131.9 | 27.7 -110.4 | 344.1
2009/9/28 17:30:44 |  5.3 126.8 | 28.6 -109.4 | 343.9
2009/9/28 17:31:04 |  5.3 121.6 | 29.6 -108.4 | 343.6
2009/9/28 17:31:24 |  5.1 116.4 | 30.5 -107.4 | 343.3
2009/9/28 17:31:44 |  4.8 111.5 | 31.4 -106.4 | 343.0
2009/9/28 17:32:04 |  4.3 106.7 | 32.3 -105.4 | 342.7
2009/9/28 17:32:24 |  3.8 102.1 | 33.2 -104.3 | 342.4
2009/9/28 17:32:44 |  3.2  97.9 | 34.0 -103.2 | 342.2
2009/9/28 17:33:04 |  2.5  93.9 | 34.9 -102.1 | 341.9
2009/9/28 17:33:24 |  1.8  90.2 | 35.8 -101.0 | 341.6
2009/9/28 17:33:44 |  1.1  86.9 | 36.6  -99.8 | 341.3
2009/9/28 17:34:04 |  0.4  83.8 | 37.4  -98.6 | 341.1

Date/Time (UTC)       Alt/Azim	  Lat/Long	Elev
=====================================================
2009/9/28 19:00:26 |  0.1 227.2 | 23.8 -137.4 | 345.3
2009/9/28 19:00:46 |  1.2 227.1 | 24.8 -136.5 | 345.1
2009/9/28 19:01:06 |  2.4 227.0 | 25.7 -135.6 | 344.8
2009/9/28 19:01:26 |  3.8 226.9 | 26.7 -134.7 | 344.5
2009/9/28 19:01:46 |  5.4 226.8 | 27.6 -133.7 | 344.2
2009/9/28 19:02:06 |  7.2 226.6 | 28.6 -132.8 | 343.9
2009/9/28 19:02:26 |  9.2 226.4 | 29.5 -131.8 | 343.6
2009/9/28 19:02:46 | 11.6 226.2 | 30.4 -130.8 | 343.3
2009/9/28 19:03:06 | 14.4 225.8 | 31.3 -129.7 | 343.0
2009/9/28 19:03:26 | 17.9 225.3 | 32.2 -128.7 | 342.8
2009/9/28 19:03:46 | 22.2 224.7 | 33.1 -127.6 | 342.5
2009/9/28 19:04:06 | 28.1 223.6 | 34.0 -126.5 | 342.2
2009/9/28 19:04:26 | 36.1 221.9 | 34.8 -125.4 | 341.9
2009/9/28 19:04:46 | 47.7 218.5 | 35.7 -124.3 | 341.6
2009/9/28 19:05:06 | 64.1 209.2 | 36.5 -123.1 | 341.4
2009/9/28 19:05:26 | 80.5 149.1 | 37.4 -121.9 | 341.1
2009/9/28 19:05:46 | 66.7  71.9 | 38.2 -120.7 | 340.8
2009/9/28 19:06:06 | 49.6  60.7 | 39.0 -119.5 | 340.6
2009/9/28 19:06:26 | 37.4  57.0 | 39.8 -118.2 | 340.3
2009/9/28 19:06:46 | 29.0  55.1 | 40.5 -116.9 | 340.1
2009/9/28 19:07:06 | 23.0  54.0 | 41.3 -115.5 | 339.8
2009/9/28 19:07:26 | 18.4  53.3 | 42.0 -114.1 | 339.6
2009/9/28 19:07:46 | 14.9  52.8 | 42.7 -112.7 | 339.3
2009/9/28 19:08:06 | 12.0  52.5 | 43.4 -111.3 | 339.1
2009/9/28 19:08:26 |  9.5  52.2 | 44.1 -109.8 | 338.8
2009/9/28 19:08:46 |  7.5  52.0 | 44.7 -108.3 | 338.6
2009/9/28 19:09:06 |  5.7  51.9 | 45.4 -106.7 | 338.4
2009/9/28 19:09:26 |  4.1  51.7 | 46.0 -105.1 | 338.1
2009/9/28 19:09:46 |  2.6  51.7 | 46.5 -103.5 | 337.9
2009/9/28 19:10:06 |  1.3  51.6 | 47.1 -101.9 | 337.7
2009/9/28 19:10:26 |  0.2  51.5 | 47.6 -100.2 | 337.5

Addendum²: As another example of what can be done, consider the command line voice synthesizer that is included in Mac OS. Using the os.system command, you can have your computer speak voice updates almost as easy as printing them. For example:

#!/usr/bin/python

# Mac OS includes a voice synthesis command called "say".   It's pretty
# easy to make a simple program like this that will use voice to announce
# the current position of the satellite.   With a little work, you could
# easily make it announce the altitude and azimuth repeatedly during a pass.

import sys
import os
import math
import ephem

iss = ephem.readtle("ISS (ZARYA)",
	"1 25544U 98067A   09270.78646569  .00012443  00000-0  87997-4 0  6860",
	"2 25544  51.6377 140.0905 0009007 135.9273 312.2213 15.74420558622113")

obs = ephem.Observer()
obs.lat = '38.0'
obs.long = '-122.0'

iss.compute(obs)
if iss.alt < 0:
	os.system('say The ISS is not currently visible.')
else:
	os.system('say "The ISS is at altitude %d, azimuth %d."' % \
		(int(math.degrees(iss.alt)+0.5), int(math.degrees(iss.az)+0.5)))

Here’s an example of the voice it outputs:
Voice synthesis example