text:
### A Pluto.jl notebook ###
# v0.17.2

# using Markdown
# using InteractiveUtils

# ╔═╡ 35da0c45-6d91-4ad9-811f-5d633684208e
begin
	using StaticArrays
	using LinearAlgebra
	using Unitful
	using Unitful.DefaultSymbols
	using SatelliteToolbox: Ellipsoid
	using Parameters
	import Proj4
	using DocStringExtensions
end

# ╔═╡ c3a5e97e-00ed-44af-9c67-fa8198900fbd
#=╠═╡ notebook_exclusive
begin
	using PlutoDevMacros
	using PlutoUtils
	using PlutoTest
end
  ╠═╡ notebook_exclusive =#

# ╔═╡ e2f50ace-ad14-4e7c-af74-abf3ca9df9fb
#=╠═╡ notebook_exclusive
md"""
# Packages
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 82c0d6f9-373d-4866-81eb-9cd2d0981310
#=╠═╡ notebook_exclusive
ToC()
  ╠═╡ notebook_exclusive =#

# ╔═╡ 2e6bcf9b-002f-4fbc-80c1-0cfd2ab1d253
#=╠═╡ notebook_exclusive
md"""
# Exports
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 1e3da0c9-2f96-4637-9093-ac7f10c1ad27
#=╠═╡ notebook_exclusive
md"""
# Helper Functions
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 48c73104-fe4c-4543-942a-6f23b0fd2547
const Point{N} = Union{Tuple{Vararg{<:Number, N}},StaticVector{N,<:Number}}

# ╔═╡ 48c73104-fe4c-4543-0001-6f23b0fd2547
const Point2D = Point{2}

# ╔═╡ 48c73104-fe4c-4543-0002-6f23b0fd2547
const Point3D = Point{3}

# ╔═╡ c4402f72-67ac-4630-a651-da81c1df71bf
"""
	ExtraOutput
Struct used inside `SatView` to use multiple dispatch to give more than the standard output.
See for example [`get_range`](@ref)
"""
struct ExtraOutput end

# ╔═╡ e796d5be-e011-45d3-ad90-58769feb5e85
#=╠═╡ notebook_exclusive
md"""
## Show/Print
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ d890aff1-dbd0-451c-bf14-bde9758c3be0
function _print_angle(io,val,displayname,last=false)
	print(io,"$displayname=")
	print(io,round(val;digits=2) * rad)
	print(io," (")
	print(io,round(rad2deg(val);digits=2) * °)
	print(io,")")
	last || print(io,", ")
end

# ╔═╡ fdbbc8d9-83d6-413e-aff8-aca18f24dfea
function _print_length(io,val,displayname,last=false)
	print(io,"$displayname=")
	mval = val < 1000 ? round(val;digits=2) * m : round(val/1000;digits=2) * km
	print(io,mval)
	last || print(io,", ")
end

# ╔═╡ cc0cae75-ba10-4a62-b0ef-22259e40a083
#=╠═╡ notebook_exclusive
md"""
## Spherical Ellipsoid
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 855553e3-491b-4e8a-a482-95855697e063
"""
	SphericalEllipsoid(r = 6371e3)
Define a spherical ellipsoid of radius `r` to be used for the various transformation in [`SatView`](@ref). Defaults to 6371km radius
"""
SphericalEllipsoid(r = 6371e3) = Ellipsoid(r,0.0)

# ╔═╡ 1f7a7093-ce8e-461f-8b91-69266de86748
#=╠═╡ notebook_exclusive
md"""
## geod_geodesic
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ f83fc65f-5f7b-498d-9ed3-0762565ad710
Proj4.geod_geodesic(e::Ellipsoid) = Proj4.geod_geodesic(e.a, e.f)

# ╔═╡ 4714c6ae-27d9-47db-b12e-126283b10606
#=╠═╡ notebook_exclusive
md"""
# EarthModel
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ f9cc2880-bab1-4be3-b67b-d43508df8d3b
begin
"""
$(TYPEDEF)
Geometrical model of the Earth ellipsoid to be used for all the view angles between satellite and points on earth.
# Fields
$(TYPEDFIELDS)

A single instance of this structure should be used for all satellites of a given simulation/constellation.\\
Changes to any of the two fields (via `setproperty!`) will trigger an automatic recomputation of the other field.\\
When called without arguments, it defaults to a spherical earth with a radius of 6371 km.

See also: [`Ellipsoid`](@ref), [`SphericalEllipsoid`](@ref), [`SatView`](@ref)
"""
@with_kw mutable struct EarthModel
	"Ellipsoid structure, used for the various point of view conversion."
	ellipsoid::Ellipsoid{Float64} = SphericalEllipsoid()
	"Extended geod structure, used for the inverse geodesic problem to compute distance and azimuth between points on earth. (Relies on GeographicLib)."
	geod::Proj4.geod_geodesic = Proj4.geod_geodesic(ellipsoid)
end
EarthModel(e::Ellipsoid) = EarthModel(;ellipsoid = e)
end

# ╔═╡ 1f2b9b74-de46-401d-8a46-0434b9f9aca1
function Base.setproperty!(value::EarthModel, name::Symbol, x)
	if name === :ellipsoid
		setfield!(value, name, x)
		setfield!(value, :geod, Proj4.geod_geodesic(x))
	else
		setfield!(value, name, x)
		setfield!(value, :ellipsoid, Ellipsoid(x.a, x.f))
	end
end

# ╔═╡ 6a5cb372-60cb-4ffc-b4f0-22e4016104e7
#=╠═╡ notebook_exclusive
md"""
# Angle Types
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ e832c1b7-8c04-4146-90f6-1628e91fea2a
const UnitfulAngleType = Union{typeof(°),typeof(rad)}

# ╔═╡ 64cf1b8b-6686-4eeb-a3cc-786300ea7c7d
const UnitfulAngleQuantity = Quantity{<:Real,<:Any,<:UnitfulAngleType}

# ╔═╡ 1d023a0c-a93a-451c-a894-1d1f6a4b78a9
#=╠═╡ notebook_exclusive
md"""
# SatViewCoordinate types
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ b8ce87d4-4768-4e8a-a16e-9e68b00b6617
abstract type SatViewCoordinate end

# ╔═╡ e3c221c6-6c4a-4b5f-93a6-30d3508ac9d2
#=╠═╡ notebook_exclusive
md"""
## LLA
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 8368ae01-ce53-449e-87dd-8dfa3f29f8f4
#=╠═╡ notebook_exclusive
md"""
Here we want to define a structure that contains useful informations and functions to perform conversions between the view from the satellite based on it's orbital position and points on ground
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ f951805e-515a-475f-893f-bb8b968e425c
#=╠═╡ notebook_exclusive
md"""
## ERA
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 86ae20a9-e69c-4d63-9119-395449e9ac09
#=╠═╡ notebook_exclusive
md"""
ERA stands for Elevation, Range and Azimuth and is used to express the position of a satellite relative to an observer in spherical coordinates.
The elevation is the angle of the pointing with respect to the local horizon of the observer, meaning the plane where the observer is located that is perpendicular to the gravity vector acting on the observe (or in an alternative definition, the plane where the observer is located that is parallel to the tangent plane to the earth ellipsoid at the given lat and lon positions of the observer.
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 9be2fd5c-4b6c-4e13-b2aa-fb7120a504b7
begin
	"""
	Elevation, Range and Azimuth for a target point on space as seen from a source point on or above the earth surface

	# Fields
	- `el::Float64`: Elevation view angle (`0 <= el <= π/2`) between source and target point [rad].
	- `r::Float64`: Range (`r >= 0`) between the source and target points [m].
	- `az::Float64`: Azimuth view angle between source and target point [rad], computed from West to North from the source point perspective. Values provided are automatically converted between -π and π

	# Constructors

		ERA(el::Real,r::Real,az::Real)
		ERA(el::UnitfulAngleQuantity,r::Real,az::Real)
		ERA(el::UnitfulAngleQuantity,r::Real,az::UnitfulAngleQuantity)
		ERA(el,r::Unitful.Length,az)

	where `UnitfulAngleQuantity` is a `Unitful.Quantity` of unit either `u"rad"` or `u"°"`.
	"""
	@with_kw_noshow struct ERA <: SatViewCoordinate
		el::Float64 # Elevation in radians
		r::Float64 # Range in meters
		az::Float64 # Azimuth in radians

		function ERA(el::Real,r::Real,az::Real)
			(isnan(el) || isnan(r) || isnan(az)) && return new(NaN,NaN,NaN)
			@assert el >= 0 && el <= π/2 "Elevation should be between 0 and π/2"
			@assert r >= 0 "Range must be positive"
			new(el,r,rem2pi(az,RoundNearest))
		end
	end

	# Define a constructor that takes combinations of real numbers and angles/lengths
	ERA(el::UnitfulAngleQuantity,r::Real,az::Real) = ERA(
		uconvert(u"rad",el) |> ustrip,
		r,
		az,
	)
	ERA(el::UnitfulAngleQuantity,r::Real,az::UnitfulAngleQuantity) = ERA(
		el,
		r,
		uconvert(u"rad",az) |> ustrip,
	)
	ERA(el,r::Unitful.Length,az) = ERA(
		el,
		uconvert(u"m",r) |> ustrip,
		az,
	)

	# Show
	function Base.show(io::IO,era::ERA)
		print(io,"ERA(")
		_print_angle(io,era.el,"el",false)
		_print_length(io,era.r,"r",false)
		_print_angle(io,era.az,"az",true)
		print(io,")")
	end
end

# ╔═╡ 16782c72-ecb1-48ec-8510-78e2e0689a10
function Base.isapprox(x1::ERA, x2::ERA)
	x1.el ≉ x2.el && return false
	# Don't care about different azimuth if elevation is 90°
	x2.el ≉ π/2 && x1.az ≉ x2.az && return false
	x1.r ≉ x2.r && return false
	return true
end

# ╔═╡ 7344190c-7989-4b55-b7be-357f7d6b7370
Base.isnan(era::ERA) = isnan(era.el)

# ╔═╡ f207d849-ebff-4e6c-95bb-50693cb7c9b6
begin
	"""
	Identify a point on or above earth using geodetic coordinates

	# Fields
	- `lat::Float64`: Latitude (`-π/2 <= lat <= π/2`) of the point [rad].
	- `lon::Float64`: Longitude of the point [rad].
	- `alt::Float64`: Altitude of the point above the reference earth ellipsoid [m].

	# Constructors
		LLA(lat::Real,lon::Real,alt::Real)
		LLA(lat::UnitfulAngleQuantity,lon::Real,alt::Real)
		LLA(lat::UnitfulAngleQuantity,lon::UnitfulAngleQuantity,alt::Real)
		LLA(lat,lon,alt::Unitful.Length)
		LLA(lat,lon) # Defaults to 0.0 altitude

	where `UnitfulAngleQuantity` is a `Unitful.Quantity` of unit either `u"rad"` or `u"°"`.
	"""
	@with_kw_noshow struct LLA <: SatViewCoordinate
		lat::Float64 # Latitude in radians
		lon::Float64 # Longitude in radians
		alt::Float64 # Altitude in meters

		function LLA(lat::Real,lon::Real,alt::Real)
			(isnan(lat) || isnan(lon) || isnan(alt)) && return new(NaN,NaN,NaN)
			l2 = rem2pi(lon,RoundNearest)
			@assert abs(lat) <= π/2 "Latitude should be between -π/2 and π/2"
			new(lat,l2,alt)
		end
	end

	# Constructor without altitude, assume it is 0
	LLA(lat,lon) = LLA(lat,lon,0.0)

	# Define a constructor that takes combinations of real numbers and angles/lengths
	LLA(lat::UnitfulAngleQuantity,lon::Real,alt::Real) = LLA(
		uconvert(u"rad",lat) |> ustrip,
		lon,
		alt)
	LLA(lat::UnitfulAngleQuantity,lon::UnitfulAngleQuantity,alt::Real) = LLA(
		lat,
		uconvert(u"rad",lon) |> ustrip,
		alt)
	LLA(lat,lon,alt::Unitful.Length) = LLA(
		lat,
		lon,
		uconvert(u"m",alt) |> ustrip)
end

# ╔═╡ 3033d5c1-d8e0-4d46-0001-7dec4ff7afbd
# Show method for LLA
function Base.show(io::IO,lla::LLA)
	print(io,"LLA(")
	_print_angle(io,lla.lat,"lat",false)
	_print_angle(io,lla.lon,"lon",false)
	_print_length(io,lla.alt,"alt",true)
	print(io,")")
end

# ╔═╡ 3033d5c1-d8e0-4d46-a4b9-7dec4ff7afbd
function Base.isapprox(x1::LLA, x2::LLA)
	x1.alt ≉ x2.alt && return false
	# Don't care about different longitude if latitude is ±90°
	abs(x1.lat) ≈ π/2 && abs(x2.lat) ≈ π/2 && return true
	# Return true if all the lat and lon are matching
	x1.lat ≈ x2.lat && (x1.lon ≈ x2.lon || abs(x1.lon) ≈ abs(x2.lon) ≈ π) && return true
	return false
end

# ╔═╡ 528beffe-0707-4661-8970-def1b1e00ea5
Base.isnan(lla::LLA) = isnan(lla.lat)

# ╔═╡ 11938cb6-46b3-0001-96c0-ef6424d1d0db
#=╠═╡ notebook_exclusive
md"""
# Inverse Geodesic Problem
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 11938cb6-46b3-0002-96c0-ef6424d1d0db
"""
    geod_inverse(geod::Proj4.geod_geodesic, lonlat1::AbstractVector{Cdouble}, lonlat2::AbstractVector{Cdouble})
	geod_inverse(geod::Proj4.geod_geodesic, lla1::LLA, lla2::LLA)
Solve the inverse geodesic problem.

# Args

- `g`       → the geod_geodesic object specifying the ellipsoid.
- `lonlat1` → point 1 (degrees), where lat ∈ [-90, 90], lon ∈ [-540, 540)
- `lonlat2` → point 2 (degrees), where lat ∈ [-90, 90], lon ∈ [-540, 540)

# Outputs

- `dist` → distance between point 1 and point 2 (meters).
- `azi1` → azimuth at point 1 (degrees) ∈ [-180, 180)
- `azi2` → (forward) azimuth at point 2 (degrees) ∈ [-180, 180)

# Remarks

If either point is at a pole, the azimuth is defined by keeping the longitude fixed,
writing lat = 90 +/- eps, and taking the limit as eps -> 0+.
"""
function geod_inverse(geod::Proj4.geod_geodesic, lonlat1::AbstractVector{Cdouble}, lonlat2::AbstractVector{Cdouble})
	dist = Ref{Cdouble}()
	azi1 = Ref{Cdouble}()
	azi2 = Ref{Cdouble}()
	ccall((:geod_inverse, Proj4.libproj), Cvoid, (Ptr{Cvoid},Cdouble,Cdouble,Cdouble,
			Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),
			pointer_from_objref(geod), lonlat1[2], lonlat1[1], lonlat2[2], lonlat2[1], dist, azi1, azi2)
	dist[], azi1[], azi2[]
