Category Archives: Python

One dimensional cellular automata code in Python..

Inspired by the KnitYak Kickstarter, I thought I would code up a simple Python program that could generate the same sort of patterns that are used in the scarves in the video. If you want to learn more about the mathematical theory of these cellular automata, google for keywords like “1D” “cellular automata” and “steve wolfram”. Without much further explanation, here is the code:

[sourcecode lang=”python”]
#!/usr/bin/env python

# __
# ___ __ __ _ _ _ / _|
# (_-</ _/ _` | ‘_| _|
# /__/\__\__,_|_| |_|
#
# Inspired by the KnitYak Kickstarter, which offers
# mathematical scarves whose patterns are algorithmically
# generated, I thought it would be good to dust off my old
# knowledge and generate the same sort of patterns using
# Python. The patterns are simply printed using ASCII,
# but the code could be adapted to output in other forms
# easily.
#
# Written by Mark VandeWettering
#

import sys
import random

WIDTH=79 # How wide is the pattern?

w = WIDTH * [0] # create the current generation
nw = WIDTH * [0] # and the next generation
w[WIDTH/2] = 1 # populate with a single one

# or alternatively, you can populate it with a random
# initial configuration. If you want to start with
# just a single one, comment the next two lines out.

#for i in range(WIDTH):
# w[i] = random.randint(0, 1)

# How wide is the neighborhood of cells that are
# examined? The traditional Wolfram 1D cellular
# automata uses a neighborhood of 3…

NEIGHBORHOOD=3

# rtab is space for the rule table. It maps all
# numbers from [0, 2**NEIGHBORHOOD) to either a 0 or 1.
rtab = (2**NEIGHBORHOOD) * [0]

# The "rule" is a number which is used to populate
# rtab. The number is in the range [0, 2**(2**NEIGHBORHOOD))
# Many rules generate uninteresting patterns, but some
# like 30 generate interesting, aperiodic patterns.
rule = 30

# This fills in the table…
for i in range(2**NEIGHBORHOOD):
if ((2**i) & rule) != 0:
rtab[i] = 1

def dump(r):
for x in r:
if x == 1:
sys.stdout.write(‘X’)
else:
sys.stdout.write(‘ ‘)
sys.stdout.write(‘\n’)

# and generates 100 lines…

for y in range(100):
dump(w)
for x in range(WIDTH):
sum = 0
for d in range(NEIGHBORHOOD):
sum = sum + (2**d) * w[(x+d+WIDTH – NEIGHBORHOOD/2) % WIDTH]
nw[x] = rtab[sum]
w, nw = nw, w
[/sourcecode]

