### A Pluto.jl notebook ###
# v0.19.35
using Markdown
using InteractiveUtils
# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
quote
local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
local el = $(esc(element))
global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
el
end
end
# ╔═╡ 448c89b6-7ae7-42ce-b019-ef9019d3e72a
begin
import Pkg
Pkg.activate("/home/hisacro/.julia/environments/v1.9/Project.toml")
using PlutoUI, DelimitedFiles, DynamicalSystems, GLMakie, Plots, OrdinaryDiffEq, ProgressMeter, LaTeXStrings, FFTW, XLSX;
#plotly(size=(450,400))
gr(size=(400,400))
using CairoMakie
end
# ╔═╡ 17342ecc-7da0-43dc-ba88-38bb85e44eca
md"""
# Abstract
Evolving the system by numerical integration. An initial check has been performed for energy deviation from the initial condition till the end of evolution.
This is followed by time series plots of $\theta_1, \theta_2$. There is also plot for verifying synchronization, details can be found on <doi>.
# Goal
Analyze times series of system variables, correlate the peaks.
# Methodology
1. Integrate the system to find the opt algorithm. After trial and error, Verner7 with the time step of $10^{-4}$ and absolute tolerance of $10^{-7}$ seems to have less deviation compared to RK4, RK6, and computationally viable on the testing machine.
2. Plot the evolved variables. Here twin axis plot is used for $\theta_1, \theta_2$.
3. Relative phases of $\theta_1, \theta_2$ are compared.
4. Check for synchronization is performed.
"""
# ╔═╡ 6d08c01e-8c5c-11ee-06e1-47df2e37668b
md"""
# Coupled Oscillator
"""
# ╔═╡ bf26074c-53a9-4df7-aee4-394dbecf38ef
TableOfContents()
# ╔═╡ 1d4cf734-84a2-4adc-ae86-b02f1601f6f2
md"""
# Problem Definition
"""
# ╔═╡ 4015e194-6e7a-4a5c-ae06-1406e6aadfed
md"""
rotor system with out of place function
"""
# ╔═╡ 3d223e16-ceb4-4cca-95a4-784da1d2dcf3
md"""
rotor system out of place function to counter trignometeric floating errors
"""
# ╔═╡ 84f64684-e2af-4ead-adab-cea67ddbaada
function rotor_1(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 = exp(-K*r)/r
pv0_r = -exp(-K*r)*(1+K*r)/(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/α^2)+1)*pv02_r*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
# ╔═╡ 9ee65bfd-94ae-495d-aff0-0ee2dc321a45
md"""
rotor system in place function
"""
# ╔═╡ 41336926-46a0-4c41-a87a-720162191d27
md"""
# Problem Solving
"""
# ╔═╡ a7d1c54b-6e6e-4331-b76c-f479130acc08
md"""
### Algorithm to Solve
"""
# ╔═╡ eb3ccf92-8fd6-4974-ac64-e9237ce25901
algo = Vern7();
# ╔═╡ 616557c2-4839-4257-b331-188c04d17ce3
diffeq = (alg = algo, reltol = 1.0e-7, abstol = 1.0e-7, atol = 1.0e-7, rtol = 1.0e-7, )
# ╔═╡ 249672c9-d9f2-4900-80bf-70b53c370e63
md"""
### Fixed Points as Initial Condition
"""
# ╔═╡ ed3c93f7-d35c-4c45-86f7-b3e1f7f30fed
md"""
The Trivial Fixed points are at $$(0,0),(0,\pm\pi),(\pm\pi,0),(\pm\pi,\pm\pi)$$
and Other 2 Critical Points
$(c1, c2)$ , $(-c1, -c2)$ for $$|1-α|\le L \le 1+α$$.
and c1, c2 take the form
$\begin{equation}
cos^{-1}\left(\frac{L^{2}+(\alpha^{2}-1)}{2L\alpha}\right),cos^{-1}\left(-\frac{L^{2}-(\alpha^{2}-1)}{2L}\right)
\end{equation}$
subsituting these directly into ODEs makes $${r\to0}$$, so we round to $(@bind s_round PlutoUI.Scrubbable(2)) significant digits
"""
# ╔═╡ 3b3b40f0-8cce-4d9d-b6f9-6c356a5393f2
md"""
# Initial Conditions
"""
# ╔═╡ 91729971-234e-4474-b96f-bb13144ffb73
begin
# initial conditions today
# energy = 14137524360429138; potential_energy = 0.14137524360429138
u0 = [-0.7165035605430603, 0.0, -2.149632692337036, 0.0]
# u0 = [-0.4441191852092743, 0.0, -2.4644575119018555, 0.0]
#e=0.5
u0 = [0.7335795847502655, -0.23913525101066202, 2.4133920290598883, -0.26152594568755844]
u0 = [-0.22337041795253754, 0.0, -2.3069841861724854, 0.0]
u0 = [-1.0157028436660767, 0.0, -2.825850486755371, 0.0]
u0 = [-0.635210394859314, 0.0, -2.9998667240142822, 0.0]
u0 = [0.7335795847502655, -0.23913525101066202, 2.4133920290598883, -0.26152594568755844]
u0 = [0.7898083,0.0,-0.6332571,0.0]
# energy = 0.5
u0 = [-1.1416815519332886, 0.0, -2.152212381362915, 0.0]
u0 = [-1.2426495552062988, 0.0, -2.6061832904815674, 0.0] #no
u0 = [-0.29311415553092957, 0.0, -2.1652581691741943, 0.0]
u0 = [-0.23045925796031952, 0.0, -2.291236639022827, 0.0] #higher_order
u0 = [-2.618450164794922, 1.2649110640673518, -3.0156140327453613, -1.0]
u0 = [0.8739768862724304, 0.0, -1.4866081476211548, 0.0]
u0=[0.08, 0.8, 0.0, -1.1]
u0 = [1.7558274269104004, 0.3162277660168379, 1.0075753927230835, 0.0] #quasi
u0 = [3.91, -0.0, 6.24, -0.0] #quasi_square
u0 = [4.14, 0.5, 1.48, -1.0] #quasi
u0 = [-2.9211299419403076, 0.3162277660168379, 2.126250743865967, 0.0] #quasi
u0 = [1.7558274269104004, 0.3162277660168379, 1.0075753927230835, 0.0] #quasi
u0 = [1.2972164154052734, 0.44721359549995787, -0.21258898079395294, 0.0] #quasi
u0 = [2.9368772506713867, 0.0, -0.7438594698905945, 0.0] #quasi
u0 = [-2.9211299419403076, 0.3162277660168379, 2.126250743865967, 0.0] #quasi
u0 = [0.804381251335144, 1.0, -0.16534698009490967, 0.0]#new_quasi
u0 = [2.9368772506713867, 0.0, -0.7438594698905945, 0.0] #quasi
u0 = [1.7558274269104004, 0.3162277660168379, 1.0075753927230835, 0.0] #quasi
u0 = [2.9368772506713867, 0.6324555320336759, -0.7438594698905945, 0.0] #quasi
u0 = [-0.18109431862831116, 0.6999999999999998, -0.02734754979610443, -1.1] #quasi in phase
u0 = [0.6220195889472961, 0.44721359549995787, 2.244325637817383, 0.0] #quasi
u0=[90.08, -0.01, 90.0, -0.91]
# u0=[c[1]+0.001,0.,c[2]-0.0002,0.]
u0=[-0.08, 0.8, 0.0, -1.1]
u0 = [-1.0426926535897931, 0.0, 2.045328181301066, 0.0]
u0 = [-1.0562073464102069, 0.0, 2.0580341015849193, 0.0] #opt
# u0=[0.01,0,pi+0.0003,0.]
u0=[-1.38, 0.0, 4.57, 0.0]
# u0 = [0.01,0.001,-0.002,0.003] #close to S1
# u0 = [pi,0.000001,-0.00000,-0.00001] #close to S2
u0 = [0.01,0.001, pi ,0.003] #close to S3
# u0 = [c1+0.01,0.001, c2-0.002 ,-0.003] #close to S4
#u0 = [3.1415927410125732, 0.0, -0.054074764251708984, 0.0]
# u0 = [4.062811374664307, 1.0, 3.044971227645874, 0.0]
end
# ╔═╡ 09ade1f2-0756-41c6-ac07-5606b09abc87
# ╠═╡ disabled = true
#=╠═╡
begin
#u0 = [θ1_init+0.000001,0.,θ2_init,0.]
# u0 = [c1+ϵ,2.05,c2+ϵ,0.0]
# u0 = [θ1_init,ω1,θ2_init,ω2]
u0=[0.08, 0.8, 0.0, -1.1]#quasi-sin_out,cos_in
# u0=[3.14,0.38,0.0,-0.0001]#quasi-out
# u0=[0.0,1.0,3.0,1.01]#both_out
# u0=[0.623,0.0,-2.474,0.0]#quasi-both_out
# u0 = [0.0,-3.2,3.14,3.0]
#u0 = [c1-ϵ,0.0,c2+ϵ,0.]
# energy 0.5
# energy 0.5 potential 0.5
u0=[0.8897242546081543, 0.0, 0.06493247300386429, 0.0]
u0 = [0.9709015488624573, 0.0, 2.5116994380950928, 0.0]
u0 = [-0.7458727955818176, 0.0, 1.0314502716064453, 0.0]
u0 = [-2.7164146900177, 0.0, 2.396667003631592, 0.0]
u0 = [-2.1372053623199463, 0.0, -1.614101529121399, 0.0] #sin_in_cos_out
u0 = [-2.9053826332092285, 0.0, 2.357250928878784, 0.0] #new
u0 = [-1.708585500717163, 0.0, -1.2012916803359985, 0.0] #sin_in_cos_out_half
u0 = [-2.40146803855896, 0.0, -1.834198236465454, 0.0]
u0 = [1.6021400690078735, 0.0, -2.24399471282959, 0.0] #sin_out_cos_out
u0 = [-0.7479982376098633, 0.0, -2.6808793544769287, 0.0] #chaos
u0 = [0.8739768862724304, 0.0, -1.4866081476211548, 0.0] #sin_out_cos_ou
# energy 0.5 potential 0.45
u0 = [0.9294304847717285, 0.3162277660168379, 2.3069841861724854, 0.0]
u0 = [-2.9211299419403076, 0.3162277660168379, 2.126250743865967, 0.0] #quasi
u0 = [-1.5608875751495361, 0.3162277660168379, 1.9605427980422974, 0.0]
u0 = [-0.9718883633613586, 0.3162277660168379, 1.0314502716064453, 0.0]
u0 = [-1.2372543811798096, 0.3162277660168379, 1.6613434553146362, 0.0]
u0 = [-1.5141088962554932, 0.3162277660168379, 1.9290480613708496, 0.0] #sin_in_cos_in
u0 = [-1.5196175575256348, 0.3162277660168379, 1.9329121112823486, 0.0]
u0 = [-0.991400420665741, 0.3162277660168379, 1.1259342432022095, 0.0]
u0 = [-0.9598749876022339, 0.3162277660168379, 0.9527135491371155, 0.0]
u0 = [-1.041481375694275, 0.3162277660168379, 0.13385231792926788, 0.0]
u0 = [-1.3936388492584229, 0.3162277660168379, -0.5537538528442383, 0.0]#quasi
u0 = [-2.574688673019409, 0.3162277660168379, -1.7388499975204468, 0.0]
u0 = [-3.0313613414764404, 0.3162277660168379, -2.004100799560547, 0.0]
u0 = [-0.9193374514579773, 0.3162277660168379, -2.291236639022827, 0.0] #choas
u0 = [-0.4845522344112396, 0.3162277660168379, -2.417215347290039, 0.0]
u0 = [-0.5905249118804932, 0.3162277660168379, -2.662233829498291, 0.0]
u0 = [-0.8897242546081543, 0.3162277660168379, -2.570263147354126, 0.0]
u0 = [-0.5400391221046448, 0.3162277660168379, -2.291236639022827, 0.0]
u0 = [0.5117882490158081, 0.3162277660168379, 2.3411126136779785, 0.0]
u0 = [2.8815248012542725, 0.3162277660168379, 1.9290480613708496, 0.0]
u0 = [2.1810054779052734, 0.3162277660168379, 1.4277540445327759, 0.0]
u0 = [1.7558274269104004, 0.3162277660168379, 1.0075753927230835, 0.0] #quasi
u0 = [1.5999176502227783, 0.3162277660168379, 0.8267349004745483, 0.0] #quasi
u0 = [1.0630860328674316, 0.3162277660168379, -0.07086299359798431, 0.0] #quasi
u0 = [0.9560887813568115, 0.3162277660168379, -0.9212188720703125, 0.0] #quasi
u0 = [1.126729130744934, 0.3162277660168379, -1.488122820854187, 0.0] #sin_out_cos_in
u0 = [1.0688594579696655, 0.3162277660168379, -0.05511565878987312, 0.0]
u0 = [1.0220575332641602, 0.3162277660168379, -0.19684164226055145, 0.0]
u0 = [0.9579108357429504, 0.3162277660168379, -0.9369662404060364, 0.0]
u0 = [1.0608805418014526, 0.3162277660168379, -1.346396803855896, 0.0]
u0 = [1.3162446022033691, 0.3162277660168379, -1.7558274269104004, 0.0] #sin_out_cos_half
u0 = [1.330649495124817, 0.3162277660168379, -1.7711701393127441, 0.0]
u0 = [1.3958609104156494, 0.3162277660168379, -1.8345640897750854, 0.0]
u0 = [2.102268695831299, 0.3162277660168379, -2.1682963371276855, 0.0]
u0 = [3.0156140327453613, 0.3162277660168379, -2.0975303649902344, 0.0]
# energy 0.5 potential 0.4
u0 = [-2.7660505771636963, 0.44721359549995787, -1.5826067924499512, 0.0] #new
u0 = [-1.8030694723129272, 0.44721359549995787, -0.722311794757843, 0.0]
u0 = [-1.2355148792266846, 0.44721359549995787, 0.559030294418335, 0.0]
u0 = [-1.7400801181793213, 0.44721359549995787, 1.7085459232330322, 0.0] #sin_out_cos_out
u0 = [-1.6455961465835571, 0.44721359549995787, 1.638767123222351, 0.0] #sin_out_cos_out
u0 = [-1.4179710149765015, 0.44721359549995787, 1.3778915405273438, 0.0]
u0 = [-2.1967527866363525, 0.44721359549995787, 1.8871625661849976, 0.0]
u0 = [-2.9629552364349365, 0.44721359549995787, 1.8345640897750854, 0.0]
u0 = [2.40146803855896, 0.44721359549995787, 1.3241767883300781, 0.0]
u0 = [2.1588151454925537, 0.44721359549995787, 1.1101869344711304, 0.0]
u0 = [1.7984302043914795, 0.44721359549995787, 0.7165035605430603, 0.0]
u0 = [1.3663203716278076, 0.44721359549995787, -0.007873665541410446, 0.0]
u0 = [1.409386157989502, 0.44721359549995787, -1.363908052444458, 0.0]
u0 = [1.42061185836792, 0.44721359549995787, 0.11810498684644699, 0.0]#quasi
u0 = [1.2972164154052734, 0.44721359549995787, -0.21258898079395294, 0.0] #quasi
u0 = [1.267660140991211, 0.44721359549995787, -1.0089129209518433, 0.0]
u0 = [1.2311201095581055, 0.44721359549995787, -0.7007562518119812, 0.0]
u0 = [1.4490694999694824, 0.44721359549995787, -1.425133466720581, 0.0]
u0 = [1.8975534439086914, 0.44721359549995787, -1.7949676513671875, 0.0]
u0 = [1.8975534439086914, 0.44721359549995787, -1.7949676513671875, 0.0]
u0 = [-0.527535617351532, 0.44721359549995787, -2.3732659816741943, 0.0] #choas
u0 = [-0.9320282340049744, 0.44721359549995787, -2.417215347290039, 0.0] #new_quasi
u0 = [-0.5736157894134521, 0.44721359549995787, -2.291236639022827, 0.0]
u0 = [-0.8973445892333984, 0.44721359549995787, -2.3069841861724854, 0.0] #choas
u0 = [-0.9311250448226929, 0.44721359549995787, -2.432962656021118, 0.0] #choas
u0 = [0.7322509288787842, 0.44721359549995787, 2.2096667289733887, 0.0] #choas
u0 = [0.5138645172119141, 0.44721359549995787, 2.417215347290039, 0.0] #new_quasi
u0 = [0.6220195889472961, 0.44721359549995787, 2.244325637817383, 0.0] #quasi
u0 = [0.6225046515464783, 0.44721359549995787, 2.24399471282959, 0.0] #choas
u0 = [0.5287531018257141, 0.44721359549995787, 2.3699734210968018, 0.0]
u0 = [0.6220195889472961, 0.44721359549995787, 2.244325637817383, 0.0] #quasi
u0 = [0.7952402234077454, 0.44721359549995787, 2.2232911586761475, 0.0] #choas
# energy 0.5 potential 0.3
u0 = [2.9368772506713867, 0.6324555320336759, 0.5802181363105774, 0.0]
u0 = [2.546109914779663, 0.6324555320336759, 0.18109431862831116, 0.0]
u0 = [2.4500322341918945, 0.6324555320336759, -0.007873665541410446, 0.0] #quasi
u0 = [2.3919014930725098, 0.6324555320336759, -0.3228203058242798, 0.0]
u0 = [2.9368772506713867, 0.6324555320336759, -0.7438594698905945, 0.0] #quasi
u0 = [0.8636139035224915, 0.6324555320336759, 2.3542261123657227, 0.0] #chaos
u0 = [-0.5653098821640015, 0.6324555320336759, -2.4644575119018555, 0.0]
u0 = [-0.8424822688102722, 0.6324555320336759, -2.5029332637786865, 0.0]
u0 = [-0.8441675305366516, 0.6324555320336759, -2.3227314949035645, 0.0]
u0 = [-0.869792640209198, 0.6324555320336759, -2.3699734210968018, 0.0]
u0 = [0.573279619216919, 0.6324555320336759, 2.40146803855896, 0.0]
u0 = [0.6489981412887573, 0.6324555320336759, 2.291236639022827, 0.0] #chaos
u0 = [-0.573279619216919, 0.6324555320336759, -2.40146803855896, 0.0]
# energy 0.5 potential 0.2
u0 = [0.7637456059455872, 0.7745966692414834, 2.324126958847046, 0.0] #chaos
u0 = [-0.796032190322876, 0.7745966692414834, -2.4802048206329346, 0.0] #chaos
u0 = [0.7794929146766663, 0.7745966692414834, 2.3322532176971436, 0.0] #chaos
u0 = [0.7007562518119812, 0.7745966692414834, 2.5216987133026123, 0.0] #chaos
# energy 0.5 potential 0.1
u0 = [-0.7165035605430603, 0.8944271909999159, -2.4695656299591064, 0.0]#chaos
u0 = [-0.6781327128410339, 0.8944271909999159, -2.40146803855896, 0.0]
u0 = [0.7007562518119812, 0.8944271909999159, 2.3763816356658936, 0.0] #chaos
# energy 1.0
# energy 1.0 potential 0.9
u0 = [1.1259342432022095, 0.44721359549995787, 3.077479839324951, 0.0]#chaos
u0 = [1.4560117721557617, 0.44721359549995787, 2.1652581691741943, 0.0]#chaos
u0 = [0.3207928538322449, 0.44721359549995787, 1.629848837852478, 0.0] #chaos
u0 = [-0.21386708319187164, 0.44721359549995787, 2.527446746826172, 0.0] #chaos
u0 = [1.2253527641296387, 0.44721359549995787, 2.9211299419403076, 0.0]
u0 = [1.3712280988693237, 0.44721359549995787, 2.1810054779052734, 0.0]
u0 = [-0.63776695728302, 0.44721359549995787, -1.529130458831787, 0.0]
# energy 1.0 potential 0.7
u0 = [0.7952402234077454, 0.7745966692414834, 2.01735782623291, 0.0] #chaos
u0 = [0.849571168422699, 0.7745966692414834, 2.0392794609069824, 0.0]
u0 = [2.55894136428833, 0.7745966692414834, 2.9486238956451416, 0.0]
u0 = [1.3936388492584229, 0.7745966692414834, 1.494260311126709, 0.0]
u0 = [-0.8096565008163452, 0.7745966692414834, 2.40146803855896, 0.0]
u0 = [-1.027974009513855, 0.7745966692414834, 2.574688673019409, 0.0]
u0 = [1.5353648662567139, 0.7745966692414834, -2.832716226577759, 0.0]
u0 = [0.048901937901973724, 0.7745966692414834, -0.8897242546081543, 0.0]
u0 = [-1.5038701295852661, 0.7745966692414834, -1.6036442518234253, 0.0]
u0 = [-2.562680244445801, 0.7745966692414834, -2.9053826332092285, 0.0]
u0 = [-0.9842082262039185, 0.7745966692414834, -2.706399440765381, 0.0]
# energy 1.0 potential 0.5
u0 = [1.8345640897750854, 1.0, -2.3329269886016846, 0.0]#chaos
u0 = [2.5939524173736572, 1.0, 1.9762901067733765, 0.0]
u0 = [1.6613434553146362, 1.0, 1.1513782739639282, 0.0]
u0 = [1.2891063690185547, 1.0, 0.7165035605430603, 0.0]
u0 = [0.804381251335144, 1.0, -0.16534698009490967, 0.0]#new_quasi
u0 = [0.7906999588012695, 1.0, -0.21258898079395294, 0.0]#new_quasi
u0 = [0.7479982376098633, 1.0, -0.4060882329940796, 0.0]
u0 = [0.7371619939804077, 1.0, -0.4802936017513275, 0.0]
u0 = [0.7266793251037598, 1.0, -0.8582295775413513, 0.0] #chaos
u0 = [1.3813059329986572, 1.0, -2.118016004562378, 0.0] #chaos
u0 = [2.9368772506713867, 1.0, -2.348701000213623, 0.0]
u0 = [0.9859780073165894, 1.0, 2.3699734210968018, 0.0] #chaos
u0 = [0.5066197514533997, 1.0, 2.291236639022827, 0.0] #chaos
u0 = [-2.6219305992126465, 1.0, 2.4089698791503906, 0.0]
u0 = [-2.2125000953674316, 1.0, 2.4088656902313232, 0.0]
u0 = [-0.7240377068519592, 1.0, 0.8109875917434692, 0.0]
u0 = [-0.8413484692573547, 1.0, 0.05511565878987312, 0.0] #new_quasi
u0 = [-0.8193594813346863, 1.0, 0.11810498684644699, 0.0]
# energy 0.3 potential 0.3
u0 = [2.968372106552124, 0.0, 0.6003814339637756, 0.0] #quasi
u0 = [2.8840129375457764, 0.0, 0.5432829260826111, 0.0]
u0 = [2.5363800525665283, 0.0, 0.16534698009490967, 0.0] #sin_in
u0 = [2.432962656021118, 0.0, -0.5056806206703186, 0.0] #sin_out
u0 = [2.391427993774414, 0.0, -0.3070729672908783, 0.0]
u0 = [2.669172763824463, 0.0, -0.7218575477600098, 0.0]
u0 = [2.9368772506713867, 0.0, -0.7438594698905945, 0.0] #quasi
u0 = [0.8761712908744812, 0.0, 2.40146803855896, 0.0] #border_quasi
u0 = [0.8582295775413513, 0.0, 2.343933582305908, 0.0] #quasi_sin_in_cos_out
u0 = [0.8636139035224915, 0.0, 2.3542261123657227, 0.0]
u0 = [0.8424822688102722, 0.0, 2.320694923400879, 0.0] #chaos
u0 = [0.7814453840255737, 0.0, 2.275489330291748, 0.0] #quasi_sin_in_cos_out
u0 = [0.7637456059455872, 0.0, 2.2694759368896484, 0.0] #quasi_sin_in_cos_out
u0 = [0.7007562518119812, 0.0, 2.268037796020508, 0.0] #chaos
u0 = [0.6489981412887573, 0.0, 2.291236639022827, 0.0] #chaos
u0 = [0.573279619216919, 0.0, 2.40146803855896, 0.0]
u0 = [0.5701944828033447, 0.0, 2.4959521293640137, 0.0]
u0 = [-0.8109875917434692, 0.0, -2.532987356185913, 0.0] #chaos_high
u0 = [-0.8551920652389526, 0.0, -2.3384788036346436, 0.0] #border_quasi
u0 = [-0.8424822688102722, 0.0, -2.320694923400879, 0.0] #chaos
u0 = [-0.8750137686729431, 0.0, -2.432962656021118, 0.0]
u0 = [-0.5668641328811646, 0.0, -2.4802048206329346, 0.0] #sin_in_cos_out
u0 = [-0.5693279504776001, 0.0, -2.417215347290039, 0.0]
u0 = [-2.7794039249420166, 0.0, -0.4576723575592041, 0.0] #sin_in_cos_out
u0 = [-0.7906999588012695, 0.0, 0.21258898079395294, 0.0]
u0 = [2.653425455093384, 0.0, -0.716033399105072, 0.0] #sin_out_cos_out
u0 = [2.7058229446411133, 0.0, 0.38580963015556335, 0.0] #sin_in_cos_out
# energy 0.3 potential 0.25
u0 = [-0.6535142660140991, 0.24494897427831774, -2.5417065620422363, -0.1]
u0 = [-0.5929373502731323, 0.24494897427831774, -2.4487099647521973, -0.1]#quasi sin_in_cos_out
u0 = [-0.5973527431488037, 0.24494897427831774, -2.4802048206329346, -0.1]
u0 = [-0.63776695728302, 0.24494897427831774, -2.5341451168060303, -0.1]
u0 = [-2.7660505771636963, 0.44721359549995787, -1.5826067924499512, 0.0]
u0 = [-0.9320282340049744, 0.44721359549995787, -2.417215347290039, 0.0] #new_quasi
u0=[0.08, 0.8, 0.0, -1.1]#quasi-sin_out,cos_in
u0 = [-2.9211299419403076, 0.3162277660168379, 2.126250743865967, 0.0] #quasi
u0 = [-1.3936388492584229, 0.3162277660168379, -0.5537538528442383, 0.0]#quasi
u0 = [-2.1372053623199463, 0.0, -1.614101529121399, 0.0]
u0 = [-1.5141088962554932, 0.3162277660168379, 1.9290480613708496, 0.0] #sin_in_cos_in
u0=[0.08, 0.8, 0.0, -1.1]
u0 = [0.16534698009490967, 1.3638181696985856, -0.010272976011037827, 0.0]
u0 = [0.21183769404888153, 1.3638181696985856, 0.0866103246808052, 0.0]
u0 = [1.0327496528625488, 1.3638181696985856, 2.24399471282959, 0.0]
u0 = [2.5690770149230957, 1.3638181696985856, 2.527446746826172, 0.0]
u0 = [-1.1755119562149048, 1.1832159566199232, -1.5196175575256348, 0.0]
u0 = [0.38580963015556335, 1.0723805294763609, 1.0666985511779785, 0.5]
#energy 1.5 potential 0.65
u0=[0.08, 0.8, 0.0, -1.1]
u0 = [0.15508675575256348, 0.6999999999999998, -0.03936832770705223, -1.1] #
u0 = [-0.16032752394676208, 0.6999999999999998, 0.02362099662423134, -1.1]
u0 = [2.6025874614715576, 0.6999999999999998, 2.55894136428833, -1.1]
u0 = [1.2923407554626465, 0.6999999999999998, 1.267660140991211, -1.1]
u0 = [2.385720729827881, 0.6999999999999998, -2.949526071548462, -1.1]
u0 = [0.7165035605430603, 0.6999999999999998, 2.0442256927490234, -1.1]
u0 = [0.7165035605430603, 0.6999999999999998, 2.0442256927490234, -1.1]
u0 = [1.4566282033920288, 0.6999999999999998, -2.6765408515930176, -1.1]
u0 = [-0.18109431862831116, 0.6999999999999998, -0.02734754979610443, -1.1] #quasi in phase
u0 = [2.4563663005828857, 1.0908712114635712, 0.007873665541410446, -1.1]
u0 = [-2.842393398284912, 1.0908712114635712, 0.7502715587615967, -1.1]
u0 = [0.7322509288787842, 1.0908712114635712, 2.2649950981140137, -1.1]
u0 = [0.16032752394676208, 0.6999999999999998, -0.02362099662423134, -1.1]
# potential 0.66
u0 = [0.002578687621280551, 0.6855654600401042, 0.1023576557636261, -1.1] #clear in phase
u0 = [-0.07086299359798431, 0.6855654600401042, -0.09567859768867493, -1.1]
u0 = [2.55894136428833, 0.6855654600401042, 2.562396764755249, -1.1]
u0 = [-0.30200788378715515, 0.6855654600401042, -2.6061832904815674, -1.1]
u0 = [-0.9419753551483154, 0.6855654600401042, -2.133763313293457, -1.1]
u0 = [-0.6220195889472961, 0.6855654600401042, -2.049452304840088, -1.1]
u0 = [-1.1056268215179443, 0.6855654600401042, -2.432962656021118, -1.1]
u0 = [0.5905249118804932, 0.6855654600401042, 2.0624496936798096, -1.1]
u0 = [-2.4959521293640137, 0.6855654600401042, 3.011141300201416, -1.1]
u0 = [2.3542261123657227, 0.6855654600401042, -2.9784982204437256, -1.1]
u0 = [0.002578687621280551, 0.6855654600401042, 0.1023576557636261, -1.1]
u0 = [3.0592548847198486, 0.6855654600401042, 3.094350576400757, -1.1]
u0 = [0.23578356206417084, 0.43588989435406705, 1.1259342432022095, -1.1]
#today
u0 = [2.23, -0.0, 5.48, -0.0] #both_out_phase
u0 = [3.91, -0.0, 6.24, -0.0] #quasi_square
u0 = [1.57, -0.0, 0.93, -0.0] #sin_in_cos_out
u0 = [0.58, -0.0, 2.5100000000000002, -0.0] #clear_sin_in_cos_out
u0 = [2.32, -0.0, 0.81, -0.0] #sharp_sin_in_cos_out
u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0 = [1.3, -0.0, 4.07, -0.0] #sin_out_cos_out
u0 = [4.14, 0.5, 1.48, -1.0] #quasi
u0 = [1.0874611139297485, 0.43588989435406705, 2.118016004562378, -1.1]
u0=[0.08, 0.8, 0.0, -1.1]
u0 = [2.32, -0.0, 0.81, -0.0] #sharp_sin_in_cos_out
u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0 = [0.58, -0.0, 2.5100000000000002, -0.0] #clear_sin_in_cos_out
u0=[0.08, 0.8, 0.0, -1.1] #sin_out_cos_in
u0=[rand(0:0.01:2*pi), .2, rand(0:0.01:2*pi), -0.5]
u0 = [2.32, -0.0, 0.81, -0.0] #sharp_sin_in_cos_out
# u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0=[1.157429,0.0,1.1882374,0.0]
u0=[1.0662,0.0,1.8030,0.0]
u0= [1.0786921977996826, 0.0, 1.8133879899978638, 0.0]
u0 = [-0.05729418992996216, 0.0, 2.3542261123657227, 0.0]
u01 = [-0.3070729672908783, 0.0, 2.731780529022217, 0.0]
u0=[0.,0.,pi,0.]
u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0 = [2.32, -0.0, 0.81, -0.0] #sharp_sin_in_cos_out
u0 = [0.723,0.0,2.42,0.0]
u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0 = [0.58, -0.0, 2.5100000000000002, -0.0] #clear_sin_in_cos_out
u0 = [1.3800000000000001, -0.0, 4.57, -0.0] #sharp_sin_out_cos_out
u0 = [2.32, -0.0, 0.81, -0.0] #sharp_sin_in_cos_out
u0=[0.08, 0.8, 0.0, -1.1]
u0=[0.722,0.,2.418,0.0]
end
╠═╡ =#
# ╔═╡ 65f37372-1ed7-4fbd-82d2-222f9bab9404
@show u0
# ╔═╡ 3d5a36ad-0ac1-4f09-9a04-78c0bb0a2595
md"""
# Parameters
### Potential order, n = ($(@bind n Scrubbable(1:1:5,default=3)))
### Exponential factor, K = $(@bind K PlutoUI.NumberField(1.5:0.1:20.0,default=2.0))
### Length between 2 rotors L = $(@bind L PlutoUI.NumberField(0.5:0.001:10., default=1.5; ))
### Additive parameter ξ = $(@bind ξ PlutoUI.NumberField(0.0:0.01:10., default=0.0; ))
#### Velocity parameters
Parameter A = $(@bind A1 PlutoUI.Slider(-1.0:0.1:20., default=1.0; show_value=true))
Parameter B = $(@bind B1 PlutoUI.Slider(0.0:0.1:20., default=0.0; show_value=true))
#### Velocities
$\dot\theta_1$ = $(@bind ω1 PlutoUI.NumberField(-10:0.05:10., default=2.05; ))
$\dot\theta_2$ = $(@bind ω2 PlutoUI.NumberField(-10.0:0.05:10., default=2.00; ))
α - rotor arm parameter = $(@bind α PlutoUI.Slider(0.2:0.01:1.0,default=1.0,show_value=true))
"""
# ╔═╡ 8d73eb5b-19b3-4b19-a594-5406bea9fc7d
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)
pr_θ1 = (α*(sin(θ1-θ2)+L*sin(θ1)))/r
pr_θ2 = (α*sin(θ2-θ1) -L*sin(θ2)) /r
du1 = (p1)/(α^2)
du2 = -pv0_r*pr_θ1
du3 = (p2)
du4 = -pv0_r*pr_θ2
return SVector{4}(du1, du2, du3, du4)
end
# ╔═╡ 5e3ee79c-36f2-433e-a1ab-6a33e82d07b5
c=round.([acos((L^2+α^2-1)/(2*L*α) ),acos((-L^2+α^2-1)/(2*L))],sigdigits=3)
# ╔═╡ 2e7c3f58-8005-4da6-9212-ef83416c3ee2
c1 = (acos( (L^2+α^2-1)/(2*L*α) ); )
# ╔═╡ 452f4196-dbae-49f3-a654-3340ea7cf19f
c2 = (acos((-L^2+α^2-1)/(2*L)); )
# ╔═╡ 938eb9e8-4de7-4ba1-848f-327e3b7591e3
md"""
θ1-init = $(@bind θ1_init Select([0, pi, -pi, c1, -c1]))
θ2-init = $(@bind θ2_init Select([0, pi, -pi, c2, -c2]))
"""
# ╔═╡ 7219ef29-1209-45cb-bb12-3dda344e0be7
k_energy(u0)=(u0[2]^2/(2*α^2))+(u0[4]^2/2)
# ╔═╡ e38fd9b7-f7cf-4212-88b2-5b6bdeed3a2d
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)+ξ
v(r) = (1-exp(-K*r^n))/r
return B1*v(t)
end
# ╔═╡ e1383092-4980-46cc-adaa-b6619306a102
p = ini(c[1],c[2])
# ╔═╡ d4e4ccd2-8cdd-4680-beb1-266c6c4d4583
params = [K,L,α,ξ,A1,B1]
# ╔═╡ a47e594d-ee55-41d1-9a5e-56677245b77b
html"<br><br>"
# ╔═╡ d81463e4-2089-4206-93c9-b66c2f3fc3ac
md"""
# ODE Integrator
"""
# ╔═╡ 75751d63-ba5f-4c20-a0e1-26684c96e05c
md"""
## Parameters for time evolution
"""
# ╔═╡ 64d3f024-8367-4a5e-8c5c-e74524f0d2de
md"""
Time of Evolution $(@bind time_end PlutoUI.Slider(10e2:10e2:10e4,default=10e2,show_value=true))
Sampling time $(@bind sampling_time PlutoUI.NumberField(0.005:0.1:10., default=0.005; ))
Transient time $(@bind transient_time PlutoUI.Slider(0:10e2:10e4,default=7000,show_value=true))
### Intermediate visual steps = $(@bind inter PlutoUI.Slider(50:50:1000,default=100;show_value=true))
time $(times = (0, time_end))
"""
# ╔═╡ 1d655439-320e-4a2e-9118-cc1145b88da2
md"""
Using Dynamical Systems
"""
# ╔═╡ 74cd42d7-031c-4727-8b05-57dee1647197
rotor_system_1 = CoupledODEs(rotor_2, u0, params; diffeq, t0 = 0.0)
# ╔═╡ 104c631c-d879-4d3b-aec6-e53963aa85e0
md"""
## Trajectory
"""
# ╔═╡ 9d24be83-f916-4cdc-82ea-4e9fd20972ea
@showtime Points_1, t = trajectory(rotor_system_1, time_end; Ttr = transient_time, Δt =sampling_time)
# ╔═╡ ae744956-f0f7-4a51-9ca2-0f6fb3940188
end_ind = lastindex(t)
# ╔═╡ 9577bcdb-e307-4dab-8571-2f7e5d61426a
md"""
## Energy Conservation
"""
# ╔═╡ fbef81d7-edf7-4b5e-8e27-09fe219cdd19
function hamil_cal(p)
# θ1,pθ1,θ2,pθ2 = p
# k, l = params
# r = (sqrt(abs(2 + l^2 - 2*l*(cos(θ1)- cos(θ2)) - 2*cos(θ2-θ1))))
# H = ((pθ1^2+pθ2^2)/2) + (exp(-k*r)/r)
# return H
θ1,pθ1,θ2,pθ2 = p
K, L, α = 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
v0=(1-exp(-K*r^n))/r
H = (((pθ1^2/α^2)+pθ2^2)/2) + v0
return H
end
# ╔═╡ 02cf1d96-5d5c-4bfd-bbec-2eb4f55f8ea6
md"""
### energy = $(hamil_cal(u0)),
### potential energy = $(hamil_cal(u0)-k_energy(u0))
"""
# ╔═╡ be066c96-8dcf-40b9-9b8d-72ea7d8371f7
function energy_check(Points)
energy = Vector{Float64}(undef,0)
energy_abs = Vector{Float64}(undef,0)
init_energy = hamil_cal(u0)
@showprogress for u in Points
push!(energy,hamil_cal(u))
end
energy_abs = energy.-init_energy
return energy, energy_abs
end
# ╔═╡ 875f041a-5750-4190-97d5-71e0114c511c
energy, energy_abs = energy_check(Points_1[1:inter:end_ind])
# ╔═╡ 04e385c3-0f3e-4827-b749-66f80aae8036
max_dev=round(maximum(abs.(energy_abs)); sigdigits=4,)
# ╔═╡ 5e47ff28-838f-407c-8ec3-4d4bbbe3e542
md"""
## Energy Deviation Plot
"""
# ╔═╡ bc4bbeec-3bc1-42dc-a639-5844e3f0070a
md"""
Maximum deviation of initial energy after evolution is
##### $(max_dev)
"""
# ╔═╡ f856999a-8bd7-43d2-8f2a-343ba8d93b8f
md"""
eaxispower = $(@bind eaxispower PlutoUI.Slider(0:1:10,default=8;show_value=true))
"""
# ╔═╡ 3a7d415d-a3fe-40fc-a600-c8a3ecc2c2b2
eyaxis1=(yformatter=y->string(round(Int, y * 10^eaxispower)))
# ╔═╡ 7b72f706-fc60-47cc-acea-1ca72c1adadf
minimum(energy)
# ╔═╡ 5fbd2359-f9b3-440e-a81b-d25f8e0eed43
function kenergy_check(Points)
kenergy = Vector{Float64}(undef,0)
@showprogress for u in Points
push!(kenergy,k_energy(u))
end
return kenergy
end
# ╔═╡ e8c746c6-15d0-49b1-9b9b-d475a82e9a35
k_energy_list = kenergy_check(Points_1[1:inter:end_ind])
# ╔═╡ a7ebc383-5702-4816-a1de-7081e185eb18
minimum(k_energy_list)
# ╔═╡ ab481719-c984-4d72-a950-1ece9cbab84e
html"<br>"
# ╔═╡ 5d0d9d19-18ac-464f-891e-898c8dd96e35
md"""
# Time Series Analysis
"""
# ╔═╡ 33a1733f-470f-49ef-b4c5-84c14dd7b2a8
md"""
Setting up time series Y-axis plot parameters to avoid cluttering.
tyaxispower = $(@bind tyaxispower PlutoUI.Slider(0:1:10,default=3;show_value=true))
tinitial = $(@bind tinitial PlutoUI.Slider(0:100:1000,default=0;show_value=true))
tsampling = $(@bind tsampling PlutoUI.Slider(0:20:1000,default=200;show_value=true))
tmargin = $(@bind tmargin PlutoUI.Slider(0:5:50,default=10;show_value=true))
tfont = $(@bind tfont PlutoUI.Slider(10:1:30,default=20;show_value=true))
"""
# ╔═╡ 27e6c24f-0d63-4e9c-a7df-878d7bebfce8
begin
Plots.plot(t[1:inter:end_ind], (energy), xlabel="Time", ylabel=L"\Delta E",legend=false, lw=3 )
# Plots.plot!(yticks=[0.4424753,0.4424784],yrotation=90)
# Plots.plot!(guidefontsize=14, tickfontsize=10)
Plots.plot!(guidefontsize=tfont+2, tickfontsize=tfont, yrotation=0, framestyle=:box, fontfamily="Computer Modern",
yaxis=eyaxis1,
# yformatter = :identity,
margin = tmargin*Plots.mm)
# Plots.plot!(yticks=(-2e-6,-0.00005))
# Plots.plot!(yticks=(0.00005,0.00005))
# Plots.plot!(ylims=(floor(max_dev;sigdigits=1),ceil(max_dev;sigdigits=1)))
Plots.annotate!([((-0.15, 1.05) , Plots.text( L" \times 10^{-%$eaxispower}", tfont-2, "Computer Modern", rotation=0), :black )])
end
# ╔═╡ 8a88ef8c-ab6a-49c8-bcad-283d1cca8b6e
tyaxis1=(formatter=y->string(round(Int, y * 10^tyaxispower)))
# ╔═╡ 33ec4de0-47bc-44e6-a505-f6c0681791ff
begin
t_inii=1500
t_finii=1800
end
# ╔═╡ 24e8ee35-1c45-4fa8-b91f-009b6097026a
md"""
### Multiplot $\theta_1$, $\theta_2$
"""
# ╔═╡ e4d6009a-d520-4cb5-ba8f-7d940f85ecc7
md"""
### $\theta_2$
"""
# ╔═╡ 263c25d2-5081-40c9-b740-a8689ea0e6e0
begin
theta2_timeseries = Points_1[:,3][1:inter:end_ind]
Plots.plot(t[1:inter:end_ind], sin.(theta2_timeseries), xlabel="Time", ylabel=L"sin\theta_2",label=false,lw=4)
# Plots.plot!(title="$(u0)")
# Plots.plot!(title="Rotor 2 time series")
# for small y value turn 90 degrees
# Plots.plot!(xticks=[7550, 7650,7750],yticks=[-2,1]*1e-8,yrotation=90)
# Plots.plot!(xlims=(7500,7600)) #S1
# Plots.plot!(yticks=[-1.00,0.00,1.00]) #S1
Plots.plot!(xlims=(7750,7950)) #S2
# Plots.plot!(yticks=[-3,3]*10e-6,) #S2
# Plots.plot!(yticks=[-.8,-0.5,-.1]) #S3
# Plots.plot!(xlims=(7500,7550)) #S3
# Plots.plot!(yticks=[0.655,0.661,0.668]) #S4
# Plots.plot!(xlims=(7500,7550)) #S4
# Plots.plot!(ylims=(-1,1))
# auto x lims
Plots.plot!(xlims=(tinitial+transient_time,tinitial+transient_time+tsampling))
Plots.plot!(guidefontsize=tfont+2, tickfontsize=tfont, yrotation=90, framestyle=:box, fontfamily="Computer Modern", yaxis=tyaxis1, formatter = :plain, margin = tmargin*Plots.mm, right_margin=(tmargin+tfont)*Plots.mm)
Plots.annotate!([((-0.25, 0.95) , Plots.text( L" \times 10^{-%$tyaxispower}", tfont-2, "Computer Modern", rotation=90), :black )])
# Plots.annotate!([(0.25,0,), ("S3",20,"Computer Modern"),:green])
end
# ╔═╡ c3d62644-72a1-4303-9573-117ac662bab8
begin
Plots.plot(t[1:inter:end_ind], sin.(Points_1[:,1][1:inter:end_ind]), xlabel="Time", ylabel=L"\sin\theta_1",label=:false,lw=4, lc=:red3, ytickfontcolor=:red3, yguidefontcolor=:red3, margin = tmargin*Plots.mm)
# Plots.plot!(title="$(u0)")
# Plots.plot!(title="Rotor 1 time series")
# Plots.plot!(xlims=(7500,7600)) #S1
# Plots.plot!(yticks=[-1.00,0.00,1.00]) #S1
Plots.plot!(twinx(),t[1:inter:end_ind], sin.(theta2_timeseries), ylabel=L"\sin\theta_2",label=false,lw=2, color=:blue3, yrotation=90, ytickfontcolor=:blue3, yguidefontcolor=:blue3, ls=:dash, )
# yticks=(0.3:0.1:0.7)
#Plots.plot!(xlims=(7750,7950)) #S2
# Plots.plot!(yticks=[-3,3]*10e-6,) #S2
# Plots.plot!(yticks=[-.8,-0.5,-.1]) #S3
# Plots.plot!(xlims=(7500,7550)) #S3
# Plots.plot!(yticks=[0.655,0.661,0.668]) #S4
# Plots.plot!(xlims=(7500,7550)) #S4
Plots.plot!(xlims=(tinitial+transient_time,tinitial+transient_time+tsampling)) #S2
Plots.plot!(guidefontsize=tfont+2, tickfontsize=tfont, yrotation=90, framestyle=:box, fontfamily="Computer Modern", yaxis=tyaxis1, formatter = :plain, )
Plots.annotate!([((-0.25, 0.95) , Plots.text( L" \times 10^{-%$tyaxispower}", tfont-2, "Computer Modern", rotation=90), :black )])
end
# ╔═╡ 9e4f5897-6b04-4596-a11a-c3f937d31056
# begin
# poin=poincaresos(rotor_system_1, (1, -3.75))
# Plots.scatter(poin[:,3],poin[:,4],marker=2)
# end
# ╔═╡ 7b620c85-cd57-409c-9919-1816257ce9ff
md"""
# Synchronization plots
To compare the experimental results, the raw data needs to be placed on
* `/tmp/co_rotation_time_series.xlsx`
* `/tmp/counter_rotation_time_series.xlsx`
"""
# ╔═╡ cd99b66f-8429-43ec-97f7-9b526e467313
md"""
## From simulated results
"""
# ╔═╡ 0cf4d4ff-8d08-408c-bc2f-792dac984ea7
begin
Plots.plot(t[1:inter:end_ind], abs.(sin.(Points_1[:,1][1:inter:end_ind]))-abs.(sin.(Points_1[:,3][1:inter:end_ind])), xlabel="Time", ylabel=L"|sin\theta_1|-|sin\theta_2|",label=false, color=:black, lw=2)
# Plots.plot!(xlims=(1060,1170))
# Plots.plot!(title="$(round.(u0,digits=5))")
Plots.plot!(guidefontsize=14, tickfontsize=12)
Plots.plot!(xguidefontsize=18,yguidefontsize=18, tickfontsize=16, fontfamily="Computer Modern",framestyle=:box)
# Plots.plot!(yrotation=90, yticks=[-6,-3,3,6]*1e-2, ylabel=false)
Plots.plot!(guidefontsize=tfont+2, tickfontsize=tfont, yrotation=90, framestyle=:box, fontfamily="Computer Modern", yaxis=tyaxis1, formatter = :plain, margin = tmargin*Plots.mm)
Plots.plot!(xlims=(tinitial+transient_time,tinitial+transient_time+tsampling))
Plots.annotate!([((-0.25, 1.00) , Plots.text( L" \times 10^{-%$tyaxispower}", tfont-2, "Computer Modern", rotation=90), :black )])
# Plots.annotate!([(0.25,0,), ("S3",20,"Computer Modern"),:green])
end
# ╔═╡ ff105624-9958-4c7b-a644-c182d5667611
md"""
## From Experimental Results
"""
# ╔═╡ 5d52d36b-0d36-46d1-afa2-b050ea71381a
begin
xf_time = XLSX.readxlsx("/tmp/co_rotation_time_series.xlsx")
# xf_time = XLSX.readxlsx("/tmp/counter_rotation_time_series.xlsx")
crp_time = xf_time["Sheet1"]
Plots.plot(float.(crp_time[:,1][2:end]), float.(crp_time[:,2][2:end]),
xlabel="t", ylabel=".",label=false, color=:black, lw=2)
Plots.plot!(xguidefontsize=18,yguidefontsize=18, tickfontsize=16, fontfamily="Computer Modern",framestyle=:box)
# Plots.plot!(yrotation=90, yticks=[-0.75,-.35,0.15,0.55], ylabel=false)
Plots.plot!(yrotation=90, yticks=[-.2,-.1,.1,.2], ylabel=false)
end
# ╔═╡ fe60650f-8ef3-4cf3-97f1-98146c39bc7f
xf_time[]
# ╔═╡ ada00ad6-2983-4acc-a1b6-0a723dde3acf
Points_1
# ╔═╡ e96744a0-035d-44f7-aaef-51eeebe3dcea
md"""
# Plots of relative phases
"""
# ╔═╡ a4e47f9f-3093-4651-ac6d-681261319a24
md"""
Plot of $\sin\theta_1$ vs $\sin\theta_2$ helps to visualuize between quasi periodic and choatic motion.
"""
# ╔═╡ 840daa90-71c1-4e82-8e8f-f4346258d12e
md"""
## From simulated data
"""
# ╔═╡ df343efa-dd6f-426e-8811-0f4177db8f5b
begin
seconds = 0:0.1:2
measurements = [8.2, 8.4, 6.3, 9.5, 9.1, 10.5, 8.6, 8.2, 10.5, 8.5, 7.2,
8.8, 9.7, 10.8, 12.5, 11.6, 12.1, 12.1, 15.1, 14.7, 13.1]
f = Makie.Figure(fonts = (; regular = "Times New Roman", bold = "Times New Roman Bold"), fontsize =tfont+5, resolution=(400,400));
# Initializing axis
ax1 = Makie.Axis(f[1, 1], xticklabelcolor= :blue3, yticklabelcolor = :blue3, xlabel=L"\sin\theta_1", xlabelcolor= :blue3,
ylabel=L"\sin\theta_2", ylabelcolor= :blue3, xticklabelfont = :regular,
yticklabelrotation=pi/2
)
ax2 = Makie.Axis(f[1, 1], xticklabelcolor = :red3, yticklabelcolor= :red3, xaxisposition = :top, yaxisposition = :right,
xlabel=L"\cos\theta_1", xlabelcolor= :red3,
ylabel=L"\cos\theta_2", ylabelcolor= :red3, yticklabelrotation=pi/2)
# plot for each axis
# can be changed to scatter plots
Makie.lines!(ax1, sin.(Points_1[:,1][1:inter:end_ind]),sin.(Points_1[:,3][1:inter:end_ind]), color = :blue3, alpha=0.6,
markersize=3.4)
Makie.lines!(ax2, cos.(Points_1[:,1][1:inter:end_ind]), cos.(Points_1[:,3][1:inter:end_ind]), color = :red3, alpha=0.7, markersize=2.4)
f
# Plots.plot(sin.(Points_1[:,1][1:inter:end_ind]), sin.(Points_1[:,3][1:inter:end_ind]), xlabel="sin(θ1)", ylabel="sin(θ2)", title="$(u0)")
# Plots.plot(sin.(Points_1[:,1][1:inter:end_ind]), sin.(Points_1[:,3][1:inter:end_ind]), xlabel=L"\sin\theta_1",ylabel=L"\sin\theta_2",
# color=:red, ls=:dash, alpha=0.5, label=:false)
# Plots.plot!(twinx(),cos.(Points_1[:,1][1:inter:end_ind]), cos.(Points_1[:,3][1:inter:end_ind]), ylabel=L"sin\theta_2",lw=2, color=:blue3, yrotation=90, ytickfontcolor=:blue3, yguidefontcolor=:blue3, ls=:dash, label=:false)
# # Plots.plot!(cos.(Points_1[:,1][1:inter:end_ind]), cos.(Points_1[:,3][1:inter:end_ind]), xlabel="sinθ₁, cosθ₁", ylabel="sinθ₂, cosθ₂",color=:black, alpha=0.3,label=false)
# # Plots.plot!(title="$(round.(u0,digits=5))")
# Plots.plot!(xguidefontsize=18,yguidefontsize=18, tickfontsize=16, fontfamily="Computer Modern",framestyle=:box)
# Plots.plot!(yrotation=90)
end
# ╔═╡ a43b12a2-c447-43f5-9311-30bf4ce6f2ee
begin
# maxs_cartp= findlocalextrema(sin.(Points_1[:,3][1:inter:end_ind]))
# maxsp = [hcat(Points_1[:,1][1:inter:end_ind][idx],Points_1[:,3][1:inter:end_ind][idx]) for idx in maxs_cartp]
end
# ╔═╡ c754247e-e624-4e9a-a303-f29be9335d10
begin
# Plots.scatter((maxsp[2][:,1]),(maxsp[2][:,2]))
# Plots.scatter!(sin.(maxsp[2][:,1]),sin.(maxsp[2][:,2]))
# ys=sin.(maxsp[1][:,2])
# xs=sin.(maxsp[1][:,1])
# zs=[r_eval([x,0,y,0]) for y in ys, x in xs]
# Plots.surface(xs,ys,zs)
end
# ╔═╡ e0cd2ac8-5924-47f4-a338-7c636dde3347
# ╠═╡ disabled = true
#=╠═╡
begin
Plots.plot(cos.(Points_1[:,1][1:inter:end_ind]), cos.(Points_1[:,3][1:inter:end_ind]), xlabel="cos(θ1)", ylabel="cos(θ2)", title="$(u0)")
# Plots.plot!(xlims=(-.69,-.681),ylims=(0.68,0.695),marker=1)
end
╠═╡ =#
# ╔═╡ 86825032-827e-4f02-ae10-f1c9d677bff2
md"""
## From experimental data
"""
# ╔═╡ a97a867f-db02-4b97-aa05-f9027f6b8ace
begin
xf = XLSX.readxlsx("/tmp/co_rotation_phase.xlsx")
xf = XLSX.readxlsx("/tmp/counter_rotation_phase.xlsx")
crp = xf["Sheet1"]
Plots.plot(float.(crp[:,1][2:end]), float.(crp[:,2][2:end]), xlabel="sinθ₁, sinθ₂", ylabel="cosθ₁, cosθ₂",color=:red, lw=:2, alpha=0.9, label=false)
Plots.plot!(float.(crp[:,3][2:end]), float.(crp[:,4][2:end]), xlabel="sinθ₁, cosθ₁", ylabel="sinθ₂, cosθ₂",color=:black, lw=:2, alpha=0.8,label=false)
Plots.plot!(xguidefontsize=18,yguidefontsize=18, tickfontsize=16, fontfamily="Computer Modern",framestyle=:box)
Plots.plot!(yrotation=90,ylabel=false)
end
# ╔═╡ 34913940-aedb-4d7c-9518-fee884b7d82d
float.(crp[:,1][2:end]),float.(crp[:,2][2,end])
# ╔═╡ acd99f44-34b1-4f57-8f81-302cc93e4ba7
md"""
# Misc
These function were written to analyze the problem but later left either due to reduncy of the similar function in another notebook or the time constrain to develop the concept further.
"""
# ╔═╡ 200aa8aa-9685-45e2-a8a2-4beadcf91be6
md"""
## Largest Lyapunov Exponent
"""
# ╔═╡ ba366079-e509-4c34-8b89-de9235030fe4
a=lyapunov(rotor_system_1,10000; Ttr=10000., Δt=0.0005)
# ╔═╡ 307429a6-db87-4629-a435-cfc203240f45
sum(a)
# ╔═╡ e2e44a89-43e0-4ce7-8887-ee25e434c00f
md"""
## Time Series & Frequency for r
"""
# ╔═╡ 6b0aca36-718b-4613-96b5-3bdd8c07feee
function r_eval( 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 r
end
# ╔═╡ 9be7453c-ad69-4188-b211-4b394a17614a
r_eval(u0)
# ╔═╡ fbc7769b-3ad4-4748-a051-a4dc03e5f5d0
function r_series(Points)
r_times = Vector{Float64}(undef,0)
for u in Points
push!(r_times,r_eval(u))
end
return r_times
end
# ╔═╡ 99da5e12-76b1-4276-93e0-8ca0600c5f2a
r_times = r_series(Points_1[1:inter:end_ind]);
# ╔═╡ 7cdee447-5522-4034-baf5-2571bcd9d8c5
Plots.plot(t[1:inter:end_ind], r_times, xlabel="t", ylabel="r")
# Plots.plot(t[1:inter:end_ind], r_times, xlabel="t", ylabel="r")
# ╔═╡ 21b71ea2-0e0c-4ab0-b5ed-8a23390fd8a6
begin
signal3=sin.(r_times)
fft_result3=fft(signal3)
frequencies3= fftfreq(length(signal3), 1 / sampling_time)
Plots.plot(frequencies3, abs.(fft_result3), xlabel="Normalized Frequency of r", ylabel="Amplitude", label="FFT Result")
end
# ╔═╡ 7b5cbc96-f4d2-4980-a4cc-b33453f78cf8
md"""
## Phase Space Visualization
"""
# ╔═╡ 14e27e6a-c419-4901-b82a-09246dc56e72
function polar2cart(Points)
y1 = α*sin.(Points[:,1])
x1 = α*cos.(Points[:,1])
y2 = sin.(Points[:,3])
x2 = cos.(Points[:,3])
((x1, y1), (L .+ x2, y2))
end
# ╔═╡ f1b0ef7d-8093-44b3-8538-c3252ef8e91f
md"""
Interactive visual stepsize $(@bind vi_step NumberField(1:20))
"""
# ╔═╡ 70a51464-516b-4032-9697-011dca24b59a
md"""
Interaction time for Visualization $(@bind phase_end PlutoUI.NumberField(1:vi_step:1000000, default=1))
"""
# ╔═╡ 5c286320-91dd-4ffe-8ea8-bc6337cba8e9
rotor1_cord, rotor2_cord = polar2cart(Points_1[1:phase_end])
# ╔═╡ a1c129bd-8050-4182-96e0-624f589d439d
begin
rotor1_x, rotor1_y = rotor1_cord[1][end], rotor1_cord[2][end]
rotor2_x, rotor2_y = rotor2_cord[1][end], rotor2_cord[2][end]
tsad1 = "Calulated end points of both rotors"
end
# ╔═╡ 647d4444-95c8-4fee-a5da-5cb0efb1cc99
phase_x2y2 = Plots.plot(rotor2_cord,
legend = false,
title = "Rotor 2 Phase Space x2 y2",
xaxis = "x",
yaxis = "y",
formatter = :plain,
widen = true,
xlims = (L-1.2, L+1.2),
ylims = (-1.2, 1.2),
aspect_ratio = 1.1
)
# ╔═╡ 27f84427-fe5f-4e18-bf12-07f82ad48860
phase_x1y1 = Plots.plot(rotor1_cord,
legend = false,
title = "Rotor 1 Phase Space x1 y1",
xaxis = "x",
yaxis = "y",
formatter = :plain,
widen = true,
xlims = (-1.2, 1.2),
ylims = (-1.2, 1.2),
aspect_ratio = 1.1
)
# ╔═╡ f05f39f1-f4e2-4e4c-a89c-3a3b107098a4
md"""
system at $(round(phase_end*sampling_time; sigdigits=2)) seconds
"""
# ╔═╡ c4300ffd-3bd5-4ff3-8579-fa7e59739d9c
begin
Plots.plot(
([0,rotor1_x],[0,rotor1_y]),
color = [:black,:red],
xlim = (-5,5),
ylim = (-2,2)
)
Plots.scatter!(
([0,rotor1_x],[0,rotor1_y]),
color = [:black,:red],
xlim = (-5,5),
ylim = (-2,2)
)
Plots.plot!(
([L,rotor2_x],[0,rotor2_y]),
color = [:black,:blue],
xlim = (-5,5),
ylim = (-2,2)
)
Plots.scatter!(
([L,rotor2_x],[0,rotor2_y]),
color = [:black,:blue],
xlim = (-5,5),
ylim = (-2,2),
title ="Combined Rotor system"
)
end
# ╔═╡ 825814d9-3e40-464c-bd46-129611efe854
rotor1_cord
# ╔═╡ Cell order:
# ╟─17342ecc-7da0-43dc-ba88-38bb85e44eca
# ╟─6d08c01e-8c5c-11ee-06e1-47df2e37668b
# ╠═448c89b6-7ae7-42ce-b019-ef9019d3e72a
# ╟─bf26074c-53a9-4df7-aee4-394dbecf38ef
# ╟─1d4cf734-84a2-4adc-ae86-b02f1601f6f2
# ╟─4015e194-6e7a-4a5c-ae06-1406e6aadfed
# ╟─3d223e16-ceb4-4cca-95a4-784da1d2dcf3
# ╟─84f64684-e2af-4ead-adab-cea67ddbaada
# ╠═8d73eb5b-19b3-4b19-a594-5406bea9fc7d
# ╟─9ee65bfd-94ae-495d-aff0-0ee2dc321a45
# ╟─41336926-46a0-4c41-a87a-720162191d27
# ╟─a7d1c54b-6e6e-4331-b76c-f479130acc08
# ╠═616557c2-4839-4257-b331-188c04d17ce3
# ╟─eb3ccf92-8fd6-4974-ac64-e9237ce25901
# ╟─249672c9-d9f2-4900-80bf-70b53c370e63
# ╟─ed3c93f7-d35c-4c45-86f7-b3e1f7f30fed
# ╟─5e3ee79c-36f2-433e-a1ab-6a33e82d07b5
# ╟─2e7c3f58-8005-4da6-9212-ef83416c3ee2
# ╟─452f4196-dbae-49f3-a654-3340ea7cf19f
# ╟─938eb9e8-4de7-4ba1-848f-327e3b7591e3
# ╠═3b3b40f0-8cce-4d9d-b6f9-6c356a5393f2
# ╠═91729971-234e-4474-b96f-bb13144ffb73
# ╟─09ade1f2-0756-41c6-ac07-5606b09abc87
# ╠═65f37372-1ed7-4fbd-82d2-222f9bab9404
# ╠═7219ef29-1209-45cb-bb12-3dda344e0be7
# ╠═02cf1d96-5d5c-4bfd-bbec-2eb4f55f8ea6
# ╠═e38fd9b7-f7cf-4212-88b2-5b6bdeed3a2d
# ╠═e1383092-4980-46cc-adaa-b6619306a102
# ╟─3d5a36ad-0ac1-4f09-9a04-78c0bb0a2595
# ╠═d4e4ccd2-8cdd-4680-beb1-266c6c4d4583
# ╟─a47e594d-ee55-41d1-9a5e-56677245b77b
# ╟─d81463e4-2089-4206-93c9-b66c2f3fc3ac
# ╟─75751d63-ba5f-4c20-a0e1-26684c96e05c
# ╟─64d3f024-8367-4a5e-8c5c-e74524f0d2de
# ╟─ae744956-f0f7-4a51-9ca2-0f6fb3940188
# ╟─1d655439-320e-4a2e-9118-cc1145b88da2
# ╠═74cd42d7-031c-4727-8b05-57dee1647197
# ╟─104c631c-d879-4d3b-aec6-e53963aa85e0
# ╟─9d24be83-f916-4cdc-82ea-4e9fd20972ea
# ╟─9577bcdb-e307-4dab-8571-2f7e5d61426a
# ╟─fbef81d7-edf7-4b5e-8e27-09fe219cdd19
# ╟─be066c96-8dcf-40b9-9b8d-72ea7d8371f7
# ╠═875f041a-5750-4190-97d5-71e0114c511c
# ╠═04e385c3-0f3e-4827-b749-66f80aae8036
# ╟─5e47ff28-838f-407c-8ec3-4d4bbbe3e542
# ╟─bc4bbeec-3bc1-42dc-a639-5844e3f0070a
# ╟─f856999a-8bd7-43d2-8f2a-343ba8d93b8f
# ╟─3a7d415d-a3fe-40fc-a600-c8a3ecc2c2b2
# ╟─27e6c24f-0d63-4e9c-a7df-878d7bebfce8
# ╠═7b72f706-fc60-47cc-acea-1ca72c1adadf
# ╠═e8c746c6-15d0-49b1-9b9b-d475a82e9a35
# ╠═a7ebc383-5702-4816-a1de-7081e185eb18
# ╟─5fbd2359-f9b3-440e-a81b-d25f8e0eed43
# ╟─ab481719-c984-4d72-a950-1ece9cbab84e
# ╟─5d0d9d19-18ac-464f-891e-898c8dd96e35
# ╟─33a1733f-470f-49ef-b4c5-84c14dd7b2a8
# ╠═8a88ef8c-ab6a-49c8-bcad-283d1cca8b6e
# ╟─33ec4de0-47bc-44e6-a505-f6c0681791ff
# ╟─24e8ee35-1c45-4fa8-b91f-009b6097026a
# ╟─c3d62644-72a1-4303-9573-117ac662bab8
# ╟─e4d6009a-d520-4cb5-ba8f-7d940f85ecc7
# ╟─263c25d2-5081-40c9-b740-a8689ea0e6e0
# ╟─9e4f5897-6b04-4596-a11a-c3f937d31056
# ╟─7b620c85-cd57-409c-9919-1816257ce9ff
# ╟─cd99b66f-8429-43ec-97f7-9b526e467313
# ╟─0cf4d4ff-8d08-408c-bc2f-792dac984ea7
# ╟─ff105624-9958-4c7b-a644-c182d5667611
# ╠═5d52d36b-0d36-46d1-afa2-b050ea71381a
# ╠═fe60650f-8ef3-4cf3-97f1-98146c39bc7f
# ╠═ada00ad6-2983-4acc-a1b6-0a723dde3acf
# ╟─e96744a0-035d-44f7-aaef-51eeebe3dcea
# ╟─a4e47f9f-3093-4651-ac6d-681261319a24
# ╟─840daa90-71c1-4e82-8e8f-f4346258d12e
# ╠═df343efa-dd6f-426e-8811-0f4177db8f5b
# ╟─a43b12a2-c447-43f5-9311-30bf4ce6f2ee
# ╟─c754247e-e624-4e9a-a303-f29be9335d10
# ╠═e0cd2ac8-5924-47f4-a338-7c636dde3347
# ╟─86825032-827e-4f02-ae10-f1c9d677bff2
# ╠═34913940-aedb-4d7c-9518-fee884b7d82d
# ╠═a97a867f-db02-4b97-aa05-f9027f6b8ace
# ╟─acd99f44-34b1-4f57-8f81-302cc93e4ba7
# ╟─200aa8aa-9685-45e2-a8a2-4beadcf91be6
# ╠═ba366079-e509-4c34-8b89-de9235030fe4
# ╠═307429a6-db87-4629-a435-cfc203240f45
# ╠═9be7453c-ad69-4188-b211-4b394a17614a
# ╟─e2e44a89-43e0-4ce7-8887-ee25e434c00f
# ╟─6b0aca36-718b-4613-96b5-3bdd8c07feee
# ╟─fbc7769b-3ad4-4748-a051-a4dc03e5f5d0
# ╟─99da5e12-76b1-4276-93e0-8ca0600c5f2a
# ╠═7cdee447-5522-4034-baf5-2571bcd9d8c5
# ╠═21b71ea2-0e0c-4ab0-b5ed-8a23390fd8a6
# ╟─7b5cbc96-f4d2-4980-a4cc-b33453f78cf8
# ╠═14e27e6a-c419-4901-b82a-09246dc56e72
# ╠═5c286320-91dd-4ffe-8ea8-bc6337cba8e9
# ╠═a1c129bd-8050-4182-96e0-624f589d439d
# ╟─647d4444-95c8-4fee-a5da-5cb0efb1cc99
# ╟─27f84427-fe5f-4e18-bf12-07f82ad48860
# ╠═f1b0ef7d-8093-44b3-8538-c3252ef8e91f
# ╠═70a51464-516b-4032-9697-011dca24b59a
# ╠═f05f39f1-f4e2-4e4c-a89c-3a3b107098a4
# ╠═c4300ffd-3bd5-4ff3-8579-fa7e59739d9c
# ╟─825814d9-3e40-464c-bd46-129611efe854