FluxTrap

Dead Simple Real Time Plotting with C/C++ and Python

2011-09-27T19:57:00.000-07:00

Concepts:

Standard input/output (Think: the input and output of your program)

Unix Pipes

Plotting with matplotlib

So, you have your code working, but you are tired of having to run a separate program to see plots? There are many simple solutions, but I am going to present what I think is the absolute simplest. What we are going to do is have your C/C++ program output the data and have python capture it and save plots in real-time! What does this entail? Well...

Your C/C++ program no longer writes to a file (using fprintf), but rather writes to the standard output (STDIN, using simply printf)

Your python script will now read data straight from STDIN using raw_input()

You will have to use a unix shell to sew all of this together

Let's look at a very simple example. Suppose we have a very important C program that outputs some very important numbers. In order for it to work with this new standard, we have it output the data to screen. Namely, it does:

 
makedata.c 
#include <stdio.h>
#include <math.h>
 
#define PI 3.14159
 
int main()
{ 
  int i,j;
  for(j=0; j<10; j++) {
    for(i=0; i<20; i++) {
      printf("%f\t",sinf(i * PI / 10.0 + j*PI/10));
    }
    printf("\n");
  }
 
  return 0;
}

Now, when we run this, we get a bunch of numbers thrown to screen!

 
$ gcc -lm -o makedata makedata.c
$ ./makedata
0.000000 0.309017 0.587785 0.809017 0.951056 1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 
0.309017 0.587785 0.809017 0.951056 1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 
0.587785 0.809017 0.951056 1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 
0.809017 0.951056 1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 
0.951056 1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 
1.000000 0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 0.951055 
0.951057 0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 0.951055 1.000000 
0.809018 0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 0.951055 1.000000 0.951059 
0.587787 0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 0.951055 1.000000 0.951059 0.809021 
0.309019 0.000003 -0.309014 -0.587783 -0.809015 -0.951055 -1.000000 -0.951058 -0.809020 -0.587789 -0.309022 -0.000005 0.309012 0.587781 0.809013 0.951055 1.000000 0.951059 0.809021 0.587792

In order to capture the data in python, we must use the raw_input() function. This function simply gets input from the user and puts it into a variable. It puts everything the user types up to when they press enter. This is why the C code is that it only prints a newline (ie: '\n') once one full line of data has been outputted to screen. If we had put a newline in the first printf statement, the python plotting program would only plot one number at a time! So, you can think of the tab (\t) as deliniating between values and the newline (\n) deliniating between different sets of data. The python code that reads this data looks like:

 
plot.py 
import numpy as np
import pylab as py
 
def plot_data(data):
  py.clf()
  py.plot(data)
  py.show()
  py.savefig("data-%.8d.png"%counter)
 
if __name__ == "__main__":
  counter = 0
  while True:
    try:
      tmp = raw_input().strip().split()
      data = np.array(tmp, dtype=np.double)
    except EOFError:
      print "Input has terminated! Exiting"
      exit()
    except ValueError:
      print "Invalid input, skipping.  Input was: %s"%tmp
      continue
 
    print "Plotting plot number %d"%counter
    plot_data(data)
    counter += 1

You can test this program by running it, typing a bunch of numbers separated by a space, then pressing enter. It will plot it, display it and save it! Then, the program will ask you again for more numbers. To exit, you type Control-D which makes the EOFError happen.

What is going on in this program is quite simple. First, "tmp" gets the long string of characters that you typed in. However, python doesn't know it contains numbers, it just looks like a bunch of random characters! Now, we use numpy and tell it to create an array out of the data. The "dtype=np.double" is us telling numpy that we are realing with valid numbers. A ValueError happens if we weren't good on our promise and the input isn't in fact all numbers.

Now for the most important part... how do we put these two things together? Unix has a very cool thing called input/output redirection. This allows us to redirect the output of one program to the input of another. So, instead of us having to type in the numbers for the python script, we can have the C/C++ program type it for us! The syntax is quite simple, all you have to do is:

 
$ ./makedata | python plot.py

And now you are done! You should have a bunch of plots coming up of sin waves with various phases. Congrats!

There is one more thing you can do to make your plots even more fancy. Sometimes, you don't want to save each figure or have to click through to see every plot, one at a time. Instead, you just want to see an animation of what is happening as it is happening! Or, you are already making an animation with many py.plot() statements, and you want it to be smoother and faster! To do this, you can to look into pylab animations. The people at scipy have a great tutorial on this issue. You can also look at a small plotting script I made which does something very similar.

The basics of this method involve: creating your plots at the beginning of your script, and saving them into variables. Then, when you get new data that you want to plot, you simply change the data in the plot with .set_data(). One thing to note is the line "py.ion()" right after I imported pylab and how I use py.draw() instead of py.show(). If you want to get started playing around with this, simply take the same code from earlier in this document, add "py.ion()" after we import pylab, delete the py.savefig() line and replace py.show() with py.draw()! This will give you a (quite slow) animation.

fortran + python = good

2010-01-06T23:00:00.000-08:00

So, after a long time working solely in FORTRAN, I decided to try to find a way to use python to get the same things done. I wanted to incorporate what I had done in FORTRAN as a library that way I wouldn't have to recode much and could instead continue developing (except now in the wonderful warmth of python)