The pattern this generates with these settings follows. Try adjusting rule to other numbers up to 255 to generate other patterns. Lots will be boring, but some will be pretty interesting.

                                                                               
                                       X                                       
                                      XXX                                      
                                     X  XX                                     
                                    XXXX XX                                    
                                   X   X  XX                                   
                                  XXX XXXX XX                                  
                                 X  X    X  XX                                 
                                XXXXXX  XXXX XX                                
                               X     XXX   X  XX                               
                              XXX   X  XX XXXX XX                              
                             X  XX XXXX X    X  XX                             
                            XXXX X    X XX  XXXX XX                            
                           X   X XX  XX  XXX   X  XX                           
                          XXX XX  XXX XXX  XX XXXX XX                          
                         X  X  XXX  X   XXX X    X  XX                         
                        XXXXXXX  XXXXX X  X XX  XXXX XX                        
                       X      XXX    X XXXX  XXX   X  XX                       
                      XXX    X  XX  XX    XXX  XX XXXX XX                      
                     X  XX  XXXX XXX XX  X  XXX X    X  XX                     
                    XXXX XXX   X   X  XXXXXX  X XX  XXXX XX                    
                   X   X   XX XXX XXXX     XXXX  XXX   X  XX                   
                  XXX XXX X X   X    XX   X   XXX  XX XXXX XX                  
                 X  X   X X XX XXX  X XX XXX X  XXX X    X  XX                 
                XXXXXX XX X  X   XXXX  X   X XXX  X XX  XXXX XX                
               X     X  X XXXXX X   XXXXX XX   XXXX  XXX   X  XX               
              XXX   XXXXX     X XX X    X  XX X   XXX  XX XXXX XX              
             X  XX X    XX   XX  X XX  XXXX X XX X  XXX X    X  XX             
            XXXX X XX  X XX X XXXX  XXX   X X  X XXX  X XX  XXXX XX            
           X   X X  XXXX  X X    XXX  XX XX XXXX   XXXX  XXX   X  XX           
          XXX XX XXX   XXXX XX  X  XXX X  X    XX X   XXX  XX XXXX XX          
         X  X  X   XX X   X  XXXXXX  X XXXXX  X X XX X  XXX X    X  XX         
        XXXXXXXXX X X XX XXXX     XXXX     XXXX X  X XXX  X XX  XXXX XX        
       X        X X X  X    XX   X   XX   X   X XXXX   XXXX  XXX   X  XX       
      XXX      XX X XXXXX  X XX XXX X XX XXX XX    XX X   XXX  XX XXXX XX      
     X  XX    X X X     XXXX  X   X X  X   X  XX  X X XX X  XXX X    X  XX     
    XXXX XX  XX X XX   X   XXXXX XX XXXXX XXXX XXXX X  X XXX  X XX  XXXX XX    
   X   X  XXX X X  XX XXX X    X  X     X    X    X XXXX   XXXX  XXX   X  XX   
  XXX XXXX  X X XXX X   X XX  XXXXXX   XXX  XXX  XX    XX X   XXX  XX XXXX XX  
 X  X    XXXX X   X XX XX  XXX     XX X  XXX  XXX XX  X X XX X  XXX X    X  XX 
XXXXXX  X   X XX XX  X  XXX  XX   X X XXX  XXX  X  XXXX X  X XXX  X XX  XXXX XX
     XXXXX XX  X  XXXXXX  XXX XX XX X   XXX  XXXXXX   X XXXX   XXXX  XXX   X   
    X    X  XXXXXX     XXX  X  X  X XX X  XXX     XX XX    XX X   XXX  XX XXX  
   XXX  XXXX     XX   X  XXXXXXXXXX  X XXX  XX   X X  XX  X X XX X  XXX X   XX 
  X  XXX   XX   X XX XXXX         XXXX   XXX XX XX XXX XXXX X  X XXX  X XX X XX
XXXXX  XX X XX XX  X    XX       X   XX X  X  X  X   X    X XXXX   XXXX  X X  X
    XXX X X  X  XXXXX  X XX     XXX X X XXXXXXXXXXX XXX  XX    XX X   XXXX XXX 
   X  X X XXXXXX    XXXX  XX   X  X X X           X   XXX XX  X X XX X   X   XX
X XXXXX X      XX  X   XXX XX XXXXX X XX         XXX X  X  XXXX X  X XX XXX X X
X     X XX    X XXXXX X  X  X     X X  XX       X  X XXXXXX   X XXXX  X   X X  
XX   XX  XX  XX     X XXXXXXXX   XX XXX XX     XXXXX      XX XX    XXXXX XX XXX
 XX X XXX XXX XX   XX        XX X X   X  XX   X    XX    X X  XX  X    X  X    
X X X   X   X  XX X XX      X X X XX XXXX XX XXX  X XX  XX XXX XXXXX  XXXXXX   
X X XX XXX XXXX X X  XX    XX X X  X    X  X   XXXX  XXX X   X     XXX     XX X
X X  X   X    X X XXX XX  X X X XXXXX  XXXXXX X   XXX  X XX XXX   X  XX   X X  
X XXXXX XXX  XX X   X  XXXX X X     XXX     X XX X  XXXX  X   XX XXXX XX XX XXX
X     X   XXX X XX XXXX   X X XX   X  XX   XX  X XXX   XXXXX X X    X  X  X    
XX   XXX X  X X  X    XX XX X  XX XXXX XX X XXXX   XX X    X X XX  XXXXXXXXX  X
 XX X  X XXXX XXXXX  X X  X XXX X    X  X X    XX X X XX  XX X  XXX        XXX 
