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

using Plots
gr( fmt=:png )

# 模型参数和初始条件
Ω = 2*pi/86400;    # 地球绕其轴旋转的角速度(rad/s)
λ = 59.93/180*pi;  # 纬度(rad)
g = 9.8195;        # 重力加速度(m/s^2)
L = 97;            # 钟摆长度(m)

initial_x = L/100; # X(m)坐标的初始值
initial_y = 0;     # Y坐标的初始值(m)
initial_xdot = 0;  # X轴速度的初始值(m/s)
initial_ydot = 0;  # Y轴速度的初始值(m/s)

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="福柯摆的负载的位置(俯视图)\收件人:ode23;RelTol=1e-5;执行时间=$(round(dt_ode23,digits=2))with",
      titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="X坐标(m)", ylabel="Y坐标(m)",
      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 )

vlen = length(x_ode23["x"].value)
half_vec = ceil(Int32, vlen/2);
plot!( x_ode23["x"].value[(half_vec-bar_w):half_vec], x_ode23["y"].value[(half_vec-bar_w):half_vec], 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="福柯摆的负载的位置(俯视图)\收件人:ode23t;RelTol=1e-5;执行时间=$(round(dt_ode23t,digits=2))with", 
          titlefont=9, guidefontsize=8, tickfontsize=8, leg=:none, xlabel="X坐标(m)", ylabel="Y坐标(m)",
          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 )

vlen = length(x_ode23t["x"].value)
half_vec = ceil(Int32, vlen/2);
plot!( x_ode23t["x"].value[(half_vec-bar_w):half_vec], x_ode23t["y"].value[(half_vec-bar_w):half_vec], lc=:yellow, lw=3 )

p2 = plot!( [initial_x], [0], mc=:black, st=:scatter, ms=8 )

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

using Serialization

# 要更新结果，请删除"结果"文件
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="麦克斯 相对功率输出",
            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="执行时间\nmodels(s)", 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)
    )