George's double pendulums

pendulums.jl

using Luxor, DynamicalSystems, Colors

function frame(scene, framenumber, datas)
    background("black")
    # datas contains two datasets
    # each dataset is an array of sarrays of 4 floats
    d1, d2 = datas
    offset = π / 2
    L = 200
    disksize = 15

    for i in 1:framenumber
        p1, p2, p3, p4 = d1[i]
        q1, q2, q3, q4 = d2[i]

        @layer begin
            translate(boxtopcenter() + (0, 100))
            setopacity(rescale(i, 1, scene.framerange.stop, 0, 1))
            sethue(HSV(10, 0.4, 0.4))
            line(O, O + polar(L, offset + p1), :stroke)
            circle(O + polar(L, offset + p1), disksize, :fill)

            @layer begin
                sethue(HSV(120, 0.4, 0.4))
                translate(O + polar(L, offset + p1))
                line(O, O + polar(L, offset + p3), :stroke)
                circle(O + polar(L, offset + p3), disksize, :fill)
            end
        end
    end
    # find the current
    p1, p2, p3, p4 = d1[framenumber]
    q1, q2, q3, q4 = d2[framenumber]

    setline(5)
    @layer begin
        translate(boxtopcenter() + (0, 100))
        sethue(HSV(10, 0.9, 0.9))
        line(O, O + polar(L, offset + p1), :stroke)
        circle(O + polar(L, offset + p1), disksize, :fill)

        sethue(HSV(120, 0.9, 0.9))
        translate(O + polar(L, offset + p1))
        line(O, O + polar(L, offset + p3), :stroke)
        circle(O + polar(L, offset + p3), disksize, :fill)
    end
end

function main()
    W, H = 800, 500
    movie = Movie(W, H, "domain")
    tempdirectory = "/tmp/temp/"
    if isdir(tempdirectory)
        rm(tempdirectory, recursive=true, force=true)
    end
    mkdir(tempdirectory)

    # double_pendulum 
    # G=10.0, L1 = 1.0, L2 = 1.0, M1 = 1.0, M2 = 1.0)
    # gravity (`G`), lengths of each rod (`L1` and `L2`) 
    # mass of each ball (`M1` and `M2`)
    M = M1 = M2 = 50
    u0 = [π / 2, 0, 0, 0.5] # initial 
    u0s = [u0 .+ [0, 0, 0.005 * i, 0] for i in 0:M-1] # initial ? 
    L1 = 200
    L2 = 200
    dp = Systems.double_pendulum(u0; L1, L2, M1, M2)
    # make trajectory - 20 time, u0s conditions, step by 0.1
    datas = trajectory.(Ref(dp), 20, u0s; Δt=0.01)

    # datas is two datasets  
    # length of animation is length of data
    frames = length(datas[1])
    animate(movie, [
            Scene(movie, (s, f) -> frame(s, f, datas), 1:frames)
        ],
        creategif=true,
        framerate=30,
        tempdirectory=tempdirectory)
    #FFMPEG.ffmpeg_exe(`-r 25 -f image2 -i $(tempdirectory)/%10d.png -c:v libx264 -pix_fmt yuv420p -y /tmp/animation.mp4`)
end

main()

Just an addendum: now the double pendulum is part of the predefined systems of the package, so one only needs to do:

using DynamicalSystems
dp1 = Systems.doublependulum(args...)

edit: forgot to say how AWESOME this is.