text:
################################################################################
#
# The LLL-Reduction Algorithm
#
# Implementration of algorithm described in :
#  Lenstra, A.K., Lenstra, H.W. & Lovász, L. Math. Ann. (1982)
#            261: 515. https://doi.org/10.1007/BF01457454
#
################################################################################

"""
  lll(b₀::AbstractMatrix{Float64})::AbstractMatrix{Float64}

Performs Lenstra-Lenstra-Lovasv reduction on the n vectors in b₀.

Parameters
----------
- b₀ : vector of n vectors of n Float64

Returns
-------
- The LLL reduction of the vectors in `b₀`
"""
function lll(b₀::Vector{Vector{Float64}})::Vector{Vector{Float64}}
    b = deepcopy(b₀) # preserve argument by working on a copy
    n = length(b)

    # initialize mu and B with zeros
    μ = zeros(n,n)
    B = zeros(n)

    # function (*) of page 521
    function ☼(l)
        abs(μ[k,l]) < 0.5 && return
        r = round(μ[k,l])
        b[k] .-= r * b[l]
        μ[k,1:l-1] .-= r * μ[l,1:l-1]
        μ[k,l] -= r
    end

    b₂ = similar(b)
    for t in 1:n
        b₂[t] = copy(b[t])
        for j in 1:t-1
            μ[t,j] = dot(b[t], b₂[j]) / B[j]
            b₂[t] .-= μ[t,j] * b₂[j]
        end
        B[t] = dot(b₂[t], b₂[t])
    end

    k = 2
    while k <= n
        ☼(k-1)
        if B[k] >= (0.75 - μ[k,k-1]^2) * B[k-1]
            foreach(☼, k-2:-1:1)
            k += 1
        else
            μₒ = μ[k,k-1]
            Bₒ = B[k] + μₒ^2 * B[k-1]
            μ[k,k-1] = μₒ * B[k-1] / Bₒ
            B[k] = B[k] * B[k-1] / Bₒ
            B[k-1] = Bₒ

            bₒ = b[k] ; b[k] = b[k-1] ; b[k-1] = bₒ
            for j in 1:k-2
                tmp = μ[k,j] ; μ[k,j] = μ[k-1,j] ; μ[k-1,j] = tmp
            end
            for t in k+1:n
                tmp = μ[k,k-1] - μₒ * μ[t,k]
                μ[t,k-1] = μ[t,k] + μ[k,k-1] * tmp
                μ[t,k]   = tmp
            end

            k = max(2, k-1)
        end
    end
    b
end


text:
using ReachabilityAnalysis, SparseArrays
Dict([
    "opts" => Dict(:T=>20.0, :alg=>GLGM06(δ=0.5)),
    "X0"   => Singleton([0.0, 5.0, 100.0, 0.0]),
    "info" => Dict(
        "name" => "projectile",
        "dim" => 4,
        "linear" => true,
        "hybrid" => false
    )
])

ext: jl

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# --------------------------------------------------------------------------- #
# Abstract type of a processing Model

abstract type Model{I <: Image} end

# --------------------------------------------------------------------------- #
# Get the image type that the model is constructed for

imgtype{I <: Image}(:: Model{I}) = I

# --------------------------------------------------------------------------- #
# Reference to the variables of a model trained by backpropagation

weights(m :: Model) = error("Not implemented")

# --------------------------------------------------------------------------- #
# Reference to other internals of a model

state(m :: Model) = error("Not implemented")

# --------------------------------------------------------------------------- #
# Forward evaluation of the model

# Low-level function needed to support automated differentiation
density(w, s, data, ::Type{Model}) = error("Not implemented")

# Convenience function built on top
function density{I <: Image}(m :: Model{I}, img :: I)
  x = reshape(img.data, (imgsize(img)..., imgchannels(I), 1))
  return density(weights(m), state(m), x, typeof(m))[:,:,1]
end

# A full image may be too large to go through the model in a single run.
# So divide the image into patches, apply the model to the patches, and
# stich the density-patches together again.
function density_patched{I <: Image}(model :: Model{I}, img :: I;
                                     patchsize = 256, overlap = 24,
                                     at = Array{Float32},
                                     callback :: Function = (i, n) -> nothing)

  s = (n, m) = imgsize(img)
  d = overlap

  # TODO: The first two can be relaxed
  @assert patchsize >= 48 "patchsize must be larger than or equal to 48"
  @assert min(n, m) > patchsize "image sizes must be larger than patchsize"
  @assert d >= 8 "overlap must be larger than or equal to 8"

  # Find parameters for suitable patch regions
  # k = number of patches in y/x direction
  # l = good length of the patches in y/x direction
  k = ceil.( Int, (s .- patchsize) ./ (patchsize .- d) ) .+ 1
  l = ceil.( Int, (s .+ (k .- 1) .* d) ./ (8 .* k) ) .* 8

  mirrorinds(a, l, n) = (a+l-1 <= n) ? (a:(a+l-1)) : ([a:n; n:-1:(n-(a+l-2-n))])

  # Pack the model, i.e. bring it to graphics memory if suitable at is given
  # TODO: These functions are only introduced in the file training.jl.
  # Should be moved?
  w = packweights(weights(model), at = at)
  s = packstate(state(model), at = at)

  # Create patch data
  patchdata = zeros(Float32, l..., imgchannels(img), prod(k))
  for i in 1:k[1], j in 1:k[2]
    y = (i == 1) ? 1 : (i-1)*(l[1] - d)
    x = (j == 1) ? 1 : (j-1)*(l[2] - d)

    # Get indices for this patch. If the patch would protrude from the
    # image, continue the patch on a mirror version after the image border
    yind = mirrorinds(y, l[1], n)
    xind = mirrorinds(x, l[2], m)

    patchdata[:,:,:,sub2ind((k[2], k[1]), j, i)] = imgdata(img)[yind,xind,:]
  end

  # Load the patch data to gpu if suitable at is given
  pd = convert(at, patchdata)

  # Since Knet does not yet support slicing in 4-d arrays, convert pd to a 2d
  # array for slicing
  sz = size(pd)
  pd = reshape(pd, (prod(sz[1:end-1]), sz[end]))

  # Bring the density patches in rectangular layout
  densities = Array{Array{Float32, 2}}(k...)
  for i in 1:k[1], j in 1:k[2]
    ind = sub2ind((k[2], k[1]), j, i)
    data = reshape(pd[:, ind:ind], (sz[1:end-1]..., 1))
    densities[i, j] = convert(at, density(w, s, data, typeof(model)))[:,:,1]

    # Give some feedback about the status via callback
    callback(ind, prod(k))
  end

  # Merge the density patches with smooth interpolation profiles
  ftrans = [ 1 - sin(r*(pi/2) / (d-1))^2 for r in 1:d ]
  rtrans = reverse(ftrans)

  sprofile = [ones(patchsize - d); ]
  eprofile = sprofile[end:-1:1]
  mprofile = sprofile .* eprofile

  dens = zeros(Float32, (k .* l .- (k.-1) .* d)...)

  for i in 1:k[1], j in 1:k[2]

    # If we are at the first or last patches in the x or y direction,
    # use start or end profiles; use middle profiles otherwise
    if i == 1
      y = 1
      yprofile = [ones(l[1] - d); ftrans]
    elseif i == k[1]
      y = (i-1)*(l[1] - d)
      yprofile = [rtrans; ones(l[1] - d)]
    else
      y = (i-1)*(l[1] - d)
      yprofile = [rtrans; ones(l[1] - 2d); ftrans]
    end

    if j == 1
      x = 1
      xprofile = [ones(l[2] - d); ftrans]
    elseif j == k[2]
      x = (j-1)*(l[2] - d)
      xprofile = [rtrans; ones(l[2] - d)]
    else
      x = (j-1)*(l[2] - d)
      xprofile = [rtrans; ones(l[2] - 2d); ftrans]
    end

    profile = yprofile .* xprofile'
    dens[y:y+l[1]-1, x:x+l[2]-1] += profile .* densities[i, j]

  end

  return dens[1:n, 1:m]
end

# --------------------------------------------------------------------------- #
# Estimate labels
# Densities are the output of the application of a neural network regressor

function label(dens :: Array{Float32, 2};
               candidates = false, level = 25, clevel = 10)

  # Find local maxima that are sufficiently high

  label = Array{Int}[]
  cands = Array{Int}[]

  m, n = size(dens)

  for i in 1:m, j in 1:n
    xsel = max(i-1, 1):min(i+1, m)
    ysel = max(j-1, 1):min(j+1, n)

    if dens[i, j] < maximum(dens[xsel, ysel])
      continue
    end

    if dens[i, j] >= level
      push!(label, [j, i])
    elseif dens[i, j] >= clevel
      push!(cands, [j, i])
    end
  end

  label = isempty(label) ? zeros(Int, 2, 0) : hcat(label...)
  cands = isempty(cands) ? zeros(Int, 2, 0) : hcat(cands...)

  if candidates
    return Label(label), Label(cands)
  else
    return Label(label)
  end
end

function label(dens :: Array{Float32, 3}; kwargs...)
  return label.([dens[:,:,i] for i in 1:size(dens, 3)]; kwargs...)
end

function label(args...; kwargs...)
  dmap = density(args...; kwargs...)
  return label(convert(Array{Float32}, dmap))
end

# --------------------------------------------------------------------------- #
# Some auxiliary functions for building neural network models

include("models/primitives.jl")

# --------------------------------------------------------------------------- #
# Class of models from (Pen, Yang, Li, et al.; 2018)
# and (Xie, Noble, Zisserman; 2015)

abstract type FCModel{I, BN} <: Model{I} end

weights(m :: FCModel) = m.weights
state(m :: FCModel) = m.bnmoments
hasbatchnorm{I, BN}(:: FCModel{I, BN}) = BN

# Different architectures of this model class

include("models/unetlike.jl")
include("models/multiscale3.jl")
include("models/fcrna.jl")

max_stars_repo_name: tscode/DCellC

max_issues_repo_name: tscode/DCellC


text:
file = "inputs/models/plant_coffee.yml"