X X XXXX    X     XXXX XXXX   X XX  XXXXX XX  X X X X  XXX X XXX  XX      X  XX
X X    XX  XXX   X   X    XX XX  XXX    X  XXXX X X XXX  X X   XXX XX    XXXX  
X XX  X XXX  XX XXX XXX  X X  XXX  XX  XXXX   X X X   XXXX XX X  X  XX  X   XXX
X  XXXX   XXX X   X   XXXX XXX  XXX XXX   XX XX X XX X   X  X XXXXXX XXXXX X   
XXX   XX X  X XX XXX X   X   XXX  X   XX X X  X X  X XX XXXXX      X     X XX X
  XX X X XXXX  X   X XX XXX X  XXXXX X X X XXXX XXXX  X     XX    XXX   XX  X  
 X X X X    XXXXX XX  X   X XXX    X X X X    X    XXXXX   X XX  X  XX X XXXXX 
XX X X XX  X    X  XXXXX XX   XX  XX X X XX  XXX  X    XX XX  XXXXXX X X     XX
 X X X  XXXXX  XXXX    X  XX X XXX X X X  XXX  XXXXX  X X  XXX     X X XX   X  
XX X XXX    XXX   XX  XXXX X X   X X X XXX  XXX    XXXX XXX  XX   XX X  XX XXX 
 X X   XX  X  XX X XXX   X X XX XX X X   XXX  XX  X   X   XXX XX X X XXX X   X 
XX XX X XXXXXX X X   XX XX X  X  X X XX X  XXX XXXXX XXX X  X  X X X   X XX XXX
 X  X X      X X XX X X  X XXXXXXX X  X XXX  X     X   X XXXXXXX X XX XX  X    
XXXXX XX    XX X  X X XXXX       X XXXX   XXXXX   XXX XX       X X  X  XXXXX   
    X  XX  X X XXXX X    XX     XX    XX X    XX X  X  XX     XX XXXXXX    XX X
X  XXXX XXXX X    X XX  X XX   X XX  X X XX  X X XXXXXX XX   X X      XX  X X X
XXX   X    X XX  XX  XXXX  XX XX  XXXX X  XXXX X      X  XX XX XX    X XXXX X  
  XX XXX  XX  XXX XXX   XXX X  XXX   X XXX   X XX    XXXX X  X  XX  XX    X XXX
XX X   XXX XXX  X   XX X  X XXX  XX XX   XX XX  XX  X   X XXXXXX XXX XX  XX   X
 X XX X  X   XXXXX X X XXXX   XXX X  XX X X  XXX XXXXX XX      X   X  XXX XX X 
XX  X XXXXX X    X X X    XX X  X XXX X X XXX  X     X  XX    XXX XXXX  X  X XX
 XXXX     X XX  XX X XX  X X XXXX   X X X   XXXXX   XXXX XX  X  X    XXXXXXX   
X   XX   XX  XXX X X  XXXX X    XX XX X XX X    XX X   X  XXXXXXXX  X      XX  
XX X XX X XXX  X X XXX   X XX  X X  X X  X XX  X X XX XXXX       XXXXX    X XXX
 X X  X X   XXXX X   XX XX  XXXX XXXX XXXX  XXXX X  X    XX     X    XX  XX    
XX XXXX XX X   X XX X X  XXX   X    X    XXX   X XXXXX  X XX   XXX  X XXX XX   
 X    X  X XX XX  X X XXX  XX XXX  XXX  X  XX XX     XXXX  XX X  XXXX   X  XX X
 XX  XXXXX  X  XXXX X   XXX X   XXX  XXXXXX X  XX   X   XXX X XXX   XX XXXX X X
  XXX    XXXXXX   X XX X  X XX X  XXX     X XXX XX XXX X  X X   XX X X    X X X
XX  XX  X     XX XX  X XXXX  X XXX  XX   XX   X  X   X XXXX XX X X X XX  XX X X
 XXX XXXXX   X X  XXXX    XXXX   XXX XX X XX XXXXXX XX    X  X X X X  XXX X X  