end

# ╔═╡ 11938cb6-46b3-0003-96c0-ef6424d1d0db
function geod_inverse(geod::Proj4.geod_geodesic, lla1::LLA, lla2::LLA)
	lonlat1 = rad2deg.(SA_F64[lla1.lon,lla1.lat])
	lonlat2 = rad2deg.(SA_F64[lla2.lon,lla2.lat])
	geod_inverse(geod,lonlat1,lonlat2)
end

# ╔═╡ 4c06d21c-ac14-4522-bf25-2e0a1ed2d6b9
begin
	export Ellipsoid, SphericalEllipsoid
	export EarthModel
	export UnitfulAngleQuantity, UnitfulAngleType, °
	export km
	export LLA, ERA
	export geod_inverse
end

# ╔═╡ 11938cb6-46b3-499b-96c0-ef6424d1d0db
#=╠═╡ notebook_exclusive
md"""
# Tests
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ d2c248b1-c48e-437b-a910-edcc59b4424f
#=╠═╡ notebook_exclusive
md"""
## LLA
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ 7b306ed5-4bda-465d-abf2-4d07cb4642c1
#=╠═╡ notebook_exclusive
@test LLA(10°,10°,1000) ≈ LLA((10+100*eps())*°,10°,1000)
  ╠═╡ notebook_exclusive =#

# ╔═╡ c45f6aac-bff3-4c98-8bd5-c91a98c9eef7
#=╠═╡ notebook_exclusive
@test LLA(90°,10°,1000) ≈ LLA(90°,130°,1000)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 1f1505a6-c6af-4fdd-9583-6e783e76de4f
#=╠═╡ notebook_exclusive
@test LLA(40°,-180°,1000) ≈ LLA(40°,180°,1000)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 980e48bd-9dd1-4195-8086-40785a7f43e1
#=╠═╡ notebook_exclusive
@test LLA(10°,10°,1000) !== LLA((10+100*eps())*°,10°,1000)
  ╠═╡ notebook_exclusive =#

# ╔═╡ fc146750-46b2-4084-a6d2-0d91c0e104e6
#=╠═╡ notebook_exclusive
LLA(1,1,NaN) |> isnan
  ╠═╡ notebook_exclusive =#

# ╔═╡ fc3e96a6-7af8-4ef0-a4ca-b33509aa512f
#=╠═╡ notebook_exclusive
md"""
## ERA
"""
  ╠═╡ notebook_exclusive =#

# ╔═╡ ed2f2fa3-be19-4e2d-ad28-36cb053ed6bd
#=╠═╡ notebook_exclusive
@test ERA(10°,1000,20°) == ERA(10°,1km,deg2rad(20)*rad)
  ╠═╡ notebook_exclusive =#

# ╔═╡ 6151e4d1-a7a6-4fa7-bc77-0e7f3b2a3cc0
#=╠═╡ notebook_exclusive
@test ERA(90°,1000,20°) ≈ ERA(90°,1km,deg2rad(90)*rad)
  ╠═╡ notebook_exclusive =#

# ╔═╡ b5449c6a-9c43-4f9c-81f7-f01277acb109
#=╠═╡ notebook_exclusive
ERA(1,1,NaN) |> isnan
  ╠═╡ notebook_exclusive =#

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
PlutoDevMacros = "a0499f29-c39b-4c5c-807c-88074221b949"
PlutoTest = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
PlutoUtils = "ed5d0301-4775-4676-b788-cf71e66ff8ed"
Proj4 = "9a7e659c-8ee8-5706-894e-f68f43bc57ea"
SatelliteToolbox = "6ac157d9-b43d-51bb-8fab-48bf53814f4a"
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"

[compat]
DocStringExtensions = "~0.8.6"
Parameters = "~0.12.3"
PlutoDevMacros = "~0.3.9"
PlutoTest = "~0.2.0"
PlutoUtils = "~0.4.13"
Proj4 = "~0.7.6"
SatelliteToolbox = "~0.9.3"
StaticArrays = "~1.2.13"
Unitful = "~1.9.2"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0-rc2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.AxisAlgorithms]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"]
git-tree-sha1 = "66771c8d21c8ff5e3a93379480a2307ac36863f7"
uuid = "13072b0f-2c55-5437-9ae7-d433b7a33950"
version = "1.0.1"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[deps.Chain]]
git-tree-sha1 = "cac464e71767e8a04ceee82a889ca56502795705"
uuid = "8be319e6-bccf-4806-a6f7-6fae938471bc"
version = "0.4.8"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f885e7e7c124f8c92650d61b9477b9ac2ee607dd"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.1"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[deps.CoordinateTransformations]]
deps = ["LinearAlgebra", "StaticArrays"]
git-tree-sha1 = "681ea870b918e7cff7111da58791d7f718067a19"
uuid = "150eb455-5306-5404-9cee-2592286d6298"
version = "0.6.2"

[[deps.Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileIO]]
deps = ["Pkg", "Requires", "UUIDs"]
git-tree-sha1 = "2db648b6712831ecb333eae76dbfd1c156ca13bb"
uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
version = "1.11.2"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Glob]]
git-tree-sha1 = "4df9f7e06108728ebf00a0a11edee4b29a482bb2"
uuid = "c27321d9-0574-5035-807b-f59d2c89b15c"
version = "1.3.0"