@testset "read_model()" begin
    model = read_model(file)

    @test all([haskey(model, i) for i in ["Metamer","Leaf"]])

    model = model["Leaf"]
    # @test typeof(model)
    @test typeof(model) <: LeafModels
    @test typeof(model.stomatal_conductance) == Medlyn{Float64}
    @test typeof(model.interception) == Translucent{Float64}
    @test typeof(model.photosynthesis) == Fvcb{Float64}

    @test model.stomatal_conductance.g0 == -0.03
    @test model.stomatal_conductance.g1 == 12.0

    @test model.interception.transparency == 0.0
    @test model.interception.optical_properties == σ(0.15, 0.9)

    @test model.photosynthesis.VcMaxRef == 200.0 # Given in the file
    @test model.photosynthesis.O₂ == 210.0 # Use default value
end;

# Test reading the meteo:

file = joinpath(dirname(dirname(pathof(PlantBiophysics))), "test", "inputs", "meteo.csv")
var_names = Dict(:temperature => :T, :relativeHumidity => :Rh, :wind => :Wind, :atmosphereCO2_ppm => :Cₐ)

@testset "read_weather()" begin
    meteo = read_weather(
        file,
        :temperature => :T,
        :relativeHumidity => (x -> x ./ 100) => :Rh,
        :wind => :Wind,
        :atmosphereCO2_ppm => :Cₐ,
        :Re_SW_f => :Ri_SW_f,
        date_format = DateFormat("yyyy/mm/dd")
    )

    @test typeof(meteo) <: Weather
    @test typeof(meteo) <: Weather
    @test NamedTuple(meteo.metadata) == (;name = "Aquiares", latitude = 15.0, altitude = 100.0, use = [:Rh, :clearness], file = file)
end;


text:
# Michelle Medina
# July 21, 2020
# Compute the cost vs. accuracy for
# Hierarchical basis vs. sparse grids

export cost1D, cost2D, cost3D, cost7D

using PyPlot
# ---------------------------------------------------------
function Base.:-(u::NodalBasis2D{T}, w::NodalBasis2D{T})::NodalBasis2D{T} where {T}
	value = NodalBasis2D(u.values - w.values)
	return value
end

using LinearAlgebra

function LinearAlgebra.norm(u::NodalBasis2D{T})::T where {T}
	return norm(u.values)
end

function cost1D(level::Int)::Int
	if level == 0
		return 2
	else
		return 2^(level-1)
	end
end

function cost2D(level::Int)::NTuple{2,Int}
    full = 0
	sparse = 0
	for l in 0:level, m in 0:level
        full = full + cost1D(l)*cost1D(m)
       	if l + m <= level
        	sparse = sparse + cost1D(l)*cost1D(m)
        end
	end
    return (full, sparse)
end

function cost3D(level::Int)::NTuple{2,Int}
    full   = 0
    sparse = 0
    for l in 0:level, m in 0:level, n in 0:level
        full = full + cost1D(l)*cost1D(m)*cost1D(n)
        if l + m + n <= level
            sparse = sparse + cost1D(l)*cost1D(m)*cost1D(n)
        end
    end
    return (full, sparse)
end

function cost7D(level::Int)::NTuple{2,Float64}
    full   = 0.0
    sparse = 0.0
    for l in 0:level, m in 0:level, n in 0:level, o in 0:level, p in 0:level, q in 0:level, r in 0:level
        full = full + 1.0 * cost1D(l)*cost1D(m)*cost1D(n)*cost1D(o)*cost1D(p)*cost1D(q)*cost1D(r)
        if l + m + n + o + p + q + r <= level
            sparse = sparse + 1.0 * cost1D(l)*cost1D(m)*cost1D(n)*cost1D(o)*cost1D(p)*cost1D(q)*cost1D(r)
        end
    end
    return (full, sparse)
end

#Cost (number of total points) vs. accurary (interpolation error)
#for both full basis and sparse basis
error_full = []
error_sparse = []
cost_full = []
cost_sparse = []

level = 3
for l in 0:level
	N = 2^l + 1

	xVal = collect(range(0.0, stop=1.0, length = N))
	yVal = collect(range(0.0, stop=1.0, length = N))
	nodalf = calcNodal2D((x,y)->exp(x)*exp(y), xVal, yVal)
	modalf = Nodal_2_H(nodalf)
	modalf_sparse = Nodal_2_H_sparse(nodalf,l)
	costF, costS = cost2D(l)

	# Compute the interpolation error
	le = 10 * (2^l+1)
	h = 1/le
	xNew = collect(range(0.0 + h/2, stop=1.0 - h/2, length = le))
	yNew = collect(range(0.0 + h/2, stop=1.0 - h/2, length = le))
	# xNew = xNew[2:end-1]
	# yNew = yNew[2:end-1]
	f_analytic = calcNodal2D((x,y)->exp(x)*exp(y), xNew, yNew)
	f_nodal_from_modalf = calcNodal2D((x,y)->evaluate(modalf, x, y), xNew, yNew)
	f_nodal_from_modalf_sparse = calcNodal2D((x,y)->evaluate(modalf_sparse, x, y), xNew, yNew)

	# error_f = sum(abs.((f_analytic - f_nodal_from_modalf).values))
	# error_s = sum(abs.((f_analytic - f_nodal_from_modalf_sparse).values))
	error_f = sqrt(sum((f_analytic.values - f_nodal_from_modalf.values).^2)/length(f_analytic.values))
	error_s = sqrt(sum((f_analytic.values - f_nodal_from_modalf_sparse.values).^2)/length(f_analytic.values))

	append!(error_full, error_f)
	append!(error_sparse, error_s)
	append!(cost_full, costF)
	append!(cost_sparse, costS)

	@show error_f, error_s, cost_full, cost_sparse

	# if l == 4
	# 	fig = figure()
	# 	subplot(1,2,1)
	# 	contourf((f_analytic - f_nodal_from_modalf).values)
	# 	colorbar()
	#
	# 	subplot(1,2,2)
	# 	contourf((f_analytic - f_nodal_from_modalf_sparse).values)
	# 	colorbar()
	# 	show()
	# end
	#
	#
end

# loglog(error_full, cost_full, "g-o", label = "Full")
# loglog(error_sparse, cost_sparse, "b-o", label = "Sparse")
# title("Cost vs Interpolation Error")
# legend()
# ylabel("Cost")
# xlabel("Interpolation Error")
# show()

# for l in 0:10
# 	@show l cost7D(l)
# end

# Heat Map of how coefficients fall off with level
# Tests how coefficients fall off with level
# n = 5
# xValue = collect(range(0.0, stop=1.0, length = 2^n + 1))
# yValue = collect(range(0.0, stop=1.0, length = 2^n + 1))
# nodalTest = calcNodal2D((x,y)->sin(pi*x)*sin(pi*y) + cos(pi*x)*cos(pi*y), xValue, yValue)
# modal = Nodal_2_H(nodalTest)
#
# A = zeros(6,6)
# for l in 0:5, m in 0:5
# 	A[l+1, m+1] = maximum(abs.(modal.levels[l+1, m+1].coefficients))
# end
#
# # use heat map
# using PyPlot
# imshow(A, vmin=0, vmax=0.5)
# title("f(x,y) = sin(pi*x)*sin(pi*y) + cos(pi*x)*cos(pi*y)")
# colorbar()
# show()


text:
# Copyright (c) Microsoft Corporation.
# Licensed under the MIT license.

module LSTMData

export LSTMInput, LSTMOutput, empty_lstm_output, load_lstm_input

struct LSTMInput
    main_params::Matrix{Float64}
    extra_params::Matrix{Float64}
    state::Matrix{Float64}
    sequence::Matrix{Float64}
end

mutable struct LSTMOutput
    objective::Float64
    gradient::Vector{Float64}
end

empty_lstm_output() = LSTMOutput(0.0, [])

function load_lstm_input(fn::AbstractString)::LSTMInput
    open(fn) do io
        line = split(readline(io), " ")
        layer_count = parse(Int, line[1])
        char_count = parse(Int, line[2])
        char_bits = parse(Int, line[3])
        readline(io)
        main_params = hcat([parse.(Float64, split(readline(io), " ")) for i=1:2*layer_count]...)
        readline(io)
        extra_params = hcat([parse.(Float64, split(readline(io), " ")) for i=1:3]...)
        readline(io)
        state = hcat([parse.(Float64, split(readline(io), " ")) for i=1:2*layer_count]...)
        readline(io)
        sequence = hcat([parse.(Float64, split(readline(io), " ")) for i=1:char_count]...)
        return LSTMInput(main_params, extra_params, state, sequence)
    end
end

end

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
# insertion.jl
# functions for inserting a point into a Delaunay tessellation

#function insertpoint()
# TODO: insertpoint! (overload push!() ?)
# INPUT: DT object, point
# OUTPUT: nothing
#
# Walk to containing tet
# Add point to DT points
# pop containing tet
# call a flip14 on the containing tet and append tets to DT tets
# loop through old tet neighbors, updating neighbor information (subroutine)
#
# Push new tets on stack
# while stack is not empty
#   pop a tet from the stack -> pabc
#   Get the neighbor opposite the new point -> abcd
#   if insphere of tet pabc and d is positive
#       Flip the tets (subroutine)
#
#

# TODO: update neighbors
# INPUT: old tet neighbors, old tet index, new tets
# OUTPUT: nothing

# TODO: flip neighboring tets


text:
# -------

type SymbolNode
    sym
    parent
    children::Array{SymbolNode}

    SymbolNode(sym, parent = nothing) = new(sym, parent, SymbolNode[])
end

type DataNode
    name
    value

    DataNode(name, value) = new(name, value)
end

global _g_records = SymbolNode(nothing, nothing)
global _data_records = {}

# -------

function getindex(sn::SymbolNode, keys...)
    haskey(sn, keys...) || error("Key not found: $(keys)")

    eval(Main, to_expr(keys...))
end

function setindex!(sn::SymbolNode, keys...)
    key_count = length(keys)
    key_count > 0 || return

    did_allocate = false

    if sn.sym == nothing
        (sn, did_allocate) = push!(sn, keys[1])
    end
    sn.sym != keys[1] && error("head symbol node does not match key $(keys[1])")

    curr_node = sn
    for i in 2:length(keys)
        # if a node higher in the tree was already a leaf
        # then this set of keys is already recorded
        if !did_allocate && length(curr_node.children) < 1
            break
        end

        (curr_node, did_allocate) = push!(curr_node, keys[i])
    end

    empty!(curr_node.children)

    nothing
end

