tthickness_distribution.jl - seaice-experiments - sea ice experiments using Granular.jl
 (HTM) git clone git://src.adamsgaard.dk/seaice-experiments
 (DIR) Log
 (DIR) Files
 (DIR) Refs
 (DIR) README
 (DIR) LICENSE
       ---
       tthickness_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