# Using GLMakie to animate crank mechanism

I am trying to animate a crank mechanism using the GLMakie package. I was following this youtube tutorial, Making animations and interactive applications in Makie.jl - YouTube , and here’s the code I wrote.

``````using GLMakie

crank = 3.3

fixed = 6.8

coupler = 2

output = 9.3

θ1 = (61/90)*pi

k1 = fixed/crank

k2 = fixed/output

k3 = ((output^2 - coupler^2 + output^2 + fixed^2)/(2 *crank* output))

P1 = [0;0]

p4x = fixed*-cos(θ1)

p4y = fixed*-sin(θ1)

function mysqrt(x)

if x < 0

0.5

else

sqrt(x)

end

end

θ2 = Array(0:0.05:2*pi)

function xycoordinates(θ2)

A = (1-k2)*cos(θ2) - k1 + k3

B = -2*sin(θ2)

C = k1 - (1 + k2)*cos(θ2 + k3)

b = mysqrt(B^2 - 4A*C)

θ4 = 2*(atan((-B - b)/2A, (-B + b)/2A))

p2x = crank*cos(θ2)

p2y = crank*sin(θ2)

p3x = output*-cos(θ4)

p3y = output*sin(θ4)

return p2x, p2y, p3x, p3y

end

function progress_for_one_step!(θ2)

step!(θ2)

u = θ2.u

return xycoordinates(u)

end

p2x, p3x, p2y, p3x = xycoordinates(θ2)

mechanism = Observable([Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)])

fig = Figure(); display(fig)

ax = Axis(fig[1,1])

lines!(ax, mechanism; linewidth=4, color = :purple)

ax.title = “crank mechanism”

function animstep!(θ2, mechanism)

p2x, p3x, p2y, p3x = progress_for_one_step!(θ2)

mechanism[] = [Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)]

end

for i in 1:1000

animstep!(θ2, mechanism)

sleep(0.001)

end
``````

Almost every line of the code runs except the line “p2x, p3x, p2y, p3x = xycoordinates(θ2)” which prevents the plot from animating. I think the problem could be linked to my “progress_for_one_step” function and how my data is presented in the program. Could someone help me in identifying and correcting the error? Thank you.

Your function `xycoordinates` seems to be intended to operate on a scalar, but here

``````p2x, p3x, p2y, p3x = xycoordinates(θ2)
``````

you are calling it with `θ2` being an array. Therefore the first error I see is

``````ERROR: MethodError: no method matching cos(::Vector{Float64})
``````
1 Like

The idea was to calculate `p2x, p3x, p2y, p3x` for all the values of `θ2`. Maybe a for loop would have worked but I don’t see how I can write the `progress_for_one_step` with it, and also the other following functions and the animations.

Then you probably have to adapt `xycoordinates` to use broadcasting? The first line should look something like

``````    A = (1-k2)*cos.(θ2) .- k1 .+ k3
``````

for example.

Edit: or the simpler solution: broadcast `xycoordinates` like

``````p2x, p3x, p2y, p3y = xycoordinates.(θ2)
``````
1 Like

This solves the method error problem. However when I print `p2x, p3x, p2y` or `p3y` it returns only a 4 element vector.

`````` p3y
(3.262944557188939, 0.49314583716287747, -8.191358441435042, -4.40359476836038)
``````

Shouldn’t each return a 126 element vector?

It’s a bit confusing cause this ` p2x, p3x, p2y, p3y = xycoordinates.(θ2)` returns the following,