function haskey(sn::SymbolNode, keys...)
    key_count = length(keys)
    key_count > 0 || return

    curr_index = 1
    curr_node = sn

    while curr_index <= key_count && curr_node != nothing && curr_node.sym == keys[curr_index]
        if curr_index == key_count || length(curr_node.children) < 1
            return true
        end

        curr_index += 1
        curr_node = find_child(curr_node, keys[curr_index])
    end

    false
end

function push!(sn::SymbolNode, child_sym)
    child = find_child(sn, child_sym)

    did_allocate = false

    if child == nothing
        did_allocate = true
        child = SymbolNode(child_sym, sn)
        push!(sn.children, child)
    end

    child, did_allocate
end

function store_data(name, value)
    node = DataNode(name, value)
    push!(_data_records, node)
    node
end

function delete!(sn::SymbolNode, keys...)
    key_count = length(keys)

    key_count == 0 && return sn

    if key_count == 1
        return sn.sym == keys[1] ? nothing : sn
    end

    curr_index = 1
    curr_node = sn

    while curr_index <= key_count && curr_node != nothing && curr_node.sym == keys[curr_index]
        if curr_index == key_count
            delete!(curr_node)
            return sn
        end

        curr_index += 1
        curr_node = find_child(curr_node, keys[curr_index])
    end

    sn
end

function delete!(sn::SymbolNode)
    parent = sn.parent

    if parent == nothing
        empty!(sn.children)
        return sn
    end

    index = find_child_index(parent, sn.sym)
    index == -1 && error("nodes parent does not contain node... something has gone horribly wrong")

    splice!(parent.children, index)

    parent
end

function show(io::IO, sn::SymbolNode, depth = 1)
    println(io, ">", sn.sym)
    for child in sn.children
        print("  "^depth)
        show(io, child, depth + 1)
    end
    println("  "^(depth - 1), "<")
end

# -------

function record(keys...)
    global _g_records
    setindex!(_g_records, keys...)
end

function fetch_records()
    build_record_tree()
end

function take_records()
    tree = build_record_tree()
    delete!(_g_records)

    for node = _data_records
        tree[node.name] = node.value
    end
    empty!(_data_records)

    tree
end

function hasrecord(keys...)
    if length(keys) > 0
        global _g_records

        sub_head = find_child(_g_records, keys[1])
        if sub_head != nothing
            return haskey(sub_head, keys...)
        end
    end

    false
end

function getrecord(keys...)
    if length(keys) > 0
        global _g_records

        sub_head = find_child(_g_records, keys[1])
        if sub_head != nothing
            return getindex(sub_head, keys...)
        end
    end

    error("Key not found: $(keys)")
end

# -------

function to_expr(syms...)
    # there might well be a better way to do this, for now we'll eval strings
    typeof(syms[1]) <: Symbol || error("first piece of key must be a Symbol")

    ex = string(syms[1])
    for sym in syms[2:end]
        symtype = typeof(sym)
        if symtype <: Symbol
            ex *= ".$sym"
        elseif symtype <: Int
            ex *= "[$sym]"
        elseif symtype <: String
            ex *= "[$(repr(sym))]"
        else
            error("unrecognized key index: $sym")
        end
    end

    parse(ex)
end

function find_child_index(sn::SymbolNode, child_sym)
    isempty(sn.children) && return -1

    for i in 1:length(sn.children)
        child_sym == sn.children[i].sym && return i
    end

    return -1
end

find_child(n::SymbolNode, k) = get(n.children, find_child_index(n, k), nothing)

function build_record_tree(sn::SymbolNode = _g_records, out = Dict(), keys = Any[])
    keys = deepcopy(keys)
    sn.sym != nothing && push!(keys, sn.sym)  # special case for _g_records

    for child in sn.children
        num_grandchildren = length(child.children)

        if num_grandchildren > 0
            out[child.sym] = Dict()
            build_record_tree(child, out[child.sym], keys)
        else
            try
                out[child.sym] = eval(Main, to_expr(keys..., child.sym))
            catch err
                out[child.sym] = err
            end
        end
    end

    out
end

# -------


text:
function metaplex_gillespie(mp::Metaplex{T1,T2}, Xi0::BitVector, Xμ0::Vector{Int}, β, D::Vector; nmax=length(Xi0), tmax=100.0, sampling_method=:tree) where {T1,T2}
    N = nv(mp.g) # number of individuals
    g = mp.g
    M = nv(mp.h) # number of populations
    h = mp.h
    #od = map(k-> k>0 ? 1/float(k) : 0, outdegree(h))
    #od = outdegree(h)
    od = outdegree(h) .> 0
    Xi = copy(Xi0) # epidemic state of each individual
    Xμ = copy(Xμ0) # location of each individual
    popcounts = zeros(Int,2,M) # number of of individuals susceptible and infected for each component
    pops_s = [Set{Int}() for i in 1:M] # Population to individuals without traversing the whole list of individuals
    pops_i = [Set{Int}() for i in 1:M]
    for k in 1:N
        if !Xi[k]
            push!(pops_s[Xμ[k]],k)
            #popcounts[1,Xμ[k]] += 1
        else
            push!(pops_i[Xμ[k]],k)
            #popcounts[2,Xμ[k]] += 1
        end
    end
    popcounts[1,:] = map(length, pops_s)
    popcounts[2,:] = map(length, pops_i)
    a = zeros(N+2*M) # aggregate reactions for each individual becoming infected and individuals migrating from each population
    for k in 1:N
        if !Xi[k]
            l = length(filter(i->Xi[i] && Xμ[i]==Xμ[k], inneighbors(g,k)))
            a[k] = β*l # probability of k recieving infection is β time its number of infected incoming neighbors in the same population
            #nreactions += l
        end
    end
    a[N+1:N+M] .= D[1].*popcounts[1,:].*od
    a[N+M+1:N+2*M] .= D[2].*popcounts[2,:].*od
    if sampling_method == :tree
        a = CategoricalTree(a)
    end
    a0 = sum(a)
    ts = zeros(nmax)
    pcs = zeros(Int,nmax,2,M)
    pcs[1,:,:] .= popcounts
    t = 0.0
    n = 1
    while n < nmax && t < tmax && a0 > 0
        τ = log(1/rand())/a0
        k = rand_categorical(a, a0)
        if k in 1:N # node k becomes infected
            Xi[k] = true
            delete!(pops_s[Xμ[k]], k)
            push!(pops_i[Xμ[k]], k)
            popcounts[1,Xμ[k]] = length(pops_s[Xμ[k]])
            popcounts[2,Xμ[k]] = length(pops_i[Xμ[k]])
            a[N+Xμ[k]] = D[1]*popcounts[1,Xμ[k]]*od[Xμ[k]]
            a[N+M+Xμ[k]] += D[2]*popcounts[2,Xμ[k]]*od[Xμ[k]]
            a[k] = 0.0
            for i in Iterators.filter(j->!Xi[j] && Xμ[j]==Xμ[k], outneighbors(g,k))
                a[i] += β
            end
        elseif k in N+1:N+M # a random susceptible node from population k-N moves to another population
            μ = k-N
            ν = rand(outneighbors(h,μ))
            i = rand(pops_s[μ])
            Xμ[i] = ν
            delete!(pops_s[μ], i)
            push!(pops_s[ν], i)
            popcounts[1,μ] = length(pops_s[μ])
            popcounts[1,ν] = length(pops_s[ν])
            a[N+μ] = D[1]*popcounts[1,μ]*od[μ]
            a[N+ν] = D[1]*popcounts[1,ν]*od[ν]
            l = length(filter(j->Xi[j] && Xμ[j]==ν, inneighbors(g,i)))
            a[i] = β*l
        else # a random infected node from k-N-M moves to another population
            μ = k-N-M
            ν = rand(outneighbors(h,μ))
            i = rand(pops_i[μ])
            Xμ[i] = ν
            delete!(pops_i[μ], i)
            push!(pops_i[ν], i)
            for j in Iterators.filter(l->!Xi[l], outneighbors(g,i))
                if Xμ[j] == μ
                    a[j] -= β
                elseif Xμ[j] == ν
                    a[j] += β
                end
            end
            popcounts[2,μ] = length(pops_i[μ])
            popcounts[2,ν] = length(pops_i[ν])
            a[N+M+μ] = D[2]*popcounts[2,μ]*od[μ]
            a[N+M+ν] = D[2]*popcounts[2,ν]*od[ν]
        end
        a0 = sum(a)
        t += τ
        ts[n+1] = t
        pcs[n+1,:,:] .= popcounts
        n += 1
    end
    return ts[1:n], pcs[1:n,:,:]
end

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
@doc raw"""
    WeightedParticleEnsemble{T}

An ensemble of `N` particles, each of dimension `n`.
"""
struct WeightedParticleEnsemble{T} <: WeightedParticleRepresentation{Vector{T}}
    positions::Matrix{T}
    weights::StatsBase.ProbabilityWeights
end

particle_dim(ens::WeightedParticleEnsemble)    = size(ens.positions, 1)

no_of_particles(ens::WeightedParticleEnsemble) = size(ens.positions, 2)

get_pos(ens::WeightedParticleEnsemble, i) = view(ens.positions, :, i)

get_weight(ens::WeightedParticleEnsemble, i) = ens.weights[i]

sum_of_weights(ens::WeightedParticleEnsemble) = sum(ens.weights)