The most successful method I tried was using f2py. What f2py does is take in a FORTRAN module (see the tvd module below, for example) and makes a python interface for it. What this means is that you can import fortran code and run it at the the speeds you would expect from compiled code!

Let me first take you through the module. You can see from the 'include "omp_lib.h"' that I'm using openMP. For those that don't know, openMP is a great way of making any code parallel on shared memory machines (ie: if you have a computer with multiple cores or hyper-threading). The lines with "!$OMP", called pragmas, are in fact calls to openMP and effect to proceeding code. In addition, I created a function called "set_num_threads" to manually change the ammount of parallelization that is used by the program (by default, the number of threads is equal to the number of CPU's on the machine). A great thing about this "!$OMP" convention is that, since ! designates a comment, you can still compile and run this code without having openMP enabled!

Now, let's move on to f2py! You can see every subroutine (and function) has lines starting with "!f2py". Similar to openMP, this is used to designate various options. For example, the "!f2py threadsafe" tells python that the function will be using various threads and that the gil shouldn't be used. As a result, we cannot take in or return in python objects (there are methods to take in python callbacks when not using threadsafe). The intent option sets if a parameter is going to be inputed or outputted. With intent, python knows how to interpret python input and what should be given back as return values.

So, let's look at the function set. It takes in, in python, u and CFL. The other parameters are optional. In addition, you can see that "dt" and "maxv" are inputs but are calculated in the function. I did this little hack so that I can return multiple values to python. So in the end, this function is defined, in the python space, as "u, dt, maxv = step(u, cfl)". The most beautiful part is that when passing u to this module from python, u is a numpy array! Note that the other functions that are called (namely, do[XYZ] are also defined in the module, however for the sake of brevity they are not listed below).

In order to use this fortran module, we need to compile it with f2py. The end result is a .so that can be imported into fortran. I compiled this particular module with `f2py -c -m tvd --f90flags="-fopenmp -lm " --opt="-O3 --ffast-math" -lgomp --fcompiler=gnu95 tvd.f90` So, let's go through this. "-c" tells f2py to look at the pragmas in the code in order to learn how each function should be handled. "-m tvd" defined the name of the resulting module. the --f90flags, as expected, dictate the flags send to the FORTRAN compiler. I have chosen to use openmp (not needed, by the way, with newer versions of GCC where openMP is built-in) and libmath. opt sets the optimization flags. -lgomp tells f2py to use openMP (this redundancy seems necessary although I have a feeling that one of my statements telling f2py to use openMP is redundant). Finally, "fcompiler" tells f2py to use gfortran.

Let's look at a little session with this module, as to understand what f2py created. (I use ipython for its tab-completion)

In [1]: import numpy

In [2]: from tvd import tvd

In [3]: u = numpy.random.random((4,100,100,100)) #create a random velocity field

In [4]: u[0,:,:,:] = 1 #set density to 1 to avoid instabilities

In [5]: %timeit output = tvd.step(u, 0.65)
10 loops, best of 3: 704 ms per loop

Amazing! It takes 704ms to evolve a 100x100x100 isothermal grid! Normal usage of this would be: `u, dt, maxv = tvd.step(u, 0.65)`. One note, I write "from tvd import tvd" because the first "tvd" represents the file and the second "tvd" represents the FORTRAN module. If instead I had just written tvd.f90 as a series of subroutines and not as a FORTRAN module, I could have simply written "import tvd". Now, I can get to doing real science instead of being swamped by the inhuman syntax of FORTRAN (a language best spelled in all caps to truly illustrate its crimes).

And now, for the code... I'm leaving out the actual computational part however I may release that in the near future (comment if you'd like to see it!).

tvd.f90:
MODULE tvd

  implicit none
  include "omp_lib.h"

  DOUBLE PRECISION :: csound=1

  CONTAINS

  SUBROUTINE set_num_threads(n)
    !f2py threadsafe
    !f2py intent(in) n
    INTEGER :: n
    CALL OMP_SET_NUM_THREADS(n)
  END SUBROUTINE set_num_threads

  SUBROUTINE step(u, n, CFL, dt, maxv)
    !f2py threadsafe
    !f2py intent(in) u
    !f2py intent(in) n
    !f2py intent(in) CFL
    !f2py intent(in), optional dt
    !f2py intent(in), optional maxv
    !f2py intent(out) :: u, dt, maxv
    !f2py depend(u) n
    INTEGER :: n, k
    DOUBLE PRECISION, DIMENSION(4,n,n,n) :: u
    DOUBLE PRECISION :: CFL, dt, maxv
  
    !$OMP PARALLEL DO shared(u) private(n,k) &
    !$OMP SCHEDULE(dynamic) REDUCTION(MAX:maxv)
    DO k=1,n
      maxv = max( maxval( abs(u(2,:,:,k)) / u(1,:,:,k) ), &
                  maxval( abs(u(3,:,:,k)) / u(1,:,:,k) ), &
                  maxval( abs(u(4,:,:,k)) / u(1,:,:,k) )  &
                )
    END DO
    !$OMP END PARALLEL DO
  
    !Calculate the timestep based on the courant condition
    if (dt .eq. 0.0) dt = CFL / (csound+maxv)
   
    !perform strang splitting using the hydro only 2nd order
    !TVD algorithm given by http://arxiv.org/abs/astro-ph/0305088
    CALL doX(u,n,dt, maxv)
     CALL doY(u,n,dt, maxv)
      CALL doZ(u,n,dt, maxv)
      CALL doZ(u,n,dt, maxv)
     CALL doY(u,n,dt, maxv)
    CALL doX(u,n,dt, maxv)
  END subroutine step
