George's double pendulums

pendulums.jl

using DynamicalSystems, Luxor

function double_pendulum(u0=rand(4);
        G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
    @inline @inbounds function eom_dp!(du, state)
    du[1] = state[2]
    gm = 0.01
    del_ = state[3] - state[1]
    den1 = (M1 + M2)*L1 - M2*L1*cos(del_)*cos(del_)
    du[2] = (M2*L1*state[2]*state[2]*sin(del_)*cos(del_) +
        M2*G*sin(state[3])*cos(del_) +
        M2*L2*state[4]*state[4]*sin(del_) -
        (M1 + M2)*G*sin(state[1]) - gm * state[2]) /den1
    du[3] = state[4]
    den2 = (L2/L1)*den1
    du[4] = (-M2*L2*state[4]*state[4]*sin(del_)*cos(del_) +
        (M1 + M2)*G*sin(state[1])*cos(del_) -
        (M1 + M2)*L1*state[2]*state[2]*sin(del_) -
        (M1 + M2)*G*sin(state[3]) - gm * state[4])/den2
    end
    return ContinuousDS(u0, eom_dp!, nothing)
end

# convert a timeseries to two trajectories
function maketrajectories(tseries, pendulumlength)
    innertrajectory = Point[]
    outertrajectory = Point[]
    @inbounds for i in 1:length(xv)
        angle1 = π/2 + xv[i][1]
        angle2 = π/2 + xv[i][3]
        pos1 = Point(pendulumlength * cos(angle1), pendulumlength * sin(angle1))
        push!(innertrajectory, pos1)
        push!(outertrajectory, pos1 + Point(pendulumlength * cos(angle2), pendulumlength * sin(angle2)))
    end
    return (innertrajectory, outertrajectory)
end

# draw a colored disk Julia-style
function drawdisk(pos, outerrad, innerrad, colpair)
    setopacity(1)
    sethue(colpair[1])
    circle(pos, outerrad, :fill)
    sethue(colpair[2])
    circle(pos, innerrad, :fill)
end

function drawrod(startpos, endpos, col, radius, linewidth)
    rodlength = norm(startpos, endpos)
    setline(linewidth)
    sethue("black")
    line(between(startpos, endpos, radius/rodlength),
         between(startpos, endpos, 1 - radius/rodlength),
         :stroke)
end

function drawtrace(pgon, fnumber, tracelength)
    traceend = max(fnumber - tracelength, 1)
    @layer begin # work backwards from end, drawing as far as traceend
        sethue(darker_blue)
        setline(3)
        @inbounds for i in fnumber:-1:traceend+1
            setopacity(rescale(i, fnumber, traceend, 1, 0))
            line(pgon[i], pgon[i-1], :stroke)
        end
    end
end

function frame(scene, framenumber)
    background("white")
    innerpos = innertrajectory[framenumber]
    outerpos = outertrajectory[framenumber]
    drawtrace(outertrajectory, framenumber, tracelength)
    drawrod(O,         innerpos, "black", markerouterradius, rodthickness)
    drawrod(innerpos,  outerpos, "black", markerouterradius, rodthickness)
    drawdisk(O,        markerouterradius, markerinnerradius, greenpair)
    drawdisk(innerpos, markerouterradius, markerinnerradius, redpair)
    drawdisk(outerpos, markerouterradius, markerinnerradius, purplepair)
end

# the Julia colors:
darker_blue = (0.251, 0.388, 0.847); lighter_blue = (0.4, 0.51, 0.878)
darker_purple = (0.584, 0.345, 0.698); lighter_purple  = (0.667, 0.475, 0.757)
darker_green  = (0.22, 0.596, 0.149); lighter_green  = (0.376, 0.678, 0.318)
darker_red  = (0.796, 0.235, 0.2); lighter_red  = (0.835, 0.388, 0.361)
bluepair = [darker_blue, lighter_blue]
purplepair = [darker_purple, lighter_purple]
greenpair = [darker_green, lighter_green]
redpair = [darker_red, lighter_red]

############################################################################
#        it don't mean a thing, if it aint got that swing
############################################################################

dp1 = double_pendulum(
    [deg2rad(-120), deg2rad(60), deg2rad(9), deg2rad(-50)],
    G=9.81, L1 = 0.85, L2 = 0.85, M1 = 1.0, M2 = 1.0)

const ΔT = 0.01
const DURATION = 15 # seconds
const rodthickness = 5
const tracelength = 600

xv = timeseries(dp1, DURATION, dt=ΔT)

pendulums = Movie(800, 600, "pendulum")
diskradius = 30
pendulumlength = 3diskradius
markerouterradius = 35
markerinnerradius = 30
innertrajectory, outertrajectory = maketrajectories(xv, pendulumlength)
animate(pendulums, [
    Scene(pendulums, frame, 1:length(xv))
    ],
    creategif=true, 
    pathname="/tmp/dynamicalpendulums-earth.gif")