Base.show(io::IO, ens::WeightedParticleEnsemble) = print(io, "Weighted particle ensemble
    # of particles: ", no_of_particles(ens),"
    particle type:  ", particle_dim(ens),"-dimensional ", eltype(ens))

function WeightedParticleEnsemble(positions::Matrix)
    N = size(positions)[2]
    return WeightedParticleEnsemble(positions, StatsBase.ProbabilityWeights(fill(1/N, N)))
end

function WeightedParticleEnsemble(model::HiddenStateModel, N::Int)
    return WeightedParticleEnsemble(hcat([initialize(model) for i in 1:N]...), StatsBase.ProbabilityWeights(fill(1/N, N)))
end

function propagate!(ens::WeightedParticleEnsemble{S}, model::HiddenStateModel{Vector{S}, ContinuousTime}, dt) where S
    propagate!(ens.positions, model, dt)
end

mean(ens::WeightedParticleEnsemble) = Statistics.mean(ens.positions, ens.weights, dims=2)
cov(ens::WeightedParticleEnsemble)  = Statistics.cov(ens.positions, ens.weights, 2, corrected=false)
var(ens::WeightedParticleEnsemble)  = Statistics.var(ens.positions, ens.weights, 2, corrected=false)

function resample!(ens::WeightedParticleEnsemble)
    N   = no_of_particles(ens)
    idx = StatsBase.sample(1:N, ens.weights, N)
    ens.positions .= view(ens.positions, :, idx)
    ens.weights   .= 1/N
    return ens.positions
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using ModelObjectsLH, Test

mo = ModelObjectsLH;

struct TestObj4 <: ModelObject
    objId :: ObjectId
end

function object_id_test()
    @testset "ObjectId" begin
        # Simplest case
        id1 = SingleId(:id1);
        o1 = ObjectId(id1);
        @test !mo.has_parent(o1)

        # Roundtrip between objectId and String
        s1 = make_string(id1);
        @test isequal(s1, "id1")
        o1a = make_object_id(s1)
        @test isequal(o1, o1a);
        @test o1 == o1a;

        # Index, no parents
        o2 = ObjectId(:id2, 2)
        println(o2);
        @test mo.own_index(o2) == [2]
        @test mo.own_name(o2) == :id2

        # Has id1 as parent
        o3 = ObjectId(:id3, 2, o1);
        println(o3)
        p3 = mo.get_parent_id(o3);
        @test mo.is_parent_of(p3, o3)
        s3 = mo.make_string(o3);
        @test isequal(s3, "id1 > id3[2]")
        o3a = make_object_id(s3)
        @test isequal(o3, o3a);
        @test o3 == o3a;

        # Make child id
        o4 = ObjectId(:id4, p3);
        obj4 = TestObj4(o4);
        childId = mo.make_child_id(obj4, :child);
        @test isequal(mo.own_id(childId),  SingleId(:child))
        @test isequal(o4, mo.get_parent_id(childId))
        @test isequal(mo.own_name(obj4), :id4)

        # Check `isequal` when "depth" of `ObjectId`s is different
        @test !isequal(obj4.objId, childId)
    end
end

function dict_object_id_test()
    @testset "Dict ObjectId" begin
        objIdV = [
            ObjectId(:id1),
            ObjectId(:id2, 1),
            make_child_id(ObjectId(:id1), :id3)];
        d = Dict{ObjectId, Int}([oId => j  for (j, oId) in enumerate(objIdV)]);
        for (id1, val) in d
            @test id1 isa ObjectId;
            @test val isa Integer;
        end
        for objId in objIdV
            @test haskey(d, objId);
        end
    end
end

function oid_broadcast_test()
    @testset "Broadcasting" begin
        sV = [ObjectId(:x), ObjectId(:y)];
        eqV = ObjectId(:x) .== sV;
        @test eqV[1]
        @test !eqV[2]
    end
end

@testset "ObjectId" begin
    object_id_test();
    dict_object_id_test();
    oid_broadcast_test();
end

# --------------


text:
const scipy_opt = PyNULL()

const igl = PyNULL()
const scipy = PyNULL()

# igl = pyimport("igl")
# scipy = pyimport("scipy")
function igl_gaussian_curvature(verts,faces)
    println("new gaussian curv")
    v = array_of_array_to_matrix(verts)
    f = array_of_array_to_matrix(faces) .-1
    K = igl.gaussian_curvature(v, f)
    m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)
    # m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)
    minv = scipy.sparse.diags(1.0 ./ m.diagonal())
    kn = minv.dot(K)
    return kn
end
function igl_mean_curvature(verts,faces)
    v = array_of_array_to_matrix(verts)
    f = array_of_array_to_matrix(faces) .-1
    l = igl.cotmatrix(v, f)
    m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_BARYCENTRIC)
    # m = igl.massmatrix(v, f, igl.MASSMATRIX_TYPE_VORONOI)

    minv = scipy.sparse.diags(1.0 ./ m.diagonal())

    hn = -minv.dot(l.dot(v))
    h = scipy.linalg.norm(hn, axis=1)
    return h
end

function igl_mean_curvature_flow(v,f,n_step = 10, delta = 0.001)
    new_v = v .- centroid(v)
    new_v = new_v / sqrt(surface_area(new_v,f))
    l = igl.cotmatrix(new_v, f .- 1)

    for i in 1:n_step
        m = igl.massmatrix(new_v,f .- 1, igl.MASSMATRIX_TYPE_BARYCENTRIC)
        s = (m - delta * l)
        b = m.dot(new_v)
        new_v = scipy.sparse.linalg.spsolve(s, b)
        new_v = new_v .- centroid(new_v)
        new_v = new_v / sqrt(surface_area(new_v,f))

    end
    return new_v
end


text:
#=##############################################################################################
Filename: vortexlattice.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this file is meant to provide convenience functions for building VortexLattice objects.
        This includes convenience constructors for common wing configurations.
=###############################################################################################

"""
    VortexLatticeSystem{T1, T2, T3, T4}

Defines a wing system for use with VortexLattice.jl.

# Fields:

* `system::VortexLattice.System` : a VortexLattice System struct
* `lifting_line_rs`
* `lifting_line_chords`
* `grids::Array{Array{Float64,2}}` : a vector of 2-D arrays, each 2-D array containing wing panel corner locations; returned by `VortexLattice.wing_to_surface_panels()`

"""
struct VortexLatticeSystem{TF}
    system::VL.System{TF}
    lifting_line_rs::Vector{Array{TF,2}}
    lifting_line_chords::Vector{Vector{TF}}
    cmacs::Vector{TF}
    airfoil_polars
end

"""
Constructor for `VortexLatticeSystem`.

# Arguments:

* `lexs::Array{Array{Float64,1}}` : array of vectors of x coordinates of the leading edge
* `leys::Array{Array{Float64,1}}` : array of vectors of y coordinates of the leading edge
* `lezs::Array{Array{Float64,1}}` : array of vectors of z coordinates of the leading edge
* `chords::Array{Array{Float64,1}}` : array of vectors of chord lengths
* `thetas::Array{Array{Float64,1}}` : array of vectors of twist angles in radians
* `phis::Array{Array{Float64,1}}` : array of vectors of dihedral angles in radians
* `nspanwisepanels::Array{Float64,1}` : vector of the number of desired spanwise panels
* `nchordwisepanels::Array{Float64,1}` : vector of the number of desired chordwise panels

# Keyword Arguments:

* `wing_cambers = [fill((xc) -> 0, length(lex)) for lex in lexs]` : vector of vectors of functions describing the camber line
* `spanwise_spacings = fill(VL.Uniform(), length(lexs))` : vector of `<: VortexLattice.AbstractSpacing` objects describing the spanwise spacing
* `chordwise_spacings = fill(VL.Uniform(), length(lexs))` : vector of `<: VortexLattice.AbstractSpacing` objects describing the chordwise spacing
* `Vref = 1.0` : velocity used to normalize
* `symmetric = true` : whether or not to mirror the geometry about the x-z plane
* `iref = 1` : index of `lexs`, `leys`, etc. to use when calculating reference properties
* `static_margin = 0.10` : static margin of the `iref`th lifting surface

# Returns:

* `vlmsystem::VortexLatticeSystem`

"""
function VortexLatticeSystem(lexs, leys, lezs, chords, thetas, phis, nspanwisepanels, nchordwisepanels;
            cambers = [fill((xc) -> 0, length(lex)) for lex in lexs],
            spanwise_spacings = fill(VL.Uniform(), length(lexs)),
            chordwise_spacings = fill(VL.Uniform(), length(lexs)),
            Vref = 1.0, symmetric = true, iref = 1, static_margin = 0.10, liftingline_x_over_c = 0.25
            ) where T <: AbstractArray

    @assert typeof(lexs[1]) <: AbstractArray "lexs must be a vector of vectors"
    @assert size(lexs) == size(leys) && size(lexs[1]) == size(leys[1])   "size of leading edge xs and ys not consistent"
    @assert size(lexs) == size(lezs) && size(lexs[1]) == size(lezs[1])   "size of leading edge xs and zs not consistent"
    @assert size(lexs) == size(chords) && size(lexs[1]) == size(chords[1]) "size of leading edge xs and chords not consistent"
    @assert size(lexs) == size(thetas) && size(lexs[1]) == size(thetas[1]) "size of leading edge xs and thetas not consistent"
    @assert size(lexs) == size(phis) && size(lexs[1]) == size(phis[1]) "size of leading edge xs and phis not consistent"

    nwings = length(lexs)

    # build grids and surfaces
    gridssurfaces = [VL.wing_to_surface_panels(lexs[i], leys[i], lezs[i], chords[i], thetas[i], phis[i], nspanwisepanels[i], nchordwisepanels[i];
        fc = cambers[i], spacing_s=spanwise_spacings[i], spacing_c=chordwise_spacings[i]) for i in 1:nwings]

    # extract grids and surfaces
    grids = [gridsurface[1] for gridsurface in gridssurfaces]
    surfaces = [gridsurface[2] for gridsurface in gridssurfaces]

    # build reference
    ref = get_reference(lexs[iref], leys[iref], chords[iref], static_margin, Vref; symmetric)

    # build freestream
    alpha = 3.0 * pi/180
    beta = 0.0
    Omegas = StaticArrays.@SVector zeros(3)
    vlmfreestream = VL.Freestream(Vref, alpha, beta, Omegas)

    # get lifting line geometry
    rs, chords = VL.lifting_line_geometry(grids)
    cmacs = [FM.trapz(r[2,:], chords[i]) / (r[2,end] - r[2,1]) for (i,r) in enumerate(rs)]
    # build system
    system = VL.steady_analysis(surfaces, ref, vlmfreestream; symmetric)
    vlmsystem = VortexLatticeSystem(system, rs, chords, cmacs, nothing)

    return vlmsystem
end

"""
    get_reference(lex, ley, chord, static_margin, Vref)

# Arguments:

* `lex`
* `ley`
* `chord`
* `static_margin`
* `Vref`

"""
function get_reference(lex, ley, chord, static_margin, Vref; symmetric = true)

    bref = ley[end] * 2
    Sref = FM.trapz(ley, chord) # area of starboard wing only
    if symmetric; Sref *= 2; end
    cref = Sref / bref
    lex_midpoint = get_midpoints(lex)
    dy = [ley[i+1] - ley[i] for i in 1:length(ley)-1]
    xmean = weightedmean(lex_midpoint, dy) # might be worth testing
    rref = [xmean + cref * (0.25 - static_margin), 0.0, 0.0]
    Vinf = Vref
    ref = VL.Reference(Sref, cref, bref, rref, Vinf)

    return ref
