Haskell.hs

oneStepH98 :: Double -> (V2 Double, V2 Double) -> (V2 Double, V2 Double)
oneStepH98 h prev = (qNew, pNew)
  where
    h2 = h / 2
    hs = pure h
    h2s = pure h2
    p2 = psPrev - h2s * nablaQ qsPrev
    qNew = qsPrev + hs * nablaP p2
    pNew = p2 - h2s * nablaQ qNew
    qsPrev = P.fst prev
    psPrev = P.snd prev
    nablaQ qs = V2 (q1 / r) (q2 / r)
      where
        q1 = qs ^. L._x
        q2 = qs ^. L._y
        r   = (q1 ^ 2 + q2 ^ 2) ** (3/2)
    nablaP ps = ps

Julia.jl

function oneStep2(h :: Float64, prev :: ((Float64, Float64), (Float64, Float64))) :: ((Float64, Float64), (Float64, Float64))
    hs = (h, h)
    h2 = h / 2
    h2s = (h2, h2)
    qsPrev = prev[1]
    psPrev = prev[2]
    function nablaQ(qs)
        q1 = qs[1]
        q2 = qs[2]
        r = (q1^2 + q2^2) ^ (3/2)
        return (q1 / r, q2 / r)
    end
    function nablaP(ps)
        return ps
    end
    p2 = psPrev - h2s * nablaQ(qsPrev)
    qNew = qsPrev + hs * nablaP(p2)
    pNew = p2 - h2s * nablaQ(qNew)
    return (qNew, pNew)
end

e = 0.6
q10 = 1 - e
q20 = 0.0
p10 = 0.0
p20 = sqrt((1 + e) / (1 - e))

h = 0.1

dummyStart = ((q10, q20), (p10, p20))

function oneStep2(h, prev)
    h2 = h / 2
    qsPrev = prev[1]
    psPrev = prev[2]
    function nablaQ(qs)
        q1 = qs[1]
        q2 = qs[2]
        r = (q1^2 + q2^2) ^ (3/2)
        return (q1 / r, q2 / r)
    end
    function nablaP(ps)
        return ps
    end
    p2 = psPrev - h2 * nablaQ(qsPrev)
    qNew = qsPrev + h * nablaP(p2)
    pNew = p2 - h2 * nablaQ(qNew)
    return (qNew, pNew)
end