Commit 8858663f authored by Javier's avatar Javier

comp_t0 cubic fits

parent 2bfbc0c6
......@@ -140,20 +140,21 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
xmax = size(Ysl, 2)
nt0 = t0_guess(t, Ysl, plat, L)
Y_aux = Matrix{uwreal}(undef, xmax, 3)
Y_aux = Matrix{uwreal}(undef, xmax, 5)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
for j = nt0-2:nt0+2
Y_aux[i, k] = uwreal(Ysl[:, i, j], id)
k = k + 1
end
end
x = t[nt0-1:nt0+1]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:3] .* x.^2 / L^3
x = t[nt0-2:nt0+2]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:5] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
t0 = x_lin_fit(par, 0.3)
@. model(x, p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
par = fit_routine(model, x, t2E, 4)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
uwerr(t0)
uwerr.(t2E)
......@@ -169,9 +170,8 @@ function comp_t0(Y::YData, plat::Vector{Int64}; L::Int64, pl::Bool=false, rw::Un
display(gcf())
figure()
plot(collect(1:xmax), value.(Y_aux[:, 2]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[2]+e[2], v[2]-e[2], alpha=0.75, color="green")
plot(collect(1:xmax), value.(Y_aux[:, 3]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[3]+e[3], v[3]-e[3], alpha=0.75, color="green")
ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$")
title(string(L"$t/a^2 = $", t[nt0]))
......@@ -188,7 +188,6 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
id = getfield.(Y, :id)
vtr = getfield.(Y, :vtr)
replica = length.(vtr)
if !all(id .== id[1])
println("IDs are not equal")
......@@ -196,25 +195,26 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
end
Ysl = isnothing(rw) ? Ysl : apply_rw.(Ysl, rw)
tmp = Ysl[1]
[tmp = cat(tmp, Ysl[k], dims=1) for k = 2:nr]
nt0 = t0_guess(t, tmp, plat, L)
xmax = size(tmp, 2)
Y_aux = Matrix{uwreal}(undef, xmax, 3)
Y_aux = Matrix{uwreal}(undef, xmax, 5)
for i = 1:xmax
k = 1
for j = nt0-1:nt0+1
for j = nt0-2:nt0+2
Y_aux[i, k] = uwreal(tmp[:, i, j], id[1], replica)
k = k + 1
end
end
x = t[nt0-1:nt0+1]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:3] .* x.^2 / L^3
par, chi2exp = lin_fit(x, t2E)
x = t[nt0-2:nt0+2]
t2E = [plat_av(Y_aux[:, j], plat) for j=1:5] .* x.^2 / L^3
t0 = x_lin_fit(par, 0.3)
@. model(x, p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
par = fit_routine(model, x, t2E, 4)
fmin(x, p) = model(x, p) .- 0.3
t0 = root_error(fmin, t[nt0], par)
if pl
uwerr(t0)
uwerr.(t2E)
......@@ -230,8 +230,8 @@ function comp_t0(Y::Vector{YData}, plat::Vector{Int64}; L::Int64, pl::Bool=false
display(gcf())
figure()
plot(collect(1:xmax), value.(Y_aux[:, 2]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[2]+e[2], v[2]-e[2], alpha=0.75, color="green")
plot(collect(1:xmax), value.(Y_aux[:, 3]) * t[nt0]^2 / L^3, "x")
fill_between(collect(plat[1]:plat[2]), v[3]+e[3], v[3]-e[3], alpha=0.75, color="green")
ylabel(L"$t^2E(x0, t)$")
xlabel(L"$x_0/a$")
title(string(L"$t/a^2 = $", t[nt0]))
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment