# PyCon 2019 Baltimore
# Stefan Karpinski
# around 3:30 in 
# https://youtu.be/kc9HwsxE1OY
# see https://juliaphysics.github.io/Measurements.jl/stable/
# 
using DifferentialEquations
using Measurements
# in julia console:
# include("julia-errorbars.jl")
#
using Plots

# 2024: DifferentialEquations and Measurements started giving error
# "ERROR: LoadError: MethodError: no method matching Float64(::Measurement{Float64})"
# which means some part is not accepting more general floats.
# Remedy: Define automatic conversion to Float64 whenever that type is expexted.
Base.Float64(x::Measurement{Float64}) = x.val

g = 9.79 ± 0.02 # Gravitational constant
L = 1.00 ± 0.01 # Length of pendulum


# initial speed and angle
u₀ = [0 ± 0, π/60 ± 0.01]
tspan = (0.0, 2π)

# Define the problem
# θ''(t) + (g/L)*θ(t) = 0
# u[2] := θ'(t) => u[2]' = θ''(t)
# insert
# u[2]' + g/L*θ(t) = 0 <=> u[2]' =  -(g/L)*θ := du[2]
function pendulum(du, u, p, t)
    θ  = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

# Pass to solvers
prop = ODEProblem(pendulum, u₀ , tspan)
sol = solve(prop, Tsit5(), reltol = 1e-6) ;

# Analytic solution
u = u₀[2] .* cos.(sqrt(g/L) .* sol.t) ;

f = plot(sol.t, u, label="analytic")
f = plot!(sol.t, sol(sol.t)[2,:] .+ .1, label="numerical (shifted up +.1) ")


# give the plot a chance to show when run from console
display(f)
readline()