thickness_distribution.jl (2026B)
1 import PyPlot 2 import ArgParse 3 import DelimitedFiles 4 import LsqFit 5 6 let 7 8 function parse_command_line() 9 s = ArgParse.ArgParseSettings() 10 ArgParse.@add_arg_table s begin 11 "--n" 12 help = "Number of bins" 13 arg_type = Int 14 default = 10 15 "tsvfile" 16 help = "input file to read" 17 required = true 18 end 19 return ArgParse.parse_args(s) 20 end 21 22 function plot_thickness_distribution(file, n) 23 data = DelimitedFiles.readdlm(file, '\t') 24 d = convert(Array{Float64, 1}, data[2:end, 1]) 25 x = convert(Array{Float64, 1}, data[2:end, 2]) 26 y = convert(Array{Float64, 1}, data[2:end, 3]) 27 N = length(d) 28 29 x_min = minimum(x) 30 x_max = maximum(x) 31 supersampling = 10 32 dx = (x_max - x_min) / (n * supersampling) 33 thickness = zeros(n * supersampling) 34 35 for ix=1:(n * supersampling) 36 x0 = x_min + (ix - 1) * dx 37 x1 = x0 + dx 38 xu = -Inf 39 xd = +Inf 40 for i=1:N 41 if (x[i] >= x0 && x[i] < x1) 42 if (xd > y[i] - d[i]) 43 xd = y[i] - d[i] 44 end 45 if (xu < y[i] + d[i]) 46 xu = y[i] + d[i] 47 end 48 end 49 end 50 if (xu != -Inf && xd != Inf) 51 thickness[ix] = xu - xd 52 end 53 end 54 55 h = LinRange(0.0, maximum(thickness)*1.2, n + 1) 56 h_counts = zeros(n) 57 for i=1:n 58 for j=1:(n * supersampling) 59 if (h[i + 1] > thickness[j] && h[i] <= thickness[j]) 60 h_counts[i] += 1 61 end 62 end 63 end 64 println(h) 65 println(h_counts) 66 g_h = h_counts / sum(h_counts) 67 println(g_h) 68 h_midpoints = (h[2:end] - h[1:end-1])*0.5 + h[1:end-1] 69 70 model(t, p) = p[1] * exp.(-t/p[2]) 71 p0 = [0.5, 0.5] 72 fit = LsqFit.curve_fit(model, h_midpoints, g_h, p0) 73 74 h_range = LinRange(0.0, maximum(h), 100) 75 PyPlot.plot(h_midpoints, g_h, marker="o", linestyle="") 76 PyPlot.plot(h_range, model(h_range, fit.param), linestyle="-") 77 PyPlot.xlabel("Thickness \$h\$ [m]") 78 PyPlot.ylabel("Thickness distribution \$g(h)\$ [m\$^{-1}\$]") 79 PyPlot.tight_layout() 80 PyPlot.savefig(file * "-itd.pdf") 81 PyPlot.clf() 82 end 83 84 function main() 85 parsed_args = parse_command_line() 86 plot_thickness_distribution(parsed_args["tsvfile"], 87 parsed_args["n"]) 88 end 89 90 main() 91 end