LissRK4.py

#!/usr/bin/env python

import sys, getopt
from math import sin, cos, pi

''' Plot a Lissajous figure, saving as PBM

    Written by PM 2Ring March 2009
    2015.03.31 Updated to use approximately constant distance
    between succesive points, using Runge-Kutta integration.


    Parametric equations of a Lissajous figure:
    x = cos(a*th), y = sin(b*th)

    Derivatives:
    dx/dth = -a*sin(a*th), dy/dth = b*cos(b*th)

    Derivative squared of arc length wrt to th
    (ds/dth)^2 = a^2*sin^2(a*th) + b^2*cos^2(b*th)

    v = ds/dt = ds/dth * dth/dt
    dth/dt = v / (ds/dth)
    th = integral(v / (ds/dth)) dt
     = integral(v / sqrt(a^2*sin^2(a*th) + b^2*cos^2(b*th))) dt
'''

class bitmap:
    ''' Basic pixel array, with PBM output '''
    def __init__(self, width, height):
        self.width, self.height = width, height

        #Pixel buffer. 0=white, 1=black.
        row = ['0']*width + ['\n']
        self.pix = [row[:] for i in xrange(height)]

    #Put pixel in buffer.
    def put(self, x, y, c=1):
        self.pix[y][x] = '01'[c]

    def savepbm(self, fname):
        with open(fname, 'w') as f:
            f.write('P1\n%d %d\n' % (self.width, self.height))
            #Technically, each row should be no more than 70 chars
            for row in self.pix:
                #Space is permitted, but not required
                f.write(''.join(row))

def liss(a, b, speed, width, height, fname):
    ''' Generate a Lissajous figure with approximately
        constant distance between successive points
        using Runge-Kutta integration
    '''
    #Pixel buffer
    pix = bitmap(width, height)

    #Origin
    ox, oy = width // 2, height // 2
    rad = min(ox, oy) - 3

    def speedfunc(theta):
        s = (a * sin(a * theta))**2 + (b * cos(b * theta))**2
        return speed / (rad * s ** 0.5)

    maxtheta = 2. * pi
    theta = 0.0
    while theta < maxtheta:
        x = rad * cos(a * theta)
        y = rad * sin(b * theta)
        pix.put(int(.5 + ox + x), int(.5 + oy + y))

        #Calculate Runge-Kutta 4 increments
        k1 = speedfunc(theta)
        k2 = speedfunc(theta + 0.5 * k1)
        k3 = speedfunc(theta + 0.5 * k2)
        k4 = speedfunc(theta + k3)

        dtheta = (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
        theta += dtheta

    pix.savepbm(fname)

def main():
    #Default args
    fname = 'liss.pbm'  #Bitmap filename.
    width = 512
    height = 512
    xfreq = 3.0
    yfreq = 2.0
    speed = 1.0

    def usage():
        s = """Plot a Lissajous figure, saving as PBM.
Usage: %s [-h] [-x x_frequency=%3f] [-y y_frequency=%3f] [-s speed=%f]
[-w width=%d] [-v vertical_height=%d] [-f filename='%s']"""
        print s % (sys.argv[0], xfreq, yfreq, speed, width, height, fname)
        raise SystemExit, 1

    try:
        opts, args = getopt.getopt(sys.argv[1:], "hx:y:s:w:v:f:")
    except getopt.GetoptError, e:
        print e.msg
        usage()

    for o, a in opts:
        if o == "-h": usage()
        elif o == "-x": xfreq = float(a)
        elif o == "-y": yfreq = float(a)
        elif o == "-s": speed = float(a)
        elif o == "-w": width = int(a)
        elif o == "-v": height = int(a)
        elif o == "-f": fname = a

    liss(xfreq, yfreq, speed, width, height, fname)

if __name__ == "__main__":
    main()