end

weightedmean(set, weights) = sum(set .* weights) / sum(weights)

get_midpoints(a) = [(a[i] + a[i+1])/2 for i in 1:length(a)-1]

"""
    simplewing_system(b, TR, AR, θroot, θtip, le_sweep, ϕ;
        xle=[0.0, 0.0], yle=[0.0, b/2], zle=[0.0, 0.0],
        cmac=b / AR, # AR = b/c => c = b/AR
        Sref=b * cmac,
        chord=[cmac * 2 / (1 + TR), cmac * 2 / (1/TR + 1)], twist=[θroot, θtip], phi=[0.0,0.0],
        static_margin=0.10,
        Vinf=1.0,
        Vref=Vinf,
        alpha=0.0,
        beta=0.0,
        Omega=[0.0, 0.0, 0.0],
        camber=fill((xc) -> 0, 2),
        npanels=50, nchordwisepanels=1,
        symmetric=true,
        spacing_s = symmetric ? VL.Cosine() : VL.Sine(),
        spacing_c=VL.Uniform(),
        mirror=!symmetric,
        kwargs...
    )

Convenience constructor for a single wing `VortexLatticeSystem`.

# Keyword Arguments:

* `wing_b`: wing span
* `wing_TR`: wing taper ratio
* `wing_AR`: wing aspect ratio
* `wing_θroot`: wing root twist (radians)
* `wing_θtip`: wing tip twist (radians)
* `xle`: leading-edge x-coordinates
* `yle`: leading-edge y-coordinates
* `zle`: leading-edge z-coordinates
* `cmac`: mean aerodynamic chord
* `Sref`: reference wing area
* `wing_chord`: chord at each wing section
* `wing_twist`: twist at each wing section
* `wing_phi`: wing section rotation about the x-axis
* `static_margin`: static margin of the wing
* `Vinf`: freestream velocity
* `Vref`: reference velocity; defaults to freestream velocity
* `alpha`: angle of attack (radians)
* `beta`: sideslip angle (radians)
* `Omega`: rotational velocity around the reference location
* `wing_camber`: camberline function for each section
* `wing_npanels`: number of spanwise panels for the wing
* `wing_nchordwisepanels`: number of chordwise panels for the wing
* `symmetric::Bool`: Flag indicating whether a mirror image across the X-Z plane should be used when calculating induced velocities. Defaults to true. From my understanding, this solves the half span and uses a mirrored aproach for induced velocities, whereas mirror will physically mirror the geometry for you and solve it all together (Tyler).
* `mirror::Bool`: mirrors the geometry over the x-z plane
* `wing_spacing_s::Int`: spanwise panel spacing scheme for the wing; defaults to Cosine (or Sine if mirror==true, which results in Cosine anyway)
* `wing_spacing_c::Int`: chordwise panel spacing scheme for the wing; defaults to Uniform()

"""
function simplewing_system(b, TR, AR, θroot, θtip, le_sweep, ϕ;
    yle = [0.0,b/2], xle = [0.0, yle[2] * tan(le_sweep)], zle = [0.0, yle[2] * tan(ϕ)],
    cmac = b/AR, # AR = b/c => c = b/AR
    Sref = b*cmac,
    chord = [cmac * 2 / (1 + TR), cmac * 2 / (1 + 1/TR)], twist=[θroot, θtip], phi=[ϕ, ϕ],
    static_margin = 0.10,
    Vref = 1.0, # in case you want the reference velocity to be different than the freestream velocity
    camber = fill((xc) -> 0, length(xle)),
    npanels = 50, nchordwisepanels=1,
    symmetric = true, # if true, we don't need to mirror geometry. Good option if parameters/geometry is symmetric and you don't care about lateral stability derivatives
    mirror = !symmetric,
    spacing_s = symmetric ? VL.Cosine() : VL.Sine(), # if mirrored, sine spacing becomes cosine spacing
    spacing_c = VL.Uniform(), # should only ever need uniform chordwise spacing
    kwargs...
)

    Vinf = 1.0
    alpha=0.0
    beta=0.0
    Omega=StaticArrays.@SVector [0.0,0.0,0.0]

    #=
    cr + ct = cmac * 2
    ct = λ * cr
    ct = cmac * 2 - cr
    λ * cr = cmac * 2 - cr
    cr = 2 * cmac / (1 + λ)

    cr = cmac * 2 - ct
    cr = ct / λ
    cmac * 2 - ct = ct / λ
    cmac * 2 = ct * (1 + 1/λ)
    ct = cmac * 2 / (1 + 1/λ)
    =#

    # Check for errors
    @assert size(xle) == size(yle)    "size of leading edge xs and ys not consistent"
    @assert size(yle) == size(zle)    "size of leading edge ys and zs not consistent"
    @assert size(xle) == size(chord)    "size of leading edge xs and chord not consistent"
    @assert size(xle) == size(twist)    "size of leading edge xs and twist not consistent"
    @assert size(xle) == size(phi)    "size of leading edge xs and phi not consistent"
    @assert yle[1] == 0.0    "yle must begin at the root (0.0)"
    @assert yle[end] == b/2    "yle must end at the tip (span/2)"
    @assert !(symmetric && mirror)    "symmetric and mirror cannot both be true"

    # construct surfaces and grids
    grid, surface = VL.wing_to_surface_panels(xle, yle, zle, chord, twist, phi, npanels, nchordwisepanels;
                        fc = camber, spacing_s, spacing_c, mirror)

    grids = [grid]
    surfaces = [surface]

    # reference parameters
    cgx = cmac * (0.25 - static_margin)
    rref = [cgx, 0.0, 0.0]
    reference = VL.Reference(Sref, cmac, b, rref, Vref)

    # initialize dummy freestream parameters
    fs = VL.Freestream(Vinf, alpha, beta, Omega)

    # perform steady state analysis
    vlmsystem = VL.steady_analysis(surfaces, reference, fs; symmetric)
    lifting_line_rs, lifting_line_chords = VL.lifting_line_geometry(grids)
    cmacs = [FM.trapz(lifting_line_r[2,:], lifting_line_chords[i]) / (lifting_line_r[end] - lifting_line_r[1]) for (i,lifting_line_r) in enumerate(lifting_line_rs)]
    wing_system = VortexLatticeSystem(vlmsystem, lifting_line_rs, lifting_line_chords, cmacs, nothing)

    return wing_system
end