X  X     XX XX XXX   XX  X   XX X  X  X X  X      X  XX  XXXXX X X XXX  X X XX 
XXXXX   X X  X   XX X XXXXX X X XXXXXXX XXXXX    XXXX XXX    X X X   XXXX X  X 
    XX XX XXXXX X X X     X X X       X     XX  X   X   XX  XX X XX X   X XXXX 
   X X  X     X X X XX   XX X XX     XXX   X XXXXX XXX X XXX X X  X XX XX    XX
X XX XXXXX   XX X X  XX X X X  XX   X  XX XX     X   X X   X X XXXX  X  XX  X X
X  X     XX X X X XXX X X X XXX XX XXXX X  XX   XXX XX XX XX X    XXXXXX XXXX  
XXXXX   X X X X X   X X X X   X  X    X XXX XX X  X  X  X  X XX  X     X    XXX
    XX XX X X X XX XX X X XX XXXXXX  XX   X  X XXXXXXXXXXXXX  XXXXX   XXX  X   
   X X  X X X X  X  X X X  X      XXX XX XXXXX             XXX    XX X  XXXXX  
  XX XXXX X X XXXXXXX X XXXXX    X  X  X     XX           X  XX  X X XXX    XX

I will experiment with larger rules: the space of possible behaviors with a neighborhood of 5 is quite a bit larger (there are 2**32 possible rules instead of just 256).

Increasing pyephem’s accuracy for satellite rise/set calculations…

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.

You reeky, motley-minded mammet!

Without further explanation:

[sourcecode lang=”python”]
#!/usr/bin/python

from random import choice
# ,, ,
# ‘ || || A Shakespearean insult generator
# \\ \\/\\ _-_, \\ \\ || =||=
# || || || ||_. || || || ||
# || || || ~ || || || || ||
# \\ \\ \\ ,-_- \\/\\ \\ \\,
#

c1 = [ "mewling", "paunchy", "pribbling", "puking", "puny", "qualling", "rank",
"reeky", "roguish", "ruttish", "saucy", "spleeny", "spongy", "surly"]
c2 = [ "idle-headed", "ill-breeding", "ill-nurtured", "knotty-pated",
"milk-livered", "motley-minded", "onion-eyed", "plume-plucked",
"pottle-deep", "pox-marked", "reeling-ripe", "rough-hewn",
"rude-growing", "rump-fed"]
c3 = [ "lewdster", "lout", "maggot-pie", "malt-worm", "mammet", "measle",
"minnow", "miscreant", "moldwarp", "mumble-news", "nut-hook",
"pidgeon-egg", "pignut", "puttock" ]

print "You %s, %s %s!" % (choice(c1), choice(c2), choice(c3))
[/sourcecode]

Inspired by an image that was posted to Facebook.

Some example python code to fetch weather forecasts…

Need to get some weather information? The website forecast.io has a nice web based service you can use up to one thousand times a day for free. I was thinking of using it for an automated sprinkler application, but just to test it, I wrote this simple chunk of python to try it out. To use it yourself, you’ll need to get your own API key and modify it to use your own latitude and longitude. It’s not that amazing, but you might find it of use.

[sourcecode lang=”python”]
#!/usr/bin/env python

# __ _
# / _|___ _ _ ___ __ __ _ __| |_
# | _/ _ \ ‘_/ -_) _/ _` (_-< _|
# |_| \___/_| \___\__\__,_/__/\__|
#
# A python program which used some publically available
# web apis to find out what the forecast will be.
#

# You’ll need an API key below… you get 1000 requests per day for free.
# Go to forecast.io and sign up.

API="PUT_YOUR_OWN_API_KEY_HERE"
URL="https://api.forecast.io/forecast/"

# Your latitude and longitude belong here, I use SF for example
LAT= 37.7833
LNG=-122.4167

directions = ["N", "NNE", "ENE", "E", "ESE", "SSE", "S", "SSW", "WSW", "W", "WNW", "NNW"]

def bearing_to_direction(bearing):
d = 360. / 12.
return directions[int((bearing+d/2)/d)]

import sys
import os
import time
import optparse
import json

import urllib2

now = time.time()
cached = False