END MODULE

SMS/Twitter Control of your server

2009-05-07T15:17:00.000-07:00

Well, I thought it was about time to post something, so let me show you this little project I've been tinkering with. It basically allows me to control my server with a cellphone. This can be applied to any computer with python installed and can do really anything (right now I have it search/download bittorrent files).

It basically works by creating a master thread on the server that constantly checks some secured twitter account for directed messages. This basically means that you can send a directed message to some secured user through twitter from your cell phone and when the computer reads this message it will act on it.

The following is the main code: it creates a main thread on the computer listening for new messages. When a new message is found, a new thread is created and the correct command is executed (the commands are defined in the configuration file which will come soon). The benefit to creating new threads is that the task can be halted later by initiating a stop command (the stop commands are hard programmed and are automatically generated for each configuration defined command).

#!/usr/bin/python
#
# Copyright 2008 Michael Gorelick 
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see:
# http://www.gnu.org/licenses/gpl.html.
from __future__ import division

'''A tool to allow for remote control of a machine with twitter'''

__author__ = 'Michael G '
__version__ = '0.1'

from sys import stdout
import threading, time, twitter, re, ConfigParser, subprocess

######### CONFIG FILE ##############
configfile = 'twitcontrol.conf'

class TwitThread(threading.Thread):
  def __init__(self,actions,twitter,command):
    threading.Thread.__init__(self)
    self._stop = threading.Event()  #Set the hook for self.stop()
    self.actions = actions
    self.twitter = twitter
    self.commandstruct = command
    self.process = None

  def run(self):
    '''Finds what action is linked to the inputed command'''
    command = self.commandstruct['text']
    for name,contents in self.actions.iteritems():
      found = contents['match'].match(command)
      if found:
        self.name = name
        print "%s: %s: Starting"%(time.ctime(),name)
        self.doCommand(found,command,contents)
        return None
      elif command.lower() == "stop %s"%name.lower():
        self.name = "killthread-%s"%name
        for thread in threading.enumerate():
          if thread.name.lower() == name.lower():
            print "%s: Killing %s"%(time.ctime(),name)
            thread.stop()
        return None

  def doCommand(self,match,command,action):
    '''Runs the command given by the config file for the given input'''
    if action["output"] == "start":
      self.message("%s: Started"%self.name)
    print "Running command: %s"%(action["command"]%match.groupdict())
    self.process = subprocess.Popen(action["command"]%match.groupdict(),\
      stdin=None, stdout=subprocess.PIPE, stderr=subprocess.PIPE, shell=True)
    #We use the following hack because Popen.wait() or Popen.interact() lock
    #the process so we cannot kill it if another thread sends a self.stop()
    while self.process.poll() == None:
      time.sleep(1)
    if action["output"] == "end":
      self.message("%s: Done"%self.name)
    if self.process.returncode != 0:
      err = self.formatstd(self.process.stderr)
      print "%s: %s: Error: %s"%(time.ctime(),self.name,err)
      if action["output"] == "output":
        self.message("%s: Error: %s"%(self.name,err))
    else:
      suc = self.formatstd(self.process.stdout)
      print "%s: %s: Success: %s"%(time.ctime(),self.name,suc)
      if action["output"] == "output":
        self.message("%s: Success: %s"%(self.name,suc))

  def formatstd(self,fd,maxsize=120):
    '''Extracts and formats data from a datastream'''
    return ", ".join([x.strip() for x in fd.readlines()])[:maxsize]

  def message(self,msg):
    '''Sends a tweet to the command sender truncated to 140 characters'''
    self.twitter.PostDirectMessage(self.commandstruct['sender_id'],msg[:140])

  def stop(self):
    '''Stops the current thread and terminates any open processes'''
    print "%s: %s: Thread Killed"%(time.ctime(),self.name)
    self.process.terminate()
    self._stop.set()

  def stopped(self):
    return self._stop.isSet()


class TwitControl:
  def __init__(self,actions,user,passwd,recieveID,timeout=5):
    self.actions = actions
    self.twitter = twitter.Api(username=user,password=passwd)
    self.recieveID = recieveID
    self.timeout = timeout

  def getMessages(self):
    '''Gets all directed messages from twitter from a given user'''
    return [x.AsDict() for x in self.twitter.GetDirectMessages() \
    if x.GetSenderId() == self.recieveID]

  def start(self):
    '''Checks for new messages and commands the threads'''
    commands = self.getMessages()
    numcommands = len(commands)
    while True:
      commands = self.getMessages()
      print ".",; stdout.flush()
      if len(commands) > numcommands:
        for command in commands[0:len(commands)-numcommands]:
          print "%s: MASTER: Found Command: %s"%(time.ctime(),command['text'])
          TwitThread(self.actions,self.twitter,command).start()
        numcommands = len(commands)
      elif len(commands) < numcommands:
        numcommands = len(commands)
      time.sleep(self.timeout)