[[deps.HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.Interpolations]]
deps = ["AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "Requires", "SharedArrays", "SparseArrays", "StaticArrays", "WoodburyMatrices"]
git-tree-sha1 = "61aa005707ea2cebf47c8d780da8dc9bc4e0c512"
uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
version = "0.13.4"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "043017e0bdeff61cfbb7afeb558ab29536bbb5ed"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.8"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OptionalData]]
git-tree-sha1 = "d047cc114023e12292533bb822b45c23cb51d310"
uuid = "fbd9d27c-2d1c-5c1c-99f2-7497d746985d"
version = "1.0.0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.PROJ_jll]]
deps = ["Artifacts", "JLLWrappers", "LibCURL_jll", "LibSSH2_jll", "Libdl", "Libtiff_jll", "MbedTLS_jll", "Pkg", "SQLite_jll", "Zlib_jll", "nghttp2_jll"]
git-tree-sha1 = "2435e91710d7f97f53ef7a4872bf1f948dc8e5f8"
uuid = "58948b4f-47e0-5654-a9ad-f609743f8632"
version = "700.202.100+0"

[[deps.Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "34c0e9ad262e5f7fc75b10a9952ca7692cfc5fbe"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.3"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoDevMacros]]
deps = ["MacroTools", "PlutoHooks"]
git-tree-sha1 = "6ab70183795e4ad00ecd406f2d05740328c9b905"
uuid = "a0499f29-c39b-4c5c-807c-88074221b949"
version = "0.3.9"

[[deps.PlutoHooks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "f297787f7d7507dada25f6769fe3f08f6b9b8b12"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.3"

[[deps.PlutoTest]]
deps = ["HypertextLiteral", "InteractiveUtils", "Markdown", "Test"]
git-tree-sha1 = "92b8ae1eee37c1b8f70d3a8fb6c3f2d81809a1c5"
uuid = "cb4044da-4d16-4ffa-a6a3-8cad7f73ebdc"
version = "0.2.0"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "1e0cb51e0ccef0afc01aab41dc51a3e7f781e8cb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.20"

[[deps.PlutoUtils]]
deps = ["Chain", "Glob", "HypertextLiteral", "InteractiveUtils", "Markdown", "PlutoDevMacros", "PlutoHooks", "PlutoTest", "PlutoUI", "PrettyTables", "Reexport", "Requires", "UUIDs"]
git-tree-sha1 = "3d3856ecfea340b4ee0c77e5c3228dd1b4478ae1"
uuid = "ed5d0301-4775-4676-b788-cf71e66ff8ed"
version = "0.4.13"

[[deps.PolynomialRoots]]
git-tree-sha1 = "5f807b5345093487f733e520a1b7395ee9324825"
uuid = "3a141323-8675-5d76-9d11-e1df1406c778"
version = "1.0.0"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "d940010be611ee9d67064fe559edbb305f8cc0eb"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.Proj4]]
deps = ["CEnum", "CoordinateTransformations", "PROJ_jll", "StaticArrays"]
git-tree-sha1 = "5f15f1c647b563e49f655fbbfd4e2ade24bd3c64"
uuid = "9a7e659c-8ee8-5706-894e-f68f43bc57ea"
version = "0.7.6"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Ratios]]
deps = ["Requires"]
git-tree-sha1 = "01d341f502250e81f6fec0afe662aa861392a3aa"
uuid = "c84ed2f1-dad5-54f0-aa8e-dbefe2724439"
version = "0.4.2"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.ReferenceFrameRotations]]
deps = ["Crayons", "LinearAlgebra", "Printf", "StaticArrays"]
git-tree-sha1 = "d526371cec370888f485756a4bf8284ab531860b"
uuid = "74f56ac7-18b3-5285-802d-d4bd4f104033"
version = "1.0.1"

[[deps.RemoteFiles]]
deps = ["Dates", "FileIO", "HTTP"]
git-tree-sha1 = "54527375d877a64c55190fb762d584f927d6d7c3"
uuid = "cbe49d4c-5af1-5b60-bb70-0a60aa018e1b"
version = "0.4.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.SQLite_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "e601efdf2ce1df61147a1d1fed83cc428af1b53a"
uuid = "76ed43ae-9a5d-5a62-8c75-30186b810ce8"
version = "3.36.1+0"

[[deps.SatelliteToolbox]]
deps = ["Crayons", "Dates", "DelimitedFiles", "Interpolations", "LinearAlgebra", "OptionalData", "Parameters", "PolynomialRoots", "PrettyTables", "Printf", "Reexport", "ReferenceFrameRotations", "RemoteFiles", "SparseArrays", "StaticArrays", "Statistics"]
git-tree-sha1 = "0a2c0f1565a51487fe58c28f528675dba1008432"
uuid = "6ac157d9-b43d-51bb-8fab-48bf53814f4a"
version = "0.9.3"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "fed34d0e71b91734bf0a7e10eb1bb05296ddbcd0"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Unitful]]
deps = ["ConstructionBase", "Dates", "LinearAlgebra", "Random"]
git-tree-sha1 = "0992ed0c3ef66b0390e5752fe60054e5ff93b908"
uuid = "1986cc42-f94f-5a68-af5c-568840ba703d"
version = "1.9.2"