if os.path.exists("WEATHER.cache"):
f = open("WEATHER.cache")
parsed = json.loads(f.read())
f.close()
if now – parsed["currently"]["time"] < 900:
cached = True

if cached:
print "::: Using cached data…"
else:
print "::: Reloading cache…"
req = urllib2.Request(URL+API+"/"+("%f,%f"%(LAT,LNG)))
response = urllib2.urlopen(req)
parsed = json.loads(response.read())
f = open("WEATHER.cache", "w")
f.write(json.dumps(parsed, indent=4, sort_keys=True))
f.close() ;

c = parsed["currently"]
print ":::", time.strftime("%F %T", time.localtime(c["time"]))
print "::: Conditions:", c["summary"]
print "::: Temperature:", ("%.1f" % c["temperature"])+u"\u00B0"
print "::: Dew Point:", ("%.1f" % c["dewPoint"])+u"\u00B0"
print "::: Humidity:", ("%4.1f%%" % (c["humidity"]*100.))
print "::: Wind:", int(round(c["windSpeed"])), "mph", bearing_to_direction(c["windBearing"])

d = parsed["daily"]["data"][0]
print "::: High:", ("%.1f" % d["temperatureMax"])+u"\u00B0"
print "::: Low:", ("%.1f" % d["temperatureMin"])+u"\u00B0"

d = parsed["hourly"]["data"]

for x in d[:12]:
print time.strftime("\t%H:%M", time.localtime(x["time"])), x["summary"], ("%.1f" % x["temperature"])+u"\u00B0"
[/sourcecode]

Mini hack o’ the day: Making IRC talk…

Okay, it’s been a while since I posted anything: I’ve been busy with travel and the holidays, and now I’m trying to get my home office/shack setup so I can pursue some other projects. It is one of those rooms that has piles of crap, some of which I haven’t seen in years, so I’m carefully working through it, tossing stuff that is useless and organizing the remainder.

As a result, I’m trying to also make the space a bit more engaging, so I’ll spend more time there. My eventual goal is to get some of my radios in here so I can listen to more shortwave and ham traffic, but I also spend a fair amount of my off hours on IRC (mostly on the #hamradio channel on irc.freenode.net using my callsign K6HX as my nickname). Instead of forcing me to sit reading a screen, I thought it would be less intrusive and allow me to get more work done if I could monitor the channel by having a voice synthesizer read the msgs that appear in the #hamradio channel. That way, I could just go about my business, but still hear the conversations.

So, that was my idea.

My first thought was to simply run pidgin (a fairly nice IRC client that runs on multiple platform) and just use the pidgin-festival plugin for voice synthesis. But it turns out it was more difficult than I had hoped: nothing I tried seem to make it work. When I cracked open the source code after an hour of flailing, I was annoyed to find a bunch of really questionable code (hence my comment on twitter earlier about chimps writing code), and I soured on the idea. I turned off the computer and went to bed.

On the drive home from work yesterday, I thought about the problem again. I knew that festival could do the speech synthesis. I knew that the curses based irc client epic5 could write a log file of each msg. I then hatched a simple idea: I’d write a little bit of Python glue that would basically act like tail -f: it would repeatedly scan for new lines at the end of the irc.log file, and then use popen to the festival voice synthesizer, spew out what needed to be read, and then continue.

So, here’s the code!

[sourcecode lang=”python”]
#!/usr/bin/env python

import sys
import os
import select
import time
import re

msgpat = re.compile("<([A-Za-z0-9_]+)> (.*)$")
actpat = re.compile("\* (.*)$")

def say(nick, s):
p = os.popen("festival –tts", "w")
p.write(”’%s says %s”’ % (nick, s))
p.close()

def act(s):
p = os.popen("festival –tts", "w")
p.write(s)
p.close()

def process(s):
m = msgpat.match(s)
if m:
say(m.group(1), m.group(2))
else:
m = actpat.match(s)
if m:
act(m.group(1))

f = open("irc.log", "r")
f.seek(0, 2)

while True:
l = f.readline()
while l:
process(l.strip())
l = f.readline()
time.sleep(1)
[/sourcecode]

It works fairly well. There are a couple of things that I need to work on. First, I think the main loop is fairly inelegant. What I want is a version of readline() that works more like a read from a pipe: blocking when there is no further data to be read. But overall, this code mostly prevents the busy wait that would result without the sleep call.

There really isn’t any reason not to use an irc client library and capture the events that you want to speak directly. I didn’t go that way because frankly I didn’t want to spend more than 10 minutes testing out the idea. I’ll probably code up that soon though.

Other things:

  • It doesn’t do good things with emoticons like :-). You could build a dictionary to map things like that to words like smile or grin. Ditto for things like “hmmm” and “hahaha”.
  • It announces the speaker before every message. If the next message is the same person, it could just say what they said without the redundant introduction.
  • It doesn’t generalize to reading my twitter feed.
  • It should also drop things that look like URLS (they aren’t interesting to hear read out loud.)

