using DelimitedFiles
using Plots
using LsqFit

@. model(N, p) = N^(2 * p[1])

function main()
    data = readdlm("analysis.csv", ',')

    ns = view(data, :, 1)
    squared_r_ees = view(data, :, 2)

    expected_ν = 0.588
    expected(N) = N^(2 * expected_ν)

    p0 = [expected_ν]
    fit = curve_fit(model, ns, squared_r_ees, p0)
    simulation_ν = fit.param[1]

    p = plot(; xlabel="N", ylabel="<Rₑₑ²>", xscale=:log10, yscale=:log10, title="10⁵ steps", dpi=300)

    plot!(p, ns, squared_r_ees; label="Simulation ν = $(round(simulation_ν; digits=3))")
    plot!(p, ns, map(N -> expected(N), ns); label="Expected ν = $expected_ν")

    savefig(p, "analysis.png")

    return nothing
end

main()