# # Making animated and interactive scientific
# # visualizations in Makie.jl
# # The process
# 0. Learn about `Observable`s
# 1. Initialize simulation in a stepping manner
# 2. Initialize the `Observable`s of the animation
# 3. Plot the `Observable`s and any other static elements
# 4. Create the "animation stepping function"
# 5. Test it
# 6. Save animations to videos
# 7. Interactive application
# Let's start!
# Load packages that will be used
using DynamicalSystems
using OrdinaryDiffEq, DifferentialEquations
using GLMakie
using DataStructures: CircularBuffer
# %% 1. Initialize simulation in a stepping manner
# u0 = [0.7,0., 2.4, 0.0]
#u0=[0.7227,1.2260, 2.419, -1.2260]
#u0=[0.7186,1.4566, 2.413, -1.1351]
#u0=[3.14,0.38,0.0,-0.0001]
#u0=[0.623,0.0,-2.474,0.0]
u0=[0.08, 0.8, 0.0, -1.1]
u0 = [0.58, -0.0, 2.5100000000000002, -0.0]
u0 = [2.32, -0.0, 0.81, -0.0]
u0 = [1.3800000000000001, -0.0, 4.57, -0.0]
u0 = [-1.380, -0.0, 4.57, -0.0]
u0 = [1.3800000000000001, 0.3, 4.57, -0.8]
u0 = [0.796469714455742, -0.10345443789882522, 2.4035673348140767, -0.1455634637437679] # k=0.5
u0 = [0.7545855293320464, -0.26795226902728103, 2.4382140255471834, -0.27874765545281494] # k=3.0
##u0 =[0.8081974727949669, 0.0747264049366238, 2.5080111959945914, 0.04154800493068017] # k=3.0
#
u0=[0.08,0.8,0.0,-1.1]
#u0 = [-1.380, -0.0, 4.57, -0.0]
#u0=[-1.38,0.,4.57,0.0]
u0=[0.001,0.01,pi,0.03]
u0 = [-0.6220195889472961, 0.0, -2.857393741607666, 0.0]
# E1, 0.1 (c1, c2)
# u0 = [0.6071736812591553, 0.0, 2.23612117767334, 0.0]
# E1, 0.1 (-c1, -c2)
# u0 = [5.344040870666504, 0.0, 3.810854196548462, 0.0]
# u0 = [5.722289562225342, 0.0, 3.9998221397399902, 0.0]
# E_en=0.1
# E1, 0.2
# u0 = [0.46397456526756287, 0.0, 2.267615795135498, 0.0]
# u0 = [0.41303494572639465, 0.0, 2.3778469562530518, 0.0]
# E_en=0.2
# # E3, 0.5
# u0 = [1.208648443222046, 0.0, 2.251868486404419, 0.0]
# u0 = [5.354092597961426, 0.0, 4.326481342315674, 0.0]
# u0 = [1.1023131608963013, 0.0, 2.105553150177002, 0.0]
# E_en=0.5
#
# # E2, 0.4
# u0 = [2.818772315979004, 0.3316624790355401, 0.4919906258583069, -0.3]
# u0 = [2.4723310470581055, 0.3316624790355401, 5.707992076873779, 0.3]
# E_en=0.4
#
# #E4, 0.7
# u0 = [3.84234881401062, 0.0, 3.1311047077178955, 0.0]
# u0 = [0.440925270318985, 0.0, 0.6072171330451965, 0.0]
# E_en =0.7
#
# #E5, 0.9
# u0 = [1.9369218349456787, 0.632455532033676, 3.331843137741089, 0.0]
# E_en=0.9
#
# #E5, 1.2
# u0 = [4.062811374664307, 1.0, 3.044971227645874, 0.0]
# E_en=1.2
#
# Table 1
# u0 = [0., 0., 0., 0.] #S1
# u0 = [pi, 0., pi, 0.] #S1
# u0 = [pi, 0., 0, 0.] #S2
# u0 = [0, 0., pi, 0.] #S3
# u0 = [0.721,0.0, 2.418, 0.0] #S4
# u0 = [-0.721,0.0, -2.418, 0.0] #S4
# Table Type I
# u0 = [0.6071736812591553, 0.0, 2.23612117767334, 0.0] # TypeI_E1
# u0 = [5.344040870666504, 0.0, 3.810854196548462, 0.0] # TypeI_E1a
# u0 = [1.1023131608963013, 0.0, 2.105553150177002, 0.0] # TypeI_E3
# Table Type II
# u0 = [2.818772315979004, 0.3316624790355401, 0.4919906258583069, -0.3] # TypeII_E2
# u0 = [2.4723310470581055, 0.3316624790355401, 5.707992076873779, 0.3] # TypeII_E2a
# u0 = [0.440925270318985, 0.0, 0.6072171330451965, 0.0] # TypeII_E4
# #E5, 1.2
# u0 = [4.062811374664307, 1.0, 3.044971227645874, 0.0]
# E_en=1.2
t0 = 8000.00
step_size = 0.05
#-------------------------------
# order of r in potential
n=3
#-------------------------------
#-------------------------------
# parameters
#
K = 2.0 #exponential factor
L = 1.5 #seperation distance
α = 1.0
ξ = 0.0
A = 1.0
B = 0.0
#-------------------------------
#params = [K, L, α, ξ]
params = [K, L, α, ξ, A, B]
function rotor_2(u, params, t)
θ1,p1,θ2,p2 = u
K,L,α,ξ,A,B = params
r=sqrt(abs(α^2+1+L^2-2*L*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)));
r=r+ξ
v0 = (1-exp(-K*r^n))/r
pv0_r = -exp(-K*r^n)*(-1-n*K*r^n+exp(K*r^n))/(r^2)
# pv02_r = -2*exp(-2*K*r)*(1+K*r)/(r^3)
pr_θ1 = (α*(sin(θ1-θ2)+L*sin(θ1)))/r
pr_θ2 = (α*sin(θ2-θ1) -L*sin(θ2)) /r
factor = ( ((1/α^2)+1)*v0*B^2 + (B*((-p1/α^2)+p2)+A) ) * pv0_r
# pr_t = (L*(α*sin(θ1)*ω1-sin(θ2)*ω2)+ α*sin(θ1-θ2)*(ω1-ω2))/r
du1 = (p1-B*v0)/(α^2)
du2 = -factor*pr_θ1
du3 = (p2+B*v0)
du4 = -factor*pr_θ2
return SVector{4}(du1, du2, du3, du4)
end
diffeq = (alg = Vern7(), reltol = 1e-10, abstol = 1e-10 , atol=1e-10,rtol=1e-10, adaptive = false, dt = step_size)
rotor_system = CoupledODEs(rotor_2, u0, params; diffeq, t0 = 0.0)
integ = rotor_system.integ
# Solve diffeq with constant step for smoother curves
function xycoords(state)
θ1 = state[1]
θ2 = state[3]
x1 = α*cos(θ1)
y1 = α*sin(θ1)
x2 = L + cos(θ2)
y2 = sin(θ2)
return x1,x2,y1,y2
end
function xytime(t)
return "$(t) sec"
end
# `integ` is an integrator. `step!(integ)` progresses the integrator
# for one step. `integ.u` is the system state at current step.
# Then `xycoords` converts the states of the integrator
# to their plottable format. So we can imagine something like
function progress_for_one_step!(integ)
step!(integ)
return xycoords(integ)
end
# to be our stepping function that returns the new data.
# If we had finite data instead of a forever-running animation,
# then the "stepping function" would simply be to progress the index `i`
# of the existing data one step forwards...
# %% 2. Initialize the `Observable`s of the animation
# Here the animated elements will be: balls and rods making the
# double pendulum, and the tail (trajectory) of the pendulum.
x1,x2,y1,y2 = xycoords(u0)
rod1 = Observable([Point2f(0, 0), Point2f(x1, y1)])
rod2 = Observable([Point2f(L, 0), Point2f(x2, y2)])
balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
# Next is the observable for the tail
tail2 = 250 # length of plotted trajectory, in units of `dt`
# The circular buffer datastructure makes making stepping-based
# animations very intuitive
traj2 = CircularBuffer{Point2f}(tail2)
fill!(traj2, Point2f(x2, y2)) # add correct values to the circular buffer
traj2 = Observable(traj2) # make it an observable
tail1 = 250
traj1 = CircularBuffer{Point2f}(tail1)
fill!(traj1, Point2f(x1, y1))
traj1 = Observable(traj1)
t_value = Observable(0.0)
i_value = Observable(u0)
##
i_value = Observable(u0)
# %% 3. Plot the `Observable`s and any other static elements
# Before plotting we need to initialie a figure
fig = Figure(); display(fig)
# in my experience it leads to cleaner code if we first initialize
# an axis and populate it accordingly.
ax = Axis(fig[1,1])
# Now we plot the observables _directly_! First the pendulum
lines!(ax, rod1; linewidth = 2, color = :black)
lines!(ax, rod2; linewidth = 2, color = :black)
scatter!(ax, balls; marker = :circle, strokewidth = 2,
strokecolor = [:red,:blue],
color = [:red,:blue], markersize = [1, 1]
)
# then its trajectory, with a nice fadeout color
c2 = to_color(:blue)
tailcol2 = [RGBAf(c2.r, c2.g, c2.b, (i/tail2)^2) for i in 1:tail2]
lines!(ax, traj2; linewidth = 1, color = tailcol2)
c1 = to_color(:red)
tailcol1 = [RGBAf(c1.r, c1.g, c1.b, (i/tail2)^2) for i in 1:tail1]
lines!(ax, traj1; linewidth = 1, color = tailcol1)
# We can also plot now any other static elements
ax.title = "Coupled Rotor at $(u0)"
ax.aspect = DataAspect()
xlims!(ax, -3, +3)
ylims!(ax, -1.5, 1.5)
# %% 4. Create the "animation stepping function"
# Using the functions of step 1, we now define a function
# that updates the observables. Makie.jl understands observable
# updates and directly reflects this on the plotted elements.
function animstep!(integ, rod1, rod2, balls, traj1, traj2)
x1,x2,y1,y2 = progress_for_one_step!(integ)
rod1[] = [Point2f(0, 0), Point2f(x1, y1)]
rod2[] = [Point2f(L, 0), Point2f(x2, y2)]
balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
t_value[]=integ.t;
i_value[]=integ.u;
push!(traj1[], Point2f(x1, y1))
push!(traj2[], Point2f(x2, y2))
traj1[] = traj1[];
traj2[] = traj2[];# <- important! Updating in-place the value of an
# `Observable` does not trigger an update!
end
# %% 5. Test it
for i in 1:100
animstep!(integ, rod1, rod2, balls, traj1, traj2)
sleep(0.001)
end
function rotor_fig(u0)
#params = [k,l,α]
params = [K, L, α, ξ, A, B]
## diffeq = (alg = Vern7(), reltol = 1e-8, abstol = 1e-8 , atol=1e-8,rtol=1e-8, adaptive = false, dt = 0.5)
##
rotor_system = CoupledODEs(rotor_2, u0, params; diffeq, t0 = 0.0)
integ = rotor_system.integ
x1,x2,y1,y2 = xycoords(u0)
rod1 = Observable([Point2f(0, 0), Point2f(x1, y1)])
rod2 = Observable([Point2f(L, 0), Point2f(x2, y2)])
balls = Observable([Point2f(x1, y1), Point2f(x2, y2)])
tail2 = 250
traj2 = CircularBuffer{Point2f}(tail2)
fill!(traj2, Point2f(x2, y2))
traj2 = Observable(traj2)
tail1 = 250
traj1 = CircularBuffer{Point2f}(tail1)
fill!(traj1, Point2f(x1, y1))
traj1 = Observable(traj1)
t_value = Observable(0.0)
i_value = Observable(u0)
fig = Figure(); display(fig)
ax = Axis(fig[1,1])
lines!(ax, rod1; linewidth = 2, color = :black)
lines!(ax, rod2; linewidth = 2, color = :black)
scatter!(ax, balls; marker = :circle, strokewidth = 2,
strokecolor = [:red,:blue],
color = [:red,:blue], markersize = [15, 15]
)
# needs tailcol function
lines!(ax, traj1; linewidth = 2, color = tailcol1)
lines!(ax, traj2; linewidth = 2, color = tailcol2)
#ax.title = "Coupled Rotor $(round.(u0; digits=3))"
#ax.title = "Coupled Rotor"
# ax.title = "Coupled Rotor"
ax.aspect = DataAspect()
ax.xgridvisible = false
ax.ygridvisible = false
# ax.backgroundcolor = "skyblue"
hidedecorations!(ax, ticks=true)
xlims!(ax, -1.5, L+1.3)
ylims!(ax, -1.5, 1.5)
return fig, integ, rod1, rod2, balls, traj1, traj2
end
## %% 5.1 changing parameters live
# l = 1.5; params=[K,L,α, ξ, A, B]; set_parameters!(rotor_system,params)
# α = 1.0; params=[K,L,α, ξ, A, B]; set_parameters!(rotor_system,params)
#c1 = acos( (L^2+α^2-1)/(2*L*α) )
#c2 = acos((-L^2+α^2-1)/(2*L))
fig, integ, rod1, rod2, balls, traj1, traj2 = rotor_fig(u0)
## to make sure integrator is changed check `rotor_system`
## %% 5.2 adding label
#function labels()
# lab1 = Label(fig[2, 1],
# @lift("time = $(round($t_value,digits=5)) units"),
# halign=:left,
# tellheight= false,
# tellwidth = false,
# valign = :top)
#
# lab2 = Label(fig[2, 1],
# "Step Size = $(step_size), k=$(k), l=$(l), α=$(α)",
# halign=:right,
# tellheight= false,
# tellwidth = false,
# valign = :top)
#end
#
#function start(n)
# for i in 1:n
# animstep!(integ, rod1, rod2, balls, traj1, traj2)
# sleep(0.001)
# end
#end
## %% 5.3
## Calculating energy
# Kinetic energies
k1_energy(p) =(p[2]^2/(2*α^2))+(p[4]^2/2)
# potential
function pot_energy(p)
θ1,pθ1,θ2,pθ2 = p
K, L, α, ξ, A, B = params
# fma(a,b,c) = a*b + c
(s1,c1) = α.*sincos(θ1)
(s2, c2) = sincos(θ2)
r = let
a = fma(s1, s2, 1)
b = fma(c1, c2, a)
c = fma(L, c1 - c2, b)
d = fma(L, L, -2*c)
(sqrt ∘ abs)(d+3.0+α^2)+ξ
end
return (1-exp(-K*r^n))/r
end
## %% 6. Save animations to videos
# epsi=0.88717005;
# fig, integ, rod1, rod2, balls, traj1, traj2 = rotor_fig([c1-epsi,0,c2+epsi,0]);
fig, integ, rod1, rod2, balls, traj1, traj2 = rotor_fig(u0)
#lab1 = Label(fig[2, 1],
# @lift("time = $(round($t_value,digits=5)) units"),
# halign=:left,
# tellheight= false,
# tellwidth = false,
# valign = :top)
#
#lab3 = Label(fig[2, 1],
# @lift("u0 = $(round.($i_value,digits=5))"),
#
# lab3 = Label(fig[2, 1],
# @lift("u0 = $(round.($i_value,digits=5))
# time = $(round($t_value,digits=5)) units"),
# halign=:left,
# tellheight= false,
# tellwidth = false,
# valign = :bottom)
#
# lab2 = Label(fig[2, 1],
# "Step Size = $(step_size), K=$(K), L=$(L)",
# halign=:right,
# tellheight= false,
# tellwidth = false,
# valign = :top)
Label(fig[2,1],"$(u0)",tellwidth=false)
function saving(name)
frames = 1:1000
record(fig, "$(name)", frames; framerate = 30) do i # i = frame number
for j in 1:5 # step 5 times per frame
animstep!(integ, rod1, rod2, balls, traj1, traj2)
end
# any other manipulation of the figure here...
#
end # for each step of this loop, a frame is recorded
end
# # %% 7. Interactive application
# # Makie.jl has tremendously strong capabilities for real-time
# # interactivity. To learn all of this takes time of course,
# # and you'll need to consult the online documentation.
# # Here we will do two interactions: 1) a play/stop button
# # 2) clicking on the screen and getting a new initial condition!
#
# fig, integ, rod, balls, traj = makefig(u0)
fig, integ, rod1, rod2, balls, traj1, traj2 = rotor_fig(u0)
# The run button is actually pretty simple, we'll add it below the plot
#
#
# lab1 = Label(fig[2, 1],
# @lift("
# time = $(round($t_value,digits=5)) units"),
# halign=:left,
# tellheight= false,
# tellwidth = false,
# valign = :top)
lab3 = Label(fig[2, 1],
@lift("u0 = $(round.($i_value,digits=5))
Time = $(round($t_value,digits=5)) units
Step Size = 0.0001,
K=$(K), L=$(L), Energy=$(round(pot_energy($i_value)+k1_energy($i_value),digits=3))
"),
halign=:center,
tellheight= false,
tellwidth = false,
valign = :top)
#lab2 = Label(fig[2, 1],
## "Step Size = $(step_size), K=$(K), L=$(L)",
# "Step Size = 0.0001,
# K = $(K), L = $(L), E = $(E_en)",
# halign=:right,
# tellheight= false,
# tellwidth = false,
# valign = :top)
#
run = Button(fig[2,1],
label = "run",
halign =:center,
tellheight=false,
tellwidth = false,
valign = :bottom)
# This button will start/stop an animation. It's actually surprisingly
# simple to do this. The magic code is:
#
isrunning = Observable(false)
on(run.clicks) do clicks; isrunning[] = !isrunning[]; end
on(run.clicks) do clicks
@async while isrunning[]
isopen(fig.scene) || break # ensures computations stop if closed window
animstep!(integ, rod1, rod2, balls, traj1, traj2)
sleep(0.00001) # or `yield()` instead
end
end
ax = content(fig[1,1])
Makie.deactivate_interaction!(ax, :rectanglezoom)
# # and we'll add a new trigger using the `select_point` function:
lpoint = select_line(ax.scene)
#
# # Now we see that we can click on the screen and the `spoint` updates!
#
# # Okay, let's do something useful when it triggers
todeg_1(z) = atan(z[2],z[1])
todeg_2(z) = atan(z[2],z[1]-L)
function ini(c1,c2)
r(θ1,θ2) = sqrt(abs(α^2+1+L^2-2*L*(α*cos(θ1)-cos(θ2))-2*α*cos(θ1-θ2)));
t=r(c1,c2)+ξ
v0(r) = exp(-K*r)/r
return B*v0(t)
end
on(lpoint) do line
θ1,θ2 = todeg_1(line[1]), todeg_2(line[2])
p = ini(θ1,θ2)
u0 = [θ1,0.0,θ2,0.0]
reinit!(integ, u0)
# Reset tail and balls to new coordinates
x1,x2,y1,y2 = xycoords(u0)
# traj[] .= fill(Point2f(x2, y2), length(traj[]))
# traj[] = traj[]
# rod[] = [Point2f(0, 0), Point2f(x1, y1), Point2f(x2, y2)]
# balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
#
rod1[] = [Point2f(0, 0), Point2f(x1, y1)]
rod2[] = [Point2f(L, 0), Point2f(x2, y2)]
balls[] = [Point2f(x1, y1), Point2f(x2, y2)]
t_value[]=integ.t;
i_value[]=integ.u;
traj1[] .= fill(Point2f(x1, y1), length(traj1[]))
traj2[] .= fill(Point2f(x2, y2), length(traj2[]))
traj1[] = traj1[];
traj2[] = traj2[];
end