How a geek tells his wife he loves her: an exercise in Python programming

I’m away from my better half this weekend, visiting my Mom and brother. I scheduled this a few weeks ago, but shortly after Carmen was granted an entry into the Nike Women’s Half Marathon on the same weekend. Rescheduling the visit with Mom would have been hard, so I decided to go anyway, and missed the opportunity to cheer her on.

But I’m a geek. On Saturday night, I hit upon a brilliant, romantic idea (or what passes for one when you are a geek). I knew I was gonna be busy during her run, so I decided to write a script that would text her ever half an hour with an encouraging message. That way, even if I got busy, she’d know I was thinking about her.

It probably would have helped had I started before 10:00pm on Saturday.

First of all, it’s pretty straightforward to send texts to someone using email. If you send an email to phonenumber@txt.att.net, the contents will get routed to the person. And using Python’s smtplib, sending email is pretty easy. Here’s a function that worked for me:

[sourcecode lang=”python”]

def emailtxt(addr, msg):
server = smtplib.SMTP("smtp.gmail.com", 587)
server.ehlo()
server.starttls()
server.ehlo()
server.login("myaccount@gmail.com", "mypassword")
server.sendmail("myaccount@gmail.com", addr, msg)
server.close()
[/sourcecode]

I decided to use my gmail account, but you can use whatever SMTP mail server you like.

So, now that I had a function to txt her, how to handle the scheduling? It turns out that Python has an interesting library called sched which can be used to implement a scheduler. Here’s the interesting bits:

[sourcecode lang=”python”]
import sched
import time

scheduler = sched.scheduler(time.time, time.sleep)

def timefor(t):
tp = time.strptime(t, "%a %b %d %H:%M:%S %Z %Y")
return time.mktime(tp)

scheduler.enterabs(timefor("Sun Oct 17 04:00:00 PDT 2010"), 1, sendit, (0,))
scheduler.enterabs(timefor("Sun Oct 17 04:30:00 PDT 2010"), 1, sendit, (1,))
scheduler.enterabs(timefor("Sun Oct 17 05:00:00 PDT 2010"), 1, sendit, (2,))
scheduler.enterabs(timefor("Sun Oct 17 05:30:00 PDT 2010"), 1, sendit, (3,))
scheduler.enterabs(timefor("Sun Oct 17 06:00:00 PDT 2010"), 1, sendit, (4,))
scheduler.enterabs(timefor("Sun Oct 17 06:30:00 PDT 2010"), 1, sendit, (5,))
…
scheduler.run()
[/sourcecode]

So, what does this do? It creates a scheduler object, which operates in real time (the time.time and time.delay functions are abstract functions that implement a time and delay function). It then uses the utility function timefor to figure out the absolute time when I wanted the event to occur, and when that time expires, it will call sendit, passing it the args in parens (which in my case is a message number in a table of nice things to send her, which I’ve omitted since they are too precious for words). When the run() method is called, the scheduler will wait until the appropriate time, and then make the calls.

I had it set to start sending messages at 4:00AM, since that is when she was going to wake up. Sadly, I screwed up the code slightly, and she didn’t get her first message until 7:30 or so. But she was successfully encourage every 30 minutes after that.

I’m very proud of her and her half-marathon performance. Neither of us are natural athletes, but she’s a constant inspiration to me, and makes me want to be stronger and better.