if __name__ == '__main__':
  configfd = ConfigParser.RawConfigParser()
  configfd.read(configfile)
  print "Reading config: %s" % configfile

  username  = configfd.get('Connection', 'username')
  password  = configfd.get('Connection', 'password')
  recieveID = configfd.getint('Connection', 'recieveID')
  if configfd.has_option('Connection', 'timeout'):
    timeout = configfd.getint('Connection', 'timeout')
  else:
    timeout = 5

  actions = {}
  for section in configfd.sections():
    if section == 'Connection':
      continue
    action = {}
    action['match']   = re.compile(configfd.get(section,'match'))
    action['command'] = configfd.get(section,'command')
    action['output']  = configfd.get(section,'output')
    actions.update({section:action})
  print "Loaded modules: %s"%", ".join(actions.keys())

  print "Listening to user: %s"%username
  control = TwitControl(actions,username,password,recieveID,timeout)
  control.start()

And the following is a sample from the configuration file:

# Copyright 2008 Michael Gorelick 
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see:
# http://www.gnu.org/licenses/gpl.html.

[Connection]
username = Secured Username
password = Password of secured username
recieveID = User ID of the command issuer

[reverseSSH]
match = reverse (?P(?:\d{1,3}\.?){4}) (?P\d+)
command = ssh -R %(PORT)s:localhost:22 %(IP)s
output = start

[Talk]
match = say (?P.+)
command = flite -t "%(text)s"
output = None

[Bittorrent]
match = bt (?P.+)
command = transmission-remote -a "`/home/fiber/projects/pyisohunt/pyisohunt.py -n 1 -s seeds -p '%%(url)s' %(search)s`"
output = output

The commands use regex in order to extract keywords. So basically, I can SMS 'd user say hello there' and my computer will initiate the command `flite -t "hello there"` once it receives it! I also have the output of each command regulated by the 'output' parameter on each command. It can either notify me when the command starts, when it ends, or it can send me the output of the command (or, of course, none of the above).

