Pkg.add(["Measures", "Serialization"])

using Plots
gr( fmt=:png )

# Параметры модели и начальные условия
Ω = 2*pi/86400;    # Угловая скорость вращения Земли вокруг оси (рад/с)
λ = 59.93/180*pi;  # Широта (рад)
g = 9.8195;        # Ускорение свободного падения (м/с^2)
L = 97;            # Длина маятника (м)

initial_x = L/100; # Начальное значение координаты x (м)
initial_y = 0;     # Начальное значение координаты y (м)
initial_xdot = 0;  # Начальное значение скорости по оси x (м/с)
initial_ydot = 0;  # Начальное значение скорости по оси y (м/с)

TotalEnergy = 1/2*(initial_xdot^2 + initial_ydot^2) + g*(L-sqrt( L^2-initial_x^2-initial_y^2));

modelName = "foucault_pendulum_model"
model = modelName in [m.name for m in engee.get_all_models()] ? engee.open( modelName ) : engee.load( "$(@__DIR__)/$modelName.engee" );

engee.set_param!( modelName, "SolverName"=>"BS3", "RelTol"=>1e-5 ) # Установим решатель, аналогичный ode23
dt_ode23 = @elapsed begin x_ode23 = engee.run( modelName ); end;

engee.set_param!( modelName, "SolverName"=>"Trapezoid", "RelTol"=>1e-5 ) # Установим решатель, аналогичный ode23t
dt_ode23t = @elapsed begin x_ode23t = engee.run( modelName ); end;

using Plots

plot( x_ode23["nE"].time, x_ode23["nE"].value, label="ode23", title="Сравнение качества решений при RelTol = 1e-5", lw=2 )
plot!( x_ode23t["nE"].time, x_ode23t["nE"].value, ylimits=(0,1.1), label="ode23t", leg=:bottomleft, lw=2 )
plot!( xlabel="Время (с)", ylabel="Относительная полная энергия",
       titlefont=11, guidefontsize=10, tickfontsize=9, legendfontsize=10 )

using Measures

bar_w = 150 # "Толщина" линий фиксирующих плоскость колебания маятника в разные моменты

lastPlane_y = maximum( x_ode23["y"].value[end-bar_w:end] )
lastPlane_x = maximum( x_ode23["x"].value[end-bar_w:end] )
plot( x_ode23["x"].value[1:2:end], x_ode23["y"].value[1:2:end],
      st=:scatter, ms=.8, markerstrokecolor=:blue, markercolor=:blue, 
      title="Положение груза маятника Фуко (вид сверху)\nРешатель: ode23; RelTol = 1e-5; Время выполнения=$(round(dt_ode23,digits=2)) с",
      titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="Координата x (м)", ylabel="Координата y (м)",
      left_margin=13mm, bottom_margin=13mm, top_margin=5mm, aspect_ratio = :equal )
scatter!( [initial_x], [0], mc=:black, ms=8 )                                                # Черный шар
plot!( x_ode23["x"].value[1:bar_w], x_ode23["y"].value[1:bar_w], lc=:black, lw=3 )           # Черная линия
plot!( x_ode23["x"].value[end-bar_w:end], x_ode23["y"].value[end-bar_w:end], lc=:red, lw=3 )
plot!( x_ode23["x"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], x_ode23["y"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], lc=:yellow, lw=3 )
annotate!( initial_x, 0.05, text("Начальное положение", :black, rotation = 90, 7, :left) )
plot!( [lastPlane_x-0.1, lastPlane_x-0.8], [lastPlane_y, lastPlane_y], lw=13, lc=:white ) # Белая линия
p1 = annotate!( lastPlane_x-0.1, lastPlane_y, text("Конечная плоскость", :red, 7, :right) )

plot( x_ode23t["x"].value[1:2:end], x_ode23t["y"].value[1:2:end], st=:scatter, ms=.8, markerstrokecolor=:blue, markercolor=:blue,
          title="Положение груза маятника Фуко (вид сверху)\nРешатель: ode23t; RelTol = 1e-5; Время выполнения=$(round(dt_ode23t,digits=2)) с", 
          titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="Координата x (м)", ylabel="Координата y (м)",
          left_margin=13mm, bottom_margin=13mm, top_margin=5mm, aspect_ratio = :equal )
plot!( x_ode23t["x"].value[1:bar_w], x_ode23t["y"].value[1:bar_w], lc=:black, lw=3 )           # Черная линия
plot!( x_ode23t["x"].value[end-bar_w:end], x_ode23t["y"].value[end-bar_w:end], lc=:red, lw=3 )
plot!( x_ode23t["x"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], x_ode23t["y"].value[Int32(ceil(end/2)-bar_w):Int32(ceil(end/2))], lc=:yellow, lw=3 )
p2 = plot!( [initial_x], [0], mc=:black, st=:scatter, ms=8 )

plot( p1, p2, size = (1000,450) )

using Serialization

# Чтобы обновить результаты, удалите файл "results"
if !isfile( "$(@__DIR__)/results.bin" )

    solvers = ["Tsit5" "DP5" "Rosenbrock23" "Trapezoid" "BS3"]
    reltols = 10 .^( collect(range(-7,-3,length=9)) )

    nE_results = zeros( length(solvers), length(reltols) )
    dt_results = zeros( length(solvers), length(reltols) )

    # Изучим все комбинации (решатель; относительная точность)
    for (i,solver) in enumerate(solvers)
        print( solver, ": ")
        for (j,reltol) in enumerate(reltols)
            print( round(reltol, sigdigits=3), ", " )
            # Установим для модели новый решатель и относительную точность
            engee.set_param!( modelName, "SolverName"=>solver, "RelTol"=>reltol )
            # Вычисляем модель и заодно измеряем время
            dt_results[i,j] = @elapsed begin sim_data = engee.run( modelName ); end;
            # Вычисляем соотношение 
            nE_results[i,j] = abs( sim_data["nE"].value[1] - sim_data["nE"].value[end] )
        end
        println()
    end
    
    serialize( "$(@__DIR__)/results.bin", (solvers, reltols, nE_results, dt_results) )
else
    (solvers, reltols, nE_results, dt_results) = deserialize( "$(@__DIR__)/results.bin" );
end

markerShapes_ = [ :circle :rect :star5 :diamond :hexagon ]
markerStrokeColors_ = [1 2 3 4 5]

plot(
      plot( reltols, nE_results', xscale=:log10, yscale=:log10,
            label=solvers, ylabel="Макс. относительная\nошибка полной энергии",
            titlefont=11, guidefontsize=8, tickfontsize=9, legendfontsize=7, left_margin=13mm, legend=:bottomleft,
            marker=markerShapes_, markerstrokecolor=markerStrokeColors_, markercolor=RGBA(1,1,1,0) ),
      plot( reltols, dt_results', xscale=:log10,
            label=solvers, ylabel="Время выполнения\nмодели (с)", xlabel="Относительная точность алгоритма численного интегрирования",
            titlefont=11, guidefontsize=8, tickfontsize=9, legendfontsize=7, left_margin=13mm, bottom_margin=13mm,
            marker=markerShapes_, markerstrokecolor=markerStrokeColors_, markercolor=RGBA(1,1,1,0) ),
      size=(800,500),
      layout=(2,1)
    )