[[deps.WoodburyMatrices]]
deps = ["LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "de67fa59e33ad156a590055375a30b23c40299d3"
uuid = "efce3f68-66dc-5838-9240-27a6d6f5f9b6"
version = "0.5.5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e2f50ace-ad14-4e7c-af74-abf3ca9df9fb
# ╠═35da0c45-6d91-4ad9-811f-5d633684208e
# ╠═c3a5e97e-00ed-44af-9c67-fa8198900fbd
# ╠═82c0d6f9-373d-4866-81eb-9cd2d0981310
# ╟─2e6bcf9b-002f-4fbc-80c1-0cfd2ab1d253
# ╠═4c06d21c-ac14-4522-bf25-2e0a1ed2d6b9
# ╟─1e3da0c9-2f96-4637-9093-ac7f10c1ad27
# ╠═48c73104-fe4c-4543-942a-6f23b0fd2547
# ╠═48c73104-fe4c-4543-0001-6f23b0fd2547
# ╠═48c73104-fe4c-4543-0002-6f23b0fd2547
# ╠═c4402f72-67ac-4630-a651-da81c1df71bf
# ╟─e796d5be-e011-45d3-ad90-58769feb5e85
# ╠═d890aff1-dbd0-451c-bf14-bde9758c3be0
# ╠═fdbbc8d9-83d6-413e-aff8-aca18f24dfea
# ╟─cc0cae75-ba10-4a62-b0ef-22259e40a083
# ╠═855553e3-491b-4e8a-a482-95855697e063
# ╟─1f7a7093-ce8e-461f-8b91-69266de86748
# ╠═f83fc65f-5f7b-498d-9ed3-0762565ad710
# ╟─4714c6ae-27d9-47db-b12e-126283b10606
# ╠═f9cc2880-bab1-4be3-b67b-d43508df8d3b
# ╠═1f2b9b74-de46-401d-8a46-0434b9f9aca1
# ╟─6a5cb372-60cb-4ffc-b4f0-22e4016104e7
# ╠═e832c1b7-8c04-4146-90f6-1628e91fea2a
# ╠═64cf1b8b-6686-4eeb-a3cc-786300ea7c7d
# ╟─1d023a0c-a93a-451c-a894-1d1f6a4b78a9
# ╠═b8ce87d4-4768-4e8a-a16e-9e68b00b6617
# ╟─e3c221c6-6c4a-4b5f-93a6-30d3508ac9d2
# ╟─8368ae01-ce53-449e-87dd-8dfa3f29f8f4
# ╠═f207d849-ebff-4e6c-95bb-50693cb7c9b6
# ╠═3033d5c1-d8e0-4d46-0001-7dec4ff7afbd
# ╠═3033d5c1-d8e0-4d46-a4b9-7dec4ff7afbd
# ╠═528beffe-0707-4661-8970-def1b1e00ea5
# ╟─f951805e-515a-475f-893f-bb8b968e425c
# ╟─86ae20a9-e69c-4d63-9119-395449e9ac09
# ╠═9be2fd5c-4b6c-4e13-b2aa-fb7120a504b7
# ╠═16782c72-ecb1-48ec-8510-78e2e0689a10
# ╠═7344190c-7989-4b55-b7be-357f7d6b7370
# ╟─11938cb6-46b3-0001-96c0-ef6424d1d0db
# ╟─11938cb6-46b3-0002-96c0-ef6424d1d0db
# ╟─11938cb6-46b3-0003-96c0-ef6424d1d0db
# ╟─11938cb6-46b3-499b-96c0-ef6424d1d0db
# ╟─d2c248b1-c48e-437b-a910-edcc59b4424f
# ╠═7b306ed5-4bda-465d-abf2-4d07cb4642c1
# ╠═c45f6aac-bff3-4c98-8bd5-c91a98c9eef7
# ╠═1f1505a6-c6af-4fdd-9583-6e783e76de4f
# ╠═980e48bd-9dd1-4195-8086-40785a7f43e1
# ╠═fc146750-46b2-4084-a6d2-0d91c0e104e6
# ╟─fc3e96a6-7af8-4ef0-a4ca-b33509aa512f
# ╠═ed2f2fa3-be19-4e2d-ad28-36cb053ed6bd
# ╠═6151e4d1-a7a6-4fa7-bc77-0e7f3b2a3cc0
# ╠═b5449c6a-9c43-4f9c-81f7-f01277acb109
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: disberd/TelecomUtils

max_issues_repo_name: disberd/TelecomUtils

max_forks_repo_name: disberd/TelecomUtils


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec <: SwaggerModel
    distinguisherMethod::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1FlowDistinguisherMethod } # spec name: distinguisherMethod
    matchingPrecedence::Any # spec type: Union{ Nothing, Int32 } # spec name: matchingPrecedence
    priorityLevelConfiguration::Any # spec type: Union{ Nothing, IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationReference } # spec name: priorityLevelConfiguration
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects} } # spec name: rules

    function IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec(;distinguisherMethod=nothing, matchingPrecedence=nothing, priorityLevelConfiguration=nothing, rules=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec, Symbol("distinguisherMethod"), distinguisherMethod)
        setfield!(o, Symbol("distinguisherMethod"), distinguisherMethod)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec, Symbol("matchingPrecedence"), matchingPrecedence)
        setfield!(o, Symbol("matchingPrecedence"), matchingPrecedence)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec, Symbol("priorityLevelConfiguration"), priorityLevelConfiguration)
        setfield!(o, Symbol("priorityLevelConfiguration"), priorityLevelConfiguration)
        validate_property(IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec

const _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec = Dict{Symbol,Symbol}(Symbol("distinguisherMethod")=>Symbol("distinguisherMethod"), Symbol("matchingPrecedence")=>Symbol("matchingPrecedence"), Symbol("priorityLevelConfiguration")=>Symbol("priorityLevelConfiguration"), Symbol("rules")=>Symbol("rules"))
const _property_types_IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec = Dict{Symbol,String}(Symbol("distinguisherMethod")=>"IoK8sApiFlowcontrolV1alpha1FlowDistinguisherMethod", Symbol("matchingPrecedence")=>"Int32", Symbol("priorityLevelConfiguration")=>"IoK8sApiFlowcontrolV1alpha1PriorityLevelConfigurationReference", Symbol("rules")=>"Vector{IoK8sApiFlowcontrolV1alpha1PolicyRulesWithSubjects}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec)
    (getproperty(o, Symbol("priorityLevelConfiguration")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1FlowSchemaSpec }, name::Symbol, val)
end

ext: jl


text:
using JLD2
using FileIO
using Statistics

# Helper functions
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

datapath = joinpath(@__DIR__, "data")
fname = joinpath(datapath, "data.jld2")
content = load(fname)
t = content["t"]
x = content["x"]
bits = content["bits"]

spb = length(x) ÷ length(bits)
ts = Iterators.partition(t, spb) |> collect
xs = Iterators.partition(x, spb) |> collect

d = 3
idx = [2, 3, 6, 7]
irange = 1:length(bits)
smats = Vector{Matrix{Float64}}(undef, length(irange))
for (m, i) in enumerate(irange)
    ti = ts[i]
    xi = xs[i]
    mat = zeros(length(idx), length(idx))
    for (k, i) in enumerate(idx)
        for (l, j) in enumerate(idx)
            ekl = abs.(getwave(xi, i, d) - getwave(xi, j, d))
            mat[k, l] = sum(ekl .^ 2) / length(ekl)
        end
    end
    smats[m] = mat .< mean(mat)
end

S0 = [
    1 0 0 1
    0 1 1 0
    0 1 1 0
    1 0 0 1
]

S1 = [
    1 1 0 0
    1 1 0 0
    0 0 1 1
    0 0 1 1
]

rbits = [argmin([norm(mat - S, 1) for S in [S0, S1]]) for mat in smats] .- 1 .|> Bool
@show ber = sum(bits .!= rbits) / length(rbits);

# Visual check
figure()
for i in irange
    ti = ts[i]
    xi = xs[i]
    ei = abs.(getwave(xi, 1, d) - getwave(xi, 8, d))
    plot(ti, ei)
    axvline.([first(ti), last(ti)], linestyle = "dashed", color = "red")
end


text:
using ITensors
using ITensorsVisualization

N = 10
s(n) = Index([QN("Sz", 0) => 1, QN("Sz", 1) => 1]; tags = "S=1/2,Site,n=$n")
l(n) = Index([QN("Sz", 0) => 10, QN("Sz", 1) => 10]; tags = "Link,l=$n")
h(n) = Index([QN("Sz", 0) => 5, QN("Sz", 1) => 5]; tags = "ham,Link,l=$n")
s⃗ = [s(n) for n in 1:N]
l⃗ = [l(n) for n in 1:N-1]
h⃗ = [h(n) for n in 1:N-1]

# Add some more indices between two of the tensors
x = Index([QN("Sz", 0) => 2]; tags = "X")
y = Index([QN("Sz", 0) => 2]; tags = "Y")

n = 2
ψn1n2 = randomITensor(l⃗[n-1], s⃗[n], s⃗[n+1], l⃗[n+1], dag(x), dag(y))
hn1 = randomITensor(dag(h⃗[n-1]), s⃗[n]', dag(s⃗[n]), h⃗[n], x, y)
hn2 = randomITensor(dag(h⃗[n]), s⃗[n+1]', dag(s⃗[n+1]), h⃗[n+1])
ELn0 = randomITensor(l⃗[n-1]', h⃗[n-1], dag(l⃗[n-1]))
ERn2 = randomITensor(l⃗[n+1]', dag(h⃗[n+1]), dag(l⃗[n+1]))

R = @visualize ELn0 * ψn1n2 * hn1 * hn2 * ERn2 pause = true
@show R ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

# Split it up into multiple contractions
R1 = @visualize ELn0 * ψn1n2 * hn1 pause = true
R2 = @visualize R1 * hn2 * ERn2 pause = true
@show R2 ≈ ELn0 * ψn1n2 * hn1 * hn2 * ERn2

max_stars_repo_name: mtfishman/ITensorsVisualization

max_issues_repo_name: mtfishman/ITensorsVisualization

max_forks_repo_name: mtfishman/ITensorsVisualization


text:
# Invertible HINT coupling layer from Kruse et al. (2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export CouplingLayerHINT

"""
    H = CouplingLayerHINT(nx, ny, n_in, n_hidden, batchsize;
        logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1) (2D)

    H = CouplingLayerHINT(nx, ny, nz, n_in, n_hidden, batchsize;
        logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1) (3D)

 Create a recursive HINT-style invertible layer based on coupling blocks.

 *Input*:

 - `nx`, `ny`, `nz`: spatial dimensions of input

 - `n_in`, `n_hidden`: number of input and hidden channels

 - `logdet`: bool to indicate whether to return the log determinant. Default is `false`.

 - `permute`: string to specify permutation. Options are `"none"`, `"lower"`, `"both"` or `"full"`.

 - `k1`, `k2`: kernel size of convolutions in residual block. `k1` is the kernel of the first and third
    operator, `k2` is the kernel size of the second operator.

 - `p1`, `p2`: padding for the first and third convolution (`p1`) and the second convolution (`p2`)

 - `s1`, `s2`: stride for the first and third convolution (`s1`) and the second convolution (`s2`)

 *Output*:

 - `H`: Recursive invertible HINT coupling layer.

 *Usage:*

 - Forward mode: `Y = H.forward(X)`

 - Inverse mode: `X = H.inverse(Y)`

 - Backward mode: `ΔX, X = H.backward(ΔY, Y)`

 *Trainable parameters:*

 - None in `H` itself

 - Trainable parameters in coupling layers `H.CL`

 See also: [`CouplingLayerBasic`](@ref), [`ResidualBlock`](@ref), [`get_params`](@ref), [`clear_grad!`](@ref)
"""
mutable struct CouplingLayerHINT <: NeuralNetLayer
    CL::AbstractArray{CouplingLayerBasic, 1}
    C::Union{Conv1x1, Nothing}
    logdet::Bool
    permute::String
    is_reversed::Bool
end

@Flux.functor CouplingLayerHINT

# Get layer depth for recursion
function get_depth(n_in)
    count = 0
    nc = n_in
    while nc > 4
        nc /= 2
        count += 1
    end
    return count +1
end

# Constructor for given coupling layer and 1 x 1 convolution
CouplingLayerHINT(CL::AbstractArray{CouplingLayerBasic, 1}, C::Union{Conv1x1, Nothing};
    logdet=false, permute="none") = CouplingLayerHINT(CL, C, logdet, permute, false)

# 2D Constructor from input dimensions
function CouplingLayerHINT(nx::Int64, ny::Int64, n_in::Int64, n_hidden::Int64, batchsize::Int64;
    logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

    # Create basic coupling layers
    n = get_depth(n_in)
    CL = Array{CouplingLayerBasic}(undef, n)
    for j=1:n
        CL[j] = CouplingLayerBasic(nx, ny, Int(n_in/2^j), n_hidden, batchsize;
            k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet)
    end

    # Permutation using 1x1 convolution
    if permute == "full" || permute == "both"
        C = Conv1x1(n_in)
    elseif permute == "lower"
        C = Conv1x1(Int(n_in/2))
    else
        C = nothing
    end

    return CouplingLayerHINT(CL, C, logdet, permute, false)
end

# 3D Constructor from input dimensions
function CouplingLayerHINT(nx::Int64, ny::Int64, nz::Int64, n_in::Int64, n_hidden::Int64, batchsize::Int64;
    logdet=false, permute="none", k1=3, k2=3, p1=1, p2=1, s1=1, s2=1)

    # Create basic coupling layers
    n = get_depth(n_in)
    CL = Array{CouplingLayerBasic}(undef, n)
    for j=1:n
        CL[j] = CouplingLayerBasic(nx, ny, nz, Int(n_in/2^j), n_hidden, batchsize;
            k1=k1, k2=k2, p1=p1, p2=p2, s1=s1, s2=s2, logdet=logdet)
    end

    # Permutation using 1x1 convolution
    if permute == "full" || permute == "both"
        C = Conv1x1(n_in)
    elseif permute == "lower"
        C = Conv1x1(Int(n_in/2))
    else
        C = nothing
    end

    return CouplingLayerHINT(CL, C, logdet, permute, false)
end

# Input is tensor X
function forward(X, H::CouplingLayerHINT; scale=1, permute=nothing, logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && ~H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    if permute == "full" || permute == "both"
        X = H.C.forward(X)
    end
    Xa, Xb = tensor_split(X)
    permute == "lower" && (Xb = H.C.forward(Xb))

    # Determine whether to continue recursion
    recursive = false
    if typeof(X) <: AbstractArray{Float32, 4} && size(X, 3) > 4
        recursive = true
    elseif typeof(X) <: AbstractArray{Float32, 5} && size(X, 4) > 4
        recursive = true
    end

    # HINT coupling
    if recursive
        # Call function recursively
        Ya, logdet1 = forward(Xa, H; scale=scale+1, permute="none")
        Y_temp, logdet2 = forward(Xb, H; scale=scale+1, permute="none")
        if logdet
            Yb, logdet3 = H.CL[scale].forward(Xa, Y_temp)[[2,3]]
        else
            Yb = H.CL[scale].forward(Xa, Y_temp)[2]
            logdet3 = 0f0
        end
        logdet_full = logdet1 + logdet2 + logdet3
    else
        # Finest layer
        Ya = copy(Xa)
        if logdet
            Yb, logdet_full = H.CL[scale].forward(Xa, Xb)[[2,3]]
        else
            Yb = H.CL[scale].forward(Xa, Xb)[2]
            logdet_full = 0f0
        end
    end

    Y = tensor_cat(Ya, Yb)
    permute == "both" && (Y = H.C.inverse(Y))
    if scale == 1
        logdet ? (return Y, logdet_full) : (return Y)
    else
        return Y, logdet_full
    end
end

# Input is tensor Y
function inverse(Y, H::CouplingLayerHINT; scale=1, permute=nothing, logdet=nothing)
    isnothing(logdet) ? logdet = (H.logdet && H.is_reversed) : logdet = logdet
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    permute == "both" && (Y = H.C.forward(Y))
    Ya, Yb = tensor_split(Y)

    # Check for recursion
    recursive = false
    if typeof(Y) <: AbstractArray{Float32, 4} && size(Y, 3) > 4
        recursive = true
    elseif typeof(Y) <: AbstractArray{Float32, 5} && size(Y, 4) > 4
        recursive = true
    end

    # Coupling layer
    if recursive
        Xa, logdet1 = inverse(Ya, H; scale=scale+1, permute="none")
        if logdet
            Y_temp, logdet2 = H.CL[scale].inverse(Xa, Yb; logdet=true)[[2,3]]
        else
            Y_temp = H.CL[scale].inverse(Xa, Yb)[2]
            logdet2 = 0f0
        end
        Xb, logdet3 = inverse(Y_temp, H; scale=scale+1, permute="none")
        logdet_full = logdet1 + logdet2 + logdet3
    else
        Xa = copy(Ya)
        if logdet
            Xb, logdet_full = H.CL[scale].inverse(Ya, Yb; logdet=true)[[2,3]]
        else
            Xb = H.CL[scale].inverse(Ya, Yb)[2]
            logdet_full = 0f0
        end
    end

    # Initial permutation
    permute == "lower" && (Xb = H.C.inverse(Xb))
    X = tensor_cat(Xa, Xb)
    if permute == "full" || permute == "both"
        X = H.C.inverse(X)
    end

    if scale == 1
        logdet ? (return X, logdet_full) : (return X)
    else
        return X, logdet_full
    end
end

# Input are two tensors ΔY, Y
function backward(ΔY, Y, H::CouplingLayerHINT; scale=1, permute=nothing)
    isnothing(permute) ? permute = H.permute : permute = permute
    permute == "both" && ((ΔY, Y) = H.C.forward((ΔY, Y)))
    Ya, Yb = tensor_split(Y)
    ΔYa, ΔYb = tensor_split(ΔY)
    recursive = false
    if typeof(Y) <: AbstractArray{Float32, 4} && size(Y, 3) > 4
        recursive = true
    elseif typeof(Y) <: AbstractArray{Float32, 5} && size(Y, 4) > 4
        recursive = true
    end
    if recursive
        ΔXa, Xa = backward(ΔYa, Ya, H; scale=scale+1, permute="none")
        ΔXa_temp, ΔXb_temp, X_temp = H.CL[scale].backward(ΔXa.*0f0, ΔYb, Xa, Yb)[[1,2,4]]
        ΔXb, Xb = backward(ΔXb_temp, X_temp, H; scale=scale+1, permute="none")
        ΔXa += ΔXa_temp
    else
        Xa = copy(Ya)
        ΔXa = copy(ΔYa)
        ΔXa_, ΔXb, Xb = H.CL[scale].backward(ΔYa.*0f0, ΔYb, Ya, Yb)[[1,2,4]]
        ΔXa += ΔXa_
    end
    permute == "lower" && ((ΔXb, Xb) = H.C.inverse((ΔXb, Xb)))
    ΔX = tensor_cat(ΔXa, ΔXb)
    X = tensor_cat(Xa, Xb)
    if permute == "full" || permute == "both"
        (ΔX, X) = H.C.inverse((ΔX, X))
    end
    return ΔX, X
end

# Input are two tensors ΔX, X
function backward_inv(ΔX, X, H::CouplingLayerHINT; scale=1, permute=nothing)
    isnothing(permute) ? permute = H.permute : permute = permute

    # Permutation
    if permute == "full" || permute == "both"
        (ΔX, X) = H.C.forward((ΔX, X))
    end
    ΔXa, ΔXb = tensor_split(ΔX)
    Xa, Xb = tensor_split(X)
    permute == "lower" && ((ΔXb, Xb) = H.C.forward((ΔXb, Xb)))

    # Check whether to continue recursion
    recursive = false
    if typeof(X) <: AbstractArray{Float32, 4} && size(X, 3) > 4
        recursive = true
    elseif typeof(X)<: AbstractArray{Float32, 5} && size(X, 4) > 4
        recursive = true
    end

    # Coupling layer backprop
    if recursive
        ΔY_temp, Y_temp = backward_inv(ΔXb, Xb, H; scale=scale+1, permute="none")
        ΔYa_temp, ΔYb, Yb = backward_inv(0f0.*ΔXa, ΔY_temp, Xa, Y_temp, H.CL[scale])[[1,2,4]]
        ΔYa, Ya = backward_inv(ΔXa+ΔYa_temp, Xa, H; scale=scale+1, permute="none")
    else
        ΔYa = copy(ΔXa)
        Ya = copy(Xa)
        ΔYa_temp, ΔYb, Yb = backward_inv(0f0.*ΔYa, ΔXb, Xa, Xb, H.CL[scale])[[1,2,4]]
        ΔYa += ΔYa_temp
    end
    ΔY = tensor_cat(ΔYa, ΔYb)
    Y = tensor_cat(Ya, Yb)
    permute == "both" && ((ΔY, Y) = H.C.inverse((ΔY, Y)))
    return ΔY, Y
end

# Clear gradients
function clear_grad!(H::CouplingLayerHINT)
    for j=1:length(H.CL)
        clear_grad!(H.CL[j])
    end
    ~isnothing(H.C) && clear_grad!(H.C)
end

# Get parameters
function get_params(H::CouplingLayerHINT)
    nlayers = length(H.CL)
    p = get_params(H.CL[1])
    if nlayers > 1
        for j=2:nlayers
            p = cat(p, get_params(H.CL[j]); dims=1)
        end
    end
    ~isnothing(H.C) && (p = cat(p, get_params(H.C); dims=1))
    return p
end

# Put parameters
function put_params!(H::CouplingLayerHINT, Params::Array{Any,1})
    nlayers = length(H.CL)
    put_params!(H.CL[1], Params[1:5])
    if nlayers > 1
        for j=2:nlayers
            put_params!(H.CL[j], Params[5*(j-1)+1:5*j])
        end
    end
    ~isnothing(H.C) && put_params!(H.C, Params[5*nlayers+1:end])
end

# Set is_reversed flag in full network tree
function tag_as_reversed!(H::CouplingLayerHINT, tag::Bool)
    H.is_reversed = tag
    nlayers = length(H.CL)
    for j=1:nlayers
        H.CL[j].is_reversed = tag
    end
    return H
end


text:
##### Beginning of file

@info("Importing the MirrorUpdater module...")

import MirrorUpdater

import TimeZones

@info("Reading config files...")

include(joinpath("config","preferences","bitbucket.jl",))
include(joinpath("config","preferences","enabled-providers.jl",))
include(joinpath("config","preferences","gitlab.jl",))
include(joinpath("config","preferences","github.jl",))
include(joinpath("config","preferences","time-zone.jl",))

include(joinpath("config","repositories","additional-repos.jl",))
include(joinpath("config","repositories",
    "do-not-push-to-these-destinations.jl",))
include(joinpath("config","repositories",
    "do-not-try-url-list.jl",))
include(joinpath("config","repositories","registries.jl",))
include(joinpath("config","repositories",
    "try-but-allow-failures-url-list.jl",))

git_hosting_providers = Any[]

if GITHUB_ENABLED
    const github_provider =
        MirrorUpdater.Hosts.GitHubHost.new_github_session(
            ;
            github_organization = GITHUB_ORGANIZATION,
            github_bot_username = GITHUB_BOT_USERNAME,
            github_bot_personal_access_token = GITHUB_BOT_PERSONAL_ACCESS_TOKEN,
            )
    push!(git_hosting_providers, github_provider)
end

if GITLAB_ENABLED
    const gitlab_provider =
        MirrorUpdater.Hosts.GitLabHost.new_gitlab_session(
            ;
            gitlab_group = GITLAB_GROUP,
            gitlab_bot_username = GITLAB_BOT_USERNAME,
            gitlab_bot_personal_access_token = GITLAB_BOT_PERSONAL_ACCESS_TOKEN,
            )
    push!(git_hosting_providers, gitlab_provider)
end

if BITBUCKET_ENABLED
    const bitbucket_provider =
        MirrorUpdater.Hosts.BitbucketHost.new_bitbucket_session(
            ;
            bitbucket_team = BITBUCKET_TEAM,
            bitbucket_bot_username = BITBUCKET_BOT_USERNAME,
            bitbucket_bot_app_password = BITBUCKET_BOT_APP_PASSWORD,
            )
    push!(git_hosting_providers, bitbucket_provider)
end

MirrorUpdater.CommandLine.run_mirror_updater_command_line!!(
    ;
    arguments = ARGS,
    git_hosting_providers = git_hosting_providers,
    registry_list = REGISTRY_LIST,
    additional_repos = ADDITIONAL_REPOS,
    do_not_try_url_list = DO_NOT_TRY_URL_LIST,
    do_not_push_to_these_destinations = DO_NOT_PUSH_TO_THESE_DESTINATIONS,
    try_but_allow_failures_url_list = TRY_BUT_ALLOW_FAILURES_URL_LIST,
    time_zone = TIME_ZONE,
    )

##### End of file


text:
### A Pluto.jl notebook ###
# v0.16.1

using Markdown
using InteractiveUtils

# ╔═╡ 35dd2eca-2446-11ec-1de6-0feeace9d4f7
md"""

#### Flash separation of a binary mixture of component A and component B

A liquid feed stream that contains 60% component A (1) and 40% component B (2) is separated in a Flash Drum operating at P = 1.21 MPa and T = 150$^{\circ}$C. The saturation pressures of pure component A and B can be modeled using the [Antoine equation](https://en.wikipedia.org/wiki/Antoine_equation):

$$\begin{eqnarray}
\ln\left(P_{1}^{sat}\right) &=& 15 -\frac{3010}{T+250}\\
\ln\left(P_{2}^{sat}\right) &=& 14 -\frac{2700}{T+205}
\end{eqnarray}$$

where $T$ denotes the temperature in units of celcius $^{\circ}$C and $P_{\star}^{sat}$ denotes the saturation pressure in units of kPa for component $\star$.

##### Assumptions

1. The liquid and vapor phases of all streams are ideal
1. The liquid and vapor phases inside the Flash Drum are ideal (Raoult's law)
1. The liquid and vapor phases inside the Flash Drum reach vapor-liquid equilibrium instantly

##### Solve:

 * Calculate the fraction of the input stream that exits the drum as liquid
 * Calculate the fraction of the input stream that exits the drum as vapor
 * Calculate the composition of the liquid ($x_{1}$ and $x_{2}$) and vapor ($y_{1}$ and $y_{2}$) in the exit streams

"""

# ╔═╡ 92ff1f65-d717-47a8-9cd3-af71c2855210
PlutoUI.LocalResource("./figs/Fig-Flash.png")

# ╔═╡ 98e43733-8711-4aa4-9081-697dd8865592
begin

	# write down some constants -
	T = 150 		# units: C
	P = 1210 		# units: kPa
	z₁ = 0.60 		# units: dimensionless
	z₂ = (1-z₁)		# units: dimensionless

	# return -
	nothing
end

# ╔═╡ 792e101b-e736-46d2-bbf3-a1197a0cc654
md"""

##### Governing equations

Let's start by writing down the governing equations. Starting from the total mol balance:

$$\dot{F} = \dot{L}+\dot{V}$$

we can divide both sides by $\dot{F}$ to give a relationship which constrains the fraction(s) of liquid and vapor in the outlet streams:

$$\hat{L}+\hat{V} = 1$$

Further, we know that the species mol balance around species $i$ is given by:

$$\dot{F}z_{i} = \dot{L}x_{i}+\dot{V}y_{i}\qquad{i=1,2}$$

However, because we assumed ideal vapor and liquid, we know that the liquid and vapor phase equlibrium compositions are related by the fugacity matching condition: $y_{i}P = x_{i}P_{i}^{sat}$. Lastly, we can sum the matching conditions to arrive at an expression which governs the pressure that is only in $x_{i}$:

$$P = x_{1}P_{1}^{sat} + x_{2}P_{2}^{sat}$$

which can be re-written (because we are in a binary mixture) as:

$$P = x_{1}P_{1}^{sat} + \left(1-x_{1}\right)P_{2}^{sat}$$

"""

# ╔═╡ 9e171e7c-953d-4a75-b6f6-cbae930dd1ab
md"""

#### Solve for $x_{1}$ and $x_{2}$:

We can re-arrange the pressure expression to solve for $x_{1}$:

$$x_{1} = \frac{P - P_{2}^{sat}}{P_{1}^{sat} -  P_{2}^{sat}}$$

"""

# ╔═╡ 8e13e386-6fbc-4bdf-8eed-de5bd35247da
# compute P1sat (units: kPa)
P1_sat = exp(15 - (3010)/(T+250))

# ╔═╡ 45a8e5fc-0929-496b-b01a-0e6f27f8b457
# compute P2sat (units: kPa)
P2_sat = exp(14 - (2700)/(T+205))

# ╔═╡ 2bd63a06-0eb6-4e4e-b2b0-78b0ac87e341
# solve for x₁:
x₁ = (P - P2_sat)/(P1_sat - P2_sat)

# ╔═╡ 6b3bcfa8-7487-4f84-abbe-185939d91a80
x₂ = 1 - x₁

# ╔═╡ 50a0b336-e14c-487a-9ba5-25a0024d9766
md"""
#### Solve for $y_{1}$ and $y_{2}$:

Once we have the liquid phase composition, we can use the fugacity matching condition to relate the liquid phase and the vapr phase becuase we have asssumed they are in equilibrium i.e.,

$$y_{i} = x_{i}\left(\frac{P_{i}^{sat}}{P}\right)\qquad{i=1,2}$$

"""

# ╔═╡ 7608d2c9-eb24-4db1-9fb3-452d52770e17
# solve for y₁:
y₁ = x₁ * (P1_sat/P)

# ╔═╡ 76f6cf32-8f87-4df0-b1db-a51ee102f709
# solve for y₂:
y₂ = x₂* (P2_sat/P)

# ╔═╡ ec5378ba-ca2c-4d31-b176-454b2cd99302
md"""

#### Solve for $\hat{L}$ and $\hat{V}$:

Now that we have the composition values in the exit streams, we can use the total mol and species mol balances to estimate the fraction
of the input feed that exited the unit as liquid $\hat{L}$ and vapur $\hat{V}$. If we divide both sides of the species mol balance by $\dot{F}$ we get the relationship:

$$z_{i} = \hat{L}x_{i}+\hat{V}y_{i}\qquad{i=1,2}$$

which we can use (along with the total mol balances) to solve for $\hat{L}$ and $\hat{V}$ i.e.,

$$\begin{eqnarray}
\hat{L}x_{1}+\hat{V}y_{1} &=& z_{1} \\
\hat{L}+\hat{V} &=& 1
\end{eqnarray}$$

This is a system of 2-equations and 2-unknowns ($\hat{L}$ and $\hat{V}$) that can be solved by substitution or by linear algebra.

###### Approach 1: Old school subsitution/elimination approach:

Let's rearrange the total mol expression to get $\hat{L}$ in terms of $\hat{V}$:

$$\hat{L} = 1 - \hat{V}$$

which can be subsituted into the species mol balance (let's use component 1):

$$z_{1} = \left(1 - \hat{V}\right)x_{1}+\hat{V}y_{1}$$

and solved for $\hat{V}$:

$$\hat{V} = \frac{z_{1} - x_{1}}{y_{1} - x_{1}}$$

Once we have $\hat{V}$, we can compute $\hat{L}$ from the total mol balance.
"""

# ╔═╡ bb56b509-5279-40f5-b898-e8c1b29dfb5e
# compute Vhat -
Vhat = (z₁ - x₁)/(y₁ - x₁)

# ╔═╡ 3f13e0a8-b58b-4ac7-899b-412154acf0cb
# what is Lhat -
Lhat = 1 - Vhat

# ╔═╡ 745177ba-a88d-496d-b0f6-2e9aca506866
md"""
###### Approach 2: Fancy linear algebra approach:

We can rewrite the 2$\times$2 system of equations into matrix-vector form and then compute the matrix inverse. In particlar,
the species and total mass balances for component 1 can be written as:

$$\begin{pmatrix}
x_{1} & y_{1} \\
1 & 1
\end{pmatrix}
\begin{pmatrix}
\hat{L} \\
\hat{V}
\end{pmatrix} =
\begin{pmatrix}
z_{1} \\
1
\end{pmatrix}$$

which is the form $\mathbf{A}\mathbf{x} = \mathbf{b}$, where the solution takes the form: $\mathbf{x} = \mathbf{A}^{-1}\mathbf{b}$.

"""

# ╔═╡ 4c0a79c3-87b4-4452-ba9e-0de5bd324ac5
# setup A -
A = [x₁ y₁ ; 1 1]

# ╔═╡ a835ed0a-a2d0-4a5b-a8a0-dffd7250e34b
# setup b -
bV = [z₁ ; 1]

# ╔═╡ 2bf7b85f-3118-4fa9-842d-f2c9f3d16dce
# compute the unknown (fractions of L and V)
x = inv(A)*bV

# ╔═╡ b091bbb6-13ce-4b1f-b36b-9964783767d4
html"""<style>
main {
    max-width: 1200px;
  	width: 80%;
	margin: auto;
	font-family: "Roboto, monospace";
}

a {
	color: blue;
	text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>
"""

# ╔═╡ 64af17f4-c632-425c-a101-9dfa5e8cbd84
# julia setup -
using PlutoUI

# ╔═╡ 79f99ade-b514-4768-9bc9-bcfa281ce02f
png_joinpathsplit__FILE__1assetsimagepng = let
    import PlutoUI
    PlutoUI.LocalResource(joinpath(split(@__FILE__, '#')[1] * ".assets", "image.png"))
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.14"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[PlutoUI]]
deps = ["Base64", "Dates", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "d1fb76655a95bf6ea4348d7197b22e889a4375f4"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.14"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═64af17f4-c632-425c-a101-9dfa5e8cbd84
# ╟─35dd2eca-2446-11ec-1de6-0feeace9d4f7
# ╟─92ff1f65-d717-47a8-9cd3-af71c2855210
# ╠═79f99ade-b514-4768-9bc9-bcfa281ce02f
# ╠═98e43733-8711-4aa4-9081-697dd8865592
# ╟─792e101b-e736-46d2-bbf3-a1197a0cc654
# ╟─9e171e7c-953d-4a75-b6f6-cbae930dd1ab
# ╠═8e13e386-6fbc-4bdf-8eed-de5bd35247da
# ╠═45a8e5fc-0929-496b-b01a-0e6f27f8b457
# ╠═2bd63a06-0eb6-4e4e-b2b0-78b0ac87e341
# ╠═6b3bcfa8-7487-4f84-abbe-185939d91a80
# ╟─50a0b336-e14c-487a-9ba5-25a0024d9766
# ╠═7608d2c9-eb24-4db1-9fb3-452d52770e17
# ╠═76f6cf32-8f87-4df0-b1db-a51ee102f709
# ╟─ec5378ba-ca2c-4d31-b176-454b2cd99302
# ╠═bb56b509-5279-40f5-b898-e8c1b29dfb5e
# ╠═3f13e0a8-b58b-4ac7-899b-412154acf0cb
# ╟─745177ba-a88d-496d-b0f6-2e9aca506866
# ╠═4c0a79c3-87b4-4452-ba9e-0de5bd324ac5
# ╠═a835ed0a-a2d0-4a5b-a8a0-dffd7250e34b
# ╠═2bf7b85f-3118-4fa9-842d-f2c9f3d16dce
# ╟─b091bbb6-13ce-4b1f-b36b-9964783767d4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Tidal averages

struct TidalAverage{T<:Quantity} <: Quantity{Float64}
    ts
    q
end

"""
Take a mask and return the indices of the rising and falling edges
"""
function edgeindices(mask)
    dm = [0;diff(mask)]
    irise = find(dm.==1)
    ifall = find(dm.==-1)
    # Drop the first ifall if it is before the first irise
    if ifall[1]<irise[1]
        ifall = ifall[2:end]
    end
    # Drop the last irise if it is after the last ifall
    if irise[end]>ifall[end]
        irise = irise[1:end-1]
    end
    irise,ifall
end

edgeindices(mask::Mask) = edgeindices(quantity(mask))

"""
Take a mask and return the times of the rising and falling edges
"""
function edgetimes(mask,ts)
    irise,ifall = edgeindices(mask)
    trise = ts[irise]
    tfall = ts[ifall]
    trise,tfall
end

"""
Integrate f sampled at times ts with a trapezoid rule
"""
function integrate(f,ts)
    dt = Dates.value.(Dates.Second.(diff(ts)))
    f1 = f[1:end-1]
    f2 = f[2:end]
    sum(dt.*(f1+f2)/2)
end

integrate(q::Quantity) = integrate(quantity(q),times(q))

"""
Average f sampled at times ts with a trapezoid rule
"""
function average(f,ts)
    fi = integrate(f,ts)
    fi./Dates.value.(Dates.Second.(ts[end]-ts[1]))
end

average(q::Quantity) = average(quantity(q),times(q))

function segment(q,mask)
    ts,qs = unzip(q)
    irise,ifall = edgeindices(mask)
    st = [ts[irise[i]:ifall[i]] for i in eachindex(irise)]
    sq = [qs[irise[i]:ifall[i]] for i in eachindex(irise)]
    st,sq
end

"""
Integrate a flux over tides given in the mask.

Returns midpoint times for the each tide and the accumulated flux over that tide
"""
function tidalbalance(f,mask)
    ts,fs = unzip(f)
    irise,ifall = edgeindices(mask)
    trise,tfall = edgetimes(mask,ts)
    tmid = trise + (tfall-trise)/2
    avg = [integrate(fs[irise[i]:ifall[i]],
                     ts[irise[i]:ifall[i]]) for i in eachindex(irise)]
    TidalAverage(tmid,avg)
end

function tidalaverage{T<:Quantity}(f::T,mask)
    ts,fs = unzip(f)
    irise,ifall = edgeindices(mask)
    trise,tfall = edgetimes(mask,ts)
    tmid = trise + (tfall-trise)/2
    avg = [average(fs[irise[i]:ifall[i]],
                   ts[irise[i]:ifall[i]]) for i in eachindex(irise)]
    TidalAverage{T}(tmid,avg)
end

"""
Compute the flood-ebb differential of f

q determines flood vs. ebb conditions
"""
function flood_ebb{T<:Quantity}(f::T,q,mask)
    ts,fs = unzip(f)
    f2 = fs[:]
    f2[quantity(q).<0] .*= -1
    tidalaverage(T(ts,f2),mask)
end


text:
using JuMP, Cbc, DelimitedFiles, CSV, DataFrames
function f()

    dat = readdlm("donor-pool2.csv", ',', '\n', comments=true)

    fr = dat[:,1]
    to = dat[:,2]
    w = dat[:,3]
    N = dat[:,4]
    N = N[1:3,]

    V = union(fr,to)       # set of all nodes
    E = collect(zip(fr,to)) # set of all edges

    W = Dict( (i,j) => k for (i,j,k) in zip(fr,to,w) )  # weights on edges

    m = Model()
    set_optimizer(m, Cbc.Optimizer)

    @variable(m, y[(i,j) in E], Bin)

    @variable(m, f_in[i in V])

    @variable(m, f_out[i in V])

    @objective(m, Max, sum( W[(i,j)]*y[(i,j)] for (i,j) in E ))

    # 1
    @constraint(m, flow_out[k in V], sum(y[(i,j)] for (i,j) in E if i==k) == f_out[k])
    @constraint(m, flow_in[k in V], sum(y[(i,j)] for (i,j) in E if j==k) == f_in[k])

    # 2
    @constraint(m, flow_ineq[k in setdiff(V,N)], f_out[k] <= f_in[k])
    @constraint(m, flow_ineq2[k in setdiff(V,N)], f_in[k] <= 1)

    # 3
    @constraint(m, flow_out_ndd[k in N], f_out[k] <= 1)

    Cycle = []
    # Note that keeping the list of subtours in Subtour is not really
    # necessary. I just want to see the subtours that are identified so
    # that I know the algorithm is working.

    # loop to iteratively identify subtours and add the violated
    # constraints to the model.
    while true # all cycles aren't under 5
        cycles_over_five = 0
        JuMP.optimize!(m)
        set = []
        for (i,j) in E
            if value.(y[(i,j)]) > 0
                set = union(set,i,j)
            end
        end
        print(length(set))
        #print(value.(y))
        # get the pairs that are a part of the solution
        searched = []
        U = []
        ExtendU = set[rand(1:end)] # choose a node at random
        push!(searched,ExtendU)
        #print("U ",ExtendU,"\n")

        # look for the next edge in our solution to
        # map all pairs to one another. capture whether
        # the pairs are part of a cycle in which we will
        # record that cycle, or if they are a part of a
        # chain in which we will ignore.
        while length(setdiff(set,searched))!=0 # while we haven't searched all edges
            push!(U,ExtendU)
            ExtendU = nothing
            for j in union(U[1],setdiff(set,searched))
                if (in((U[end],j),E) && JuMP.value(y[(U[end],j)]) > 0)
                    ExtendU = j
                    print(length(searched),"\n")
                    push!(searched,ExtendU)
                    break
                end
            end
            # if we've found the end of a cycle
            if U[1]==ExtendU
                print("Cycle Found","\n")
                # we've found a cycle, lets add it to our cycle list and create
                # a constraint if necessary
                #print("hi")
                if length(U) > 5
                    push!(Cycle,U)
                    cycles_over_five = cycles_over_five + 1
                end
                U = []
                rn = setdiff(set,searched)
                ExtendU = rn[rand(1:end)]
                push!(searched,ExtendU)
            end
            # if we've found the end of a chain
            if isnothing(ExtendU)
                print("Chain Found","\n")
                U = []
                rn = setdiff(set,searched)
                ExtendU = rn[rand(1:end)] # don't pick the same one twice..
                push!(searched,ExtendU)
                # we've found a chain!
            end
        end
        for z in 1:length(Cycle)
            print("constraint added")
            @constraint(m, sum(y[(i,j)] for i in Cycle[z], j in Cycle[z] if in((i,j),E)) <= length(Cycle[z]) - 1)
        end
        if cycles_over_five==0
            JuMP.optimize!(m)
            for (i,j) in E
                if value.(y[(i,j)]) > 0
                    println(i, " to ", j)
                end
            end
            break
        end
    end
    return(JuMP.value.(y),E)
end

function find_chains(y,E,soln,ndd)
    Chain = []
    for i in ndd
        stop=false
        searched = []
        U = []
        ExtendU = i # start at ndd to find chain
        push!(searched,ExtendU)

        while stop==false
        # convert solution to a usable format for health professionals
            push!(U,ExtendU)
            ExtendU = nothing
            for j in union(U[1],setdiff(soln,searched))
                if (in((U[end],j),E) && y[(U[end],j)] > 0)
                    ExtendU = j
                    print(length(searched),"\n")
                    push!(searched,ExtendU)
                    break
                end
            end
            if isnothing(ExtendU)
                print("Chain Found","\n")
                push!(Chain,U)
                U = []
                rn = setdiff(soln,searched)
                ExtendU = rn[rand(1:end)] # don't pick the same one twice..
                push!(searched,ExtendU)
                stop=true
                # we've found a chain!
            end
        end
    end
    return(Chain)
end
# find cycles
function find_cycles(foundchains,y,E,soln)
    fcarray = []
    Cycle = []
    for i in 1:length(foundchains)
        union!(fcarray,foundchains[i])
    end
    tosearch = setdiff(soln,fcarray) # search only pairs a part of cycles
    searched = []
    U = []
    ExtendU = tosearch[rand(1:end)] # start at ndd to find chain
    push!(searched,ExtendU)
    while length(setdiff(tosearch,searched))!=0 # while we haven't searched all edges
        push!(U,ExtendU)
        ExtendU = nothing
        for j in tosearch
            if (in((U[end],j),E) && y[(U[end],j)] > 0)
                ExtendU = j
                print(length(searched),"\n")
                push!(searched,ExtendU)
                break
            end
        end
        if U[1]==ExtendU
            # found a cycle
            push!(U,ExtendU)
            push!(Cycle,U)
            U = []
            rn = setdiff(tosearch,searched)
            ExtendU = rn[rand(1:end)]
            push!(searched,ExtendU)
        end
    end
    push!(U,ExtendU)
    push!(Cycle,U)
    return(Cycle)
end
y, E = f()
soln = []
for (i,j) in E
    if y[(i,j)] > 0
        union!(soln,i,j)
    end
end
ndd = [13,978,999] # enter the NDDs here
foundchains = find_chains(y,E,soln,ndd)
Cycle = find_cycles(foundchains,y,E,soln)

CSV.write("cycles.csv",DataFrame(Cycle),delim=',')
CSV.write("chains.csv",DataFrame(foundchains),delim=',')

max_stars_repo_name: trouze/UMNCsci

max_issues_repo_name: trouze/UMNCsci

max_forks_repo_name: trouze/UMNCsci


text:
export SparseGP, GPmodel, GPODE
export samplegp, gpmeanvar

# https://github.com/SciML/DiffEqFlux.jl/blob/c59971fd4d3ee84aff39f88b7073d7e8cf51c34c/src/neural_de.jl#L38

abstract type TrajType end
struct SampleTraj <: TrajType end
struct MeanTraj <: TrajType end

abstract type SparseGPMethod end
struct FITC <: SparseGPMethod end
struct VLB <: SparseGPMethod end
struct Vanilla <: SparseGPMethod end

###
# Struct that contains everything needed for prediction
###
struct SparseGP{K, T<:Real, N, A<: NTuple{N, Array{<:Array{<:Real,1},1}}, M<:SparseGPMethod}
    kernel::K
    σ_n::T
    inP::A # inducing points
    mean::Function
    trafo::Function
    method::M
    # type? FITC, SOR, PITC
end

#ToDo: Update σ_n in this object

function zeromean(n)
    return x -> fill(0, n)
end
identitytrafo(x) = x

function SparseGP(kernel, Z, U; σ_n = 1e-6, mean = zeromean(length(U[1])), trafo = identitytrafo, method = Vanilla())
    indP = (trafo.(Z), U)
    N = length(indP)
    SparseGP{typeof(kernel), typeof(σ_n), N, typeof(indP), typeof(method)}(kernel, σ_n, indP, mean, trafo, method)
end
function SparseGP(kernel, Z, X, Y; σ_n = 1e-6, mean = zeromean(length(Y[1])), trafo = identitytrafo, method = FITC())
    indP = (trafo.(Z), trafo.(X), Y)
    N = length(indP)
    SparseGP{typeof(kernel), typeof(σ_n), N, typeof(indP), typeof(method)}(kernel, σ_n, indP, mean, trafo, method)
end

function (sgp::SparseGP)(x::T) where T <: Real
    Z = sgp.inP[1]
    ker = sgp.kernel
    Kx = kernelmatrix(ker, [[x]], Z)
end

function (sgp::SparseGP)(x::Array{T,1}) where T <: Real
    Z = sgp.inP[1]
    ker = sgp.kernel
    Kx = kernelmatrix(ker, [x], Z)
end

###
# GP model, that contains the sparse GP object with all necessary data, and
struct GPmodel{SGP <: SparseGP, T<:Real}
    sgp::SGP
    KinvU::Array{T,2}
    Σ::Array{T,2}
end

# maybe with data
function GPmodel(sgp::SparseGP)
    KiU, Σ = computeKinvU(sgp)
    GPmodel(sgp, KiU, Σ)
end

function (gpm::GPmodel)(x, traj::MeanTraj = MeanTraj())
    Kx = gpm.sgp(gpm.sgp.trafo(x))
    μ = gpm.sgp.mean
    return (μ(x) .+ (Kx * gpm.KinvU)[:])
end

function (gpm::GPmodel)(x, traj::SampleTraj)
    Kx = gpm.sgp(gpm.sgp.trafo(x))
    μ = gpm.sgp.mean
    m =  (μ(x) .+ (Kx * gpm.KinvU))[:]

    Σ = gpm.Σ

    std = sqrt.(diag(Kx * (Σ \ Kx')))
    # var = diag(Kx*Kuu*Kx')
    s = randn(length(x)) .* std
    return (m .+ s)
end

# this can almost certainly be removed. might be interesting for comparison purposes
function (gpm::GPmodel)(xv::MS) where MS  <: Array{<:Measurement{<:Real}, 1}
    x = getfield.(xv, :val)
    Sigma = diagm(getfield.(xv, err))

    Kx = gpm.sgp(gpm.sgp.trafo(x))
    μ = gpm.sgp.mean
    m =  (μ(x) .+ (Kx * gpm.KinvU))[:]

    Σ = gpm.Σ

    var = sqrt.(diag(Kx * (Σ \ Kx')))
    # var = diag(Kx*Kuu*Kx')

    return (m .± var)
end

####
# functions to facilitate efficient computation, as per Q-C&R
function computeKinvU(sgp::SparseGP)
    return _computeKinvU(sgp, sgp.inP, sgp.method)
end

function _computeKinvU(sgp::SparseGP, indP::NTuple{2, Array{<:Array{<:Real,1},1}}, method::M) where M<:SparseGPMethod
    Z = indP[1]
    U = indP[2]
    σ_n = sgp.σ_n
    μ = sgp.mean
    U = U .- μ.(Z)
#     vU = reduce(vcat, U)
    # ToDo: might have to make output(?) dimensions more explicit
    vU = reshape(reduce(vcat, U), :, length(Z)*length(Z[1]))
    vU = permutedims(vU)
    ker = sgp.kernel
    K = kernelmatrix(ker, Z) + σ_n * I
    KinvU = K \ vU
    return (KinvU, K)
end

function _computeKinvU(sgp::SparseGP, indP::NTuple{3, Array{<:Array{<:Real,1},1}}, method::M) where M <: SparseGPMethod
    Z = indP[1]
    X = indP[2]
    Y = indP[3]
    μ = sgp.mean
    Y = Y .- μ.(X)
    ker = sgp.kernel

    Kff = kernelmatrix(ker, X)
    Kfu = kernelmatrix(ker, X, Z)
    Kuu = kernelmatrix(ker, Z)

    Qff = Kfu * ( Kuu \ Kfu' )

    noise = sgp.σ_n
    Λ = _computelambda(Kff, Qff, noise, method)

    Σ = _computesigma(Kuu, Kfu, Λ)

    vY = reshape(reduce(vcat, Y), :, length(X)*length(X[1]))
    vY = permutedims(vY)
    KinvU = Σ \ (Kfu' * (Λ \ vY))
    return (KinvU, Σ)
end

function _computesigma(Kuu, Kfu, Λ)
    return Kuu + Kfu' * (Λ \ Kfu)
end

function _computelambda(Kff, Qff, noise, method::VLB)
    return noise * I
end

function _computelambda(Kff, Qff, noise, method::FITC)
    return Diagonal(diag( Kff - Qff) .+ noise)
end

#####
# complete GPODE construct
#####
basic_tgrad(u,p,t) = zero(u)
struct GPODE{M<:GPmodel,T,A,K,TT<: TrajType} #<: NeuralDELayer
    model::M
    # p::P, parameters, maybe one day
    tspan::T
    args::A
    kwargs::K
    trajtype::TT

    function GPODE(model::GPM,tspan,args...; sample = false, kwargs...) where GPM <: GPmodel
        tt = sample ? SampleTraj() : MeanTraj()
        new{typeof(model),typeof(tspan),typeof(args),typeof(kwargs), typeof(tt)}(
            model,tspan,args,kwargs,tt)
    end
end

function GPODE(sgp::SGP,tspan,args...;sample::Bool=false, kwargs...) where SGP <: SparseGP
    gpm = GPmodel(sgp)
    return GPODE(gpm,tspan,args..., sample = sample, kwargs...)
end

function (gp::GPODE)(x0)
    dudt_(u,p,t) = gp.model(u, gp.trajtype)
    ff = ODEFunction{false}(dudt_,tgrad=basic_tgrad)
    prob = ODEProblem{false}(ff,x0,getfield(gp,:tspan))
    solve(prob,gp.args...;gp.kwargs...)
end

function samplegp(x, gp, n, args = gp.args; kwargs...)
    #     dudt_(u,p,t) = gp.model(u, GaussianProcessODEs.SampleTraj())
    #     ff = ODEFunction{false}(dudt_,tgrad=GaussianProcessODEs.basic_tgrad)
    #     prob = ODEProblem{false}(ff,x,getfield(gp,:tspan))
    #     ensemble = EnsembleProblem(prob)
    #     solve(ensemble, gp.args..., EnsembleThreads(); trajectories = n, gp.kwargs...)

        sample_gpODE = GPODE(gp.model, gp.tspan,
            args...; sample = true, kwargs...);
        t = @elapsed sols = ThreadsX.collect(sample_gpODE(x) for i in 1:n)
        EnsembleSolution(sols, t, true)
end

using DifferentialEquations.EnsembleAnalysis
function gpmeanvar(samplesols, ts)
    meanvar = timepoint_meanvar(samplesols,ts)
    mtrajs = meanvar[1]
    stdtrajs = sqrt.(meanvar[2])
    return mtrajs, stdtrajs
end


text:
"""
    TypicalSection <: UnsteadyModel

Typical section structural model with state variables ``h, \\theta, \\dot{h},
\\dot{\\theta}``, inputs ``\\mathcal{L}, \\mathcal{M}``, and parameters ``k_h,
k_\\theta, m, S_\\theta, I_\\theta``
"""
struct TypicalSection <: UnsteadyModel end

"""
    TypicalSection()

Initialize an object of type [`TypicalSection`](@ref)
"""
TypicalSection()

# --- Traits --- #

number_of_states(::Type{TypicalSection}) = 4
number_of_inputs(::Type{TypicalSection}) = 2
number_of_parameters(::Type{TypicalSection}) = 5

inplaceness(::Type{TypicalSection}) = OutOfPlace()

rate_jacobian_type(::Type{TypicalSection}) = Constant()
state_jacobian_type(::Type{TypicalSection}) = Constant()
input_jacobian_type(::Type{TypicalSection}) = Invariant()
parameter_jacobian_type(::Type{TypicalSection}) = Linear()
time_gradient_type(::Type{TypicalSection}) = Zeros()

# --- Methods --- #

function get_residual(::TypicalSection, dx, x, y, p, t)
    # extract rates
    dh, dθ, dhdot, dθdot = dx
    # extract states
    h, θ, hdot, θdot = x
    # extract inputs
    L, M = y
    # extract parameters
    kh, kθ, m, Sθ, Iθ = p
    # calculate state rates
    return section_residual(dh, dθ, dhdot, dθdot, h, θ, hdot, θdot, L, M,
        kh, kθ, m, Sθ, Iθ)
end

# --- Performance Overloads --- #

function get_rate_jacobian(::TypicalSection, p)
    # extract parameters
    kh, kθ, m, Sθ, Iθ = p
    # return jacobian
    return section_rate_jacobian(m, Sθ, Iθ)
end

function get_state_jacobian(::TypicalSection, p)
    # extract parameters
    kh, kθ, m, Sθ, Iθ = p
    # return jacobian
    return section_state_jacobian(kh, kθ)
end

function get_input_jacobian(::TypicalSection)
    # return jacobian
    return section_input_jacobian()
end

function get_parameter_jacobian(::TypicalSection, dx, x, y, p, t)
    # extract rates
    dh, dθ, dhdot, dθdot = dx
    # extract states
    h, θ, hdot, θdot = x
    # extract inputs
    L, M = y
    # extract structural parameters
    kh, kθ, m, Sθ, Iθ = p
    # return jacobian
    return section_parameter_jacobian(h, θ, dhdot, dθdot)
end

# --- Convenience Methods --- #

function set_states!(x, model::TypicalSection; h, theta, hdot, thetadot)

    x[1] = h
    x[2] = theta
    x[3] = hdot
    x[4] = thetadot

    return x
end

function set_inputs!(y, model::TypicalSection; L, M)

    y[1] = L
    y[2] = M

    return y
end

function set_parameters!(p, model::TypicalSection; kh, ktheta, m, Stheta, Itheta)

    p[1] = kh
    p[2] = ktheta
    p[3] = m
    p[4] = Stheta
    p[5] = Itheta

    return p
end

function separate_states(model::TypicalSection, x)

    return (h = x[1], theta = x[2], hdot = x[3], thetadot = x[4])
end

function separate_inputs(model::TypicalSection, y)

    return (L = y[1], M = y[2])
end

function separate_parameters(model::TypicalSection, p)

    return (kh = p[1], ktheta = p[2], m = p[3], Stheta = p[4], Itheta = p[5])
end

# --- Plotting --- #

@recipe function f(model::TypicalSection, dx, x, y, p, t)

    framestyle --> :origin
    grid --> false
    xticks --> false
    xlims --> (-1.0, 1.0)
    ylims --> (-1.5, 1.5)
    label --> @sprintf("t = %6.3f", t)

    h = x[1]
    θ = x[2]

    xplot = [-0.5*cos(θ),  0.5*cos(θ)]
    yplot = [ 0.5*sin(θ)-h, -0.5*sin(θ)-h]

    xplot, yplot
end

# --- Internal Methods for this Model --- #

function section_residual(dh, dθ, dhdot, dθdot, h, θ, hdot, θdot, L, M, kh, kθ, m, Sθ, Iθ)
    r1 = dh - hdot
    r2 = dθ - θdot
    r3 = m*dhdot + Sθ*dθdot + kh*h + L
    r4 = Sθ*dhdot + Iθ*dθdot + kθ*θ - M
    return SVector(r1, r2, r3, r4)
end

section_rate_jacobian(m, Sθ, Iθ) = @SMatrix [1 0 0 0; 0 1 0 0; 0 0 m Sθ; 0 0 Sθ Iθ]

section_state_jacobian(kh, kθ) = @SMatrix [0 0 -1 0; 0 0 0 -1; kh 0 0 0; 0 kθ 0 0]

section_input_jacobian() = @SMatrix [0 0; 0 0; 1 0; 0 -1]

function section_parameter_jacobian(h, θ, dhdot, dθdot)
    return @SMatrix [0 0 0 0 0; 0 0 0 0 0; h 0 dhdot dθdot 0; 0 θ 0 dhdot dθdot]
end

# --- Internal Methods for Couplings --- #

# section velocities for coupling with 2D aerodynamic models
function section_velocities(U, θ, hdot, θdot)

    u = U
    v = U*θ + hdot
    ω = θdot

    return SVector(u, v, ω)
end

# section accelerations for coupling with 2D aerodynamic models
function section_accelerations(dhdot, dθdot)
    udot = 0
    vdot = dhdot
    ωdot = dθdot
    return SVector(udot, vdot, ωdot)
end

# derivatives
section_velocities_U(θ) = SVector(1, θ, 0)
section_velocities_θ(U) = SVector(0, U, 0)
section_velocities_hdot() = SVector(0, 1, 0)
section_velocities_θdot() = SVector(0, 0, 1)
section_accelerations_dhdot() = SVector(0, 1, 0)
section_accelerations_dθdot() = SVector(0, 0, 1)