pendulum logo

pendulum.jl

using Luxor

using ODE # I should be using OrdinaryDiff, but that's for another day

function pendulum(t, y)
   Y = [6 * (2 * y[3] - 3 * cos(y[1] - y[2]) * y[4])/(16 - 9 * cos(y[1] - y[2])^2);
        6 * (8 * y[4] - 3 * cos(y[1] - y[2]) * y[3])/(16 - 9 * cos(y[1] - y[2])^2)]
   [Y[1];
    Y[2];
    - (Y[1] * Y[2] * sin(y[1] - y[2]) + 3 * sin(y[1]))/2;
    - (sin(y[2]) - Y[1] * Y[2] * sin(y[1] - y[2]))/2;]
end

initial = [-π/6, π/2, 0.5, 0.9]

T, xv = ode23(pendulum, initial, 0.0:0.03:300, points=:specified);

"""
    draw a pendulum between pos1 and pos2, in colorname, with diskradius 
"""
function drawpendulum(pos1, pos2, colorname, diskrad)
    rodlength = distance(pos1, pos2)
    setlinecap("round")
    setline(3)
    sethue("gold3")
    line(between(pos1, pos2, diskrad/rodlength), between(pos1, pos2, 1 - diskrad/rodlength), :stroke)

    sethue(colorname)
    setline(3)
    circle(pos2, diskrad, :fill)
    @layer begin
        setopacity(0.9)
        circle(pos2, diskrad, :stroke)
    end
end

# for trails, don't worry about drawing pos1 or rod, just do last disk
function drawpendulumtrail(pos1, pos2, colorname, diskrad)
    rodlength = distance(pos1, pos2)
    sethue("white")
    circle(pos2, diskrad, :fill)
end

moviewidth = 128 # square!
pendulums = Movie(moviewidth, moviewidth, "pendulum")

function backdrop(scene, framenumber)
    background("white")
    sethue("grey5")
    bw = blend(O, 20, Point(0, 0), 80,  "grey30", "grey5")
    setblend(bw)
    squircle(O, moviewidth ÷ 2 - 8, moviewidth ÷ 2 - 8, :fill)
end

function frame(scene, framenumber)
    translate(0, -10)
    diskradius = moviewidth/20
    pendulumlength = 28

    # draw the colored disk at the middle
    sethue(cols[1])
    setline(3)
    circle(O, diskradius, :fill)
    setopacity(0.9)
    circle(O, diskradius, :stroke)

    # draw the trails, by going through all previous positions
    # inefficient, but whatever
    for prev in 1:framenumber
        angle1 = π/2 + xv[prev][1]
        angle2 = π/2 + xv[prev][2]
        pos1 = Point(pendulumlength * cos(angle1), pendulumlength * sin(angle1))
        pos2 = pos1 + Point(pendulumlength * cos(angle2), pendulumlength * sin(angle2))
        setopacity(0.7)
        drawpendulumtrail(pos1, pos2, cols[3], 1)
    end

    # draw the pendulums already
    angle1 = π/2 + xv[framenumber][1]
    angle2 = π/2 + xv[framenumber][2]

    # start and end points for pendulum rods 1 and 2
    pos1 = Point(pendulumlength * cos(angle1), pendulumlength * sin(angle1))
    pos2 = pos1 + Point(pendulumlength * cos(angle2), pendulumlength * sin(angle2))
    setopacity(1)
    drawpendulum(O,    pos1, cols[2], diskradius)
    drawpendulum(pos1, pos2, cols[3], diskradius)
end

cols = ["forestgreen", "mediumorchid3", "brown3", "royalblue3"]

totaltime = 2000 # number of frames

animate(pendulums, [
    Scene(pendulums, backdrop,  1:totaltime),
    Scene(pendulums, frame,  1:totaltime)
    ],
    creategif=true,
    pathname="/tmp/pendulums.gif")

20 frames looks like this