"""
    simpleairplanesystem(;
        # wing geometry
        wing_b=2.0, wing_TR=0.8, wing_AR=8.0, wing_θroot=0.0, wing_θtip=0.0,
        wing_xle=[0.0,0.0], wing_yle=[0.0, wing_b/2], wing_zle=[0.0,0.0],
        wing_cmac=wing_b/wing_AR, # AR = b/c => c = b/AR
        wing_chord=[wing_cmac * 2 / (1 + wing_TR), wing_cmac * 2 / (1/wing_TR + 1)],
        wing_theta=[wing_θroot, wing_θtip],
        wing_phi=[1.0,1.0] * 3pi/180.0,
        wing_camber=fill((xc) -> 0, length(wing_xle)),
        # tail geometry
        tail_b=0.4, tail_TR=0.7, tail_AR=6.0, tail_θroot=0.0, tail_θtip=0.0,
        tail_incidence=3pi/180, tail_position=[1.5,0.0,0.1], tail_cmac=tail_b/tail_AR,
        tail_xle=[tail_position[1],tail_position[1]],
        tail_yle=[0.0,tail_b/2],
        tail_zle=[tail_position[3],tail_position[3]],
        tail_chord=[tail_cmac*2/(1+tail_TR),tail_cmac*2/(1/tail_TR+1)],
        tail_twist=[tail_θroot, tail_θtip] .- tail_incidence,
        tail_phi=[0.0,0.0],
        tail_camber=fill((xc) -> 0, length(tail_xle)),
        # vertical tail geometry
        vertical_tail_b=0.4, vertical_tail_TR=0.7, vertical_tail_AR=6.0, vertical_tail_θroot=0.0, vertical_tail_θtip=0.0,
        vertical_tail_cmac=vertical_tail_b/vertical_tail_AR,
        vertical_tail_xle=[tail_position[1], tail_position[1]],
        vertical_tail_yle=[0.0,0.0],
        vertical_tail_zle=[tail_position[3], tail_position[3] + vertical_tail_b],
        vertical_tail_chord=[cmac_vertical_tail * 2 / (1 + vertical_tail_TR), cmac_vertical_tail * 2 / (1/vertical_tail_TR + 1)],
        vertical_tail_theta=[vertical_tail_θroot, vertical_tail_θtip],
        vertical_tail_phi=[0.0, 0.0],
        vertical_tail_camber=fill((xc) -> 0, length(vertical_tail_xle)),
        # other parameters/reference values
        Sref=wing_b*cmac,
        static_margin=0.10,
        Vinf=1.0,
        Vref=Vinf,
        alpha=0.0,
        beta=0.0,
        Omega=@SVector [0.0,0.0,0.0]
        symmetric=[true, true, false],
        mirror=[false, false, false],
        # panel setup
        wing_npanels=50, wing_nchordwisepanels=3,
        wing_spacing_s= mirror[1] ? VL.Sine() : VL.Cosine(), # if mirrored, Sine becomes Cosine
        wing_spacing_c=VL.Uniform(),
        tail_npanels=20, tail_nchordwisepanels=1,
        tail_spacing_s= mirror[2] ? VL.Sine() : VL.Cosine(),
        tail_spacing_c=VL.Uniform(),
        vertical_tail_npanels=20, vertical_tail_nchordwisepanels=1,
        vertical_tail_spacing_s= mirror[3] ? VL.Sine() : VL.Cosine(),
        vertical_tail_spacing_c=VL.Uniform(),
        kwargs...)

Convenience constructor for a simple airplane `VortexLatticeSystem`.

# Keyword Arguments:

* `wing_b`: wing span
* `wing_TR`: wing taper ratio
* `wing_AR`: wing aspect ratio
* `wing_θroot`: wing root twist
* `wing_θtip`: wing tip twist
* `wing_xle`: x-coordinates of wing leading edge
* `wing_yle`: y-coordinates of wing leading edge
* `wing_zle`: z-coordinates of wing leading edge
* `wing_cmac`: wing mean aerodynamic chord
* `wing_chord`: chord at each wing section
* `wing_theta`: twist (radians) at each wing section
* `wing_phi`: wing dihedral angle
* `wing_camber`: camberline function for each section

* `tail_b`: tail span
* `tail_TR`: tail taper ratio
* `tail_AR`: tail aspect ratio
* `tail_θroot`: tail root twist
* `tail_θtip`: tail tip twist
* `tail_incidence`: tail incidence angle (subtracts from twist)
* `tail_position`: position of the y=0 leading edge point of the tail
* `tail_cmac`: tail mean aerodynamic chord
* `tail_xle`: x-coordinates of tail leading edge
* `tail_yle`: y-coordinates of tail leading edge
* `tail_zle`: z-coordinates of tail leading edge
* `tail_chord`: chord at each tail section
* `tail_theta`: twist (radians) at each tail section
* `tail_phi`: tail dihedral angle
* `tail_camber`: camberline function for each section

* `vertical_tail_b`: vertical tail span
* `vertical_tail_TR`: vertical tail taper ratio
* `vertical_tail_AR`: vertical tail aspect ratio
* `vertical_tail_θroot`: vertical tail root twist
* `vertical_tail_θtip`: vertical tail tip twist
* `vertical_tail_cmac`: vertical tail mean aerodynamic chord
* `vertical_tail_xle`: x-coordinates of vertical tail leading edge
* `vertical_tail_yle`: y-coordinates of vertical tail leading edge
* `vertical_tail_zle`: z-coordinates of vertical tail leading edge
* `vertical_tail_chord`: chord at each vertical tail section
* `vertical_tail_theta`: twist (radians) at each vertical tail section
* `vertical_tail_phi`: vertical tail dihedral angle
* `vertical_tail_camber`: camberline function for each section

* `Sref`: wing reference area
* `static_margin`: wing static margin
* `Vinf`: freestream velocity
* `Vref`: reference velocity; defaults to freestream velocity
* `alpha`: wing angle of attack (radians)
* `beta`: wing slipstream angle (radians)
* `Omega`: rotational velocity around the reference location
* `symmetric`: uses a mirror image when calculating induced velocities
* `mirror`: adds mirrored image to the geometry and solves without symmetry

* `wing_npanels`: number of spanwise panels for the wing
* `wing_nchordwisepanels`: number of chordwise panels for the wing
* `wing_spacing_s`: spanwise panel spacing scheme for the wing; defaults to Cosine()
* `wing_spacing_c`: chordwise panel spacing scheme for the wing; defaults to Uniform()
* `tail_npanels`: number of spanwise panels for the tail
* `tail_nchordwisepanels`: number of chordwise panels for the tail
* `tail_spacing_s`: spanwise panel spacing scheme for the tail; defaults to Cosine()
* `tail_spacing_c`: chordwise panel spacing scheme for the tail
* `vertical_tail_npanels`: number of spanwise panels for the vertical tail
* `vertical_tail_nchordwisepanels`: number of chordwise panels for the vertical tail
* `vertical_tail_spacing_s`: spanwise panel spacing scheme for the vertical tail
* `vertical_tail_spacing_c`: chordwise panel spacing scheme for the vertical tail

"""
function simpleairplanesystem(;
            # wing geometry
            wing_b=2.0, wing_TR=0.8, wing_AR=8.0, wing_θroot=0.0, wing_θtip=0.0,
            wing_xle=[0.0,0.0], wing_yle=[0.0, wing_b/2], wing_zle=[0.0,0.0],
            wing_cmac=wing_b/wing_AR, # AR = b/c => c = b/AR
            wing_chord=[wing_cmac * 2 / (1 + wing_TR), wing_cmac * 2 / (1/wing_TR + 1)],
            wing_theta=[wing_θroot, wing_θtip],
            wing_phi=[1.0,1.0] * 3pi/180.0,
            wing_camber=fill((xc) -> 0, length(wing_xle)),
            # tail geometry
            tail_b=0.4, tail_TR=0.7, tail_AR=6.0, tail_θroot=0.0, tail_θtip=0.0,
            tail_incidence=3pi/180, tail_position=[1.5,0.0,0.1], tail_cmac=tail_b/tail_AR,
            tail_xle=[tail_position[1],tail_position[1]],
            tail_yle=[0.0,tail_b/2],
            tail_zle=[tail_position[3],tail_position[3]],
            tail_chord=[tail_cmac*2/(1+tail_TR),tail_cmac*2/(1/tail_TR+1)],
            tail_twist=[tail_θroot, tail_θtip] .- tail_incidence,
            tail_phi=[0.0,0.0],
            tail_camber=fill((xc) -> 0, length(tail_xle)),
            # vertical tail geometry
            vertical_tail_b=0.4, vertical_tail_TR=0.7, vertical_tail_AR=6.0, vertical_tail_θroot=0.0, vertical_tail_θtip=0.0,
            vertical_tail_cmac=vertical_tail_b/vertical_tail_AR,
            vertical_tail_xle=[tail_position[1], tail_position[1]],
            vertical_tail_yle=[0.0,0.0],
            vertical_tail_zle=[tail_position[3], tail_position[3] + vertical_tail_b],
            vertical_tail_chord=[cmac_vertical_tail * 2 / (1 + vertical_tail_TR), cmac_vertical_tail * 2 / (1/vertical_tail_TR + 1)],
            vertical_tail_theta=[vertical_tail_θroot, vertical_tail_θtip],
            vertical_tail_phi=[0.0, 0.0],
            vertical_tail_camber=fill((xc) -> 0, length(vertical_tail_xle)),
            # other parameters/reference values
            Sref=wing_b*cmac,
            static_margin=0.10,
            Vinf=1.0,
            Vref=Vinf,
            alpha=0.0,
            beta=0.0,
            Omega=(StaticArrays.@SVector [0.0,0.0,0.0]),
            symmetric=[true, true, false],
            mirror=[false, false, false],
            # panel setup
            wing_npanels=50, wing_nchordwisepanels=3,
            wing_spacing_s= mirror[1] ? VL.Sine() : VL.Cosine(), # if mirrored, Sine becomes Cosine
            wing_spacing_c=VL.Uniform(),
            tail_npanels=20, tail_nchordwisepanels=1,
            tail_spacing_s= mirror[2] ? VL.Sine() : VL.Cosine(),
            tail_spacing_c=VL.Uniform(),
            vertical_tail_npanels=20, vertical_tail_nchordwisepanels=1,
            vertical_tail_spacing_s= mirror[3] ? VL.Sine() : VL.Cosine(),
            vertical_tail_spacing_c=VL.Uniform(),
            kwargs...)

    # Check for errors
    @assert size(wing_xle) == size(wing_yle)    "size of leading edge xs and ys not consistent for the wing"
    @assert size(wing_yle) == size(wing_zle)    "size of leading edge ys and zs not consistent for the wing"
    @assert size(wing_xle) == size(wing_chord)    "size of leading edge xs and wing_chord not consistent"
    @assert size(wing_xle) == size(wing_twist)    "size of leading edge xs and wing_twist not consistent"
    @assert size(wing_xle) == size(wing_phi)    "size of leading edge xs and wing_phi not consistent"
    @assert wing_yle[1] == 0.0    "wing yle must begin at the root (0.0)"
    @assert wing_yle[end] == wing_b/2    "wing yle must end at the tip (span/2)"

    @assert size(tail_xle) == size(tail_yle)    "size of leading edge xs and ys not consistent for the tail"
    @assert size(tail_yle) == size(tail_zle)    "size of leading edge ys and zs not consistent for the tail"
    @assert size(tail_xle) == size(tail_chord)    "size of leading edge xs and tail_chord not consistent"
    @assert size(tail_xle) == size(tail_twist)    "size of leading edge xs and tail_twist not consistent"
    @assert size(tail_xle) == size(tail_phi)    "size of leading edge xs and tail_phi not consistent"
    @assert tail_yle[1] == 0.0    "tail yle must begin at the root (0.0)"
    @assert tail_yle[end] == tail_b/2    "tail yle must end at the tip (span/2)"

    @assert size(vertical_tail_xle) == size(vertical_tail_yle)    "size of leading edge xs and ys not consistent for the vertical tail"
    @assert size(vertical_tail_yle) == size(vertical_tail_zle)    "size of leading edge ys and zs not consistent for the vertical tail"
    @assert size(vertical_tail_xle) == size(vertical_tail_chord)    "size of leading edge xs and vertical_tail_chord not consistent"
    @assert size(vertical_tail_xle) == size(vertical_tail_twist)    "size of leading edge xs and vertical_tail_twist not consistent"
    @assert size(vertical_tail_xle) == size(vertical_tail_phi)    "size of leading edge xs and vertical_tail_phi not consistent"
    @assert vertical_tail_yle[1] == 0.0    "vertical tail yle must begin at the root (0.0)"

    @assert !(symmetric[1] && mirror[1])    "symmetric and mirror cannot both be true for the wing"
    @assert !(symmetric[1] && mirror[2])    "symmetric and mirror cannot both be true for the tail"
    @assert !(symmetric[1] && mirror[3])    "symmetric and mirror cannot both be true for the vertical tail"

    # construct wing surfaces and grids
    wing_grid, wing_surface = VL.wing_to_surface_panels(wing_xle, wing_yle, wing_zle, wing_chord, wing_theta, wing_phi, wing_npanels, wing_nchordwisepanels;
                                    fc = wing_camber, spacing_s=wing_spacing_s, spacing_c=wing_spacing_c, mirror=mirror[1])

    # construct tail surfaces and grids
    tail_grid, tail_surface = VL.wing_to_surface_panels(tail_xle, tail_yle, tail_zle, tail_chord, tail_theta, tail_phi, tail_npanels, tail_nchordwisepanels;
                                    fc = tail_camber, spacing_s=tail_spacing_s, spacing_c=tail_spacing_c, mirror=mirror[2])

    # construct vertical tail surfaces and grids
    vertical_tail_grid, vertical_tail_surface = VL.wing_to_surface_panels(vertical_tail_xle, vertical_tail_yle, vertical_tail_zle, vertical_tail_chord,
                                                    vertical_tail_theta, vertical_tail_phi, vertical_tail_npanels, vertical_tail_nchordwisepanels;
                                                    fc = vertical_tail_camber, spacing_s=vertical_tail_spacing_s, spacing_c=vertical_tail_spacing_c, mirror=mirror[3])

    # assemble system
    grids = [wing_grid, tail_grid, vertical_tail_grid]
    surfaces = [wing_surface, tail_surface, vertical_tail_surface]

    # reference parameters
    cgx = cmac * (0.25 - static_margin)
    rref = [cgx, 0.0, 0.0]
    reference = VL.Reference(Sref, cmac, wing_b, rref, Vref)

    # initialize freestream parameters
    fs = VL.Freestream(Vinf, alpha, beta, Omega)

    # perform steady state analysis
    vlmsystem = VL.steady_analysis(surfaces, reference, fs; symmetric)
    wing_system = VortexLatticeSystem(vlmsystem, grids)

    return wing_system