It’s too bad I’m such a geek.

Passing the Amateur Extra test by guessing…

The Amateur Extra test is 50 questions, multiple choice, with 4 answers per question. A passing grade is 35 or more. A few minutes of programming this morning, even before I had any coffee yielded that the exact probability of passing was:

      4677523340461106447
------------------------------
158456325028528675187087900672

or about 1 in 33.9 billion.

This wasn’t that interesting of a question, but to solve it, I hacked up a quick but limited implementation of rational arithmetic in Python. I was wondering if there was a better way to implement this in Python so overloading would “just work”. I didn’t know how, and the problem was simple enough, so I didn’t try. Here’s my solution. 

#!/usr/bin/env python

def gcd(a, b):
    if (a < b):
        a, b = b, a
    while b != 0:
        a, b = b, a%b
    return a 

class Rational:
    def __init__(self, a, b):
        self.a = a 
        self.b = b 
    def __str__(self):
        return "[%d / %d]" % (self.a, self.b)
    def pow(self, p):
        return Rational(pow(self.a, p), pow(self.b, p))
    def mult(self, r):
        tmpa = self.a * r.a ;
        tmpb = self.b * r.b ;
        d = gcd(tmpa, tmpb)
        return Rational(tmpa//d, tmpb//d)
    def imult(self, i):
        tmpa = self.a * i ;
        tmpb = self.b ;
        d = gcd(tmpa, tmpb)
        return Rational(tmpa//d, tmpb//d)
    def add(self, r):
        tmpa = self.a * r.b + r.a * self.b
        tmpb = self.b * r.b 
        d = gcd(tmpa, tmpb)
        return Rational(tmpa//d, tmpb//d)

p = Rational(1, 4)
q = Rational(3, 4)

def fact(n):
        c = 1 
        while n > 1:
                c = c * n 
                n = n - 1 
        return c 

def comb(a, b):
        return fact(a)/(fact(b)*fact(a-b))

total = Rational(0, 1)

for t in range(35, 51):
        x = p.pow(t).mult(q.pow(50-t)).imult(comb(50, t))
        total = total.add(x)

print "Exact probability is",  total
print "Only about 1 in", total.b // total.a, "will pass"

Python Gadget of the Day – Pydot

Part of a directory treepydot is an interface to the GraphViz suite of programs for drawing abstract graphs and networks. Nifty. I’ve had need of such a thing quite a few times, and it only took me about two minutes to convert a part of my home directory hierarchy into the picture on the right. Fun stuff.

Unix HistoryNot that impressed? How ’bout this slightly nicer diagram of the history of Unix variants? Pretty cool.

pypod

My ipod is on the fritz, so this isn’t immediately useful to me, but check out pypod: a Python library and script for manipulating the song database on your ipod. Neat. When mine comes back from the shop, I’ll have to try it out.

TinyP2P

Ed Felton of freedom-to-tinker has released a tiny 15 line Python program called TinyP2P which allows you to create a simple (if not secure or scaleable) file sharing network. Get the code here. It’s cute, and might not be bad for tiny bits of file sharing.

Addendum: Actually trying to run it, I got 

localhost - - [15/Dec/2004 15:53:33] "POST /RPC2 HTTP/1.0" 200 -
Traceback (most recent call last):
  File "../tinyp2p.py", line 14, in ?
    for url in pxy(ar[3]).f(pw(ar[3]),0,[]):
  File "/u0/markv/my-python/lib/python2.3/xmlrpclib.py", line 1029, in __call__
    return self.__send(self.__name, args)
  File "/u0/markv/my-python/lib/python2.3/xmlrpclib.py", line 1316, in __request    verbose=self.__verbose
  File "/u0/markv/my-python/lib/python2.3/xmlrpclib.py", line 1080, in request
    return self._parse_response(h.getfile(), sock)
  File "/u0/markv/my-python/lib/python2.3/xmlrpclib.py", line 1219, in _parse_response
    return u.close()
  File "/u0/markv/my-python/lib/python2.3/xmlrpclib.py", line 742, in close
    raise Fault(**self._stack[0])
xmlrpclib.Fault: