 |
tthickness_distribution.jl - seaice-experiments - sea ice experiments using Gra… |
 |
git clone git://src.adamsgaard.dk/seaice-experiments |
 |
Log |
|
Files |
|
Refs |
|
README |
|
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… |
|
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 |