end

function VL.Freestream(freestream::Freestream)
    VL.Freestream(freestream.vinf, freestream.alpha, freestream.beta, freestream.Omega)
end

# outputs = [Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Outputs(0.24285036855335476 - 0.0im, 0.09374476338448777 - 0.0im, 0.0, 0.0, 1.2308487178526557, 0.0, 0.2772045889709332, 0.3183884946733561, 4.534990314410568 - 0.0im, 1.5050759240380058, 0.13763120546656513, 1.4099532010335012, 0.5442681845675164, 0.9835773158795993, 0.9917546651665418), Outputs(0.7217628407058008 - 0.0im, 0.21697762721530664 - 0.0im, 0.0, 0.0, 1.7326960734792525, 0.0, 0.2589657900820015, 0.245706809640596, 6.7673827545271354 - 0.0im, 1.6376114138887694, 0.05416136893272035, 1.569136116792526, 0.47171648663229143, 0.9996516983041511, 0.9998258339851752), Outputs(1.2990599048174767 - 0.0im, 0.31355858702643263 - 0.0im, 0.0, 0.0, 2.0129587976546834, 0.0, 0.22688706695390065, 0.1431274468578776, 8.953184950774602 - 0.0im, 1.3892205399175588, 0.01383177005724867, 1.3505053059805476, 0.32597614162717065, 0.9989877735121309, 0.999493758615896), Outputs(1.7619807902797593 - 0.0im, 0.3720701211173133 - 0.0im, 0.0, 0.0, 2.09694555206726, 0.0, 0.19036789138116791, 0.07972442835883256, 11.104190789025736 - 0.0im, 1.09077738226134, 0.019353098469995278, 1.067410153478332, 0.2254005419795043, 0.9960145868560558, 0.9980053040220056), Outputs(2.2711845655764025 - 0.0im, 0.4405002858849476 - 0.0im, 0.0, 0.0, 2.173383411501406, 0.0, 0.1657314413795813, 0.0556670648944034, 13.265213925277973 - 0.0im, 0.9558687456440128, 0.02470688289088402, 0.9386953965963786, 0.18206164167668523, 0.9863118571582342, 0.9931323462450683), Outputs(2.4372150640530514 - 0.0im, 0.4554901469388099 - 0.0im, 0.0, 0.0, 2.1629611788994447, 0.0, 0.15314916243164545, 0.04256438043828989, 14.343403366255124 - 0.0im, 0.8800488657666018, 0.02782677887567497, 0.8655033880045069, 0.16175358145153387, 0.9771513649187485, 0.9885096686015512), Outputs(2.6371721454533454 - 0.0im, 0.47293752842306824 - 0.0im, 0.0, 0.0, 2.1681319604225773, 0.0, 0.14389224181868238, 0.0387234776796615, 15.424864001541529 - 0.0im, 0.8574901315613467, 0.028785248170534933, 0.8445005855017522, 0.15144859631086924, 0.9608451485828865, 0.9802270903126921), Outputs(2.9035116834001253 - 0.0im, 0.5070499699969895 - 0.0im, 0.0, 0.0, 2.1978736826245666, 0.0, 0.13838725719185638, 0.03708290981589751, 16.51325595934029 - 0.0im, 0.8476745782368323, 0.029258614602765128, 0.8355344940565491, 0.14591218715080329, 0.9309406549944107, 0.9648526597332935), Outputs(2.915378995839232 - 0.0im, 0.4988635277901568 - 0.0im, 0.0, 0.0, 2.132452708206054, 0.0, 0.12893626553987916, 0.029308601308119464, 17.59133845621036 - 0.0im, 0.7995946884337315, 0.03243082128604525, 0.7887875033803128, 0.13497295451969854, 0.888830318342529, 0.9427779793474862), Outputs(2.672487776120149 - 0.0im, 0.4658323109725899 - 0.0im, 0.0, 0.0, 1.9807514705340548, 0.0, 0.11703583739096063, 0.014783645024415093, 18.662919158779463 - 0.0im, 0.7077670820565516, 0.039347704185111965, 0.6983307468035365, 0.121723672045743, 0.8261271529063575, 0.908915371696594), Outputs(2.6753732822512957 - 0.0im, 0.4714559098077924 - 0.0im, 0.0, 0.0, 1.925896978460363, 0.0, 0.11765311844886973, 0.013955301300062325, 19.763948184544113 - 0.0im, 0.702698187395941, 0.039939768326604105, 0.6931521258065394, 0.1221476899224832, 0.689154508910361, 0.8301533044627125), Outputs(2.397999681058518 - 0.0im, 0.4549131149780635 - 0.0im, 0.0, 0.0, 1.774720023483803, 0.0, 0.1305556999113481, 0.01385081781673042, 20.895071940836615 - 0.0im, 0.7020550192002749, 0.04000541332594471, 0.6908722123615864, 0.1310620817257423, 0.4256466033071096, 0.652415974135451), Outputs(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2} where {TF1,TF2}}, sym::Symbol)
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2}} where TF1, sym::Symbol) where TF2
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2}} where TF2, sym::Symbol) where TF1
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,Float64} where TF1}, sym::Symbol)
#     return getfield.(obj, sym)
# end

# outputs2 = Outputs{Float64,Float64}[Outputs{Float64,Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Outputs{Float64,Float64}(0.24285036855335476, 0.09374476338448777, 0.0, 0.0, 1.2308487178526557, 0.0, 0.2772045889709332, 0.3183884946733561, 4.534990314410568, 1.5050759240380058, 0.13763120546656513, 1.4099532010335012, 0.5442681845675164, 0.9835773158795993, 0.9917546651665418), Outputs{Float64,Float64}(0.7217628407058008, 0.21697762721530664, 0.0, 0.0, 1.7326960734792525, 0.0, 0.2589657900820015, 0.245706809640596, 6.7673827545271354, 1.6376114138887694, 0.05416136893272035, 1.569136116792526, 0.47171648663229143, 0.9996516983041511, 0.9998258339851752), Outputs{Float64,Float64}(1.2990599048174767, 0.31355858702643263, 0.0, 0.0, 2.0129587976546834, 0.0, 0.22688706695390065, 0.1431274468578776, 8.953184950774602, 1.3892205399175588, 0.01383177005724867, 1.3505053059805476, 0.32597614162717065, 0.9989877735121309, 0.999493758615896), Outputs{Float64,Float64}(1.7619807902797593, 0.3720701211173133, 0.0, 0.0, 2.09694555206726, 0.0, 0.19036789138116791, 0.07972442835883256, 11.104190789025736, 1.09077738226134, 0.019353098469995278, 1.067410153478332, 0.2254005419795043, 0.9960145868560558, 0.9980053040220056), Outputs{Float64,Float64}(2.2711845655764025, 0.4405002858849476, 0.0, 0.0, 2.173383411501406, 0.0, 0.1657314413795813, 0.0556670648944034, 13.265213925277973, 0.9558687456440128, 0.02470688289088402, 0.9386953965963786, 0.18206164167668523, 0.9863118571582342, 0.9931323462450683), Outputs{Float64,Float64}(2.4372150640530514, 0.4554901469388099, 0.0, 0.0, 2.1629611788994447, 0.0, 0.15314916243164545, 0.04256438043828989, 14.343403366255124, 0.8800488657666018, 0.02782677887567497, 0.8655033880045069, 0.16175358145153387, 0.9771513649187485, 0.9885096686015512), Outputs{Float64,Float64}(2.6371721454533454, 0.47293752842306824, 0.0, 0.0, 2.1681319604225773, 0.0, 0.14389224181868238, 0.0387234776796615, 15.424864001541529, 0.8574901315613467, 0.028785248170534933, 0.8445005855017522, 0.15144859631086924, 0.9608451485828865, 0.9802270903126921), Outputs{Float64,Float64}(2.9035116834001253, 0.5070499699969895, 0.0, 0.0, 2.1978736826245666, 0.0, 0.13838725719185638, 0.03708290981589751, 16.51325595934029, 0.8476745782368323, 0.029258614602765128, 0.8355344940565491, 0.14591218715080329, 0.9309406549944107, 0.9648526597332935), Outputs{Float64,Float64}(2.915378995839232, 0.4988635277901568, 0.0, 0.0, 2.132452708206054, 0.0, 0.12893626553987916, 0.029308601308119464, 17.59133845621036, 0.7995946884337315, 0.03243082128604525, 0.7887875033803128, 0.13497295451969854, 0.888830318342529, 0.9427779793474862), Outputs{Float64,Float64}(2.672487776120149, 0.4658323109725899, 0.0, 0.0, 1.9807514705340548, 0.0, 0.11703583739096063, 0.014783645024415093, 18.662919158779463, 0.7077670820565516, 0.039347704185111965, 0.6983307468035365, 0.121723672045743, 0.8261271529063575, 0.908915371696594), Outputs{Float64,Float64}(2.6753732822512957, 0.4714559098077924, 0.0, 0.0, 1.925896978460363, 0.0, 0.11765311844886973, 0.013955301300062325, 19.763948184544113, 0.702698187395941, 0.039939768326604105, 0.6931521258065394, 0.1221476899224832, 0.689154508910361, 0.8301533044627125), Outputs{Float64,Float64}(2.397999681058518, 0.4549131149780635, 0.0, 0.0, 1.774720023483803, 0.0, 0.1305556999113481, 0.01385081781673042, 20.895071940836615, 0.7020550192002749, 0.04000541332594471, 0.6908722123615864, 0.1310620817257423, 0.4256466033071096, 0.652415974135451), Outputs{Float64,Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]

