Double pendulum

Inspired by Mark M's gravity simulation, I went ahead and made a simple double-pendulum simulation, that might be fun to play with. Here it is:

size(440,440)
speed(30)
 
from Numeric import *
 
def setup():
    global lineset,z
    global L1,L2,M1,M2,G
    global CX,CY
    global STEP
    L1,L2 = 100,100 #arm lengths (pixels)
    M1,M2 = 1.0,1.0 #masses
    a1 = 3.0 #initial angle of first arm
    a2 = 2.0 #initial angle of second arm 
    G = 150 #Gravity constant
    STEP = 0.1 #step size for RK4-integrator
    CX,CY = 220,220  #where the pendulum is attached
    z = a1,0.0,a2,0.0 #initial state.
    lineset = [] #keeps track of old points
    #build initial set of old points:
    for i in range(10):
        z = nextstep(z)
        lineset.append(points(z))
 
def draw_fireflies():
    """
    Fancy drawing of the current
    and old states of the pendulum
    """
    global lineset,z
    z = nextstep(z) #next state
    lineset = lineset[1:] + [points(z),] #update set of old states
    rect(0,0,500,500) #black background
    m = n = float(len(lineset)+1) #used to caluclate alpha
    strokewidth(20)
    for x1,y1,x2,y2 in lineset: 
        m-=1.0 # states drawn weaker the older they are
        nostroke()
        fill(1,1,1,1-m/n) #older=more alpha
        oval(x1-10,y1-10,20,20)
        oval(x2-10,y2-10,20,20)
        
def draw_classic():
    """
    Classic drawing of the pendulum state
    """
    global z,CX,CY
    strokewidth(2)
    stroke(0,0,0,1)
    nofill()
    z = nextstep(z)
    x1,y1,x2,y2 = points(z)
    oval(CX-10,CY-10,20,20)
    oval(x1-10,y1-10,20,20)
    oval(x2-10,y2-10,20,20)
    line(CX,CY,x1,y1)
    line(x1,y1,x2,y2)
 
draw = draw_classic
#draw = draw_fireflies
 
def derivs(z):
    """
    The double pendulum equations of motion.
    gives derivatives for angles and angular
    velocities, given the current angle and 
    angular velocities.
    """
    global M1,M2,L1,L2,G
    a1,w1,a2,w2 = z
    da1 =  w1
    da2 =  w2
    dw1 =  (M2*L1*w1**2*sin(a2-a1)*cos(a2-a1) + M2*G*sin(a2)*cos(a2-a1) + M2*L2*w2**2*sin(a2-a1) - (M1+M2)*G*sin(a1))/((M1+M2)*L1-M2*L1*cos(a2-a1)**2)
    dw2 =  (-M2*L2*w2**2*sin(a2-a1)*cos(a2-a1) + (M1+M2)*(G*sin(a1)*cos(a2-a1) - L1*w1**2*sin(a2-a1)-G*sin(a2)))/((M1+M2)*L2-M2*L2*cos(a2-a1)**2)
    return array((da1,dw1,da2,dw2))
 
def nextstep(z): 
    """
    4th order Runge-Kutta to give next state
    """
    global STEP
    k1 = derivs(z)
    k2 = derivs(z+STEP/2*k1)
    k3 = derivs(z+STEP/2*k2)
    k4 = derivs(z+STEP*k3)
    return z + STEP/6.0*(k1 + 2*k2 + 2*k3 + k4)
 
def points(z):
    """
    From state of angles, calculate 
    coordinates of arm positions
    """
    global L1,L2,CX,CY
    a1,w1,a2,w2 = z
    x1 = L1*sin(a1)+CX
    y1 = L1*cos(a1)+CY
    x2 = L2*sin(a2) + x1
    y2 = L2*cos(a2) + y1
    return x1,y1,x2,y2

Nice work. I see you're using Runge-Kutta. I was looking for a comprehensible example to add to a physics library. Would you mind if I use your code as a basis? Let me see if I can post some code of what I'm working on...

Go ahead and use whatever you feel like! Just glad if it's useful. Looking forward to seeing some more code! :)