Now, the interesting part is the bittorrent section. Here I send the keywords to another little program I made (really just a bunch of regex and rss tricks... I'll include it at the end). This program returns the URL of the highest seeded match to the keywords from IsoHunt.com. This is sent to transmission and my server starts downloading the torrent! (I use Transmission for bittorrent on my server with the Clutch web UI) So essentially, if a friend recommends a movie (that is in the public domain, of course ;) when I am out, I can tell my server to start downloading it so it will be ready when I am home.

Also, I can set my server to start a reverse-ssh connection to any IP. This is great because I can keep my router locked down but still connect to my server from the outside. Note: normally this sets my server to set up a reverse connection to an intermediary server so that everyone can be behind a firewall, but I'd like to keep the location of this server unknown for now.

Finally, the following is the source for pyisohunt.py. It was a true hack so don't expect clean code!

#!/bin/env python
#
# Copyright 2008 Michael Gorelick 
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see:
# http://www.gnu.org/licenses/gpl.html.

import feedparser, re, getopt
from sys import argv

class pyisohunt:
  def __init__(self, search,cat):
    self.search = search
    self.category = cat
    self.isohunt = feedparser.parse( \
      "http://isohunt.com/js/rss/%s?iht=%d"%(search.replace(" ","+"),cat))
    self.results = self.parse(self.isohunt.entries)

  def parse(self,entries):
    results = []
    findstats = re.compile(u"Size: (?P[0-9.]+) MB, in (?P[0-9]+) files")
    findhealth = re.compile(u"Seeds: (?P[0-9]+) \ \; \| \ \; Leechers: (?P[0-9]+) \ \; \| \ \; Downloads: (?P[0-9]+)\<")
    fieldtype = {"size":float, "files":int, "seeds":int, "leechers":int, "downloads":int, "title":unicode,"url":str} 
    for item in entries:
      tmp = {}
      tmp['title'] = item.title[:item.title.rfind('[')]
      tmp['url'] = item.enclosures[0].href
      tmp.update(findstats.search(item.summary_detail.value).groupdict())
      tmp.update(findhealth.search(item.summary_detail.value).groupdict())
      results.append(dict([(title,fieldtype[title](value)) for title,value in tmp.iteritems()]))
    return results

  def sort(self,key="seeds"):
    return sorted(self.results,lambda x, y: y[key] - x[key])

if __name__ == "__main__":
  categories = {"all":-1,"video":1,"tv":3,"audio":2,"musicvideo":10,\
                "game":4,"app":5,"pic":6,"anime":7,"comic":8,"book":9,\
                "misc":0,"unclassified":11}
  sortfields = ["size","files","seeds","leechers","downlads"]
  help = """%s [-c category] [-s sort] [-n maxresults] search
  Search through IsoHunt torrents on the commandline.  By Michael G.
  -c  Category name. Can be: """ + ", ".join(categories.keys()) + """
  -s  Sort field.  Can be: """ + ", ".join(sortfields) + """
  -p  Print string (standard python printf notation... try with '-p ""' to see fields)
  -n  Max number of results (Note: sort is applied before results are truncated\n"""

  #Default valus
  category = categories["all"]
  sortfield = "seeds"
  maxresults = None
  printf = """Title: %(title)s
  URL: %(url)s
  Size: %(size)0.2f MB
  Files: %(files)d
  Seeds: %(seeds)d
  Leechers: %(leechers)d
  Downloads: %(downloads)d
  """

  try:
    opts, args = getopt.getopt(argv[1:], "hc:s:n:p:")
  except Exception:
    print "Usage: " + help%argv[0]
    exit(1)
  errors = []
  if args == []:
    errors.append("Must have search terms!")
  else:
    search = " ".join(['"%s"'%x for x in args])
  for o, a in opts:
    if o == "-c":
      try:
        category = categories[a]
      except KeyError:
        errors.append("Invalid category")
    elif o == "-p":
      try:
        printf = str(a)
        test = a%{'files': 1, 'title': "1", 'url':"1",u'downloads': 1, u'leechers': 1, u'seeds': 1, u'size':1}
        del test
      except KeyError, e:
        errors.append("Invalid printf string: Invalid Key: %s"%",".join(e))
      except ValueError:
        errors.append("Invalid print string")
    elif o == "-s":
      if not a in sortfields:
        errors.append("Invalid sort field")
      else:
        sortfield = a
    elif o == "-n":
      try:
        maxresults = int(a)
      except ValueError:
        errors.append("Invalid max number of results (must be integer)")
    elif o == "-h":
      print help%argv[0]
      exit(0)
  if errors != []:
    print "Usage: " + help%argv[0] + "\n".join(errors)
    exit(1)

  results = pyisohunt(search,category)
  for result in results.sort(sortfield)[:maxresults]:
    if printf:
      print printf%result
    else:
      print result

Well, I don't really know what else to say about this... ask questions and I'll be sure to give you more details than you ever wanted!

eeePC 901 Antenna Mod

2008-12-22T14:41:00.000-08:00

Well, after the treachery of exams and then the confusing transition into vacation I finally think it's time to update! Well, I've done quite a lot: everything seems ready for the final version of the 3D whiteboard along with the construction of the 3D display. Furthermore everything is in motion to begin the fusion device first semester of the next academic year (damn you Darryl for leaving to Singapore!). But, in the realm of the actual rather than the theoretical, I have modded my eeePC 901 to support external antennas (for only $20).

I love my eeePC, it has changed the way I work. I carry it around with me everywhere scarcely ever noticing its weight allowing me to do work almost anywhere! It's powerful enough for testing, enough hard drive space for my tools, runs of my simulation and even some music. The only downside when I first got it: the antennas weren't powerful enough to pick up those elusive wifi signals. Although, after cracking the magic device open I found that the internal antennas were hooked up using U.FL connection. And so the hacking began!

[picture of wireless module to come in very near future]

Above you can see the actual wifi module from the expansion bay (this is the removable region in the center of the bottom of the laptop). This led me to buy a U.FL to RP-SMA cabel from hyperlinktech (http://www.hyperlinktech.com/item.aspx?id=1123) as seen below:

I chose RP-SMA for the outside connection because I already had an RP SMA antenna from an old router and it seems to be the most popular kind (so getting new antennas shouldn't be too much of a problem). I also bought a 9dB 2.4GHz antenna from the same store. 2.4GHz is the frequency of wireless internet (802.11a/b/n). In order to make everything fit into the eee, I had to sand the stopper almost all the way down until it looks like:

In order to make this hack a nice clean one, I decided to completely gut out the computer so I could put everything back where I wanted to. There are MANY screws so make sure to keep everything organized. I normally take a sheet of paper and put the screws on different parts of it while marking in marker where those screws should be (and how many of them I took out!). I also took out the screen because I was going to machine the hole, but I later realized it was easy enough to do at home. Here's a beauty shot of everything dismantled:

I decided to have the RP SMA connection come out of the "lock" region of the case (right behind the eth0 port). So, I had to make sure that the wire didn't block or bend the motherboard (scary thought! one crack and you have transformed your eeepc into a nice paper weight!). So, I drilled the hole for the connection making sure that it's center was .5mm below the stand of the motherboard (circled in red). To make the hole i used a standard dremel tool outfitted with a 1/8" drill bit.

With that completed I snaked the wire underneath where the motherboard sits. I now plugged the cable to the wifi module, unplugging one of the internal antennas (I recommend doing this before remounting everything as it takes a good amount of force. I even detached the wifi module to make sure there'd be no damage to the motherboard). Once everything is hooked up, it should look like this:

This has been a tremendous gain! The 3dB antenna fits nicely in the carrying case (and so does the eee even with the added width) and the 9dB fits nicely in my backpack. Below is a comparison between only internal antennas, external mod without any antenna and finally external mod with 3dB antenna using `iwlist wifi0 scanning`


#before
fiber@mercury:~$ cat wireless.before
 Cell 01 - Address: 00:1E:58:40:6D:8F
           ESSID:"telepathy"
           Mode:Managed
           Channel:11
           Quality:100/100  Signal level:-31 dBm  Noise level:-81 dBm
#after without antenna
fiber@mercury:~$ cat wireless.afternoext
 Cell 01 - Address: 00:1E:58:40:6D:8F
           ESSID:"telepathy"
           Mode:Managed
           Channel:11
           Quality:65/100  Signal level:-64 dBm  Noise level:-81 dBm
#after with 3dB antenna
fiber@mercury:~$ cat wireless.after3db
 Cell 01 - Address: 00:1E:58:40:6D:8F
           ESSID:"telepathy"
           Mode:Managed
           Channel:11
           Quality:100/100  Signal level:-29 dBm  Noise level:-81 dBm

30% gain with a cheap antenna I just had lying around. I'll post some stats on the 9dB antenna once I am in a controllable enough area (ie: home). Remember, dBm work on the log base 10 scale, so every unit change in dB is a factor of 10 change in strength.

So tell me how this mod works for you or if there are any points that I didn't clearly explain! Happy Hacking!

new project

2008-11-18T09:58:00.000-08:00

Well, it's taking longer than I expected for the parts for my most recent project to come in. I'm waiting for some wireless pigtails (U.FL to RP-SMA) so that my EEE PC has an external antenna (an easy hack, but one I haven't seen any documentation online for).

In the meantime, here's an old video I made a while ago playing with non-Newtonian fluids. This is really easy to make (just saturate water with corn starch) and very fun! Non-Newtonian fluids are fluids whose viscosity is dependant on the applied strain. So, the strain from the speaker constantly hitting the fluid causes it to become much more viscous to the point of appearing solid. This causes all sorts of fun! In this video I played a 100hz sin wave through the speaker. I did this later as a demonstration at the Dragon Academy and I have to say it was a great success, I recommend it for any teachers or demonstrators. Hope everyone enjoys!

obama fever

2008-11-04T21:32:00.000-08:00

OBAMA WON! Truly a day to remember. So, in the spirit of science funding not being slashed to death I am going to release one of my first python scripts that downloads all images off of the NASA Astronomy Picture of the Day (apod).

I tried keeping the code commented so that it could also serve as a good learning tool. One thing I'd like to include is some way of saving the description of the image. This could be done by saving it into a text file with the same filename as the image, but I suspect there is a better way (using ImageMagick to embed the information?). If you have an idea or preference, post a comment! So, here it goes:


#!/bin/env python
#
# Copyright 2008 Michael Gorelick 
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see:
# http://www.gnu.org/licenses/gpl.html.

import re, urllib, sys, os

### Check/Do Help
help = """NASA Image Getter - Michael G. (GPL 2008)
./%s [-h] [imagepath]
imagepath  - path to save images
-h, --help - this help"""%sys.argv[0]

if "-h" in sys.argv or "--help" in sys.argv:
exit(help)
elif len(sys.argv) == 1:
print "Not enough program arguments"
exit(help)

### CONFIG
allowedChars = '[-a-zA-Z0-9*:\'.,&!?\(\)//\n"+;_ ]'
allowedFiles = '[jpg|jpeg|png|gif|mov|avi|mpeg|mpg]'
url = "http://antwrp.gsfc.nasa.gov/apod/archivepix.html"
pbarwidth = 80
try:
imagePath = sys.argv[1]
os.makedirs(imagePath)
except Exception, e:
#errno == 17 means the directory already exists.
if e.errno != 17:
 exit(e)
print "Saving images to %s"%imagePath

##############################################
## This section looks for links to image pages
print "Searching for image pages in", url
html = urllib.urlopen(url).read()
pagePattern = re.compile('([0-9]{6}).html">(' + allowedChars + '*)')
pages = pagePattern.findall(html)
totalpages = len(pages)
print "Found", totalpages, "image pages."

#############################################
## Now we go to the pages and extract a list
##  of images and download
try:
alreadyHave = open(imagePath + "/tracker.log").readlines()
except IOError, e:
alreadyHave = ""
totalhave = len(alreadyHave)
try:
tracker = open(imagePath + "/tracker.log", 'a+')
except IOError, e:
exit(e)

base = url[0:url.rfind('/')]
imagePattern = re.compile('<a href="image/(' + allowedChars + '*.' + \
allowedFiles + ')"', re.IGNORECASE)
sanitizePattern = re.compile('[//\n\r\t]')
imageGet = urllib.URLopener()

for i, page in enumerate(pages):
try:
 if page[0]+"\n" not in alreadyHave:
   #Extract information from image page
   title = sanitizePattern.sub("", page[1])
   currentpage = base + "/ap" + page[0] + ".html"
   print "Getting " + title
   content = urllib.urlopen(currentpage).read()
   image = imagePattern.findall(content)

   #Parse data
   filename = imagePath + "/" + (title + \
     image[0][image[0].rfind('.'):]).replace(" ","_")
   print "\tFilename: " + filename
   url = base + "/image/" + image[0]
   print "\tURL: " + url

   print "\n[" + "="*(pbarwidth*(i-totalhave)/(totalpages-totalhave)) \
     +"@"+"="*(pbarwidth*(i-totalpages)/(totalhave-totalpages)-1)+"]",
   sys.stdout.flush()

   #Download
   imageGet.retrieve(url, filename)
   tracker.write(page[0] + "\n")

   print "\r" + " "*(pbarwidth+2),"\r",
   sys.stdout.flush()
except IndexError, e:
 print "\tImage not found\n"
except IOError, e:
 print "\t404 Error\n"

tracker.close()
print "Done"

Distractions

2008-11-03T09:06:00.000-08:00

Well, exams are finally done but I am still contending with all the work I had to put off. In the works is a post on the topic of conscience and the release of my mac remote code. For now, I will distract everyone with an old youtube video of mine: exploding flint (well, actually ferrocium)

The video pretty much explains itself. It is a nice little demo on the properties of ferrocium. In the end, this stability (or lack thereof) is why it proves to be so useful in lighters. Enjoy!

failed telescope night

2008-10-26T00:49:00.000-07:00

So, after a while spent getting my telescope aligned (and, of course, staring at M41 and M42 for a while), the clouds came in before I could get the camera ready. When I was done cursing the clouds and the cold I had been sitting in, I decided to take pictures of Toronto to gauge the camera.

I was using my 8" Celestron (f#=10, focal length=2032 mm). It's a great telescope although it takes a magical touch to get aligned.

This is the main picture of interest: the symbol of BMO on top of First Canadian Place (which is 2580m away from my balcony). From my calculations, the symbol is 6m tall which means I'm capturing about .1332° or 2.32e-3rad of the sky! (I got this number by doing )

By the way, the picture of the moon I took before was with my astroscan (f#=4.1, focal length=445cm). It's a great telescope to bring around and it's really quick to get set up and start observing, although the optics aren't that great which leads to a lot of artifacts. Also, I like the tactility that comes with manually aiming the beast.

While I was first testing the camera I took the following picture from 37m away:
Since my camera is 9.5cm tall and it takes up 46.1% of the height and 34.6% of the width of the image we can get... (5.563e-3, 7.413e-3)radians = (.318,.424)degrees=(height,width) degrees!

So, if we take the differences in f# into account, both measurements agree! *phew*, well at least now I can say tonight wasn't a waste.

oscilloscope!

2008-10-25T16:46:00.000-07:00

So, I was at active surplus (with trusted, albeit slightly slow, sidekick Darryl) looking for parts for the fusor or the 3D display when I ran into this amazing old school oscilloscope. After quickly brainstorming some reason I could justify buying this beautiful piece of science equipment (or history depending on your view) I got it! The justification is actually quite good: instead of buying a Geiger counter, we can use activated aluminum and a charged plate (neutron hits the aluminum, shoots of an electron which gets captured by the plate which turns into a signal that the oscilloscope can see) to verify fusion is happening.

Anyways, here's a little video of the piece in action: hooked up to my speakers playing Roads by Portishead. Sorry about the video quality... also the cameras refresh rate doesn't do justice to the display.

As for the fusor route, we are going to get florescent transformers (~20,000V at a couple miliamps) for power, and hopefully e-bay will have the correct vacuum pump (it's hard to find ones that go down to 5microns). As for the vacuum container... we might just have to get one custom made. This is one piece that I don't really want to try building myself since it's the only part of the project that could seriously maim us.

Oh, and stay tuned... tomorrow I'm going to post a video I did a while ago about getting a cool explosion using ferrocium which you can get out of a common lighter!

panoramas

2008-10-24T07:56:00.001-07:00

So, a little bit ago I got a new camera (well, I guess it's no longer new) so I thought, since I couldn't afford any great lenses, I'd get acquainted with good post-processing. First thing to attend to was mending a mosaic of pictures together into one seamless picture. So, bellow is my evolution in trying to get a good panorama:

Balcony View:
This one didn't look too bad. I was lucky that lighting conditions didn't change too much as I took all the photos.

Night Panorama:
I was rushed. As a result I didn't get enough pictures to get good perspective. Also, you can see the lighting conditions changed drastically as I took the pictures.

Dawn Panorama:
That's not change, that's more of the same! I took these pictures fast enough to have some resemblance to homogeneous lighting, even though it's homogeneously dark.

Toronto Through a fishbowl:
This one turned out great! I took more than enough pictures and did it quickly. You can see a few artifacts from mending (check the train tracks), but all in all it turned out great. In fact, there are parts of the image that were cut out of all the photos that I had to reconstruct completely using photoshop (special prize to whoever can guess where). If you want to see a full resolution (17mb, 5989x4058) version, go here.

Sunset:
So, I was doing work and I look out the window and I see a rainbow AND the moon. So, I quickly try to get as many pictures as I can for a panorama. Didn't turn out exactly how I wanted, but it is still quite cool!

For all these images, I loaded them into photoshop and used the automated function "photomerge". Once the images were merged I blended the edges and applied any sort of photoshop magic I could to make the image seem smooth. This was hard in some cases (ie: in the last image, sunrise, you can see a drastic color change in the sky).

These skills are actually quite useful for many things, especially the telescope camera. I used basically the same technique to make the moon picture from the telescope camera (although I probably should have used a more stack oriented program... that's what the future is for!).

So, feel free to download all these images, just remember they are all licenced under the Creative Commons (Attribution-Noncommercial-Share Alike 2.5 Canada).

Projects to come

2008-10-23T07:47:00.001-07:00

I also want to use this blog as a place to write down the projects I want to try doing on top of talking about the ones I've already completed. So, here's a list of what's coming up and what I would like to come up:

In the works:

3D Whiteboard (almost done!)
3D Volumetric Display (anyone want to donate an optical bench or some mirrors? :) (with help from Darryl)
Fusor (anyone know how to make a cheap vacuum chamber?) (together with Darryl)
python Hobbyist Observer Tool (phot?, almost done... but not ready for the public yet)
pyremote (python tool set for the mac remote... basically done, I'm just making some modules ftw)

To come:

Ferrofluid (I'm doing this one for the kids at a local highschool)
Hologram of Saturn (well, I have the camera, now I just need some clear nights!) (Helping Darryl)
Micron sized tractor beam! (too cool to say no)
Tape as an x-ray source (again, vacuum chamber?)

Materials Wishlist:

Vacuum Pump
Vacuum Chamber! (or at least the materials to make one)
NeHe Laser
Heavy Water (harder to get in Canada than I thought... united nuclear has some if someone wouldn't mind bringing me some back from the US)
Microscope
Oscilloscope
Wave Generator
LCD Projector (3d (whiteboard|display))
That paper thing... what's it called?... money?

On top of those personal projects, I'm also working with Chris Matzner on protostellar outflow-driven turbulence. To put shortly, we are basically seeing if jet's shot out by young stars, protostars, can cause a feedback effect and, by the creation of super-sonic turbulence, slow down star formation rates. I'll put a post about that sooner rather than later... I think I will start a one post per day rule to get this thing going.

Telescope Camera HOWTO v1

2008-10-22T17:55:00.000-07:00

Alright, well... after many failed attempts to start a blog I will, instead of starting and ending with an introductory post, get right into it. This will be kind of odd as I've already completed this project in full, but soon i'm going to be doing it again, so stay tuned for version two of this howto!

Anyways, I've always wanted a camera for my telescope but they were always too damn expensive! So, I did what any hacker would do: made my own. The good thing about something like this is that it can be used for any kind of optical system... in fact I will be making one for a microscope soon (and be using properly milled pieces).

So, for this project you'll need:

A webcam (any kind will do, I got a cheap $8 logitec express from canada computers)
Some sort of metal case. It should have at least an internal volume of 2.5x5x3.5cm (might be different depending on what webcam you get)
Computer fan with heat sink (see image)
Camera film tube.

(Total cost: ~$20)

And now to begin the fun! Find out how to open the webcam and do so. For my webcam there was a screw covered by a small plastic cap that held the whole unit together. You can always just smash it apart, but you don't want to destroy the delicate CCD or the circuitry!

Once open, disconnect the mic (you might want to keep it for later, it became really useful in my assembly) and take the circuitry out of the casing. Make sure the pins on the board for the mic are stowed away, some cameras are already quite compact (in the image below I had to fold the pins down, making sure not to touch anything else). Now take the lens off of the CCD (it unscrews). Be careful because you don't want the CCD to get dirty! You'll see that there is still a plastic piece where the lens used to be that is screwed onto the board (I'll call this the lens holder). Screw this off and break the glass off of it! This glass is an infrared filter and we don't want it. Make sure that everything is clean and screw this piece back on. Now you have a visible+infrared camera, but we want more!

Also, a good tip is to take the lens and glue something long to it... then you can use this as a lens cap so the CCD doesn't get dirty when you store it. Save this step for last so you can see what would best serve your needs.

Now, look at the main cable that connects the USB to the circuit board... you need to identify the ground and power. Voltmeters are useful for this (just plug in the USB cable and check where the voltage flows). Remember where this is for future reference.

Let's take care of the casing for now. First, decide what you are going to use as the front and back of the case. Screw the fan/heat sink securely to the back of the case. You want to make sure that it is secure so that you don't have any excess vibrations. Next, drill a hole in the front the size of the circle of the lens holder from before. This is where the CCD is going! Also, figure out how you want the wires to come out of the case... I drilled a hole in the side and had all the wires (USB+power to the fan) coming out of it. And finally, take that film canister and cut the base off so it becomes just a tube. Glue the lip side down onto the case so that it is centered on the hole for the CCD. I recommend using epoxy or something that gives a bit... more rigid glues (super glue and the like) seem to chip off after extensive use with things like this. And now the case is ready... almost done!

Now you want to get all the wires through the hole in the case you prepared. Solder the wires for the fan onto the correct parts of the USB wire (I found it useful to branch off the wires coming off the CCD board using the wires from the microphone so I had more room to do my work). Before reattaching the USB wires to the CCD board, apply some glue to the outside of the lens holder (being extra careful not to get any on the ccd... this is how i ruined my first CCD) and get it securely fitted into the hole in the case. Now, the final step: attach the usb cable and close the case up!

You are finally ready to take pictures! Use whatever webcam program you prefer. I personally use a script I made in python that uses openCV to control the webcam and find the FWHM (full width half maximum) of a selected region so you know everything is in focus. In a later post I will release this code! For now I'll leave you with a picture of the moon I took while testing this unit. I used gqcam to take the pictures and photoshop to put them together. They were taken off my 4.5" Astroscan telescope (not the best... in fact most distortion comes from the antique optics in the giant red beast, but when it clears up I'll take/post some pictures with my 8" Celestron!).

Hope everyone enjoyed this and please post pictures you take with the camera here! Oh, and if you have any question, just post a comment.