# struct Outputs{TF1,TF2}
#     Np::TF1
#     Tp::TF1
#     a::TF2
#     ap::TF2
#     u::TF2
#     v::TF2
#     phi::TF2
#     alpha::TF2
#     W::TF1
#     cl::TF2
#     cd::TF2
#     cn::TF2
#     ct::TF2
#     F::TF2
#     G::TF2
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2} where {TF1,TF2}}, sym::Symbol)
#     println("Sherlock! Method 1")
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2}}, sym::Symbol) where {TF1,TF2}
#     println("Sherlock! Method 2")
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2} where TF2}, sym::Symbol) where TF1
#     println("Sherlock! Method 3")
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs{TF1,TF2} where TF1}, sym::Symbol) where TF2
#     println("Sherlock! Method 4")
#     return getfield.(obj, sym)
# end

# function Base.getproperty(obj::Vector{Outputs}, sym::Symbol)
#     println("Sherlock! Using catchall method")
#     return getfield.(obj, sym)
# end

# outputs = Outputs[Outputs{Float64,Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0), Outputs{Complex{Float64},Float64}(0.24285036855335476 - 0.0im, 0.09374476338448777 - 0.0im, 0.0, 0.0, 1.2308487178526557, 0.0, 0.2772045889709332, 0.3183884946733561, 4.534990314410568 - 0.0im, 1.5050759240380058, 0.13763120546656513, 1.4099532010335012, 0.5442681845675164, 0.9835773158795993, 0.9917546651665418), Outputs{Complex{Float64},Float64}(0.7217628407058008 - 0.0im, 0.21697762721530664 - 0.0im, 0.0, 0.0, 1.7326960734792525, 0.0, 0.2589657900820015, 0.245706809640596, 6.7673827545271354 - 0.0im, 1.6376114138887694, 0.05416136893272035, 1.569136116792526, 0.47171648663229143, 0.9996516983041511, 0.9998258339851752), Outputs{Complex{Float64},Float64}(1.2990599048174767 - 0.0im, 0.31355858702643263 - 0.0im, 0.0, 0.0, 2.0129587976546834, 0.0, 0.22688706695390065, 0.1431274468578776, 8.953184950774602 - 0.0im, 1.3892205399175588, 0.01383177005724867, 1.3505053059805476, 0.32597614162717065, 0.9989877735121309, 0.999493758615896), Outputs{Complex{Float64},Float64}(1.7619807902797593 - 0.0im, 0.3720701211173133 - 0.0im, 0.0, 0.0, 2.09694555206726, 0.0, 0.19036789138116791, 0.07972442835883256, 11.104190789025736 - 0.0im, 1.09077738226134, 0.019353098469995278, 1.067410153478332, 0.2254005419795043, 0.9960145868560558, 0.9980053040220056), Outputs{Complex{Float64},Float64}(2.2711845655764025 - 0.0im, 0.4405002858849476 - 0.0im, 0.0, 0.0, 2.173383411501406, 0.0, 0.1657314413795813, 0.0556670648944034, 13.265213925277973 - 0.0im, 0.9558687456440128, 0.02470688289088402, 0.9386953965963786, 0.18206164167668523, 0.9863118571582342, 0.9931323462450683), Outputs{Complex{Float64},Float64}(2.4372150640530514 - 0.0im, 0.4554901469388099 - 0.0im, 0.0, 0.0, 2.1629611788994447, 0.0, 0.15314916243164545, 0.04256438043828989, 14.343403366255124 - 0.0im, 0.8800488657666018, 0.02782677887567497, 0.8655033880045069, 0.16175358145153387, 0.9771513649187485, 0.9885096686015512), Outputs{Complex{Float64},Float64}(2.6371721454533454 - 0.0im, 0.47293752842306824 - 0.0im, 0.0, 0.0, 2.1681319604225773, 0.0, 0.14389224181868238, 0.0387234776796615, 15.424864001541529 - 0.0im, 0.8574901315613467, 0.028785248170534933, 0.8445005855017522, 0.15144859631086924, 0.9608451485828865, 0.9802270903126921), Outputs{Complex{Float64},Float64}(2.9035116834001253 - 0.0im, 0.5070499699969895 - 0.0im, 0.0, 0.0, 2.1978736826245666, 0.0, 0.13838725719185638, 0.03708290981589751, 16.51325595934029 - 0.0im, 0.8476745782368323, 0.029258614602765128, 0.8355344940565491, 0.14591218715080329, 0.9309406549944107, 0.9648526597332935), Outputs{Complex{Float64},Float64}(2.915378995839232 - 0.0im, 0.4988635277901568 - 0.0im, 0.0, 0.0, 2.132452708206054, 0.0, 0.12893626553987916, 0.029308601308119464, 17.59133845621036 - 0.0im, 0.7995946884337315, 0.03243082128604525, 0.7887875033803128, 0.13497295451969854, 0.888830318342529, 0.9427779793474862), Outputs{Complex{Float64},Float64}(2.672487776120149 - 0.0im, 0.4658323109725899 - 0.0im, 0.0, 0.0, 1.9807514705340548, 0.0, 0.11703583739096063, 0.014783645024415093, 18.662919158779463 - 0.0im, 0.7077670820565516, 0.039347704185111965, 0.6983307468035365, 0.121723672045743, 0.8261271529063575, 0.908915371696594), Outputs{Complex{Float64},Float64}(2.6753732822512957 - 0.0im, 0.4714559098077924 - 0.0im, 0.0, 0.0, 1.925896978460363, 0.0, 0.11765311844886973, 0.013955301300062325, 19.763948184544113 - 0.0im, 0.702698187395941, 0.039939768326604105, 0.6931521258065394, 0.1221476899224832, 0.689154508910361, 0.8301533044627125), Outputs{Complex{Float64},Float64}(2.397999681058518 - 0.0im, 0.4549131149780635 - 0.0im, 0.0, 0.0, 1.774720023483803, 0.0, 0.1305556999113481, 0.01385081781673042, 20.895071940836615 - 0.0im, 0.7020550192002749, 0.04000541332594471, 0.6908722123615864, 0.1310620817257423, 0.4256466033071096, 0.652415974135451), Outputs{Float64,Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
module PressureFieldContact

using DocStringExtensions
using StaticArrays
using Rotations: Quat, Rotation, SPQuat, RotMatrix
using ForwardDiff: Dual
using RigidBodyDynamics
using RigidBodyDynamics.Spatial: vector_to_skew_symmetric, vector_to_skew_symmetric_squared
using GeometryTypes: HomogenousMesh, Face, Point
using ColorTypes: RGBA
using MeshCatMechanisms
using DocStringExtensions
using CoordinateTransformations: Translation

include(joinpath("math_kernel", "MathKernel.jl"))
using .MathKernel

include(joinpath("clip", "Clip.jl"))
using .Clip

include(joinpath("obb", "Binary_BB_Trees.jl"))
using .Binary_BB_Trees

include(joinpath("radau", "Radau.jl"))
using .Radau

include(joinpath("geometry", "Geometry.jl"))
using .Geometry

using LinearAlgebra
using GenericLinearAlgebra

include("structs.jl")
include("body_inertia.jl")
# include("matrix_transform.jl")
include("mechanism_scenario.jl")
include("extensions.jl")
include("contact_algorithms_non_friction.jl")
include("contact_algorithms_friction.jl")
include("contact_algorithms_normal.jl")
include("vis_meshcat.jl")
include("example_integrator.jl")
include("utility.jl")

export
    # structs.jl
    MeshInertiaInfo,
    ContactProperties,
    # eTree,
    InertiaProperties,
    get_c_prop,
    get_ind_tri,
    get_ind_tet,
    get_ϵ,
    get_Ē,
    MeshID,
    MeshDict,
    MeshCacheDict,
    MeshCache,

    # mesh_inertia.jl
    makeInertiaInfo,
    makeInertiaTensor,
    centroidVolumeCombo,
    equiv_volume,

    # mesh_body_utility.jl
    newBodyFromInertia,
    outputJointTransform_ParentChild,

    # mechanism_scenario
    BristleID,
    Bristle,
    Regularized,
    ContactInstructions,
    TractionCache,
    calc_p_dA,
    spatialStiffness,
    TypedElasticBodyBodyCache,
    TypedMechanismScenario,
    DiscreteControl,
    MechanismScenario,
    finalize!,
    get_mechanism,
    num_partials,
    num_x,
    type_dual,
    get_state,
    set_state_spq!,
    addMesh!,
    add_body_contact!,
    add_contact!,
    add_body!,
    add_body_from_inertia!,  # this needs to be here
    add_friction_regularize!,
    add_friction_bristle!,

    # extensions.jl
    principal_value!,

    # contact_algorithms_non_friction.jl
    calcXd!,
    calcXd,
    refreshJacobians!,
    normal_wrench,
    forceAllElasticIntersections!,
    calcTriTetIntersections!,
    refreshBodyBodyCache!,
    addGeneralizedForcesThirdLaw!,

    # contact_algorithms_friction.jl

    # contact_algorithms_normal.jl
    normal_wrench,

    # vis_meshcat.jl
    set_body_mesh_visual!,
    set_mesh_visual!,
    HomogenousMesh_32,
    play_recorded_data,

    # example_integrator.jl
    integrate_scenario_radau,

    # utility.jl
    num_partials,
    type_dual,
    fill_with_nothing!,
    Radau_for_MechanismScenario,
    get_bristle_d0,
    get_bristle_d1,
    as_static_vector,
    find_mesh,
    find_mesh_id,

    # obb/Binary_BB_Trees.jl
    transform!,
    eMesh_half_plane,
    eMesh_sphere,
    eMesh_box,
    eMesh_cylinder,
    as_tet_eMesh,
    as_tri_eMesh,
    n_tri,
    n_tet,
    n_point,
    get_tri,
    get_tet,
    get_point,
    create_swept_mesh,
    f_swept_triv,
    f_swept_circle,

    # math_kernel/MathKernel.jl
    make_pd_gains,
    basic_dh

end