``````126-element Vector{NTuple{4, Float64}}:
(3.3, 0.0, -5.694607616035193e-16, -9.3)
(3.2958758593033886, 0.16493125859323848, -3.575547441280739, -8.585188436787558)
(3.283513745417485, 0.3294502749345329, -6.406050177980947, -6.741848494084563)
(3.262944557188939, 0.49314583716287747, -8.191358441435042, -4.40359476836038)
(3.234219706876097, 0.655608791623702, -9.059692957512048, -2.10046745168941)
(3.1974109916451274, 0.8164330655399257, -9.299493352799107, -0.09707410182553187)
(3.1526104141144997, 0.9752166819824206, -9.17138658677134, 1.541320237974818)
(3.09992995239635, 1.1315627646029895, -8.854311906471532, 2.844496557023141)
(3.0395012802095205, 1.2850805296185466, -8.455865787554524, 3.8714769510956564)
(2.9714754377638335, 1.4353862625670597, -8.035875190729774, 4.6813149775478395)
(2.89602245423823, 1.5821042773938698, -7.625577377449992, 5.32358616540795)
(2.8133309227963688, 1.7248678554711754, -7.240040434941139, 5.837106689141233)
(2.723607529201938, 1.863320162203617, -6.885502387781394, 6.251388395217235)
⋮
(2.7545521899692273, -1.817262290572204, 7.001109501653403, 6.121639138813921)
(2.8419349766330186, -1.6773209557475526, 7.3664180702933395, 5.676784724970269)
(2.9222144059063537, -1.5331871920653966, 7.76130731451458, 5.123681173693625)
(2.9951898210247645, -1.3852212588715342, 8.177158413287465, 4.429907480298212)
(3.060678821455318, -1.2337929939397787, 8.594493631737825, 3.5531224597553637)
(3.1185177188029507, -1.079280889070144, 8.974161142141433, 2.440170443817147)
(3.1685619459462075, -0.9220711440564553, 9.242755549872747, 1.0302765867943209)
(3.210686418379762, -0.7625567013879915, 9.271268734381575, -0.730462904519232)
(3.2447858468605295, -0.6011362640979135, 8.852712283402793, -2.849471043420044)
(3.2707750005759477, -0.4382132992133052, 7.7056524917506914, -5.2070067866650955)
(3.2885889201766174, -0.2741950292977381, 5.566920703166823, -7.449791532966048)
(3.298183080140848, -0.1094914146069375, 2.4258061282479093, -8.97805461267389)
``````

I see, sorry. You probably want this

``````coords = xycoordinates.(θ2)
p2x = map(t -> t, coords)
p3x = map(t -> t, coords)
p2y = map(t -> t, coords)
p3y = map(t -> t, coords)
@show p2x p3x p2y p3y
``````

then?

Edit: also found this nice solution

``````unzip(a) = map(x->getfield.(a, x), fieldnames(eltype(a)))
p2x, p3x, p2y, p3y = unzip(xycoordinates.(θ2))
``````
1 Like

Thanks I got the difference between the two. I thought the second method would work but all hit an error in this line
`mechanism = Observable([Point2f0(0,0), Point2f0(p2x,p2y), Point2f0(p3x,p3y), Point2f0(p4x,p4y)])`
the error is:

``````ERROR: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float32
Closest candidates are:
convert(::Type{T}, ::ColorTypes.Gray24) where T<:Real at C:\Users\ogwen\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:114
convert(::Type{T}, ::ColorTypes.Gray) where T<:Real at C:\Users\ogwen\.julia\packages\ColorTypes\6m8P7\src\conversions.jl:113
convert(::Type{T}, ::Static.StaticFloat64{N}) where {N, T<:AbstractFloat} at C:\Users\ogwen\.julia\packages\Static\A5kpr\src\float.jl:26
...
Stacktrace:
 macro expansion
@ C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\util.jl:11 [inlined]
 convert_ntuple
@ C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\util.jl:7 [inlined]
 Point{2, Float32}(x::Tuple{Vector{Float64}, Vector{Float64}})
@ GeometryBasics C:\Users\ogwen\.julia\packages\GeometryBasics\WMp6v\src\fixed_arrays.jl:19
 Point{2, Float32}(::Vector{Float64}, ::Vararg{Vector{Float64}})
@ StaticArrays C:\Users\ogwen\.julia\packages\StaticArrays\OWJK7\src\convert.jl:4
 top-level scope
@ REPL:1
``````

It used to run past this point though, don’t know what changed.

I noticed that `Point2f0` isn’t available for me. `p2...` and `p3...` are vectors here (as they should?) and you probably need some broadcasting again

``````mechanism = Observable([Point2f(0,0), Point2f.(p2x,p2y), Point2f.(p3x,p3y), Point2f(p4x,p4y)])
``````

?

1 Like

Thanks, that got it through. There however seems to be some conversion issue with GLMakie. This line is meant to plot the lines but generates an error.

`````` lines!(ax, mechanism; linewidth=4, color = :purple)
ERROR: `Makie.convert_arguments` for the plot type Lines{Tuple{Vector{AbstractVector}}} and its conversion trait PointBased() was unsuccessful.

The signature that could not be converted was:
::Vector{AbstractVector}

