Skip to content

Commit

Permalink
more cumulant tweaks
Browse files Browse the repository at this point in the history
  • Loading branch information
joott committed Jul 19, 2022
1 parent 84685d9 commit d808fcc
Showing 1 changed file with 12 additions and 11 deletions.
23 changes: 12 additions & 11 deletions cumulant_analysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@ cd("/share/tmschaef/jkott/modelB/KZ/cumulants")

using Plots
using DelimitedFiles
using Statistics

const L = 32
const z = 3.906e0
Expand All @@ -19,28 +20,28 @@ const t_e = (m²e - m_b)/m_a
const t_c_steps = trunc(Int, t_c/Δt)

skip = 10
maxt = (trunc(Int, t_e / Δt)+1)÷skip+1
maxt = (trunc(Int, t_e / Δt)+1)÷skip+2

id_max = 16
series_max = 16
run_max = 32
n_files = id_max * series_max * run_max

M4 = zeros(Float32, maxt, n_files)
M2 = zeros(Float32, maxt, n_files)
M4 = zeros(Float64, div(maxt,10)+1, n_files)
M2 = zeros(Float64, div(maxt,10)+1, n_files)

function cumulant_run(fn, i)
df = readdlm(fn, ' ')
M = df[:,3]
M = df[1:10:end,3]

M4[:,i] .= M.^4
M2[:,i] .= M.^2
M4[:,i] .= (M*(L^3/2)^(1/4)).^4
M2[:,i] .= (M*(L^3/2)^(1/4)).^2
end

function collect_files()
Threads.@threads for id in 1:id_max
for series in 1:series_max, run in 1:run_max
cumulant_run("sum_L_$L"*"_id_$id"*"_series_$series"*"_run_$run.dat", run + (series-1)*series_max + (id-1)*series_max*id_max)
cumulant_run("sum_L_$L"*"_id_$id"*"_series_$series"*"_run_$run.dat", run + (series-1)*run_max + (id-1)*series_max*run_max)
end
end
end
Expand All @@ -50,19 +51,19 @@ collect_files()
M4_err = std(M4, dims=2) / sqrt(n_files)
M2_err = std(M2, dims=2) / sqrt(n_files)

M4 .= sum(M4, dims=2)/n_files
M2 .= sum(M2, dims=2)/n_files
M4 .= sum(M4, dims=2)
M2 .= sum(M2, dims=2)

c4 = M4 .- 3 * M2.^2
c4_err = sqrt.(M4_err.^2 .+ (6 * M2 .* M2_err).^2)

output_file = open("data_cumulant_$L.jl","w")

write(output_file, "c4 = ")
show(output_file, C_tot)
show(output_file, c4)
write(output_file, "\n\n")

write(output_file, "c4_err = ")
show(output_file, c4_err)

close(output_file)
close(output_file)

0 comments on commit d808fcc

Please sign in to comment.