read_ms added

src/juobs.jl

@@ -7,7 +7,7 @@ include("juobs_reader.jl")
 include("juobs_tools.jl")
 include("juobs_obs.jl")
 
-export read_mesons, read_ms1
+export read_mesons, read_ms1, read_ms
 export get_matrix, uwgevp_tot, energies, uwdot
 export corr_obs, plat_av, lin_fit, x_lin_fit, y_lin_fit, fit_routine
 export meff, dec_const_pcvc 

src/juobs_reader.jl

@@ -213,4 +213,60 @@ function read_ms1(path::String; v::String="1.2")
     [W[k, :] = prod(mean(exp.(.-r_data[k]), dims=2), dims=1) for k = 1:nrw]
     close(data)
     return W
-end
\ No newline at end of file
+end
+
+@doc raw"""
+read_ms(path::String)
+
+Reads openQCD ms dat files at a given path. This method return: 
+    t(t): flow time values
+    Wsl(x0, t, icfg): the time-slice sums of the densities of the Wilson plaquette action
+    Ysl(x0, t, icfg): the time-slice sums of the densities of the Yang-Mills action 
+    Qsl(x0, t, icfg): the time-slice sums of the densities of the topological charge
+
+Examples:
+```@example    
+t, W, Y, Q = read_ms(path)
+```
+"""
+function read_ms(path::String)  
+    data = open(path, "r")
+
+    dn = read(data, Int32)
+    nn = read(data, Int32)
+    tvals = read(data, Int32)
+    eps = read(data, Float64)
+
+    fsize = filesize(path)
+    datsize=4 + 3*8*(nn + 1) * tvals # measurement size of each cnfg
+
+    ncfg = Int32((fsize - 3*4 - 8) / datsize)
+
+    vcfg = Vector{Int32}(undef, ncfg)
+    
+    # x0, t, cfg 
+    Wsl = Array{Float64}(undef, tvals, nn + 1, ncfg) 
+    Ysl = Array{Float64}(undef, tvals, nn + 1, ncfg)
+    Qsl = Array{Float64}(undef, tvals, nn + 1, ncfg)
+
+    for icfg = 1:ncfg
+        vcfg[icfg] = read(data, Int32)
+        for iobs = 1:3
+            for inn = 0:nn
+                tmp = Vector{Float64}(undef, tvals)
+                read!(data, tmp)
+                if iobs == 1
+                    Wsl[:, inn + 1, icfg] = tmp
+                elseif iobs == 2
+                    Ysl[:, inn + 1, icfg] = tmp
+                elseif iobs == 3
+                    Qsl[:, inn + 1, icfg] = tmp
+                end
+
+            end
+        end
+    end
+    close(data)
+    t = Float64.(0:nn) .* dn .* eps
+    return (t, Wsl, Ysl, Qsl)
+end