A few years ago, I created my own Python implementation of the Plan13 satellite prediction code written by James Miller (G3RUH). The Plan13 algorithm isn’t very complicated: you can easily run it on processors like the Arduino (in fact, I used it for my ANGST satellite tracker) But somehow, I managed to misplace the source code for the Python version, probably on a hard drive for a computer that died, and so when I wanted to do some ISS calculations, I decided that I’d go ahead and use PyEphem, a much more complete package that can do all sorts of astronomical calculations, and which includes an implementation of the fairly standard SGP4 model.
It’s fairly simple to write some code to predict ISS passes from any position on earth. Here’s a quick example, with the TLEs for the ISS hard coded, as well as an observer position:
[sourcecode lang=”python”]
#!/usr/bin/env python
from math import degrees
import ephem
# create an observer
obs = ephem.Observer()
obs.lat = ‘38.0’
obs.lon = ‘-122’
obs.elevation = 100.
# create the iss object…
iss = ephem.readtle("ISS (ZARYA) ",
"1 25544U 98067A 15067.48441091 .00017347 00000-0 26069-3 0 9998",
"2 25544 51.6448 227.0292 0008846 83.1031 17.2346 15.54929647932364")
# figure out the next pass…
# only interested in rt (rise time), tt (transit time) and st (set time)
#
rt, razi, tt, televation, st, sazi = obs.next_pass(iss)
print "Rise Time: ", rt
print "Transit Time:", tt
print "Set Time: ", st
[/sourcecode]
But I was somewhat chagrined that if I ran this little snippet multiple times, I got different answers. Not by a lot, but by several seconds:
pi@blueberrypi ~ $ !. ./iss Rise Time: 2015/3/8 22:30:12 Transit Time: 2015/3/8 21:00:52 Set Time: 2015/3/8 21:04:09 pi@blueberrypi ~ $ !! ./iss Rise Time: 2015/3/8 22:30:12 Transit Time: 2015/3/8 21:00:53 Set Time: 2015/3/8 21:04:10 pi@blueberrypi ~ $ !! ./iss Rise Time: 2015/3/8 22:30:13 Transit Time: 2015/3/8 21:00:53 Set Time: 2015/3/8 21:04:10 pi@blueberrypi ~ $ !! ./iss Rise Time: 2015/3/8 22:30:03 Transit Time: 2015/3/8 21:00:49 Set Time: 2015/3/8 21:04:11
Digging in the code, it’s not hard to see what’s going on. The code is stepping forward by intervals of about 1 minute, trying to catch the satellite as it peaks over the horizon. When it does, it then uses twenty iterations of the secant method to find the point where the satellite crosses zero altitude. But there is an early out: if the interval drops below ten seconds, it goes ahead and exists. Similarly, it only locates the time of the transit to within 15 seconds. But these can be fixed with some quick, and probably undisciplined changes to the code will make things behave better.
Both edits are in the file riset_cir.c which is part of libastro-3.7.5 (your version might change, but will probably be the same). Near the top you’ll find the declaration for TMACC
[sourcecode lang=”C”]
#define TMACC (10./3600./24.0) /* convergence accuracy, days */
[/sourcecode]
Change this to 0.5 seconds instead.
[sourcecode lang=”C”]
#define TMACC (0.5/3600./24.0) /* convergence accuracy, days */
[/sourcecode]
You might want to increase MAXPASSES in the find_0alt function, but I had no difficulty with leaving it at 20.
Similarly, we need to change the error in find_transit:
[sourcecode lang=”C”]
static int
find_transit (double dt, Now *np, Obj *op)
{
#define MAXLOOPS 10
#define MAXERR (0.25/60.) /* hours */
[/sourcecode]
The number of iterations is pretty low (just 10) and the MAXERR was 15 seconds. I decided to add some iterations to ensure convergence, and to set the MAXERR to be just one second.
[sourcecode lang=”C”]
static int
find_transit (double dt, Now *np, Obj *op)
{
#define MAXLOOPS 20
#define MAXERR (1./3600.) /* hours */
[/sourcecode]
With these changes in place, you get times which are accurate to around one second. Note: I do not mean to imply that these numbers are somehow absolutely better. It is known that the SGP4 model can only locate satellites to errors of around 1km per day at best, so trying to converge them to insane levels of accuracy isn’t meaningful. But for my application, I want to count down to transits, and having the estimate jump around by a few seconds made the count down look funny. Of course, I could have just computed the times once, and used that, but I still think it’s kind of odd that the estimates vary by a human discernable amount depending solely on what time you decide to try to compute from. This additional accuracy makes the program a very tiny bit slower, but is worth it for my application.
Addendum: A few more lines of code make the code a bit more useful. This now outputs the next pass expressed in the localtime rather than UTC, and gives you a bit of a countdown.
[sourcecode lang=”python”]
#!/usr/bin/env python
from math import degrees
import ephem
import datetime
# create an observer
obs = ephem.Observer()
obs.lat = ‘38.0’
obs.lon = ‘-122’
obs.elevation = 100.
# create the iss object…
iss = ephem.readtle("ISS (ZARYA) ",
"1 25544U 98067A 15067.48441091 .00017347 00000-0 26069-3 0 9998",
"2 25544 51.6448 227.0292 0008846 83.1031 17.2346 15.54929647932364")
# figure out the next pass…
# only interested in rt (rise time), tt (transit time) and st (set time)
#
n = ephem.now()
rt, razi, tt, televation, st, sazi = obs.next_pass(iss)
print "%d" % round((rt – n) * 3600 * 24), "seconds until the next pass."
rt = ephem.localtime(rt).strftime("%x %X")
tt = ephem.localtime(tt).strftime("%x %X")
st = ephem.localtime(st).strftime("%x %X")
print "Rise Time: ", rt, "Azimuth: %5.1f" % degrees(razi)
print "Transit Time:", tt, "Elevation: %5.1f" % degrees(televation)
print "Set Time: ", st, "Azimuth: %5.1f" % degrees(sazi)
[/sourcecode]
Addendum2: A few more lines of code fetches a current set of orbital elements from celestrak.com. A future revision of this will probably cache a local copy of these orbital elements so you don’t unnecessarily hammer their server.
[sourcecode lang=”python”]
#!/usr/bin/env python
from math import degrees
import ephem
import datetime
import urllib2
req = urllib2.Request("http://celestrak.com/NORAD/elements/amateur.txt")
response = urllib2.urlopen(req)
data = response.read().splitlines()
objs = []
for idx in range(0, len(data), 3):
objs.append(ephem.readtle(data[idx], data[idx+1], data[idx+2]))
# create an observer
obs = ephem.Observer()
obs.lat = ‘38.0’
obs.lon = ‘-122’
obs.elevation = 100.
# create the iss object…
iss = ephem.readtle("ISS (ZARYA) ",
"1 25544U 98067A 15067.48441091 .00017347 00000-0 26069-3 0 9998",
"2 25544 51.6448 227.0292 0008846 83.1031 17.2346 15.54929647932364")
# figure out the next pass…
# only interested in rt (rise time), tt (transit time) and st (set time)
#
tlist = []
for obj in objs:
r = obs.next_pass(obj)
tlist.append((r[0], obj, r))
tlist.sort()
for t, obj, r in tlist:
n = ephem.now()
rt, razi, tt, televation, st, sazi = r
print obj.name
print "%d" % round((rt-n) * 3600 * 24), "seconds until the next pass."
rt = ephem.localtime(rt).strftime("%x %X")
tt = ephem.localtime(tt).strftime("%x %X")
st = ephem.localtime(st).strftime("%x %X")
print "Rise Time: ", rt, "Azimuth: %5.1f" % degrees(razi)
print "Transit Time:", tt, "Elevation: %5.1f" % degrees(televation)
print "Set Time: ", st, "Azimuth: %5.1f" % degrees(sazi)
print
[/sourcecode]
The output is now a list of all the amateur satellites that it fetched, sorted in order of the passes which are soonest appear first.
BEESAT-2 410 seconds until the next pass. Rise Time: 03/08/15 18:33:03 Azimuth: 272.2 Transit Time: 03/08/15 18:37:54 Elevation: 8.1 Set Time: 03/08/15 18:42:45 Azimuth: 9.4 CUBESAT XI-V (CO-58) 526 seconds until the next pass. Rise Time: 03/08/15 18:34:59 Azimuth: 31.0 Transit Time: 03/08/15 18:41:08 Elevation: 16.4 Set Time: 03/08/15 18:47:16 Azimuth: 153.1 XIWANG-1 (HOPE-1) 1228 seconds until the next pass. Rise Time: 03/08/15 18:46:41 Azimuth: 91.6 Transit Time: 03/08/15 18:53:13 Elevation: 8.5 Set Time: 03/08/15 18:59:45 Azimuth: 6.5 OSCAR 7 (AO-7) 1251 seconds until the next pass. Rise Time: 03/08/15 18:47:04 Azimuth: 229.2 Transit Time: 03/08/15 18:54:42 Elevation: 8.5 Set Time: 03/08/15 19:02:20 Azimuth: 317.9 AAUSAT3 2080 seconds until the next pass. Rise Time: 03/08/15 19:00:54 Azimuth: 29.7 Transit Time: 03/08/15 19:07:39 Elevation: 19.3 Set Time: 03/08/15 19:14:25 Azimuth: 157.0 PCSAT (NO-44) 2515 seconds until the next pass. Rise Time: 03/08/15 19:08:08 Azimuth: 173.2 Transit Time: 03/08/15 19:15:24 Elevation: 24.3 Set Time: 03/08/15 19:22:39 Azimuth: 41.1 CUTE-1.7+APD II (CO-65) 2940 seconds until the next pass. Rise Time: 03/08/15 19:15:14 Azimuth: 81.8 Transit Time: 03/08/15 19:17:56 Elevation: 1.9 Set Time: 03/08/15 19:20:37 Azimuth: 31.8 JAS-2 (FO-29) 2956 seconds until the next pass. Rise Time: 03/08/15 19:15:30 Azimuth: 255.5 Transit Time: 03/08/15 19:19:08 Elevation: 2.2 Set Time: 03/08/15 19:22:47 Azimuth: 307.5 M-CUBED & EXP-1 PRIME 3635 seconds until the next pass. Rise Time: 03/08/15 19:26:49 Azimuth: 354.1 Transit Time: 03/08/15 19:32:20 Elevation: 12.4 Set Time: 03/08/15 19:37:52 Azimuth: 247.0 MOZHAYETS 4 (RS-22) 3940 seconds until the next pass. Rise Time: 03/08/15 19:31:54 Azimuth: 347.2 Transit Time: 03/08/15 19:37:10 Elevation: 9.9 Set Time: 03/08/15 19:42:26 Azimuth: 248.6 CUBESAT XI-IV (CO-57) 4134 seconds until the next pass. Rise Time: 03/08/15 19:35:08 Azimuth: 194.6 Transit Time: 03/08/15 19:42:31 Elevation: 27.0 Set Time: 03/08/15 19:49:53 Azimuth: 334.6 HAMSAT (VO-52) 4222 seconds until the next pass. Rise Time: 03/08/15 19:36:36 Azimuth: 235.6 Transit Time: 03/08/15 19:40:29 Elevation: 4.5 Set Time: 03/08/15 19:44:22 Azimuth: 310.9 PRISM (HITOMI) 4468 seconds until the next pass. Rise Time: 03/08/15 19:40:41 Azimuth: 13.2 Transit Time: 03/08/15 19:47:14 Elevation: 82.2 Set Time: 03/08/15 19:53:46 Azimuth: 192.3 CUTE-1 (CO-55) 4562 seconds until the next pass. Rise Time: 03/08/15 19:42:16 Azimuth: 198.2 Transit Time: 03/08/15 19:49:30 Elevation: 23.6 Set Time: 03/08/15 19:56:44 Azimuth: 332.9 STRAND-1 4714 seconds until the next pass. Rise Time: 03/08/15 19:44:47 Azimuth: 18.5 Transit Time: 03/08/15 19:52:18 Elevation: 49.9 Set Time: 03/08/15 19:59:48 Azimuth: 181.1 BEESAT-3 4722 seconds until the next pass. Rise Time: 03/08/15 19:44:55 Azimuth: 342.9 Transit Time: 03/08/15 19:44:56 Elevation: -0.0 Set Time: 03/08/15 19:44:56 Azimuth: 342.9 YUBILEINY (RS-30) 5115 seconds until the next pass. Rise Time: 03/08/15 19:51:28 Azimuth: 279.5 Transit Time: 03/08/15 19:57:25 Elevation: 3.7 Set Time: 03/08/15 20:03:21 Azimuth: 342.8 TRITON-1 5411 seconds until the next pass. Rise Time: 03/08/15 19:56:24 Azimuth: 92.7 Transit Time: 03/08/15 20:00:09 Elevation: 3.9 Set Time: 03/08/15 20:03:54 Azimuth: 25.1 UOSAT 2 (UO-11) 6509 seconds until the next pass. Rise Time: 03/08/15 20:14:43 Azimuth: 95.0 Transit Time: 03/08/15 20:18:33 Elevation: 4.5 Set Time: 03/08/15 20:22:23 Azimuth: 22.8 SEEDS II (CO-66) 6557 seconds until the next pass. Rise Time: 03/08/15 20:15:30 Azimuth: 127.3 Transit Time: 03/08/15 20:21:07 Elevation: 16.4 Set Time: 03/08/15 20:26:43 Azimuth: 6.7 UWE-3 6694 seconds until the next pass. Rise Time: 03/08/15 20:17:48 Azimuth: 76.8 Transit Time: 03/08/15 20:20:03 Elevation: 1.3 Set Time: 03/08/15 20:22:18 Azimuth: 35.9 RADIO ROSTO (RS-15) 7326 seconds until the next pass. Rise Time: 03/08/15 20:28:20 Azimuth: 128.2 Transit Time: 03/08/15 20:36:14 Elevation: 6.0 Set Time: 03/08/15 20:44:09 Azimuth: 56.8 GOMX 1 7721 seconds until the next pass. Rise Time: 03/08/15 20:34:55 Azimuth: 139.6 Transit Time: 03/08/15 20:41:44 Elevation: 28.3 Set Time: 03/08/15 20:48:33 Azimuth: 0.7 CUBEBUG-2 (LO-74) 9556 seconds until the next pass. Rise Time: 03/08/15 21:05:29 Azimuth: 126.2 Transit Time: 03/08/15 21:11:12 Elevation: 15.9 Set Time: 03/08/15 21:16:56 Azimuth: 7.5 VELOX-PII 9914 seconds until the next pass. Rise Time: 03/08/15 21:11:28 Azimuth: 125.2 Transit Time: 03/08/15 21:17:05 Elevation: 15.3 Set Time: 03/08/15 21:22:43 Azimuth: 8.0 FUNCUBE-1 (AO-73) 9936 seconds until the next pass. Rise Time: 03/08/15 21:11:49 Azimuth: 118.8 Transit Time: 03/08/15 21:17:04 Elevation: 11.9 Set Time: 03/08/15 21:22:20 Azimuth: 11.1 ZACUBE-1 (TSHEPISOSAT) 10153 seconds until the next pass. Rise Time: 03/08/15 21:15:26 Azimuth: 121.9 Transit Time: 03/08/15 21:20:52 Elevation: 13.4 Set Time: 03/08/15 21:26:17 Azimuth: 9.6 DELFI-C3 (DO-64) 11107 seconds until the next pass. Rise Time: 03/08/15 21:31:21 Azimuth: 142.9 Transit Time: 03/08/15 21:37:20 Elevation: 29.2 Set Time: 03/08/15 21:43:19 Azimuth: 359.6 HUMSAT-D 12799 seconds until the next pass. Rise Time: 03/08/15 21:59:32 Azimuth: 133.5 Transit Time: 03/08/15 22:05:21 Elevation: 20.4 Set Time: 03/08/15 22:11:09 Azimuth: 4.0 EAGLE 2 13942 seconds until the next pass. Rise Time: 03/08/15 22:18:36 Azimuth: 136.2 Transit Time: 03/08/15 22:24:17 Elevation: 22.0 Set Time: 03/08/15 22:29:58 Azimuth: 3.0 CUBEBUG-1 (CAPITAN BETO) 15152 seconds until the next pass. Rise Time: 03/08/15 22:38:45 Azimuth: 129.1 Transit Time: 03/08/15 22:44:45 Elevation: 18.7 Set Time: 03/08/15 22:50:44 Azimuth: 4.6 SPROUT 15611 seconds until the next pass. Rise Time: 03/08/15 22:46:24 Azimuth: 92.6 Transit Time: 03/08/15 22:50:05 Elevation: 4.0 Set Time: 03/08/15 22:53:46 Azimuth: 24.0 JUGNU 16144 seconds until the next pass. Rise Time: 03/08/15 22:55:18 Azimuth: 204.1 Transit Time: 03/08/15 23:00:40 Elevation: 6.1 Set Time: 03/08/15 23:06:02 Azimuth: 124.2 DUCHIFAT-1 16575 seconds until the next pass. Rise Time: 03/08/15 23:02:29 Azimuth: 35.0 Transit Time: 03/08/15 23:07:53 Elevation: 12.7 Set Time: 03/08/15 23:13:17 Azimuth: 148.7 SRMSAT 17209 seconds until the next pass. Rise Time: 03/08/15 23:13:03 Azimuth: 176.4 Transit Time: 03/08/15 23:15:40 Elevation: 1.1 Set Time: 03/08/15 23:18:17 Azimuth: 139.8 SWISSCUBE 18892 seconds until the next pass. Rise Time: 03/08/15 23:41:05 Azimuth: 73.7 Transit Time: 03/08/15 23:43:37 Elevation: 1.4 Set Time: 03/08/15 23:46:09 Azimuth: 31.4 SAUDISAT 1C (SO-50) 19177 seconds until the next pass. Rise Time: 03/08/15 23:45:51 Azimuth: 170.5 Transit Time: 03/08/15 23:52:00 Elevation: 17.7 Set Time: 03/08/15 23:58:10 Azimuth: 47.9 BEESAT 19303 seconds until the next pass. Rise Time: 03/08/15 23:47:57 Azimuth: 77.8 Transit Time: 03/08/15 23:50:51 Elevation: 2.0 Set Time: 03/08/15 23:53:45 Azimuth: 28.8 TISAT 1 19890 seconds until the next pass. Rise Time: 03/08/15 23:57:44 Azimuth: 141.2 Transit Time: 03/09/15 00:03:57 Elevation: 28.7 Set Time: 03/09/15 00:10:10 Azimuth: 359.5 ITUPSAT 1 21056 seconds until the next pass. Rise Time: 03/09/15 00:17:09 Azimuth: 104.3 Transit Time: 03/09/15 00:22:11 Elevation: 8.1 Set Time: 03/09/15 00:27:12 Azimuth: 14.2 KKS-1 (KISEKI) 24703 seconds until the next pass. Rise Time: 03/09/15 01:17:56 Azimuth: 96.4 Transit Time: 03/09/15 01:22:04 Elevation: 5.2 Set Time: 03/09/15 01:26:12 Azimuth: 20.4 LUSAT (LO-19) 25871 seconds until the next pass. Rise Time: 03/09/15 01:37:24 Azimuth: 29.2 Transit Time: 03/09/15 01:44:13 Elevation: 20.0 Set Time: 03/09/15 01:51:01 Azimuth: 157.9 ITAMSAT (IO-26) 26776 seconds until the next pass. Rise Time: 03/09/15 01:52:30 Azimuth: 42.8 Transit Time: 03/09/15 01:57:50 Elevation: 7.6 Set Time: 03/09/15 02:03:10 Azimuth: 133.0 TECHSAT 1B (GO-32) 27579 seconds until the next pass. Rise Time: 03/09/15 02:05:52 Azimuth: 41.3 Transit Time: 03/09/15 02:11:26 Elevation: 8.3 Set Time: 03/09/15 02:17:00 Azimuth: 134.4 HORYU 2 31680 seconds until the next pass. Rise Time: 03/09/15 03:14:13 Azimuth: 41.7 Transit Time: 03/09/15 03:19:13 Elevation: 8.3 Set Time: 03/09/15 03:24:13 Azimuth: 137.6 ISS (ZARYA) 38130 seconds until the next pass. Rise Time: 03/09/15 05:01:43 Azimuth: 115.2 Transit Time: 03/09/15 04:15:25 Elevation: -74.9 Set Time: 03/09/15 05:01:43 Azimuth: 115.2
Addendum3: Hmm. If you look carefully, the last entry in the test output is bogus: the time of maximum elevation does not fall between the rise and set time, and the elevation is negative. Staring at the code some more, I’m now worried about the logic of the root finder. I’ll have to ponder it some more.
Hi Mark, I’d like to reproduce this for study but I’m not having any luck trying different times near the date of your post. Please add printing the value of ephem.now so the next time there is a bogus result we can look into it. Thanks.