### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 8efe4e7f-a63c-4a73-bea3-48c02b292d12
begin
    import Pkg
    # careful: this is _not_ a reproducible environment
    # activate the global environment
    Pkg.activate()

    using Plots, PlutoUI
	import HDF5
	import Pda
	using Statistics
end

# ╔═╡ e5f8844e-7a41-4ac6-a843-9904309f9fa8
file = "U10.000singlewallDT1.000.h5"

# ╔═╡ ad67db84-c555-438d-a4d1-aae9f81d3981
@bind frame PlutoUI.Slider(0:98, show_value=true)

# ╔═╡ f1c78cf0-a784-4273-a6c6-090f64cfed36
begin
	x = Pda.get_h5data(file, "config/$frame/x")
	theta = Pda.get_h5data(file, "config/$frame/theta")
end

# ╔═╡ be3457aa-a2a0-43e6-b240-fa0e076d396c
function compute_inf_linedensity(file; frame::Int=80, )
		right = Pda.get_h5attr(file, "right")
		threshold = right/2.0
		Lx = right-threshold
		xin = [xi for xi in x if xi > threshold]
		ninf = length(xin) / Lx
	return ninf
end

# ╔═╡ c6cd5b84-d00d-4898-a769-e9ebd819ac93
begin
	U0 = Pda.get_h5attr(file, "U0")
	tauR = Pda.get_h5attr(file, "tauR")
	ell = U0*tauR
	DT = Pda.get_h5attr(file, "DT")
	delta = sqrt(DT*tauR)
	# Lx = 20 # box length in x 
	# # bulk line density, take only the latter half of the box
	# ninf = sum( x .>= Lx/2)/(Lx/2) 
	# # number of particles on the wall 
	# Nw = 2*ninf*ell
end

# ╔═╡ 09a21a89-e9a0-4236-88f8-d954d4165a55
sum( x .<= eps(Float64) ) /(compute_inf_linedensity(file)*ell)

# ╔═╡ ef496cd8-2a47-4c1e-9aac-c4fdbfd98ff6
begin
	qx_wall = Float64[]
	qy_wall = Float64[]
	theta_wall = Float64[]
	N = convert(Int, Pda.get_h5attr(file, "N") )
	for i in 1:N
		if x[i]  <= eps(Float64)
			push!(qx_wall, cos(theta[i]) )
			push!(qy_wall, sin(theta[i]) )
			push!(theta_wall, mod(theta[i], 2pi))
		end
	end
	h1 = histogram(qx_wall, label = "qx")
	h2 = histogram(qy_wall, label = "qy")
	plot(h1, h2, framestyle=:box, layout = (2,1) )
end

# ╔═╡ 65f836ea-1839-4c63-889c-9ad062ccf6b5
begin
	histogram(theta_wall, bins=50, framestyle=:box, xticks = pi*(0:1/2:2))
	let
		theta = LinRange(pi/2, 3pi/2, 100)
		lee = [	
		0.47883923959186386 0.09804957410222426;
		0.5198535717065695 0.2995920825021341;
		0.5576166579070441 0.5216095918830788;
		0.5936994405206624 0.6869355109518869;
		0.6330218896210373 0.8176132098108629;
		0.6706910537986635 0.9246701632248963;
		0.7116140370966255 1.014401719930961;
		0.7476234833362387 1.0899639172173865;
		0.7476234833362387 1.0899639172173865;
		0.788528454191463 1.1576482440119886;
		0.8277930061544658 1.2174598467269784;
		0.8277930061544658 1.2174598467269784;
		0.8654081330038772 1.2583751104066245;
		0.8654081330038772 1.2583751104066245;
		0.9013789811519084 1.286693243597059;
		0.9406126546416456 1.308709595035256;
		0.9765719233622026 1.3228545089968935;
		1.0239549405878874 1.3196675932352986;
		1.0549903794256732 1.3070447306850326;
		1.0925540421529758 1.2849679089033574;
		1.1301151316741729 1.2597414828486164;
		1.174194152260465 1.212462680470202;
		1.2084589647605855 1.1525931806178398;
		1.2492545743563286 1.0864193258071344;
		1.2884084784567968 1.0107979447802835;
		1.2884084784567968 1.0107979447802835;
		1.325912957443674 0.9162802247180892;
		1.325912957443674 0.9162802247180892;
		1.3634058570030767 0.8075892854270975;
		1.3634058570030767 0.8075892854270975;
		1.4008846039288998 0.6815755226342425;
		1.4008846039288998 0.6815755226342425;
		1.4416017307384266 0.5193387374950225;
		1.487170637659772 0.29568080922185613;
		1.5196843834051366 0.092505601548041;
		1.5196843834051366 0.092505601548041];
		#scatter!(lee[:,1]*pi, lee[:,2]*210, )
		plot!(theta, 40*cos.(theta).^2 )
		# This simple theory doesnt work!???? why
	end
end

# ╔═╡ f9faa467-6f14-46e6-84b5-fb9762214e00
let

scatter(lee[:,1]*pi, lee[:,2], )
theta = LinRange(0.5*pi, 1.5*pi,200)
plot!(theta, cos.(theta).^2*1.35)
end

# ╔═╡ e47394f4-c336-4937-aa0c-aa174bf409b4
 sum( x .> sqrt(eps(Float64)))

# ╔═╡ cffe8f9b-b7ab-448b-bf6c-c2123bdc2616
begin
	dx = Pda.get_h5data(file, "dx/m")
	dxtotal =mean( dx*N) # average over frames
	dt = Pda.get_h5attr(file, "dt")
	ninf = compute_inf_linedensity(file)
	# 1st scale by ninf kbt
	# 2nd scale by ninf ksts, for DT=0 this should be unity 
	dxtotal/(dt*ninf*DT), dxtotal/(dt*ninf*0.5*U0^2*tauR)
end

# ╔═╡ eaea4043-ad2e-4194-b2d7-52a21012ee4c
begin
 	1 + 0.5*(ell/delta)^2
end

# ╔═╡ 089cef24-5c87-4a59-b23e-cb7b952cd3d3
let
	y = [xi for xi in x if xi > 0.]
	
end

# ╔═╡ Cell order:
# ╠═8efe4e7f-a63c-4a73-bea3-48c02b292d12
# ╠═e5f8844e-7a41-4ac6-a843-9904309f9fa8
# ╠═ad67db84-c555-438d-a4d1-aae9f81d3981
# ╠═f1c78cf0-a784-4273-a6c6-090f64cfed36
# ╟─be3457aa-a2a0-43e6-b240-fa0e076d396c
# ╟─c6cd5b84-d00d-4898-a769-e9ebd819ac93
# ╠═09a21a89-e9a0-4236-88f8-d954d4165a55
# ╠═ef496cd8-2a47-4c1e-9aac-c4fdbfd98ff6
# ╠═65f836ea-1839-4c63-889c-9ad062ccf6b5
# ╠═f9faa467-6f14-46e6-84b5-fb9762214e00
# ╠═e47394f4-c336-4937-aa0c-aa174bf409b4
# ╠═cffe8f9b-b7ab-448b-bf6c-c2123bdc2616
# ╠═eaea4043-ad2e-4194-b2d7-52a21012ee4c
# ╠═089cef24-5c87-4a59-b23e-cb7b952cd3d3