Makie needs to convert all plot input arguments to types that can be consumed by the backends (typically Arrays with Float32 elements).
You can define a method for `Makie.convert_arguments` (a type recipe) for these types or their supertypes to make this set of arguments convertible (See http://makie.juliaplots.org/stable/documentation/recipes/index.html).

Alternatively, you can define `Makie.convert_single_argument` for single arguments which have types that are unknown to Makie but which can be converted to known types and fed back to the conversion pipeline.
``````

According to this error, I think `p2...` and `p3...` should be arrays of `Float32` type.

We can fix this with splatting

``````mechanism = Observable([Point2f(0,0), Point2f.(p2x,p2y)..., Point2f.(p3x,p3y)..., Point2f(p4x,p4y)])
``````
1 Like

I just looked up the script for the demo: I believe we are doing it wrong…

1 Like

I suspected so. My interpretation must have been wrong. Especially in the implementation of the `xycoordinates` and `progress_for_one_step` functions.

In the first step, we’re supposed to initialize the simulation in a stepping manner, but since his problem is different from what we’re solving I don’t quite understand how he achieves that. I can’t draw the boundary between when he is just solving the diff eqn, or when he is initializing it in a stepping manner.
In our case `θ2` is the known variable from which we’re analytically solving for `θ4` and hence `p2...` and `p3...`

OK, the following program does something

``````using GLMakie

crank = 3.3
fixed = 6.8
coupler = 2
output = 9.3
θ1 = (61/90)*pi
k1 = fixed/crank
k2 = fixed/output
k3 = ((output^2 - coupler^2 + output^2 + fixed^2)/(2 *crank* output))
P1 = [0;0]
p4x = fixed*-cos(θ1)
p4y = fixed*-sin(θ1)

function mysqrt(x)
if x < 0
0.5
else
sqrt(x)
end
end

θ2 = Array(0:0.05:2*pi)

function xycoordinates(θ2)
A = (1-k2)*cos(θ2) - k1 + k3
B = -2*sin(θ2)
C = k1 - (1 + k2)*cos(θ2 + k3)
b = mysqrt(B^2 - 4A*C)
θ4 = 2*(atan((-B - b)/2A, (-B + b)/2A))
p2x = crank*cos(θ2)
p2y = crank*sin(θ2)
p3x = output*-cos(θ4)
p3y = output*sin(θ4)
return p2x, p2y, p3x, p3y
end

fig = Figure(); display(fig)
ax = Axis(fig[1,1])
ax.title = "crank mechanism"

p2x, p3x, p2y, p3y = xycoordinates(θ2)
mechanism = Observable([Point2f(0,0), Point2f(p2x,p2y), Point2f(p3x,p3y), Point2f(p4x,p4y)])
lines!(ax, mechanism; linewidth=4, color = :purple)
for i in 2:length(θ2)
global p2x, p3x, p2y, p3x = xycoordinates(θ2[i])
mechanism[] = [Point2f(0,0), Point2f(p2x,p2y), Point2f(p3x,p3y), Point2f(p4x,p4y)]
sleep(0.001)
end
``````

Edit: fixed some mistakes…

1 Like

Yes!!!
This is exactly what I wanted. Thank you so much, I’m deeply grateful for your help. Could you mind explaining the logic behind how you did it?

Unfortunately there was not much logic to it, only experience and experimentation…

2 Likes

@goerch

I am trying the code but it only show axis. I was hoping to see real crank there…

Maybe the code or the package installed in my Julia is the one at fault?

Hi, I just wanted to add that this approach can also be used to record videos of moving objects using GLMakie. I have been searching for a way to do that for several hours. And I made the following adaptation of your code to record a video of a moving circle. I did not find any other way to do this with GLMakie, as the old plots would always stick around. Another advantage of this approach is that it is possible to have parts of the plot that are not moving only be plotted once, and then have a moving part on top.

``````using GLMakie

# GLMakie
fig = Figure(resolution = (1920, 1080))
ax1 = fig[1, 1] = Axis(fig,
aspect=1.0,
limits = (-2, 2, -2, 2))

# What to plot
frames = 100 # frames
t_range = range(0.0, 10.0, length=frames) # time vector
theta = range(0.0 ,2 * pi ,length=150) # angle

# circle function
# define circle (parametric)
a = radius * cos.( theta )
b = radius * sin.( theta )
return a, b
end

# set everything up
a, b = circleFun(theta,radius) # call function
mechanism = Observable([Point2f(a[i],b[i]) for i in eachindex(a)]) # make mechanism
lines!(ax1, mechanism, linewidth=5, color=:red) # plot initial
frame_iterator = 1:1:frames # frame iterator

# record video of moving circle
record(fig, "movingCircle.mp4", frame_iterator; framerate = 12) do i