### 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
# ╔═╡ fd93fe8f-09a1-4efe-bb18-93afacccacd3
begin
import Pkg
Pkg.activate("/home/hisacro/.julia/environments/v1.9/Project.toml")
using Markdown
using InteractiveUtils
using ImplicitPlots, HypertextLiteral
using ImplicitPlots.Contour
using PlutoUI,OrdinaryDiffEq
using ForwardDiff, LinearAlgebra
using DelimitedFiles, DynamicalSystems, OrdinaryDiffEq, ProgressLogging, FFTW
using Plots ; gr(size=(500,500))
using IntervalRootFinding, IntervalArithmetic, StaticArrays,IntervalConstraintProgramming
#using MDBM
using LaTeXStrings
end
# ╔═╡ ce05288d-e162-4e7e-9b46-a00fa0e646c0
TableOfContents()
# ╔═╡ e58a6df6-47a6-4d3f-bd1d-091524f88f9b
md"""
# Abstract
Introduction of an alternative system parameter, $r_{ex}$ capturing the dynamics of Libration motioin. $r_{ex}$ is the maxima of distance between rotors at a given energy level.
Energy regimes are suggested based on the features of the system such that,
* $E_1$ - Start of libration Type I.
* $E_2$ - Start of libration Type II.
* $E_3$ - Merged wells in Type I.
* $E_4$ - Merged wells in Type II.
* $E_5$ - Uncoupled motion.
# Goal
As documented in the bif\_area.jl, there is an inherent issue with the area of the contour as parameter. Intent of introducing a new system parameter, $r_{ex}$ is to simplify the study while capturing all the dynamics.
# Methodology
1. Fix an energy, E and set kinetic energy as 0. Solve for $V(r)=E$.
2. List all such $r_{ex}$ values. It is found for this system, there are atmost two values of $r_{ex}$ for each energy levels.
3. Associate values of $r_{ex}$ with the libration types I & II.
4. Loop the steps 1-3 for energies 0 to 1.5
5. Report the maximum value of $V$ and corresponding $r_{ex}$ as $r_{end}$
6. Plot $r_{ex}$ as a function of energy, distinguish the solution curves of Type I, Type II.
7. List the energy levels where the contours merge, these are $V(0),V(L+2),V(2-L),V(L),V(r_{end})$.
8. Mark the state configurations, S1 to S4 on the plot. These are functions of $(\theta_1,\theta_2)$.
9. Distinguish the energy regime with vertical ranges.
"""
# ╔═╡ 557a0674-05c0-4ddb-8470-38caceda731b
md"""
# System Initialization
"""
# ╔═╡ 24a878b3-d44e-4d0d-9bae-8b3390b2fa1f
md"""
## Setting up parameters
Potential order, n = ($(@bind n Scrubbable(1:1:5,default=3)))
Exponential factor, K = $(@bind K PlutoUI.NumberField(0.5:0.1:20.0,default=2.0))
Total Energy, $E$ = $(@bind energy PlutoUI.NumberField(0.0:0.1:10., default=0.5; ))
Potential Energy $P_{energy}$ = $(@bind p_energy PlutoUI.NumberField(0.0:0.01:10., default=0.5; ))
momenta $p_2$ = $(@bind p2 PlutoUI.NumberField(-10.0:0.1:10.,default=0.0; ))
Pivot lengths, L = $(@bind L PlutoUI.Slider(0.25:0.01:5., default=1.500; show_value=true ))
"""
# ╔═╡ 8b506079-b9b3-46b8-9de0-44369c1edbf4
md"""
# Plot of r_ex as a function of energy
"""
# ╔═╡ 0a065bda-583d-4875-8af1-d79fcdb4d5e5
md"""
Font size = $(@bind tfont PlutoUI.Slider(10:1:30,default=15;show_value=true))
Plot margin = $(@bind phase_margin PlutoUI.Slider(-50:5:50,default=0;show_value=true))
Regime Labels = $(@bind regime_y PlutoUI.Slider(-5:0.05:5,default=0.5;show_value=true))
"""
# ╔═╡ bfa94741-62ee-439f-be98-e1795dfbb950
md"""
### Pivot lengths, L = $(L)
"""
# ╔═╡ 025b18be-2c5a-46e0-a902-60dbb2b8860e
# ╠═╡ disabled = true
#=╠═╡
md"""
## (Retired) Outline
1. Finding the smallest of r values when energy is fixed.
2. For each of the r and energy find the corresponding (θ1, θ2)
3. Store argmax(sin.(θ1)) and the corresponding θ2 in a list
4. Plot above
"""
╠═╡ =#
# ╔═╡ 91c0f7e2-f6bc-4b27-bbc8-b52f245afc3f
# ╠═╡ disabled = true
#=╠═╡
md"""
## (Retired) Find the list of (θ1, θ2) for Type I
"""
╠═╡ =#
# ╔═╡ 3e4a9935-49b9-4490-b344-66e53e8e2e08
# ╠═╡ disabled = true
#=╠═╡
# type 1
rindex = findfirst(item -> item == p_energy, pr_list_1)
# type 2
# rindex = findfirst(item -> item == p_energy, pr_list_2)
╠═╡ =#
# ╔═╡ 9f58e0a0-ff90-40f9-bf3e-ffbca5d510d2
# ╠═╡ disabled = true
#=╠═╡
length(r_list_2)
╠═╡ =#
# ╔═╡ cd2465ee-1a31-4ed3-a3bf-a472c8a31535
# ╠═╡ disabled = true
#=╠═╡
# pr_list_2[5793]
╠═╡ =#
# ╔═╡ 903c1a55-8805-4c7d-af11-20c8697f8a3b
# ╠═╡ disabled = true
#=╠═╡
begin
nd1_list = Vector{Float64}(undef,0)
pnd1_list = Vector{Float64}(undef,0)
nd2_list = Vector{Float64}(undef,0)
pnd2_list = Vector{Float64}(undef,0)
# @progress for i in 4700:100:length(r_list_2)
# @progress for i in 0.447:0.01:0.497
@progress for i in 0.290:0.001:0.301
# rindex = findfirst(item -> item == i, pr_list_2)
# mandelbrot(x,y)= r_list_2[i]-sqrt((2+L^2-2*L*(cos(x)-cos(y))-2*cos(x-y)))
# mandelbrot(x,y)= r_list_1[rindex]-sqrt((2+L^2-2*L*(cos(x)-cos(y))-2*cos(x-y)))
# mandelbrot(x,y)= find_r(p_energy)[1] - sqrt((2+L^2-2*L*(cos(x)-cos(y))-2*cos(x-y)))
mandelbrot(x,y)= find_r(i)[2] - sqrt((2+L^2-2*L*(cos(x)-cos(y))-2*cos(x-y)))
# mandelbrot(x,y)= r_list_2[rindex] - sqrt((2+L^2-2*L*(cos(x)-cos(y))-2*cos(x-y)))
# type I
# Mandelbrotmdbm=MDBM_Problem(mandelbrot,[-pi:0.1:pi,0:0.1:2*pi])
# type II
Mandelbrotmdbm=MDBM_Problem(mandelbrot,[0:0.1:2*pi,-pi:0.1:1*pi])
MDBM.solve!(Mandelbrotmdbm,3)
theta_energy_x, theta_energy=getinterpolatedsolution(Mandelbrotmdbm)
nd1 = argmax(cot.(theta_energy_x))
nd2 = argmin(cot.(theta_energy_x))
# push!(pnd1_list, pr_list_1[i])
# push!(nd1_list,abs(theta_energy[nd1]-theta_energy[nd2]))
# push!(pnd2_list, pr_list_2[i])
push!(pnd2_list, i)
push!(nd2_list,abs(theta_energy[nd1]-theta_energy[nd2]))
# push!(nd1_list, argmax(cos.(theta_energy_x)))
# push!(nd2_list, argmin(cos.(theta_energy_x)))
# push!(nd3_list, argmax(cos.(theta_energy_x)))
# push!(nd4_list, argmin(cos.(theta_energy_x)))
end
end
╠═╡ =#
# ╔═╡ f52bd23a-6548-4dcb-871f-a69d91892c00
# ╠═╡ disabled = true
#=╠═╡
begin
nd2_list = [1.3770916187552225, 1.7758461799611789, 2.086784169126725, 2.3446304601268153, 2.573522975507724, 2.7738991419420675, 2.9546556281699115, 3.120043345333971, 3.2703312622650795, 3.4129925717356393, 3.544661659569203, 3.6677729585649246, 3.7851833498611303, 3.8930731828536373, 4.001938502180938, 4.105862524438064, 4.19951669494214, 4.293057454246856, 4.382567304844319, 4.471926235862469, 4.557121392023524, 4.6379049179277585, 4.720330327621814, 4.8019914847125795, 4.8805894789627935, 4.960530214520084, 5.038035227252008, 5.115399131170232, 5.19497259007017, 5.2742361607558434, 5.356843963591441]
pnd2_list = [0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6]
end
╠═╡ =#
# ╔═╡ 2da79af3-5867-453f-9f48-dca7485892dd
#=╠═╡
@show nd2_list
╠═╡ =#
# ╔═╡ 0439cac0-4f21-4e52-869e-d5d772cb3636
#=╠═╡
begin
for i1 in nd2_list
if i1 < 0.001
samindy = findfirst(item -> item == i1, nd2_list)
println(pnd2_list[samindy])
end
end
end
╠═╡ =#
# ╔═╡ bbb7b365-d744-4d90-b7cb-8c9f553397a6
#=╠═╡
@show pnd2_list
╠═╡ =#
# ╔═╡ 69f6ce2b-1e2f-4a52-bca7-f8d77184d9f1
#=╠═╡
begin
scatter(pnd2_list,nd2_list, label=false)
end
╠═╡ =#
# ╔═╡ 9e9b036e-422a-4ef7-ae90-058c6a5b0f1e
# ╠═╡ disabled = true
#=╠═╡
begin
scatter_data = Tuple{Float64, Float64}[]
# Push tuples into the array
push!(scatter_data, (1.0, 2.0))
push!(scatter_data, (1.0, 3.0))
# Check the result
println(scatter_data)
end
╠═╡ =#
# ╔═╡ 62113f97-6068-46c4-bfaf-9653b233b17d
# ╠═╡ disabled = true
#=╠═╡
nd1_list
╠═╡ =#
# ╔═╡ e54f0ca2-43e0-412e-a344-18dac0cf6316
# ╠═╡ disabled = true
#=╠═╡
begin
nd1 = argmax(cot.(theta_energy_x))
nd2 = argmin(cot.(theta_energy_x))
nd3 = argmax(cot.(theta_energy))
nd4 = argmin(cot.(theta_energy))
end
╠═╡ =#
# ╔═╡ e935699b-9e41-4214-a120-684f32d598e1
#=╠═╡
nd2
╠═╡ =#
# ╔═╡ 40181b4b-4ac1-494b-99b0-3efc71a9783b
#=╠═╡
begin
Plots.scatter(theta_energy_x,theta_energy, markersize=0.5, label=false)
Plots.scatter!([theta_energy_x[nd1]],[theta_energy[nd1] ],color=:white, label=true)
Plots.scatter!([theta_energy_x[nd2]],[theta_energy[nd2] ],color=:black,label=true)
Plots.scatter!([theta_energy_x[nd3]],[theta_energy[nd3] ],color=:cyan, label=true)
Plots.scatter!([theta_energy_x[nd4]],[theta_energy[nd4] ],color=:red, label=true, title="E = $(p_energy)")
Plots.plot!([theta_energy_x[nd1],theta_energy_x[nd2]], [theta_energy[nd1],theta_energy[nd2]], color=:blue, lw=3)
# type I
Plots.plot!(ylim=(0,2*pi),xlim=(-pi,pi))
# type II
# Plots.plot!(xlim=(0-0.1,2*pi),ylim=(-pi,pi))
end
╠═╡ =#
# ╔═╡ ce731a78-bf66-11ef-0fdd-8142bce085c0
md"""
# Appendix
"""
# ╔═╡ 223c8bde-eb57-47c3-9605-8e69ae07b6c4
md"""
## Finding r for given p_energy
"""
# ╔═╡ b9ef0687-55aa-4d74-a2a8-77a189540ede
md"""
accuracy = $(@bind acc PlutoUI.Slider(0.000001:0.00001:0.1, default=0.00001; show_value=true))
"""
# ╔═╡ fa0180f5-b21f-4f3c-82ee-d5e0a4029ddb
begin
function find_r(pr)
potential(r1) = -pr+((1-exp(-K*r1^n))/r1);
# newer syntax
# rts_r=IntervalRootFinding.roots(potential, 0..L+2.0);
# r_values = sort([mid(root_region(u)) for u ∈ rts_r if u.status ==
# :unique && mid(root_region(u))>0])
# 3 roots for odd power but 1 being < 0
# midpoints = mid.(interval.(rts_r));
rts_r=IntervalRootFinding.roots(potential, 0..L+2.0,Newton, 1e-10);
midpoints = mid.(interval.(rts_r));
r_values = sort([mid(interval(u)) for u ∈ rts_r if u.status == :unique && mid(interval(u))>0])
return r_values
end
end
# ╔═╡ 3e23043e-463d-4bc1-9b67-f759c30756c8
find_r(p_energy)
# ╔═╡ ae09d270-cb4e-4cef-be6d-346ce46a8dfc
find_r(0.7)
# ╔═╡ 27d101aa-29fd-4f2d-8fb4-3c114ea3c2a3
md"""
## Calculating r_end
"""
# ╔═╡ 84aa6ae4-2bb8-4c34-9823-8443fb464689
function rmax()
# for newer version of root finding.
# dvdr(r)=(exp(K*r^n)-K*n*(r^n)-1)
# dvdr(r)=(exp(-K*r^n)*(exp(K*r^n)-K*n*(r^n)-1))/r^2
# rts_vmax=roots(dvdr, 0.5..1.5)
# return [ mid.(root_region(rt)) for rt in rts_vmax ][1]
# Implicit differentiation of the potential, V(r)
dvdr(r)=(exp(-K*r^n)*(exp(K*r^n)-K*n*(r^n)-1))/r^2
# Uses IntervalRootingFinding.jl to search for roots
rts_vmax=roots(dvdr, 0.5..1.5, Newton, 1e-8);
# Returns the maxima of r
return mid.(interval.(rts_vmax))[1];
end
# ╔═╡ e197012c-6b6c-412e-a80f-9ec0bdf7a19d
vend = rmax();
# ╔═╡ 858303a4-d196-47c7-9642-58948f6e2732
r_end = vend;
# ╔═╡ 901f0475-7e9f-47b6-bdeb-80f14a9e73ce
V(r)=(1-exp(-K*r^n))/r
# ╔═╡ c25774a3-fe67-47aa-85a7-5dccb551fe0d
begin
# for bifn of libration types
pr_list_1 = Vector{Float64}(undef,0)
r_list_1 = Vector{Float64}(undef,0)
pr_list_2 = Vector{Float64}(undef,0)
r_list_2 = Vector{Float64}(undef,0)
p_range = 0.0001:0.001:1.0
for pr in p_range
# IntervalBox(-pi..pi, 1)
potential(r1) = -pr+((1-exp(-K*r1^n))/r1);
rts_r=roots(potential, 0..L+2.0, Newton, 1e-8)
# 3 roots for odd power but 1 being < 0
midpoints = mid.(interval.(rts_r));
r_values = sort([mid(interval(u)) for u ∈ rts_r if u.status ==
:unique && mid(interval(u))>0])
if length(r_values) > 0
push!(pr_list_1,pr)
push!(r_list_1,r_values[1])
end
if length(r_values) > 1
push!(pr_list_2,pr)
push!(r_list_2,r_values[2])
end
end
Plots.plot(grid = false)
Plots.plot!(pr_list_1,r_list_1,lw=3, color=:blue3)
Plots.plot!(pr_list_2,r_list_2,lw=3, color=:red2)
Plots.vline!([V(L),V(L+2),V(2-L), V(r_end)], lw=:2, ls=:dash, color=:gray)
## Regimes Annotate
Plots.annotate!([((0+V(L+2))/2,2+regime_y,(L"E_1",tfont,"Computer Modern",:black))])
Plots.annotate!([((V(2-L)+V(L+2))/2,2+regime_y,(L"E_2",tfont,"Computer Modern",:black))])
Plots.annotate!([((V(2-L)+V(L))/2,2+regime_y,(L"E_3",tfont,"Computer Modern",:black))])
Plots.annotate!([((V(r_end)+V(L))/2,2+regime_y,(L"E_4",tfont,"Computer Modern",:black))])
Plots.annotate!([((V(r_end)+1.0)/2,2+regime_y,(L"E_5",tfont,"Computer Modern",:black))])
Plots.annotate!([((V(2-L)+V(L))/2,1+regime_y,("II",tfont,"Computer Modern",:red3))])
Plots.annotate!([((V(2-L)+V(L))/2,-0.1+regime_y,("I",tfont,"Computer Modern",:blue3))])
Plots.scatter!([V(L+2)],[L+2], color="black", label="V(L1)", markersize=5 )
# Plots.annotate!([(V(L+2)+0.05,L+2-1.15,(L"E_2",20,"Computer Modern",:red2))])
# Plots.annotate!([(V(L+2)+0.1,L+2-0.15,("V(L2)",18,"Computer Modern",:black))])
Plots.annotate!([(V(L+2)+0.1,L+2-0.15,(L"V(S2)",tfont,"Computer Modern",:black))])
Plots.scatter!([V(abs(2-L))],[abs(2-L)], color="black", markersize=5 )
# Plots.annotate!([(V(2-L)+0.05,2-L+0.01,(L"E_3",20,"Computer Modern",:blue3))])
#Plots.annotate!([(V(2-L)+0.05,2-L+0.25,("V(L3)",18,"Computer Modern",:black))])
Plots.annotate!([(V(2-L)+0.05,2-L+0.25,(L"V(S3)",tfont,"Computer Modern",:black))])
Plots.scatter!([V(L)] ,[L], color="black", markersize=5 )
# Plots.annotate!([(V(L)+0.1,L+0.00,("E₄",20,"Computer Modern",:red2))])
# Plots.annotate!([(V(L)+0.05,L+0.15,(L"V(L4)",18,"Computer Modern",:black))])
Plots.annotate!([(V(L)+0.05,L+0.15,(L"V(S1)",tfont,"Computer Modern",:black))])
Plots.annotate!([(0.1,0.4,(L"V(S4)",tfont,"Computer Modern",:black))])
Plots.scatter!([V(vend)] ,[vend], color="black" , markersize=5)
#Plots.annotate!([(0.9,vend-0.5,(L"V(r_{end})",18,"Computer Modern",:black))])
Plots.scatter!([0.0] ,[0.0], color="black" , markersize=5)
# Plots.annotate!([(0.20,0.15,("E₁",20,"Computer Modern",:blue3))])
Plots.plot!(markersize=5, xlabel="Energy", ylabel=L"r_{ex}",framestyle=:box, fontfamily="Computer Modern", legend=false, xlims=[-0.01,1], xticks=(0.0:0.2:1.2))
Plots.plot!(
# xticks=([ V(2-L), V(vend)],[ L"V(2-L)", L"V(r_{end})"]),
xticks=([V(L+2), V(2-L), V(L), V(vend)],[ round(V(L+2),sigdigits=2), round(V(2-L),sigdigits=2), round(V(L),sigdigits=2), round(V(r_end),sigdigits=2)]),
tickfontsize=tfont, guidefontsize=tfont,
grid=false,
margin=phase_margin*Plots.mm
)
Plots.plot!(twiny(),[],[],
xticks=([V(L+2), V(2-L), V(L), V(r_end)],[L"V(L+2)\quad", L"\quad \quad V(2-L)", L"V(L)", L"V(r_{end})"]),
# xticks=([V(L+2), V(L)],[L"V(L+2)", L"V(L)" ]),
markersize=5, tickfontsize=tfont-2, guidefontsize=tfont-2,framestyle=:box, fontfamily="Computer Modern", legend=false, xlims=[-0.01,1],
)
Plots.vline!([1], lw=:3, color=:black, ylabelfontsize=tfont+5)
end
# ╔═╡ 8ad38c27-3627-4770-9383-ff456ee20992
md"""
-------------
#### r corresponding to maximum of potential, $r_{end}$ = $(round(r_end,digits=4)),
#### $V(r_{end})$ = $(round(V(r_end),digits=4)).
------------
S1-S4 are stationary configurations found in rotor\_transcend\_fixed.jl.
* S1 - (0,0), (pi,pi)
* S2 - (pi,0)
* S3 - (0,pi)
* S4 - (c1,c2), (-c1,c2)
------------
"""
# ╔═╡ 9fdded08-f0b8-4f53-9859-2d8148ee9524
V(0.001);
# ╔═╡ fdfa8522-5fd1-4e2a-bed6-ce92f39b9771
find_r(V(r_end-0.01))
# ╔═╡ Cell order:
# ╠═ce05288d-e162-4e7e-9b46-a00fa0e646c0
# ╟─e58a6df6-47a6-4d3f-bd1d-091524f88f9b
# ╟─557a0674-05c0-4ddb-8470-38caceda731b
# ╟─fd93fe8f-09a1-4efe-bb18-93afacccacd3
# ╟─24a878b3-d44e-4d0d-9bae-8b3390b2fa1f
# ╟─8b506079-b9b3-46b8-9de0-44369c1edbf4
# ╟─0a065bda-583d-4875-8af1-d79fcdb4d5e5
# ╟─bfa94741-62ee-439f-be98-e1795dfbb950
# ╟─c25774a3-fe67-47aa-85a7-5dccb551fe0d
# ╟─8ad38c27-3627-4770-9383-ff456ee20992
# ╟─9fdded08-f0b8-4f53-9859-2d8148ee9524
# ╟─e197012c-6b6c-412e-a80f-9ec0bdf7a19d
# ╟─858303a4-d196-47c7-9642-58948f6e2732
# ╠═025b18be-2c5a-46e0-a902-60dbb2b8860e
# ╠═91c0f7e2-f6bc-4b27-bbc8-b52f245afc3f
# ╠═3e4a9935-49b9-4490-b344-66e53e8e2e08
# ╠═9f58e0a0-ff90-40f9-bf3e-ffbca5d510d2
# ╠═cd2465ee-1a31-4ed3-a3bf-a472c8a31535
# ╠═903c1a55-8805-4c7d-af11-20c8697f8a3b
# ╠═f52bd23a-6548-4dcb-871f-a69d91892c00
# ╠═2da79af3-5867-453f-9f48-dca7485892dd
# ╠═0439cac0-4f21-4e52-869e-d5d772cb3636
# ╠═bbb7b365-d744-4d90-b7cb-8c9f553397a6
# ╠═69f6ce2b-1e2f-4a52-bca7-f8d77184d9f1
# ╠═9e9b036e-422a-4ef7-ae90-058c6a5b0f1e
# ╠═62113f97-6068-46c4-bfaf-9653b233b17d
# ╠═e54f0ca2-43e0-412e-a344-18dac0cf6316
# ╠═e935699b-9e41-4214-a120-684f32d598e1
# ╠═40181b4b-4ac1-494b-99b0-3efc71a9783b
# ╟─ce731a78-bf66-11ef-0fdd-8142bce085c0
# ╟─223c8bde-eb57-47c3-9605-8e69ae07b6c4
# ╟─b9ef0687-55aa-4d74-a2a8-77a189540ede
# ╠═fa0180f5-b21f-4f3c-82ee-d5e0a4029ddb
# ╠═3e23043e-463d-4bc1-9b67-f759c30756c8
# ╠═ae09d270-cb4e-4cef-be6d-346ce46a8dfc
# ╠═fdfa8522-5fd1-4e2a-bed6-ce92f39b9771
# ╟─27d101aa-29fd-4f2d-8fb4-3c114ea3c2a3
# ╠═84aa6ae4-2bb8-4c34-9823-8443fb464689
# ╟─901f0475-7e9f-47b6-bdeb-80f14a9e73ce