text:
module PiCraft

using Sockets
using GZip

include("blocks.jl")
include("turtle.jl")
include("drawing.jl")
include("NBT.jl")
include("buildTools.jl")

export World, Block, connectToWorld, mc_send, getBlock, setBlock, setBlocks, getHeight, getPlayerIds
export setting, saveWorld, restoreWorld, post, getChatEvents, getTile, setTile, getPos, setPos, pollBlockHits
export clearEvents, camera
export turtle, move, yaw, pitch, roll
export drawLine
export parseNBT, importSchematic, exportSchematic
export buildModel, copyModel, cutModel, flip, rotate, rotate!

mutable struct World
    s::Sockets.TCPSocket
    function World()
        try
            s = Sockets.connect("localhost", 4711)
            new(s)
        catch
            @warn("Unable to connect to Minecraft World. Use `connectToWorld()` after starting Minecraft")
            new(TCPSocket())
        end
    end
end

global const minecraftWorld = Ref{World}()

"""
    connectToWorld(address = "localhost", port = 4711)

Connect to the Minecraft API.
"""
function connectToWorld(address = "localhost", port = 4711)
    try
        minecraftWorld[].s = connect(address, port)
    catch
        error("Unable to connect to the Minecraft World.")
    end
    @info("Successfully connected to the Minecraft World.")
end

"""
    mc_send(cmd, output=true)

Communicate with the Minecraft API.
"""
function mc_send(cmd, output=true)
    if minecraftWorld[].s.status == Base.StatusInit || minecraftWorld[].s.status == Base.StatusUninit
        error("Connection to Minecraft World is not initialised. Use `PiCraft.connectToWorld()` first.")
    end
    write(minecraftWorld[].s, cmd*"\n")
    if output
        s = readline(minecraftWorld[].s)
        if occursin("Fail", s)
            println("Return: $s")
            error("Error from Minecraft in command $cmd")
        end
        return split(strip(s), ',')
    end
end

"""
    getBlock(pos::Tuple{Real, Real, Real})

Get the Block information from the specified coordinates.
"""
function getBlock(pos::Tuple{Real, Real, Real})
    pos = round.(Int, pos)
    res = PiCraft.mc_send("world.getBlockWithData($(pos[1]),$(pos[2]),$(pos[3]))", true)
    return Block(parse(Int, res[1]), (parse(Int, res[2])))
end

"""
    setBlock(pos::Tuple{Real, Real, Real}, block::Block)

Place the specified `Block` at the given coordinates.
"""
function setBlock(pos::Tuple{Real, Real, Real}, block::Block)
    pos = round.(Int, pos)
    PiCraft.mc_send("world.setBlock($(pos[1]),$(pos[2]),$(pos[3]),$(block.id),$(block.data))", false)
end

"""
    setBlocks(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real}, block::Block)

Set an entire region to the specified block type defined by the corners `p1` and `p2`.
"""
function setBlocks(p1::Tuple{Real, Real, Real}, p2::Tuple{Real, Real, Real}, block::Block)
    p1 = round.(Int, p1)
    p2 = round.(Int, p2)
    PiCraft.mc_send("world.setBlocks($(p1[1]),$(p1[2]),$(p1[3]),$(p2[1]),$(p2[2]),$(p2[3]),$(block.id),$(block.data))", false)
end

"""
    getHeight(x::Int, z::Int)

Get the height of the world at the specified `x` and `z` coordinates.
"""
function getHeight(x::Int, z::Int)
    h = PiCraft.mc_send("world.getHeight($x,$z)", true)
    return parse(Int, h[1])
end

"""
    getPlayerIds()

Return an array of all Player Id's connected to the server.
"""
function getPlayerIds()
    return parse.(Int,PiCraft.mc_send("world.getPlayerIds()"))
end

setting(name, status) = PiCraft.mc_send("world.setting($(name),$(status))", false)

"Save the World's progress"
saveWorld() = PiCraft.mc_send("world.checkpoint.save()", false)

"Restore the world to the previous savepoint."
restoreWorld() = PiCraft.mc_send("world.checkpoint.restore()", false)

"""
    post(m...)

Post a message to chat
"""
function post(m...)
    s = ""
    for x in m
        s *= string(x)
    end
    PiCraft.mc_send("chat.post($(s)", false)
end

function getChatEvents()
    ChatEventsRaw = []
    for j in PiCraft.mc_send("events.chat.posts()", true)
        for s in split(j, "|")
            length(s) > 0 && push!(ChatEventsRaw, s)
        end
    end
    ChatEvents = Array{Tuple{Int, String}, 1}(undef, Int(length(ChatEventsRaw)/2))
    for i in 1:length(ChatEvents)
        ChatEvents[i] = (parse(Int, ChatEventsRaw[2*i - 1]), ChatEventsRaw[2*i])
    end
    return ChatEvents
end

"""
    getTile()
Return the tile's coordinates on which the player is standing.
"""
function getTile()
    pos = PiCraft.mc_send("player.getTile()", true)
    return tuple(parse.(Int,pos)...)
end

"""
    setTile(pos::Tuple{Real, Real, Real})

Teleport the player on top of the specified tile
"""
function setTile(pos::Tuple{Real, Real, Real})
    pos = round.(Int, pos)
    PiCraft.mc_send("player.setTile($(pos[1]),$(pos[2]),$(pos[3]))", false)
end

"""
    getPos()

Return the player's coordinates.
"""
function getPos()
    pos = PiCraft.mc_send("player.getPos()", true)
    return tuple(parse.(Float64, pos)...)
end

"""
    setPos(pos::Tuple{Real, Real, Real})

Teleport the player to the specified coordinates
"""
function setPos(pos::Tuple{Real, Real, Real})
    PiCraft.mc_send("player.setPos($(pos[1]),$(pos[2]),$(pos[3]))", false)
end

#Entity Commands
"""
    getTile(entityId::Int)

Get the tile on which the specified `entity` is.
"""
function getTile(entityId::Int)
    pos = PiCraft.mc_send("player.getTile($(entityId))", true)
    return tuple(parse.(Int, pos)...)
end

"""
    setTile(entityId::Int)

Teleport the `entity` on the specified tile.
"""
function setTile(entityId::Int, pos::Tuple{Real, Real, Real})
    pos = round.(Int,pos)
    PiCraft.mc_send("player.setTile($(entityId),$(pos[1]),$(pos[2]),$(pos[3]))", false)
end

"""
    getPos(entityId::Int)

Get the position of the specified `entity`.
"""
function getPos(entityId::Int)
    pos = PiCraft.mc_send("player.getPos($(entityId))", true)
    return tuple(parse.(Float64, pos)...)
end

"""
    setPos(entityId::Int, pos::Tuple{Real, Real, Real})

Set the position of the specified `entity`.
"""
function setPos(entityId::Int, pos::Tuple{Real, Real, Real})
    PiCraft.mc_send("player.setPos($(entityId),$(pos[1]),$(pos[2]),$(pos[3]))", false)
end

"""
    pollBlockHits()

Returns an array of all the events which have occurred since the last time the function was called.
Each event is described with a tuple `((x, y, z), face, entityId)`. `x`, `y` and `z` are the coordinates of the block.
`face` is the block's face number which was hit and `entityId` identifies the player who hit the block using a sword.
"""
function pollBlockHits()
    rawEvents = []
    for s in PiCraft.mc_send("events.block.hits()", true)
        for j in split(s, "|")
            length(j) > 0 && push!(rawEvents, parse(Int,j))
        end
    end
    Events = []
    for i in 0:(Int(length(rawEvents)/5) - 1)
        push!(Events, ( (rawEvents[1 + 5*i], rawEvents[2 + 5*i], rawEvents[3 + 5*i]) , rawEvents[4 + 5*i], rawEvents[5 + 5*i]))
    end
    return Events
end

"Clear all events in the buffer"
clearEvents() = PiCraft.mc_send("events.clear()", false)

# Camera functions, Only supported by Minecraft: Pi Edition
module camera
    using PiCraft

    setNormal(entityId::Int) = PiCraft.mc_send("camera.mode.setNormal($(entityId))", false)
    setFixed() = PiCraft.mc_send("camera.mode.setFixed()", false)
    setFollow(entityId::Int) = PiCraft.mc_send("camera.mode.setFollow($(entityId))", false)
    setPos(pos::Tuple{Real, Real, Real}) = PiCraft.mc_send("camera.setPos($(pos[1]),$(pos[2]),$(pos[3]))", false)

end

end # module


text:
using BlockDiagonalMatrices
using LinearAlgebra: Diagonal, I
using Random
using Test

@testset "base_maths.jl" begin
    rng = MersenneTwister(123456)
    N1, N2, N3 = 3, 4, 5
    N = N1 + N2 + N3
    b1 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N3, N3)])
    b2 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N3, N3), rand(rng, N2, N2)])
    b3 = BlockDiagonal([rand(rng, N1, N1), rand(rng, N2, N2), rand(rng, N2, N2)])
    A = rand(rng, N, N + N1)
    B = rand(rng, N + N1, N + N2)
    A′, B′ = A', B'
    a = rand(rng, N)
    b = rand(rng, N + N1)

    @testset "Addition" begin
        @testset "BlockDiagonal + BlockDiagonal" begin
            @test b1 + b1 isa BlockDiagonal
            @test Matrix(b1 + b1) == Matrix(b1) + Matrix(b1)
            @test_throws DimensionMismatch b1 + b3
        end

        @testset "BlockDiagonal + Matrix" begin
            @test b1 + Matrix(b1) isa Matrix
            @test b1 + Matrix(b1) == b1 + b1
            @test_throws DimensionMismatch b1 + Matrix(b3)

            # Matrix + BlockDiagonal
            @test Matrix(b1) + b1 isa Matrix
            @test Matrix(b1) + b1 == b1 + b1
            @test_throws DimensionMismatch Matrix(b1) + b3

            # If the AbstractMatrix is diagonal, we should return a BlockDiagonal.
            # Test the StridedMatrix method.
            @test Matrix(Diagonal(b1)) + b1 isa BlockDiagonal  # StridedMatrix
            @test b1 + Matrix(Diagonal(b1)) isa BlockDiagonal

            # Test the AbstractMatrix method.
            @test Matrix(Diagonal(b1))' + b1 isa BlockDiagonal
            @test b1 + Matrix(Diagonal(b1))' isa BlockDiagonal
        end

        @testset "BlockDiagonal + Diagonal" begin
            D = Diagonal(randn(rng, N))
            D′ = Diagonal(randn(rng, N + N1))

            @test b1 + D isa BlockDiagonal
            @test b1 + D == Matrix(b1) + D
            @test_throws DimensionMismatch b1 + D′

            # Diagonal + BlockDiagonal
            @test D + b1 isa BlockDiagonal
            @test D + b1 == D + Matrix(b1)
            @test_throws DimensionMismatch D′ + b1
        end

        @testset "BlockDiagonal + UniformScaling" begin
            @test b1 + 5I isa BlockDiagonal
            @test b1 + 5I == Matrix(b1) + 5I

            # UniformScaling + BlockDiagonal
            @test 5I + b1 isa BlockDiagonal
            @test 5I + b1 == 5I + Matrix(b1)
        end
    end  # Addition

    @testset "Subtraction" begin
        @testset "BlockDiagonal - BlockDiagonal" begin
            @test b1 - b1 isa BlockDiagonal
            @test Matrix(b1 - b1) == Matrix(b1) - Matrix(b1)
            @test_throws DimensionMismatch b1 - b3
        end

        @testset "BlockDiagonal - Matrix" begin
            @test b1 - Matrix(b1) isa Matrix
            @test b1 - Matrix(b1) == b1 - b1
            @test_throws DimensionMismatch b1 - Matrix(b3)

            # Matrix - BlockDiagonal
            @test Matrix(b1) - b1 isa Matrix
            @test Matrix(b1) - b1 == b1 - b1
            @test_throws DimensionMismatch Matrix(b1) - b3

            # If the AbstractMatrix is diagonal, we should return a BlockDiagonal.
            # Test the StridedMatrix method.
            @test Matrix(Diagonal(b1)) - b1 isa BlockDiagonal  # StridedMatrix
            @test b1 - Matrix(Diagonal(b1)) isa BlockDiagonal

            # Test the AbstractMatrix method.
            @test Matrix(Diagonal(b1))' - b1 isa BlockDiagonal
            @test b1 - Matrix(Diagonal(b1))' isa BlockDiagonal
        end

        @testset "BlockDiagonal - Diagonal" begin
            D = Diagonal(randn(rng, N))
            D′ = Diagonal(randn(rng, N + N1))

            @test b1 - D isa BlockDiagonal
            @test b1 - D == Matrix(b1) - D
            @test_throws DimensionMismatch b1 - D′

            # Diagonal + BlockDiagonal
            @test D - b1 isa BlockDiagonal
            @test D - b1 == D - Matrix(b1)
            @test_throws DimensionMismatch D′ - b1
        end

        @testset "BlockDiagonal - UniformScaling" begin
            @test b1 - 5I isa BlockDiagonal
            @test b1 - 5I == Matrix(b1) - 5I

            # UniformScaling - BlockDiagonal
            @test 5I - b1 isa BlockDiagonal
            @test 5I - b1 == 5I - Matrix(b1)
        end
    end  # Subtraction

    @testset "Opposite" begin
        @test -b1 isa BlockDiagonal
        @test -b1 == -Matrix(b1)
    end

    @testset "Multiplication" begin

        @testset "BlockDiagonal * BlockDiagonal" begin
            @test b1 * b1 isa BlockDiagonal
            @test Matrix(b1 * b1) ≈ Matrix(b1) * Matrix(b1)
            @test_throws DimensionMismatch b3 * b1
        end

        @testset "BlockDiagonal * Number" begin
            @test b1 * 2 ≈ Matrix(b1) * 2 ≈ 2 * b1
            @test b1 / 2 ≈ Matrix(b1) / 2
            @test b1 * complex(2, 1) ≈ Matrix(b1) * complex(2, 1) ≈ complex(2, 1) * b1
            @test b1 / complex(2, 1) ≈ Matrix(b1) / complex(2, 1)
        end

        @testset "BlockDiagonal * Vector" begin
            @test b1 * a isa Vector
            @test b1 * a ≈ Matrix(b1) * a
            @test_throws DimensionMismatch b1 * b
        end

        @testset "BlockDiagonal * Matrix" begin
            @test b1 * A isa Matrix
            @test b1 * A ≈ Matrix(b1) * A
            @test_throws DimensionMismatch b1 * B

            # Matrix * BlockDiagonal
            @test A′ * b1 isa Matrix
            @test A′ * b1 ≈ A′ * Matrix(b1)
            @test_throws DimensionMismatch A * b1

            # degenerate cases
            m = rand(0, 0)
            @test m * BlockDiagonal([m]) == m * m == m
            m = rand(5, 0)
            @test m' * BlockDiagonal([m]) == m' * m == rand(0, 0)
            @test m * BlockDiagonal([m']) == m * m' == zeros(5, 5)
        end

        @testset "BlockDiagonal * Diagonal" begin
            D = Diagonal(randn(rng, N))
            D′ = Diagonal(randn(rng, N + N1))

            @test b1 * D isa BlockDiagonal
            @test b1 * D ≈ Matrix(b1) * D
            @test_throws DimensionMismatch D′ * b1

            # Diagonal * BlockDiagonal
            @test D * b1 isa BlockDiagonal
            @test D * b1 ≈ D * Matrix(b1)
            @test_throws DimensionMismatch D′ * b1
        end

        @testset "Non-Square BlockDiagonal * Non-Square BlockDiagonal" begin
    	    b4 = BlockDiagonal([ones(2, 4), 2 * ones(3, 2)])
            b5 = BlockDiagonal([3 * ones(2, 2), 2 * ones(4, 1)])

            @test b4 * b5 isa Array
            @test b4 * b5 == [6 * ones(2, 2) 4 * ones(2, 1); zeros(3, 2) 8 * ones(3, 1)]
            # Dimension check
            @test sum(size.(b4.blocks, 1)) == size(b4 * b5, 1)
            @test sum(size.(b5.blocks, 2)) == size(b4 * b5, 2)
        end
    end  # Multiplication
end

ext: jl


text:
using Gridap
using GridapGmsh

model = GmshDiscreteModel("demo/demo.msh")

order = 1
diritags = ["boundary1", "boundary2"]
fespace = ConformingFESpace(Float64,model,order,diritags)

ufun1(x) = 0.0
ufun2(x) = 1.0
V = TestFESpace(fespace)
U = TrialFESpace(fespace,[ufun1,ufun2])

trian = Triangulation(model)
quad = CellQuadrature(trian,order=2)

a(v,u) = inner(∇(v), ∇(u))
t_Ω = LinearFETerm(a,trian,quad)

assem = SparseMatrixAssembler(V,U)

op = LinearFEOperator(V,U,assem,t_Ω)

ls = LUSolver()
solver = LinearFESolver(ls)

uh = solve(solver,op)
writevtk(trian,"demo",cellfields=["uh"=>uh])


text:
function local_solution(n, m, supp::Vector{Vector{Vector{UInt16}}}, coe; nb=0, numeq=0,
    startpoint=[], QUIET=false)
    model = Model(optimizer_with_attributes(Ipopt.Optimizer))
    set_optimizer_attribute(model, MOI.Silent(), QUIET)
    if QUIET == true
        set_optimizer_attribute(model, "print_level", 0)
    end
    if isempty(startpoint)
        @variable(model, x[1:n], start = 0)
    else
        @variable(model, x[i=1:n], start = startpoint[i])
    end
    @NLobjective(model, Min, sum(coe[1][j]*prod(x[supp[1][j][k]] for k=1:length(supp[1][j])) for j=1:length(supp[1])))
    for i = 1:nb
        @NLconstraint(model, x[i]^2-1==0)
    end
    for i in 1:m-numeq
        @NLconstraint(model, sum(coe[i+1][j]*prod(x[supp[i+1][j][k]] for k=1:length(supp[i+1][j])) for j=1:length(supp[i+1]))>=0)
    end
    for i in m-numeq+1:m
        @NLconstraint(model, sum(coe[i+1][j]*prod(x[supp[i+1][j][k]] for k=1:length(supp[i+1][j])) for j=1:length(supp[i+1]))==0)
    end
    optimize!(model)
    status = termination_status(model)
    objv = objective_value(model)
    if QUIET == false
        println("optimum = $objv")
    end
    return objv,value.(x),status
end

function local_solution(n, m, supp::Vector{Array{UInt8, 2}}, coe; nb=0, numeq=0,
    startpoint=[], QUIET=false)
    model = Model(optimizer_with_attributes(Ipopt.Optimizer))
    set_optimizer_attribute(model, MOI.Silent(), QUIET)
    if QUIET == true
        set_optimizer_attribute(model, "print_level", 0)
    end
    if isempty(startpoint)
        @variable(model, x[1:n], start = 0)
    else
        @variable(model, x[i=1:n], start = startpoint[i])
    end
    @NLobjective(model, Min, sum(coe[1][j]*prod(x[k]^supp[1][k,j] for k=1:n) for j=1:size(supp[1], 2)))
    for i = 1:nb
        @NLconstraint(model, x[i]^2-1==0)
    end
    for i in 1:m-numeq
        @NLconstraint(model, sum(coe[i+1][j]*prod(x[k]^supp[i+1][k,j] for k=1:n) for j=1:size(supp[i+1], 2))>=0)
    end
    for i in m-numeq+1:m
        @NLconstraint(model, sum(coe[i+1][j]*prod(x[k]^supp[i+1][k,j] for k=1:n) for j=1:size(supp[i+1], 2))==0)
    end
    optimize!(model)
    status = termination_status(model)
    objv = objective_value(model)
    if QUIET == false
        println("optimum = $objv")
    end
    return objv,value.(x),status
end

"""
    ref_sol,upper_bound,rel_gap = refine_sol(opt, sol, data, QUIET=false, tol=1e-4)

Refine the obtained solution by a local solver.
Return the refined solution, the upper bound given by the local solver and the
relative optimality gap.
"""
function refine_sol(opt, sol, data::upop_data; QUIET=false, tol=1e-4)
    n = data.n
    nb = data.nb
    supp = data.supp
    coe = data.coe
    for i = 1:n
        if sol[i] == 0
            sol[i] = 1e-10
        end
    end
    upper_bound,rsol,status = local_solution(n, 0, [supp], [coe], nb=nb, numeq=0, startpoint=sol, QUIET=QUIET)
    if status == MOI.LOCALLY_SOLVED
        gap = abs(upper_bound)>1 ? abs((opt-upper_bound)/upper_bound) : abs(opt-upper_bound)
        if gap < tol
            println("Global optimality certified!")
        end
    else
        rsol,upper_bound,gap = sol,nothing,nothing
        println("The local solver failed!")
    end
    return rsol,upper_bound,gap
end

function refine_sol(opt, sol, data::Union{cpop_data,mcpop_data}; QUIET=false, tol=1e-4)
    n = data.n
    nb = data.nb
    numeq = data.numeq
    if typeof(data) == cpop_data && !isempty(data.gb)
        m = length(data.pop)-1
        supp = Vector{Array{UInt8,2}}(undef, m+1)
        coe = Vector{Vector{Float64}}(undef, m+1)
        supp[2:m+1-numeq] = data.supp[2:end]
        coe[2:m+1-numeq] = data.coe[2:end]
        for k in [1; [k for k=m+2-numeq:m+1]]
            mon = monomials(data.pop[k])
            coe[k] = coefficients(data.pop[k])
            supp[k] = zeros(UInt8, n, length(mon))
            for i = 1:length(mon), j = 1:n
                @inbounds supp[k][j,i] = MultivariatePolynomials.degree(mon[i], data.x[j])
            end
        end
    else
        m = data.m
        supp = data.supp
        coe = data.coe
    end
    for i = 1:n
        if sol[i] == 0
            sol[i] = 1e-10
        end
    end
    upper_bound,rsol,status = local_solution(n, m, supp, coe, nb=nb, numeq=numeq, startpoint=sol, QUIET=QUIET)
    if status == MOI.LOCALLY_SOLVED
        gap = abs(upper_bound)>1 ? abs((opt-upper_bound)/upper_bound) : abs(opt-upper_bound)
        if gap < tol
            println("Global optimality certified!")
        end
    else
        rsol,upper_bound,gap = sol,nothing,nothing
        println("The local solver failed!")
    end
    return rsol,upper_bound,gap
end


text:
function createelementdf(elements)

  #=
    mutable struct Element
      length::Float64
      id::Float64                 # Inside diameter
      od::Float64                 # Outside diameter
      npind::Int
    end
  =#

  element_types = [Float64, Float64, Float64, Int]
  element_names = [:length, :id, :od, :npind]

  elementdf = DataFrame(element_types, element_names, 0)
  for i in 1:size(elements, 1)
      append!(elementdf,
        DataFrame(
          length=elements[i].length,
          id=elements[i].id,
          od=elements[i].od,
          npind=elements[i].npind
        )
      )
  end

  elementdf

end


text:
using Test

push!(LOAD_PATH, joinpath(@__DIR__, ".."))
using Insolation

using Dates
using Statistics
using Roots
using Optim

using CLIMAParameters
using CLIMAParameters.Planet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

FT = Float32
@testset "Orbital Params" begin
    include("test_orbit_param.jl")
end
@testset "Types" begin
    include("test_types.jl")
end
@testset "Zenith Angle" begin
    include("test_zenith_angle.jl")
end
@testset "Insolation" begin
    include("test_insolation.jl")
end
@testset "Equinox Date" begin
    include("test_perihelion.jl")
    include("test_equinox.jl")
end

FT = Float64
@testset "Orbital Params" begin
    include("test_orbit_param.jl")
end
@testset "Types" begin
    include("test_types.jl")
end
@testset "Zenith Angle" begin
    include("test_zenith_angle.jl")
end
@testset "Insolation" begin
    include("test_insolation.jl")
end
@testset "Equinox Date" begin
    include("test_perihelion.jl")
    include("test_equinox.jl")
end


text:
"""
    covfn(x1, x2, gp_params)
Calculate GP covariance between x1 and x2 based on a squared exponential covariance function

Returns GP covariance matrix
"""
function covfn(x1, x2, gp_params::SVGP_params)
    rho = exp.(gp_params.log_rho)
    R   = exp.(-Distances.pairwise(SqEuclidean(), x1 ./ rho, x2 ./ rho, dims=1))

    return exp(gp_params.log_kappa[1]) * R
end

# function covfn(x1, x2, params::SVGP_params)
#     rho = exp.(params.log_rho)
#     x2    = x2 ./ rho
#     R     = exp.(-pairwise(haversine(), x1', x2')) # WORK ON BEING ABLE TO USE THIS KERNEL FOR EARTH SURFACE
#     return exp(params.log_kappa[1]) * R .+ (2.0 * params.mean_factor)^2
# end

"""
    covfn_lonlat(x1, x2, gp_params)
Calculate GP covariance between x1 and x2 based on the Haversine distance on a sphere

Warning: Currently does not work due to lack of auto-differentiation through Haversine distance call

Returns GP covariance matrix
"""
function covfn_lonlat(x1, x2, gp_params::SVGP_params)
    rho = exp.(gp_params.log_rho)
    R   = exp.(-Distances.pairwise(Haversine(), x1 ./ rho, x2 ./ rho, dims=1))

    return exp(gp_params.log_kappa[1]) * R
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
using StaticArrays
using Random

Random.seed!(10)

# data
T = Float64
n = 25
A = rand(SMatrix{n,n,T,n^2})
M = A * A'
b = rand(SVector{n,T})

# Basic inverse
x0 = inv(M) * b
@benchmark $x0 = inv($M) * $b

# Basic \
x1 = M \ b
@benchmark $x1 = $M \ $b
@show norm(x1 - x0)

# StaticArrays QR
qrM = StaticArrays.qr(M)
@benchmark $qrM = StaticArrays.qr($M)
x2 = qrM.R \ (qrM.Q' * b)
@benchmark $x2 = $qrM.R \ ($qrM.Q' * $b)
@show norm(x2 - x0)

# Custom QR
sol = SDMGSSolver(M)
factorize!(sol, M)
@benchmark factorize!(sol, M)
qr_solve!(sol, b)
@benchmark qr_solve!(sol, b)
x3 = sol.xs
@show norm(x3 - x0)


text:
using CubicEoS
using CubicEoSDatabase: Data

C₁C₅ = load(
    BrusilovskyEoSMixture,
    # names=("methane",),
    names = ("methane", "n-pentane"),
    component_dbs = (Data.martinez(), Data.brusilovsky_comp()),
    mix_eos_db = Data.brusilovsky_mix()
)

#= pure CH₄ =#
# V = 1e-6
# N = [0.01]
# RT = 140 * CubicEoS.GAS_CONSTANT_SI

#= C₁C₅ mixture
[0.547413, 0.452587] mole fractions from (Mikyska, 2013, Example 2)
V = 1e-6 from old vtflash
(5e3 * V) is (5 kmol m⁻³ * V) just an ordinary value
T = 370 Kelvins at 5 kmol m⁻³ should correspond to two-phase
T = 440 Kelvins at 5 kmol m⁻³ should correspond to single-phase
=#

V = 1e-6
N = (5000 * V) .* [0.547413, 0.452587]
RT = 300 * CubicEoS.GAS_CONSTANT_SI

issinglephase, tries = vt_stability(C₁C₅, N, V, RT)
dump(tries)

dump(CubicEoS.__vt_flash_init_conc_choose(tries))

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
# This file was generated, do not modify it. # hide
ŷ = MLJ.predict(LogisticModel, X)
residuals = [1 - pdf(ŷ[i], y[i,1]) for i in 1:nrow(y)]
r = report(LogisticModel)

k = collect(keys(fp.fitted_params_given_machine))[3]
println("\n Coefficients:  ", fp.fitted_params_given_machine[k].coef)
println("\n y \n ", y[1:5,1])
println("\n ŷ \n ", ŷ[1:5])
println("\n residuals \n ", residuals[1:5])
println("\n Standard Error per Coefficient \n", r.report_given_machine[k].stderror)


text:
#!/usr/bin/env julia
using Gtk.ShortNames, Gtk

win = GtkWindow("Sample App")
v = GtkBox(:v)
push!(win, v)

l = GtkLabel("You clicked 0 times.")
r = GtkButton("Reset")
b = GtkButton("Click")
push!(v, l)
push!(v, r)
push!(v, b)
set_gtk_property!(v, :expand, l, true)

niter = 0
function click(widget)
    global niter += 1
    set_gtk_property!(l, :label, "You clicked $niter times.")
    return nothing
end

function resetcount(widget)
    global niter = 0
    set_gtk_property!(l, :label, "You clicked 0 times.")
    return nothing
end

signal_connect(click, b, "clicked")
signal_connect(resetcount, r, "clicked")

showall(win)

if !isinteractive()
    c = Condition()
    signal_connect(win, :destroy) do widget
        notify(c)
    end
    wait(c)
end

max_stars_repo_name: goropikari/JuliaGUIAppExample

max_issues_repo_name: goropikari/JuliaGUIAppExample

max_forks_repo_name: goropikari/JuliaGUIAppExample


text:
# https://adventofcode.com/2020/day/11

# The following rules are applied to every seat simultaneously:
#
#   * If a seat is empty (L) and there are no occupied seats adjacent to it, the seat
#     becomes occupied.
#   * If a seat is occupied (#) and four or more seats adjacent to it are also occupied,
#     the seat becomes empty.
#   * Otherwise, the seat's state does not change.
#
# Floor (.) never changes; seats don't move, and nobody sits on the floor.

function decode_seat(char::Char)
    if char == 'L'
        return :Empty
    elseif char == '.'
        return :Floor
    elseif char == '#'
        return :Occupied
    else
        error("invalid input: $(char)")
    end
end

function read_layout(string::AbstractString)
    rows = strip.(split(string, '\n', keepempty=false))
    layout = Array{Symbol,2}(undef, length(rows), length(rows[1]))

    for (i, row) in enumerate(rows)
        for (j, char) in enumerate(row)
            layout[i, j] = decode_seat(char)
        end
    end
    return layout
end

@assert read_layout("L.#\n..L\n##.") == [:Empty :Floor :Occupied; :Floor :Floor :Empty; :Occupied :Occupied :Floor]

function count_adjecent_occupied_seats(layout, x, y)
    n, k = size(layout)
    occupied_seats = 0
    for i = (x - 1):(x + 1)
        if i < 1 || i > n
            continue
        end
        for j = (y - 1):(y + 1)
            if i == x && j == y || j < 1 || j > k
                continue
            end
            if layout[i, j] == :Occupied
                occupied_seats += 1
            end
        end
    end
    return occupied_seats
end

@assert count_adjecent_occupied_seats(read_layout("...\n...\n..."), 2, 2) == 0
@assert count_adjecent_occupied_seats(read_layout(".#.\n..#\n#.."), 2, 2) == 3
@assert count_adjecent_occupied_seats(read_layout(".#.\n.LL\n#.."), 2, 2) == 2
@assert count_adjecent_occupied_seats(read_layout("#####\n#...#\n#...#\n#...#\n#####"), 3, 3) == 0
@assert count_adjecent_occupied_seats(read_layout(".#.\n.#.\n..#"), 1, 1) == 2
@assert count_adjecent_occupied_seats(read_layout(".#.\n.##\n..."), 1, 3) == 3

function new_state(layout, x, y)
    occupied_seats = count_adjecent_occupied_seats(layout, x, y)
    if layout[x, y] == :Empty && occupied_seats == 0
        return :Occupied
    elseif layout[x, y] == :Occupied && occupied_seats >= 4
        return :Empty
    end
    return layout[x, y]
end

@assert new_state(read_layout("...\n...\n..."), 2, 2) == :Floor
@assert new_state(read_layout("...\n.L.\n..."), 2, 2) == :Occupied
@assert new_state(read_layout("...\n.L#\n..."), 2, 2) == :Empty
@assert new_state(read_layout("...\n.#.\n..."), 2, 2) == :Occupied
@assert new_state(read_layout(".#.\n###\n.#."), 2, 2) == :Empty
@assert new_state(read_layout(".#.\n##L\n.#."), 2, 2) == :Occupied
@assert new_state(read_layout("#.#\n.#.\n#.#"), 2, 2) == :Empty
@assert new_state(read_layout("###\n###\n###"), 2, 2) == :Empty
@assert new_state(read_layout("L..\n...\n..."), 1, 1) == :Occupied
@assert new_state(read_layout("##.\n##.\n..."), 1, 1) == :Occupied

function apply_rules(layout)
    n, k = size(layout)
    updated = Array{Symbol,2}(undef, n, k)
    for i = 1:n
        for j = 1:k
            updated[i, j] = new_state(layout, i, j)
        end
    end
    return updated
end

example1 = split(read("data/day-11-example1.txt", String), "\n\n", keepempty=false)

@assert apply_rules(read_layout(example1[1])) == read_layout(example1[2])
@assert apply_rules(read_layout(example1[2])) == read_layout(example1[3])
@assert apply_rules(read_layout(example1[3])) == read_layout(example1[4])
@assert apply_rules(read_layout(example1[4])) == read_layout(example1[5])
@assert apply_rules(read_layout(example1[5])) == read_layout(example1[6])
# Equilibrium?!
@assert apply_rules(read_layout(example1[6])) == read_layout(example1[6])

# How many seats end up occupied?

function part1(string)
    layout = read_layout(string)
    while true
        updated_layout = apply_rules(layout)
        if layout == updated_layout
            break
        end
        layout = updated_layout
    end
    return sum(layout .== :Occupied)
end

@assert part1(example1[1]) == 37

# Now, instead of considering just the eight immediately adjacent seats, consider the first
# seat in each of those eight directions. (...)
# Also, people seem to be more tolerant than you expected: it now takes five or more visible
# occupied seats for an occupied seat to become empty (rather than four or more from the
# previous rules). The other rules still apply (...)

function first_visible_seat(layout, x, y, dx, dy)
    n, k = size(layout)
    while true
        x += dx
        y += dy
        try
            if layout[x, y] in (:Empty, :Occupied)
                return layout[x, y]
            end
        catch BoundsError
            # we bumped the wall
            break
        end
    end
    return :Floor
end

@assert first_visible_seat(read_layout("L..\n...\n..."), 1, 1, +1, 0) == :Floor
@assert first_visible_seat(read_layout("L..\n...\n..."), 1, 1, 0, +1) == :Floor
@assert first_visible_seat(read_layout("L..\n...\n..."), 1, 1, +1, +1) == :Floor
@assert first_visible_seat(read_layout("L#.\n...\n..."), 1, 1, 0, +1) == :Occupied
@assert first_visible_seat(read_layout("LL#\n...\n..."), 1, 1, 0, +1) == :Empty
@assert first_visible_seat(read_layout("L#L\n...\n..."), 1, 1, 0, +1) == :Occupied
@assert first_visible_seat(read_layout("L.#\n...\n..."), 1, 1, 0, +1) == :Occupied
@assert first_visible_seat(read_layout("L..\n#..\n..."), 1, 1, +1, 0) == :Occupied
@assert first_visible_seat(read_layout("L..\n#..\nL.."), 1, 1, +1, 0) == :Occupied
@assert first_visible_seat(read_layout("L..\n...\n#.."), 1, 1, +1, 0) == :Occupied
@assert first_visible_seat(read_layout("L..\n.#.\n..."), 1, 1, +1, +1) == :Occupied
@assert first_visible_seat(read_layout("L..\n.#.\n..L"), 1, 1, +1, +1) == :Occupied
@assert first_visible_seat(read_layout("L..\n...\n..#"), 1, 1, +1, +1) == :Occupied
@assert first_visible_seat(read_layout("L..\n.L.\n..#"), 1, 1, +1, +1) == :Empty
@assert first_visible_seat(read_layout("#..\n...\n..L"), 3, 3, -1, -1) == :Occupied

function count_first_visible_seats(layout, x, y)
    n, k = size(layout)
    occupied_seats = 0

    for dx in (-1, 0, +1)
        for dy in (-1, 0, +1)
            dx == dy == 0 && continue
            if first_visible_seat(layout, x, y, dx, dy) == :Occupied
                occupied_seats += 1
            end
        end
    end
    return occupied_seats
end

@assert count_first_visible_seats(read_layout("#.#..\n.....\n.#L..\n.....\n...#."), 3, 3) == 3
@assert count_first_visible_seats(read_layout("L...#\n.....\n#....\n.#...\n....."), 1, 1) == 2

no_occupied_seats_example = "
.##.##.
#.#.#.#
##...##
...L...
##...##
#.#.#.#
.##.##.
"

@assert count_first_visible_seats(read_layout(no_occupied_seats_example), 4, 4) == 0
@assert count_first_visible_seats(read_layout("#.##.##.##\n#######.##"), 1, 3) == 5

function new_state_2(layout, x, y)
    occupied_seats = count_first_visible_seats(layout, x, y)
    if layout[x, y] == :Empty && occupied_seats == 0
        return :Occupied
    elseif layout[x, y] == :Occupied && occupied_seats >= 5
        return :Empty
    end
    return layout[x, y]
end

function apply_rules_2(layout)
    n, k = size(layout)
    updated = Array{Symbol,2}(undef, n, k)
    for i = 1:n
        for j = 1:k
            updated[i, j] = new_state_2(layout, i, j)
        end
    end
    return updated
end

example2 = split(read("data/day-11-example2.txt", String), "\n\n", keepempty=false)

@assert apply_rules_2(read_layout(example2[1])) == read_layout(example2[2])
@assert apply_rules_2(read_layout(example2[2])) == read_layout(example2[3])
@assert apply_rules_2(read_layout(example2[3])) == read_layout(example2[4])
@assert apply_rules_2(read_layout(example2[4])) == read_layout(example2[5])
@assert apply_rules_2(read_layout(example2[5])) == read_layout(example2[6])
# Equilibrium?!
@assert apply_rules_2(apply_rules_2(read_layout(example2[6]))) == apply_rules_2(read_layout(example2[6]))

function part2(string)
    layout = read_layout(string)
    while true
        updated_layout = apply_rules_2(layout)
        if layout == updated_layout
            break
        end
        layout = updated_layout
    end
    return sum(layout .== :Occupied)
end

@assert part2(example2[1]) == 26

test = read("data/day-11.txt", String)
println("Part 1: $(result1 = part1(test))")
println("Part 2: $(result2 = part2(test))")

@assert result1 == 2238
@assert result2 == 2013

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
issquare(A) = size(A, 1) == size(A, 2) && ndims(A) == 2

perf_hvcat(A, B) =  [A B; B A]

function perf_hvcat_setind(A, B)
    @assert issquare(A) && size(A) == size(B)
    n = size(A, 1)
    C = Matrix{Float64}(2n, 2n)
    C[1:n, 1:n] = A
    C[1:n, (n+1):end] = B
    C[(n+1):end, 1:n] = B
    C[(n+1):end, (n+1):end] = A
end

perf_hcat(A, B) = [A B B A]

function perf_hcat_setind(A, B)
    @assert issquare(A) && size(A) == size(B)
    n = size(A, 1)
    C = Matrix{Float64}(n, 4n)
    C[:, 1:n] = A
    C[:, (n+1):2n] = B
    C[:, (2n+1):3n] = B
    C[:, (3n+1):end] = A
    return C
end

perf_vcat(A, B) = [A; B; B; A]

function perf_vcat_setind(A, B)
    @assert issquare(A) && size(A) == size(B)
    n = size(A, 1)
    C = Matrix{Float64}(4n, n)
    C[1:n, :] = A
    C[(n+1):2n, :] = B
    C[(2n+1):3n, :] = B
    C[(3n+1):4n, :] = A
    return C
end

function perf_catnd(n)
    A = samerand(1, n, n, 1)
    B = samerand(1, n, n)
    return cat(3, A, B, B, A)
end

function perf_catnd_setind(n)
    A = samerand(1, n, n, 1)
    B = samerand(1, n, n)
    C = Array{Float64}(1, n, 4n, 1)
    C[1, :, 1:n, 1] = A
    C[1, :, (n+1):2n, 1] = B
    C[1, :, (2n+1):3n, 1] = B
    C[1, :, (3n+1):4n, 1] = A
    return C
end


text:
function section{T<:AbstractFloat}(x::T, b::T, n::Int)
    fld(n * (x + b/2), b)
end

function density{T<:AbstractFloat}(file_out::IO, C::AbstractMatrix{T}, n::Int, bw::T = zero(Float64))
    num_atoms = size(C, 2)

    b::T = 2 * maxabs(C)
    println(b)

    if (bw != zero(T))
        bw < b && error("Your box specification is not big enough to contain all atoms.")

        b = bw
    else
        b *= (1 + 1/n)
    end

    println(file_out, "density")
    println(file_out, "$b   "^3, "90.00000   "^3)
    println(file_out, "$n   "^3)

    indices, counts = Int[], Int[]
    @inbounds for i in 1:num_atoms
        xs = section(C[1, i], b, n)
        ys = section(C[2, i], b, n)
        zs = section(C[3, i], b, n)
        index = zs + n*ys + n*n*xs + 1
        push!(indices, index)
        push!(counts, 1)
    end
    ρ = sparsevec(indices, counts)
    full_ρ = full(ρ)

    remaining_zeros = n^3 - length(full_ρ)

    print_joined(file_out, full_ρ, '\n')
    println(file_out, "\n0" ^ remaining_zeros)
end

function density2{T<:AbstractFloat}(file_out::IO, C::AbstractMatrix{T}, n::Int, bw::T = zero(Float64))
    num_atoms = size(C, 2)

    b::T = 2 * maxabs(C)
    println(b)

    if (bw != zero(T))
        bw < b && error("Your box specification is not big enough to contain all atoms.")

        b = bw
    else
        b *= (1 + 1/n)
    end

    println(file_out, "density")
    println(file_out, "$b   "^3, "90.00000   "^3)
    println(file_out, "$n   "^3)

    indices, counts = Int[], Int[]
    @inbounds for i in 1:num_atoms
        xs = section(C[1, i], b, n)
        ys = section(C[2, i], b, n)
        zs = section(C[3, i], b, n)
        index = zs + n*ys + n*n*xs + 1
        push!(indices, index)
        push!(counts, 1)
    end
    ρ = sparsevec(indices, counts)
    full_ρ = full(ρ)

    remaining_zeros = n^3 - length(full_ρ)

    print_joined(file_out, full_ρ, '\n')
    println(file_out, "\n0" ^ remaining_zeros)
end

max_stars_repo_name: NickMcNutt/Alpine

max_issues_repo_name: NickMcNutt/Alpine

max_forks_repo_name: NickMcNutt/Alpine


text:
#typealias Vec3 Vector3{Float64}
#typealias Vec2 Vector2{Float64}
#typealias Mat2 Matrix2x2{Float64}
#typealias Mat3 Matrix3x3{Float64}

import Base: Func, AddFun

immutable SubFun <: Func{2} end
call(::SubFun, x, y) = x - y

abstract AbstractGaussPoint

immutable GaussPoint2 <: AbstractGaussPoint
    local_coords::Point2
    weight::Float64
end

immutable GaussPoint3 <: AbstractGaussPoint
    local_coords::Point3
    weight::Float64
end

getweight(gp::AbstractGaussPoint) = gp.weight

# Lint does not refognize enums
@lintpragma("Ignore use of undeclared variable DofType")
@lintpragma("Ignore use of undeclared variable Du")
@lintpragma("Ignore use of undeclared variable Dv")
@lintpragma("Ignore use of undeclared variable Dw")
@lintpragma("Ignore use of undeclared variable Gu1")
@lintpragma("Ignore use of undeclared variable Gv1")
@lintpragma("Ignore use of undeclared variable Gu2")
@lintpragma("Ignore use of undeclared variable Gv2")

@enum DofType Du Dv Dw Gu1 Gv1 Gu2 Gv2 Gu3 Gv3 K1 K2

immutable Dof
    eq_n::Int
    dof_type::DofType # Not here
    active::Bool
end
@inline isactive(dof::Dof) = dof.active

immutable DofVals
    free_dof_values::Vector{Float64}
    presc_dof_values::Vector{Float64}
end
DofVals() = DofVals(zeros(0), zeros(0))

function get_value(dv::DofVals, dof::Dof)
    if isactive(dof)
        return dv.free_dof_values[dof.eq_n]
    else
        return dv.presc_dof_values[dof.eq_n]
    end
end

function mod_value(dv::DofVals, dof::Dof, f::Func, val)
    if isactive(dof)
        dv.free_dof_values[dof.eq_n] = f(dv.free_dof_values[dof.eq_n], val)
    else
        dv.presc_dof_values[dof.eq_n] =  f(dv.presc_dof_values[dof.eq_n] , val)
    end
end

type TimeStepper{T <: Range}
    timesteps::T
    t::Float64
    dt::Float64
end

#########
# Nodes #
#########
abstract AbstractFENode

immutable FENode2 <: AbstractFENode
    n::Int
    coords::Point2
    dofs::Vector{Dof}
end
@inline get_dof(node::AbstractFENode, i::Int) = node.dofs[i]
@inline get_dofs(node::AbstractFENode) = node.dofs

#TODO: Dispatch on dimension
FENode2(n::Int, c::Point2) = FENode2(n, c, Array(Dof, 0))
# Check size here?
#FENode2(n::Int, c::Vector{Float64}) = FENode2(n, Point2(c[1], c[2]), Array(Dof, 0))
function FENode2{T <: Number}(n::Int, c::Vector{T})
    c_v = convert(Vector{Float64}, c)
    FENode2(n, Point2(c[1], c[2]), Array(Dof, 0))
end

# TODO: Move
@inline get_coord(node::FENode2) = Point3(node.coords[1], node.coords[2], 0.0)
@inline get_displacement(node::FENode2, dv::DofVals) = Point3(get_value(dv, node.dofs[1]),
                                                              get_value(dv, node.dofs[2]),
                                                              0.0)

immutable FENode3 <: AbstractFENode
    n::Int
    coords::Point3
    dofs::Vector{Dof}
end

FENode3(n::Int, c::Point3) = FENode2(n, c, Array(Dof, 0))
# Check size here?
#FENode3(n::Int, c::Vector{Float64}) = FENode2(n, Point3(c[1], c[2], c[3]), Array(Dof, 0))
#FENode3(n::Int, c::Vector{Int}) = FENode2(n, convert(Vector{Float64}, c), Array(Dof, 0))
function FENode3{T <: Number}(n::Int, c::Vector{T})
    c_v = convert(Vector{Float64}, c)
    FENode2(n, Point3(c[1], c[2], c[3]), Array(Dof, 0))
end

# TODO: Add more sets
# Only constant loads right now
#=
immutable DirichletBC
    value::Float64
    dof_types::Vector{DofType}
    node_set::NodeSet
end
=#

immutable DirichletBC{f <: FastAnonymous.Fun}
    func::f
    dof_types::Vector{DofType}
    node_set::NodeSet
end
function DirichletBC(f::String, dof_types::Vector{DofType}, node_set::NodeSet)
    prog = string("ff = @anon (x,y,z,t) -> ", f)
    eval(parse(prog))
    DirichletBC(ff, dof_types,  node_set)
end

immutable NodeLoad{f <: FastAnonymous.Fun}
    func::f
    dof_types::Vector{DofType}
    node_set::NodeSet
end

function NodeLoad(f::String, dof_types::Vector{DofType}, node_set::NodeSet)
    prog = string("ff = @anon (x,y,z,t) -> ", f)
    eval(parse(prog))
    NodeLoad(ff, dof_types,  node_set)
end

function evaluate(bc::Union(DirichletBC, NodeLoad), node::FENode2, t::Number)
    return bc.func(node.coords.x, node.coords.y, 0.0, t)
end

immutable EdgeLoad
    value::Float64
    dof_types::Vector{DofType}
    edge_set::EdgeSet
end

# TODO Add more loads.

function fill_from_start!{T}(v1::AbstractArray{T}, v2::AbstractArray{T})
    if length(v2) > length(v1)
        error("Array to be filled is shorter than array to fill with")
    end
    for i in 1:length(v2)
        v1[i] = v2[i]
    end
end

function assemble!{T <: Real}(A::Array, X::AbstractArray, I::AbstractVector{T})
   count = 1
   for i in I
        A[count] = X[i]
        count += 1
   end
end


text:
"""
    A = JopProjection(u::AbstractArray)

`A*m` is the projection of vector `m` onto the vector `u`.
"""
function JopProjection(u::AbstractArray{T,N}) where {T,N}
    u_normed = zeros(eltype(u), size(u))
    for idx in CartesianIndices(size(u)[1:end-1])
        u_normed[idx,:] .= u[idx,:] ./ norm(u[idx,:])
    end
    n = size(u)
    dom = JetSpace(T, n)
    rng = JetSpace(T, n[1:end-1])
    JopLn(dom = dom, rng = rng, df! = JopProjection_df!, df′! = JopProjection_df′!, s = (u=u_normed,))
end

export JopProjection

function JopProjection_df!(d::AbstractArray, m::AbstractArray; u, kwargs...)
    for idx in CartesianIndices(size(d))
        d[idx] = dot(u[idx,:], m[idx,:])
    end
    d
end

function JopProjection_df′!(m::AbstractArray, d::AbstractArray; u, kwargs...)
    for idx in CartesianIndices(size(d))
        for i = 1:size(m,ndims(m))
            m[idx,i] = u[idx,i]*d[idx]
        end
    end
    m
end


text:
#Constants

const accl_g = 9.80665

# ---------------------------------------------------------------------------------------------
# Theodorsen solver input
type TheoDef
    alpha_amp :: Float64
    h_amp :: Float64
    alpha_mean :: Float64
    alpha_zl :: Float64
    k :: Float64
    phi :: Float64
    pvt :: Float64
end

type DelVortDef
    flag :: Int8
    limit :: Int16
    dist :: Float64
end

type TheoDefwFlap
    alpha_amp :: Float64
    h_amp :: Float64
    alpha_mean :: Float64
    alpha_zl :: Float64
    k :: Float64
    phi :: Float64
    pvt :: Float64
    beta_amp :: Float64
    xf :: Float64
    psi :: Float64
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# Data structures used to represent the airfoil kinematics
type KinemPar
    alpha :: Float64
    h :: Float64
    alphadot :: Float64
    hdot :: Float64
    u :: Float64
    udot :: Float64
end

# Added by Laura Merchant 2016
type KinemParwFlap
    alpha :: Float64
    h :: Float64
    alphadot :: Float64
    hdot :: Float64
    u :: Float64
    udot :: Float64
    n :: Float64
    ndot :: Float64
end

#For airfoil with 2DOF in pitch and plunge
type KinemPar2DOF
    alpha :: Float64
    h :: Float64
    alphadot :: Float64
    hdot :: Float64
    u :: Float64
    udot :: Float64
    alphaddot :: Float64
    hddot :: Float64
    alpha_pr :: Float64
    alpha_pr2 :: Float64
    alpha_pr3 :: Float64
    h_pr :: Float64
    h_pr2 :: Float64
    h_pr3 :: Float64
    alphadot_pr :: Float64
    alphadot_pr2 :: Float64
    alphadot_pr3 :: Float64
    hdot_pr :: Float64
    hdot_pr2 :: Float64
    hdot_pr3 :: Float64
    alphaddot_pr :: Float64
    alphaddot_pr2 :: Float64
    alphaddot_pr3 :: Float64
    hddot_pr :: Float64
    hddot_pr2 :: Float64
    hddot_pr3 :: Float64
end

type KinemPar2DFree
    alpha :: Float64
    h :: Float64
    alphadot :: Float64
    hdot :: Float64
    u :: Float64
    udot :: Float64
    alphaddot :: Float64
    hddot :: Float64
    alpha_pr :: Float64
    h_pr :: Float64
    alphadot_pr :: Float64
    alphadot_pr2 :: Float64
    alphadot_pr3 :: Float64
    hdot_pr :: Float64
    hdot_pr2 :: Float64
    hdot_pr3 :: Float64
    alphaddot_pr :: Float64
    alphaddot_pr2 :: Float64
    alphaddot_pr3 :: Float64
    hddot_pr :: Float64
    hddot_pr2 :: Float64
    hddot_pr3 :: Float64
    u_pr :: Float64
    udot_pr :: Float64
    udot_pr2 ::Float64
    udot_pr3 :: Float64
end

type TwoDOFPar
    x_alpha :: Float64
    r_alpha :: Float64
    kappa :: Float64
    w_alpha :: Float64
    w_h :: Float64
    w_alphadot :: Float64
    w_hdot :: Float64
    cubic_h_1 :: Float64
    cubic_h_3 :: Float64
    cubic_alpha_1 :: Float64
    cubic_alpha_3 :: Float64
end

type TwoDFreePar
    r_g :: Float64
    x_g :: Float64
    kappa :: Float64
end

# ---------------------------------------------------------------------------------------------
# Various motion kinematics definitions
abstract MotionDef

immutable EldUpDef <: MotionDef
    amp :: Float64
    K :: Float64
    a :: Float64
end

immutable EldUptstartDef <: MotionDef
    amp :: Float64
    K :: Float64
    a :: Float64
    tstart :: Float64
end

immutable EldRampReturnDef <: MotionDef
    amp :: Float64
    K :: Float64
    a :: Float64
end

immutable ConstDef <: MotionDef
    amp :: Float64
end

function (eld::EldUpDef)(t)
    sm = pi*pi*eld.K/(2*(eld.amp)*(1 - eld.a))
    t1 = 1.
    t2 = t1 + ((eld.amp)/(2*eld.K))
    ((eld.K/sm)*log(cosh(sm*(t - t1))/cosh(sm*(t - t2))))+(eld.amp/2)
end

function (eld::EldUptstartDef)(t)
    sm = pi*pi*eld.K/(2*(eld.amp)*(1 - eld.a))
    t1 = eld.tstart
    t2 = t1 + ((eld.amp)/(2*eld.K))
    ((eld.K/sm)*log(cosh(sm*(t - t1))/cosh(sm*(t - t2))))+(eld.amp/2)
end

function (cons::ConstDef)(t)
    cons.amp
end

function (eld::EldRampReturnDef)(tt)
    fr = eld.K/(pi*abs(eld.amp));
    t1 = 1.
    t2 = t1 + (1./(2*pi*fr));
    t3 = t2 + ((1/(4*fr)) - (1/(2*pi*fr)));
    t4 = t3 + (1./(2*pi*fr));
    t5 = t4+1.;

    nstep = round(Int,t5/0.015) + 1
    g = zeros(nstep)
    t = zeros(nstep)

    for i = 1:nstep
        t[i] = (i-1.)*0.015
        g[i] = log((cosh(eld.a*(t[i] - t1))*cosh(eld.a*(t[i] - t4)))/(cosh(eld.a*(t[i] - t2))*cosh(eld.a*(t[i] - t3))))
    end
    maxg = maximum(g);

    gg = log((cosh(eld.a*(tt - t1))*cosh(eld.a*(tt - t4)))/(cosh(eld.a*(tt - t2))*cosh(eld.a*(tt - t3))))

    return eld.amp*gg/(maxg);

end

immutable SinDef <: MotionDef
  mean :: Float64
  amp :: Float64
  k :: Float64
  phi :: Float64
end

immutable CosDef <: MotionDef
  mean :: Float64
  amp :: Float64
  k :: Float64
  phi :: Float64
end

function (kin::SinDef)(t)
  (kin.mean) + (kin.amp)*sin(2*kin.k*t + kin.phi)
end

function (kin::CosDef)(t)
  (kin.mean) + (kin.amp)*cos(2*kin.k*t + kin.phi)
end
# ---------------------------------------------------------------------------------------------

# Added by Laura Merchant 2016
type KinemDefwFlap
    alpha :: MotionDef
    h :: MotionDef
    u :: MotionDef
    n :: MotionDef
end

type KinemDef
    alpha :: MotionDef
    h :: MotionDef
    u :: MotionDef
end

type KinemDef3D
    alpha :: MotionDef
    h :: MotionDef
    u :: MotionDef
    vartype :: String
    vary ::Int8
    add :: Array{Float64,2}

    function KinemDef3D(alpha :: MotionDef, h::MotionDef, u::MotionDef, vartype = "Constant", vary = 0, add = zeros(1,1))
        new(alpha, h, u, vartype, vary, add)
    end
end

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

# ---------------------------------------------------------------------------------------------
#Definition of a 2D vortex blob
type TwoDVort
    x :: Float64
    z :: Float64
    s :: Float64
    vc :: Float64
    vx :: Float64
    vz :: Float64
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
#Definition of a 3D vortex blob
type ThreeDVort
    x :: Vector{Float64}
    s :: Vector{Float64}
    vc :: Float64
    v :: Vector{Float64}
    ds :: Vector{Float64}
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
#Definition of flowfield in 2D
immutable TwoDFlowField
    velX :: MotionDef
    velZ :: MotionDef
    u :: Vector{Float64}
    w :: Vector{Float64}
    tev :: Vector{TwoDVort}
    lev :: Vector{TwoDVort}
    extv :: Vector{TwoDVort}
    function TwoDFlowField(velX = ConstDef(0.), velZ = ConstDef(0.))
        u = [0;]
        w = [0;]
        tev = TwoDVort[]
        lev = TwoDVort[]
        extv = TwoDVort[]
        new(velX, velZ, u, w, tev, lev, extv)
    end
end

#Definition of flowfield in 3D
immutable ThreeDFieldSimple
    f2d :: Vector{TwoDFlowField}
    function ThreeDFieldSimple()
        f2d = TwoDFlowField[]
    new(f2d)
    end
end

# immutable ThreeDFlowField
#     velX :: MotionDef
#     velY :: MotionDef
#     velZ :: MotionDef
#     u :: Vector{Float64}
#     v :: Vector{Float64}
#     w :: Vector{Float64}
#     tev :: Vector{ThreeDVort}
#     lev :: Vector{ThreeDVort}
#     extv :: Vector{ThreeDVort}
#     function ThreeDFlowField(velX = ConstDef(0.), velY = ConstDef(0.), velZ = ConstDef(0.))
#         u = [0;]
#         v = [0;]
#         w = [0;]
#         tev = ThreeDVort[]
#         lev = ThreeDVort[]
#         extv = ThreeDVort[]
#         new(velX, velY, velZ, u, v, w, tev, lev, extv)
#     end
# end

immutable TwoDFlowData
    tev :: Vector{TwoDVort}
    lev :: Vector{TwoDVort}
    extv :: Vector{TwoDVort}
    bv :: Vector{TwoDVort}
end

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

# ---------------------------------------------------------------------------------------------
# BEGIN TwoDSurfwFlap
# Added by Laura Merchant 2016
immutable TwoDSurfwFlap
    c :: Float64
    uref :: Float64
    coord_file :: String
    pvt :: Float64
    ndiv :: Int8
    naterm :: Int8
    dynamics_type :: String
    kindef :: KinemDefwFlap
    cam_af :: Vector{Float64}
    cam :: Vector{Float64}
    cam_slope :: Vector{Float64}
    cam_tder :: Vector{Float64}
    theta :: Vector{Float64}
    x :: Vector{Float64}
    x_b :: Vector{Float64}
    kinem :: KinemParwFlap
    bnd_x :: Vector{Float64}
    bnd_z :: Vector{Float64}
    uind :: Vector{Float64}
    wind :: Vector{Float64}
    downwash :: Vector{Float64}
    a0 :: Vector{Float64}
    aterm :: Vector{Float64}
    a0dot :: Vector{Float64}
    adot :: Vector{Float64}
    a0prev :: Vector{Float64}
    aprev :: Vector{Float64}
    bv :: Vector{TwoDVort}
    bv_prev :: Vector{TwoDVort}
    lespcrit :: Vector{Float64}
    levflag :: Vector{Int8}

    function TwoDSurfwFlap(c, uref, coord_file, pvt, ndiv, naterm, dynamics_type, kindef, x_b=1., lespcrit=zeros(1))
        x_b[1] = x_b[1]*c #Dimensional value
        theta = zeros(ndiv)
        x = zeros(ndiv)
        cam = zeros(ndiv)
        cam_af = zeros(ndiv)
        cam_slope = zeros(ndiv)
        cam_tder = zeros(ndiv)
        bnd_x = zeros(ndiv)
        bnd_z = zeros(ndiv)
        kinem = KinemParwFlap(0, 0, 0, 0, 0, 0, 0, 0)

        dtheta = pi/(ndiv-1)
        for ib = 1:ndiv
            theta[ib] = real(ib-1.)*dtheta
            x[ib] = c/2.*(1-cos(theta[ib]))
        end
        if (coord_file != "FlatPlate")
            cam, cam_slope = camber_calc(x, coord_file)
        end

        # Calculates the initial condition for airfoil pitch
        # Various options for prescribing the airfoil pitch kinematics
        if (typeof(kindef.alpha) == EldUpDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == EldUptstartDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == EldRampReturnDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == ConstDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = 0.
        elseif (typeof(kindef.alpha) == SinDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == CosDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        end
        # ---------------------------------------------------------------------------------------------

        # Calculates the initial condition for airfoil plunge
        # Various options for prescribing the airfoil plunge kinematics
        if (typeof(kindef.h) == EldUpDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldUptstartDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldUpIntDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldUpInttstartDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldRampReturnDef)
            kinem.h= kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == ConstDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = 0.
        elseif (typeof(kindef.h) == SinDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == CosDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        end
        # ---------------------------------------------------------------------------------------------

        # Calculates the initial condition for forward velocity of airfoil
        # Various options for prescribing the airfoil forward motion
        if (typeof(kindef.u) == EldUpDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/c
        elseif (typeof(kindef.u) == EldRampReturnDef)
            kinem.u, kinem.udot = kindef.u(0.)
            kinem.u = kinem.u*uref
            kinem.udot = kinem.udot*uref*uref/c
        elseif (typeof(kindef.u) == ConstDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = 0.
        end
        # ---------------------------------------------------------------------------------------------

        # ---------------------------------------------------------------------------------------------
        # Calculates the initial condition for flap deflection beta (n)
        # Various options for prescribing the flap kinematics
        if (typeof(kindef.n) == EldUpDef)
            kinem.n = kindef.n(0.)
            kinem.ndot = ForwardDiff.derivative(kindef.n,0.)*uref/c
        elseif (typeof(kindef.n) == EldRampReturnDef)
            kinem.n = kindef.n(0.)
            kinem.n = ForwardDiff.derivative(kindef.n,0.)*uref/c
        elseif (typeof(kindef.n) == ConstDef)
            kinem.n = kindef.n(0.)
            kinem.n = 0.
        elseif (typeof(kindef.n) == SinDef)
            kinem.n = kindef.n(0.)
            kinem.ndot = ForwardDiff.derivative(kindef.n,0.)*uref/c
        elseif (typeof(kindef.n) == CosDef)
            kinem.n = kindef.n(0.)
            kinem.ndot = ForwardDiff.derivative(kindef.n,0.)*uref/c
        end
# ---------------------------------------------------------------------------------------------

# Defines the aerofoil as a single flat plate
# Calls the calculation of the camber and it's derivatives (spatial and temporal)
#Deflection of flap is modelled as camber variation
if (coord_file != "FlatPlate")
    (cam_af, cam_slope) = camber_calc(x, coord_file)
end

for i = 1:ndiv
    if x[i] < x_b[1]
        cam[i] = cam_af[i];
    else
        cam[i] = cam_af[i] + (x_b[1] - x[i])*tan(kinem.n);
        cam_tder[i] = (x_b[1] - x[i])*kinem.ndot*sec(kinem.n)*sec(kinem.n);
    end
end
# ---------------------------------------------------------------------------------------------

# Populates the arrays for bnd_x, bnd_z based on initial conditions
for i = 1:ndiv
    bnd_x[i] = -((c - pvt*c)+((pvt*c - x[i])*cos(kinem.alpha))) + (cam[i]*sin(kinem.alpha))
    bnd_z[i] = kinem.h + ((pvt*c - x[i])*sin(kinem.alpha))+(cam[i]*cos(kinem.alpha))
end
# ---------------------------------------------------------------------------------------------

# Defines the arrays for various parameters such that they can be later populated
uind = zeros(ndiv)
wind = zeros(ndiv)
downwash = zeros(ndiv)
a0 = zeros(1)
a0dot = zeros(1)
aterm = zeros(naterm)
adot = zeros(3)
a0prev = zeros(1)
aprev = zeros(3)
bv = TwoDVort[]
bv_prev = TwoDVort[]
for i = 1:ndiv-1
    push!(bv,TwoDVort(0,0,0,0.02*c,0,0))
    push!(bv_prev,TwoDVort(0,0,0,0.02*c,0,0))
end
# ---------------------------------------------------------------------------------------------
levflag = [0]

new(c, uref, coord_file, pvt, ndiv, naterm, dynamics_type, kindef, cam_af, cam, cam_slope, cam_tder, theta, x, x_b, kinem, bnd_x, bnd_z, uind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv, bv_prev, lespcrit,levflag)
end
end
# ---------------------------------------------------------------------------------------------
# END TwoDSurfwFlap

immutable TwoDSurf
    c :: Float64
    uref :: Float64
    coord_file :: String
    pvt :: Float64
    ndiv :: Int8
    naterm :: Int8
    kindef :: KinemDef
    cam :: Vector{Float64}
    cam_slope :: Vector{Float64}
    theta :: Vector{Float64}
    x :: Vector{Float64}
    kinem :: KinemPar
    bnd_x :: Vector{Float64}
    bnd_z :: Vector{Float64}
    uind :: Vector{Float64}
    wind :: Vector{Float64}
    downwash :: Vector{Float64}
    a0 :: Vector{Float64}
    aterm :: Vector{Float64}
    a0dot :: Vector{Float64}
    adot :: Vector{Float64}
    a0prev :: Vector{Float64}
    aprev :: Vector{Float64}
    bv :: Vector{TwoDVort}
    lespcrit :: Vector{Float64}
    levflag :: Vector{Int8}

    function TwoDSurf(coord_file, pvt, kindef,lespcrit=zeros(1), c=1., uref=1., ndiv=70, naterm=35)
        theta = zeros(ndiv)
        x = zeros(ndiv)
        cam = zeros(ndiv)
        cam_slope = zeros(ndiv)
        bnd_x = zeros(ndiv)
        bnd_z = zeros(ndiv)
        kinem = KinemPar(0, 0, 0, 0, 0, 0)

        dtheta = pi/(ndiv-1)
        for ib = 1:ndiv
            theta[ib] = real(ib-1.)*dtheta
            x[ib] = c/2.*(1-cos(theta[ib]))
        end
        if (coord_file != "FlatPlate")
            cam, cam_slope = camber_calc(x, coord_file)
        end

        if (typeof(kindef.alpha) == EldUpDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == EldRampReturnDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == ConstDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = 0.
        elseif (typeof(kindef.alpha) == SinDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        elseif (typeof(kindef.alpha) == CosDef)
            kinem.alpha = kindef.alpha(0.)
            kinem.alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/c
        end

        if (typeof(kindef.h) == EldUpDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
          elseif (typeof(kindef.h) == EldUpIntDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == EldRampReturnDef)

            kinem.h= kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == ConstDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = 0.
        elseif (typeof(kindef.h) == SinDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        elseif (typeof(kindef.h) == CosDef)
            kinem.h = kindef.h(0.)*c
            kinem.hdot = ForwardDiff.derivative(kindef.h,0.)*uref
        end

        if (typeof(kindef.u) == EldUpDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/c
        elseif (typeof(kindef.u) == EldRampReturnDef)
            kinem.u, kinem.udot = kindef.u(0.)
            kinem.u = kinem.u*uref
            kinem.udot = kinem.udot*uref*uref/c
        elseif (typeof(kindef.u) == ConstDef)
            kinem.u = kindef.u(0.)*uref
            kinem.udot = 0.
        end

        for i = 1:ndiv
            bnd_x[i] = -((c - pvt*c)+((pvt*c - x[i])*cos(kinem.alpha))) + (cam[i]*sin(kinem.alpha))
            bnd_z[i] = kinem.h + ((pvt*c - x[i])*sin(kinem.alpha))+(cam[i]*cos(kinem.alpha))
        end
        uind = zeros(ndiv)
        wind = zeros(ndiv)
        downwash = zeros(ndiv)
        a0 = zeros(1)
        a0dot = zeros(1)
        aterm = zeros(naterm)
        adot = zeros(3)
        a0prev = zeros(1)
        aprev = zeros(3)
        bv = TwoDVort[]
        for i = 1:ndiv-1
            push!(bv,TwoDVort(0,0,0,0.02*c,0,0))
        end
        levflag = [0]
        new(c, uref, coord_file, pvt, ndiv, naterm, kindef, cam, cam_slope, theta, x, kinem, bnd_x, bnd_z, uind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv,lespcrit,levflag)
    end
end

immutable TwoDSurf_2DOF
    c :: Float64
    uref :: Float64
    coord_file :: String
    pvt :: Float64
    ndiv :: Int8
    naterm :: Int8
    cam :: Vector{Float64}
    cam_slope :: Vector{Float64}
    theta :: Vector{Float64}
    x :: Vector{Float64}
    strpar :: TwoDOFPar
    kinem :: KinemPar2DOF
    bnd_x :: Vector{Float64}
    bnd_z :: Vector{Float64}
    uind :: Vector{Float64}
    wind :: Vector{Float64}
    downwash :: Vector{Float64}
    a0 :: Vector{Float64}
    aterm :: Vector{Float64}
    a0dot :: Vector{Float64}
    adot :: Vector{Float64}
    a0prev :: Vector{Float64}
    aprev :: Vector{Float64}
    bv :: Vector{TwoDVort}
    lespcrit :: Vector{Float64}
    levflag :: Vector{Int8}

    function TwoDSurf_2DOF(c, uref, coord_file, pvt, ndiv, naterm, strpar, kinem ,lespcrit=zeros(1))
        theta = zeros(ndiv)
        x = zeros(ndiv)
        cam = zeros(ndiv)
        cam_slope = zeros(ndiv)
        bnd_x = zeros(ndiv)
        bnd_z = zeros(ndiv)

        dtheta = pi/(ndiv-1)
        for ib = 1:ndiv
            theta[ib] = real(ib-1.)*dtheta
            x[ib] = c/2.*(1-cos(theta[ib]))
        end
        if (coord_file != "FlatPlate")
            cam, cam_slope = camber_calc(x, coord_file)
        end

        for i = 1:ndiv
            bnd_x[i] = -((c - pvt*c)+((pvt*c - x[i])*cos(kinem.alpha))) + (cam[i]*sin(kinem.alpha))
            bnd_z[i] = kinem.h + ((pvt*c - x[i])*sin(kinem.alpha))+(cam[i]*cos(kinem.alpha))
        end
        uind = zeros(ndiv)
        wind = zeros(ndiv)
        downwash = zeros(ndiv)
        a0 = zeros(1)
        a0dot = zeros(1)
        aterm = zeros(naterm)
        adot = zeros(3)
        a0prev = zeros(1)
        aprev = zeros(3)
        bv = TwoDVort[]
        for i = 1:ndiv-1
            push!(bv,TwoDVort(0,0,0,0.02*c,0,0))
        end
        levflag = [0]
        new(c, uref, coord_file, pvt, ndiv, naterm, cam, cam_slope, theta, x, strpar, kinem, bnd_x, bnd_z, uind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv,lespcrit,levflag)
    end
end

immutable TwoDFreeSurf
    c :: Float64
    uref :: Float64
    coord_file :: String
    pvt :: Float64
    ndiv :: Int8
    naterm :: Int8
    cam :: Vector{Float64}
    cam_slope :: Vector{Float64}
    theta :: Vector{Float64}
    x :: Vector{Float64}
    strpar :: TwoDFreePar
    kinem :: KinemPar2DFree
    bnd_x :: Vector{Float64}
    bnd_z :: Vector{Float64}
    uind :: Vector{Float64}
    wind :: Vector{Float64}
    downwash :: Vector{Float64}
    a0 :: Vector{Float64}
    aterm :: Vector{Float64}
    a0dot :: Vector{Float64}
    adot :: Vector{Float64}
    a0prev :: Vector{Float64}
    aprev :: Vector{Float64}
    bv :: Vector{TwoDVort}
    lespcrit :: Vector{Float64}
    levflag :: Vector{Int8}

    function TwoDFreeSurf(c, uref, coord_file, pvt, ndiv, naterm, strpar, kinem ,lespcrit=zeros(1))
        theta = zeros(ndiv)
        x = zeros(ndiv)
        cam = zeros(ndiv)
        cam_slope = zeros(ndiv)
        bnd_x = zeros(ndiv)
        bnd_z = zeros(ndiv)

        dtheta = pi/(ndiv-1)
        for ib = 1:ndiv
            theta[ib] = real(ib-1.)*dtheta
            x[ib] = c/2.*(1-cos(theta[ib]))
        end
        if (coord_file != "FlatPlate")
            cam, cam_slope = camber_calc(x, coord_file)
        end

        for i = 1:ndiv
            bnd_x[i] = -((c - pvt*c)+((pvt*c - x[i])*cos(kinem.alpha))) + (cam[i]*sin(kinem.alpha))
            bnd_z[i] = kinem.h + ((pvt*c - x[i])*sin(kinem.alpha))+(cam[i]*cos(kinem.alpha))
        end
        uind = zeros(ndiv)
        wind = zeros(ndiv)
        downwash = zeros(ndiv)
        a0 = zeros(1)
        a0dot = zeros(1)
        aterm = zeros(naterm)
        adot = zeros(3)
        a0prev = zeros(1)
        aprev = zeros(3)
        bv = TwoDVort[]
        for i = 1:ndiv-1
            push!(bv,TwoDVort(0,0,0,0.02*c,0,0))
        end
        levflag = [0]
        new(c, uref, coord_file, pvt, ndiv, naterm, cam, cam_slope, theta, x, strpar, kinem, bnd_x, bnd_z, uind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv,lespcrit,levflag)
    end
end

#----------------------------------------------------------------------------------------------
# Definition of a 3D surface
immutable patch
    x :: Float64
    y :: Float64
    z :: Float64
    pvt :: Float64
    coord_file :: String
    c :: Float64
    twist :: Float64
    lc :: Float64
    nspan :: Int8
end

# immutable ThreeDSurf
#     cref :: Float64
#     bref :: Float64
#     sref :: Float64
#     uref :: Float64
#     patchdata :: Vector{patch}
#     ndiv :: Int8
#     nspan :: Int8
#     naterm :: Int8
#     kindef :: KinemDef3D
#     c :: Vector{Float64}
#     pvt :: Vector{Float64}
#     cam :: Array{Float64,2}
#     cam_slope :: Array{Float64,2}
#     theta :: Vector{Float64}
#     psi :: Vector{Float64}
#     x :: Array{Float64,2}
#     xle :: Vector{Float64}
#     yle :: Vector{Float64}
#     zle :: Vector{Float64}
#     kinem :: Vector{KinemPar}
#     bnd_x :: Array{Float64,2}
#     bnd_z :: Array{Float64,2}
#     uind :: Array{Float64,2}
#     vind :: Array{Float64,2}
#     wind :: Array{Float64,2}
#     downwash :: Array{Float64,2}
#     a0 :: Vector{Float64}
#     aterm :: Array{Float64,2}
#     a0dot :: Vector{Float64}
#     adot :: Array{Float64,2}
#     a0prev :: Vector{Float64}
#     aprev :: Array{Float64,2}
#     bv :: Array{ThreeDVort,2}
#     lespcrit :: Vector{Float64}
#     levflag :: Vector{Int8}

#     function ThreeDSurf(cref, bref, sref, patchdata, kindef, uref=1., ndiv=70, naterm=35)
#         nspan = 0
#         for i = 1:length(patchdata)-1
#             nspan += patchdata[i].nspan
#         end

#         lespcrit = zeros(nspan)

#         psi = zeros(nspan)
#         cam = zeros(nspan,ndiv)
#         cam_slope = zeros(nspan,ndiv)
#         cam1 = zeros(ndiv)
#         cam_slope1 = zeros(ndiv)
#         cam2 = zeros(ndiv)
#         cam_slope2 = zeros(ndiv)
#         bnd_x = zeros(nspan,ndiv)
#         bnd_z = zeros(nspan,ndiv)
#         c = zeros(nspan)
#         pvt = zeros(nspan)
# 	yle = zeros(nspan)
# 	xle = zeros(nspan)
# 	zle = zeros(nspan)
#        	twist = zeros(nspan)

#         kinem = KinemPar[]

# 	theta = zeros(ndiv)
#         x = zeros(nspan,ndiv)

#         dtheta = pi/(ndiv-1)
#         for ib = 1:ndiv
#             theta[ib] = real(ib-1.)*dtheta
# 	    for is = 1:nspan
#   	        x[is,ib] = c[is]/2.*(1-cos(theta[ib]))
# 	    end
# 	end

# 	nspan = 0
#         for i = 1:length(patchdata) - 1
#             if i != 1
#                 nspan += 1
#                 psi[nspan] = (patchdata[i].y - patchdata[1].y)*pi/(patchdata[length(patchdata)].y - patchdata[1].y)
#                 yle[nspan] = patchdata[i].y
#                 xle[nspan] = patchdata[i].x
#                 zle[nspan] = patchdata[i].z
#                 c[nspan] = patchdata[i].c
#                 twist[nspan] = patchdata[i].twist
#                 if (patchdata[i].coord_file != "FlatPlate")
#                     cam[nspan,:], cam_slope[nspan,:] = camber_calc(0.5*(1-cos(theta)), patchdata[i].coord_file)
#                 end
#                 lespcrit[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].lc, patchdata[i+1].lc, yle[nspan])
#                 pvt[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].pvt, patchdata[i+1].pvt, yle[nspan])
#             end
#             startpsi = (patchdata[i].y - patchdata[1].y)*pi/(patchdata[length(patchdata)].y - patchdata[1].y)
#             lenpsi = (patchdata[i+1].y - patchdata[i].y)*pi/(patchdata[length(patchdata)].y - patchdata[1].y)
#             dpsi =  lenpsi/(patchdata[i].nspan + 1.)
# #            divspan = (patchdata[i+1].y - patchdata[i].y)/patchdata[i].nspan
#             for j = 1:patchdata[i].nspan
#                 nspan += 1
#                 psi[nspan] = startpsi + real(j)*dpsi
#                 #                yle[nspan] = patchdata[i].y + (i-1)*divspan
#                 yle[nspan] = -bref*cos(psi[nspan])/2.
#                 xle[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].x, patchdata[i+1].x, yle[nspan])
#                 zle[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].z, patchdata[i+1].z, yle[nspan])
#                 c[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].c, patchdata[i+1].c, yle[nspan])
#                 twist[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].twist, patchdata[i+1].twist, yle[nspan])
#                 if (patchdata[i].coord_file != "FlatPlate")
#                     cam1, cam_slope1 = camber_calc(0.5*(1-cos(theta)), patchdata[i].coord_file)
#                 end
#                 if (patchdata[i+1].coord_file != "FlatPlate")
#                     cam2, cam_slope2 = camber_calc(0.5*(1-cos(theta)), patchdata[i+1].coord_file)
#                 end
#                 for k = 1:ndiv
#                     cam[nspan,k] = interp(patchdata[i].y, patchdata[i+1].y, cam1[k], cam2[k], yle[nspan])*c[nspan]
#                     cam_slope[nspan,k] = interp(patchdata[i].y, patchdata[i+1].y, cam_slope1[k], cam_slope2[k], yle[nspan])*c[nspan]
#                 end
#                 lespcrit[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].lc, patchdata[i+1].lc, yle[nspan])
#                 pvt[nspan] = interp(patchdata[i].y, patchdata[i+1].y, patchdata[i].pvt, patchdata[i+1].pvt, yle[nspan])
#             end
#         end

#         inkim = KinemPar(0, 0, 0, 0, 0, 0)
#         for i = 1:nspan
#             push!(kinem, inkim)
#         end

#         if (kindef.vartype == "Constant")
#             i = 1
#             if (typeof(kindef.alpha) == EldUpDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == EldRampReturnDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == ConstDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = 0.
#             elseif (typeof(kindef.alpha) == SinDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == CosDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             end

#             if (typeof(kindef.h) == EldUpDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == EldUpIntDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == EldRampReturnDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == ConstDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = 0.
#             elseif (typeof(kindef.h) == SinDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == CosDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             end

#             if (typeof(kindef.u) == EldUpDef)
#                 kinem[i].u = kindef.u(0.)*uref
#                 kinem[i].udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/cref
#             elseif (typeof(kindef.u) == EldRampReturnDef)
#                 kinem[i].u, kinem[i].udot = kindef.u(0.)
#                 kinem[i].u = kinem[i].u*uref
#                 kinem[i].udot = kinem[i].udot*uref*uref/cref
#             elseif (typeof(kindef.u) == ConstDef)
#                 kinem[i].u = kindef.u(0.)*uref
#                 kinem[i].udot = 0.
#             end

#             for i = 2:nspan
#                 kinem[i] = kinem[1]
#             end
#         elseif (typeof(kindef.vartype) == "Read")
#             # Provide the max amplitude as the input, and the read array goes from 0 to 1 as a function of yle
#             i = nspan
#             if (typeof(kindef.alpha) == EldUpDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == EldRampReturnDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == ConstDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = 0.
#             elseif (typeof(kindef.alpha) == SinDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             elseif (typeof(kindef.alpha) == CosDef)
#                 kinem[i].alpha = kindef.alpha(0.)
#                 kinem[i].alphadot = ForwardDiff.derivative(kindef.alpha,0.)*uref/cref
#             end

#             if (typeof(kindef.h) == EldUpDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == EldUpIntDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == EldRampReturnDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == ConstDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = 0.
#             elseif (typeof(kindef.h) == SinDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             elseif (typeof(kindef.h) == CosDef)
#                 kinem[i].h = kindef.h(0.)*cref
#                 kinem[i].hdot = ForwardDiff.derivative(kindef.h,0.)*uref
#             end

#             if (typeof(kindef.u) == EldUpDef)
#                 kinem[i].u = kindef.u(0.)*uref
#                 kinem[i].udot = ForwardDiff.derivative(kindef.u,0.)*uref*uref/cref
#             elseif (typeof(kindef.u) == EldRampReturnDef)
#                 kinem[i].u, kinem[i].udot = kindef.u(0.)
#                 kinem[i].u = kinem[i].u*uref
#                 kinem[i].udot = kinem[i].udot*uref*uref/cref
#             elseif (typeof(kindef.u) == ConstDef)
#                 kinem[i].u = kindef.u(0.)*uref
#                 kinem[i].udot = 0.
#             end

# addspl = Spline1D(kindef.add[:,1], kindef.add[:,2])

# if (kindef.vary == 1)

#     for i = 1:nspan
#         kinem.h[i] = kinem.h[nspan]*evaluate(addspl,yle[i])
#     end
# elseif (kindef.vary == 2)
#     addspl = Spline1D(kinem.add[:,1], kinem.add[:,2])
#     for i = 1:nspan
#         kinem.alpha[i] = kinem.alpha[nspan]*evaluate(addspl,yle[i])
#     end
# end

# end

# for j = 1:nspan
#     for i = 1:ndiv
#         bnd_x[j,i] = -((c[j] - pvt[j]*c[j])+((pvt[j]*c[j] - x[j,i])*cos(kinem[j].alpha))) + (cam[j,i]*sin(kinem[j].alpha))
#         bnd_z[j,i] = kinem[j].h + ((pvt[j]*c[j] - x[j,i])*sin(kinem[j].alpha))+(cam[j,i]*cos(kinem[j].alpha))
#     end
# end

# uind = zeros(nspan,ndiv)
# vind = zeros(nspan,ndiv)
# wind = zeros(nspan,ndiv)
# downwash = zeros(nspan,ndiv)
# a0 = zeros(nspan)
# a0dot = zeros(nspan)
# aterm = zeros(nspan,naterm)
# adot = zeros(nspan,3)
# a0prev = zeros(nspan)
# aprev = zeros(nspan,3)
# bv = Array(ThreeDVort,nspan,ndiv-1)

# for j = 1:nspan
#     for i = 1:ndiv-1
#         bv[j,i] = ThreeDVort([0; 0; 0;], [0; 0; 0;], 0.02*cref, [0; 0; 0;], [0; 0; 0;])
#     end
# end
# levflag = zeros(nspan)

# new(cref, bref, sref, uref, patchdata, ndiv, nspan, naterm, nbterm, kindef, c, pvt, cam, cam_slope, theta, psi, x, xle, yle, zle, kinem, bnd_x, bnd_z, uind, vind, wind, downwash, a0, aterm, a0dot, adot, a0prev, aprev, bv,lespcrit,levflag)
# end
# end

immutable ThreeDSurfSimple
    cref :: Float64
    AR :: Float64
    uref :: Float64
    pvt :: Float64
    lespcrit :: Vector{Float64}
    coord_file :: String
    ndiv :: Int8
    nspan :: Int8
    naterm :: Int8
    kindef :: KinemDef
    psi :: Vector{Float64}
    yle :: Vector{Float64}
    s2d :: Vector{TwoDSurf}
    a03d :: Vector{Float64}
    bc :: Vector{Float64}
    nshed :: Vector{Float64}

    function ThreeDSurfSimple(AR, kindef, coord_file, pvt, lespcrit = [10.;], nspan = 20, cref = 1., uref=1., ndiv=70, naterm=35)

        bref = AR*cref

        psi = zeros(nspan)
        yle = zeros(nspan)

        s2d = TwoDSurf[]

	for i = 1:nspan
            psi[i] = real(i)*(pi/2)/nspan
            yle[i] = -bref*cos(psi[i])/2.
        end

        for i = 1:nspan
            push!(s2d, TwoDSurf(coord_file, pvt,  kindef, lespcrit, cref, uref, ndiv, naterm))
        end

        a03d = zeros(nspan)
        bc = zeros(nspan)
        nshed = [0.;]

new(cref, AR, uref, pvt, lespcrit, coord_file,  ndiv, nspan, naterm, kindef, psi, yle, s2d, a03d, bc, nshed)
end
end

# ---------------------------------------------------------------------------------------------
# BEGIN EldUpIntDef
immutable EldUpIntDef <: MotionDef
    amp :: Float64
    K :: Float64
    a :: Float64
end

immutable EldUpInttstartDef <: MotionDef
    amp :: Float64
    K :: Float64
    a :: Float64
    tstart :: Float64
end

function (eld::EldUpIntDef)(t)

    dt = 0.015
    nsteps = t/dt + 1
    nsteps = round(Int,nsteps)
    dt = t/(nsteps-1)
    sm = pi*pi*eld.K/(2*eld.amp*(1 - eld.a))
    t1 = 1.
    t2 = t1 + ((eld.amp)/(2*eld.K))

    prev_h = 0
    amp = 0
    for i = 1:nsteps
      tmpt = (i-1)*dt
      if (eld.amp == 0.)
      	 hdot = 0.
      else
         hdot = ((eld.K/sm)*log(cosh(sm*(tmpt - t1))/cosh(sm*(tmpt - t2))))+(eld.amp/2.)
      end
      amp = prev_h + hdot*dt
      prev_h = amp
    end
    if (nsteps == 1)
      amp = 0.
    end
    amp
end

function (eld::EldUpInttstartDef)(t)

    dt = 0.015
    nsteps = t/dt + 1
    nsteps = round(Int,nsteps)
    dt = t/(nsteps-1)
    sm = pi*pi*eld.K/(2*eld.amp*(1 - eld.a))
    t1 = eld.tstart
    t2 = t1 + ((eld.amp)/(2*eld.K))

    prev_h = 0
    amp = 0
    for i = 1:nsteps
      tmpt = (i-1)*dt
      if (eld.amp == 0.)
      	 hdot = 0.
      else
         hdot = ((eld.K/sm)*log(cosh(sm*(tmpt - t1))/cosh(sm*(tmpt - t2))))+(eld.amp/2.)
      end
      amp = prev_h + hdot*dt
      prev_h = amp
    end
    if (nsteps == 1)
      amp = 0.
    end
    amp
end

# END EldUpIntDef
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# BEGIN KelvinCondition
immutable KelvinCondition
    surf :: TwoDSurf
    field :: TwoDFlowField
end

immutable KelvinConditionwFlap
    surf :: TwoDSurfwFlap
    field :: TwoDFlowField
end

immutable KelvinCondition2DOF
    surf :: TwoDSurf_2DOF
    field :: TwoDFlowField
    kelv_enf :: Float64
end

immutable KelvinCondition2DFree
    surf :: TwoDFreeSurf
    field :: TwoDFlowField
    kelv_enf :: Float64
end

immutable KelvinConditionLLTldvm
    surf :: ThreeDSurfSimple
    field :: ThreeDFieldSimple
end

function (kelv::KelvinCondition)(tev_iter::Array{Float64})
    #Update the TEV strength
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = tev_iter[1]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf, [kelv.field.u[1],kelv.field.w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val = val + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val = val + kelv.field.lev[iv].s
    end

    #Add kelv_enforced if necessary - merging will be better
    # val is the value of Gam_b + sum Gam_tev + Gam_lev which will equal zero
    # if the condition is satified

    return val
end

function (kelv::KelvinCondition2DOF)(tev_iter::Array{Float64})
    #Update the TEV strength
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = tev_iter[1]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val = val + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val = val + kelv.field.lev[iv].s
    end
    val = val + kelv.kelv_enf

    #Add kelv_enforced if necessary - merging will be better
    # val is the value of Gam_b + sum Gam_tev + Gam_lev which will equal zero
    # if the condition is satified

    return val
end

function (kelv::KelvinCondition2DFree)(tev_iter::Array{Float64})
    #Update the TEV strength
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = tev_iter[1]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val = val + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val = val + kelv.field.lev[iv].s
    end
    val = val + kelv.kelv_enf

    #Add kelv_enforced if necessary - merging will be better
    # val is the value of Gam_b + sum Gam_tev + Gam_lev which will equal zero
    # if the condition is satified

    return val
end

function (kelv::KelvinConditionwFlap)(tev_iter::Array{Float64})
    #Update the TEV strength
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = tev_iter[1]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val = val + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val = val + kelv.field.lev[iv].s
    end

    #Add kelv_enforced if necessary - merging will be better
    # val is the value of Gam_b + sum Gam_tev + Gam_lev which will equal zero
    # if the condition is satified

    return val
end

function (kelv::KelvinConditionLLTldvm)(tev_iter::Array{Float64})
    val = zeros(2*kelv.surf.nspan)
    lhs = zeros(kelv.surf.nspan,kelv.surf.nspan)
    rhs = zeros(kelv.surf.nspan)

    #Assume symmetry condition for now
    for i = 1:kelv.surf.nspan
        nlev = length(kelv.field.f2d[i].lev)
        ntev = length(kelv.field.f2d[i].tev)
        kelv.field.f2d[i].tev[ntev].s = tev_iter[i]

        #Update incduced velocities on airfoil
        update_indbound(kelv.surf.s2d[i], kelv.field.f2d[i])

        #Calculate downwash
        update_downwash(kelv.surf.s2d[i], [kelv.field.f2d[i].u[1], kelv.field.f2d[i].w[1]])

        #Calculate first two fourier coefficients
        update_a0anda1(kelv.surf.s2d[i])

        kelv.surf.bc[i] = kelv.surf.s2d[i].a0[1] + 0.5*kelv.surf.s2d[i].aterm[1]

        kelv.surf.a03d[i] = 0

        for n = 1:kelv.surf.nspan
            nn = 2*n - 1
            kelv.surf.a03d[i] = kelv.surf.a03d[i] - real(nn)*tev_iter[n+kelv.surf.nspan]*sin(nn*kelv.surf.psi[i])/sin(kelv.surf.psi[i])
        end

        val[i] = kelv.surf.s2d[i].uref*kelv.surf.s2d[i].c*pi*(kelv.surf.bc[i] + kelv.surf.a03d[i])

        for iv = 1:ntev
            val[i] = val[i] + kelv.field.f2d[i].tev[iv].s
        end
        for iv = 1:nlev
            val[i] = val[i] + kelv.field.f2d[i].lev[iv].s
        end
    end

    for i = 1:kelv.surf.nspan
        for n = 1:kelv.surf.nspan
            nn = 2*n - 1
            lhs[i,n] = sin(nn*kelv.surf.psi[i])*(sin(kelv.surf.psi[i]) + (nn*pi/(2*kelv.surf.AR)))
        end
        rhs[i] = pi*sin(kelv.surf.psi[i])*kelv.surf.bc[i]/(2*kelv.surf.AR)
    end

    val[kelv.surf.nspan+1:2*kelv.surf.nspan] = lhs*tev_iter[kelv.surf.nspan+1:2*kelv.surf.nspan] - rhs

    return val
end

# END KelvinCondition
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# BEGIN KelvinKutta
immutable KelvinKutta
    surf :: TwoDSurf
    field :: TwoDFlowField
end

immutable KelvinKutta2DOF
    surf :: TwoDSurf_2DOF
    field :: TwoDFlowField
    kelv_enf :: Float64
end

immutable KelvinKutta2DFree
    surf :: TwoDFreeSurf
    field :: TwoDFlowField
    kelv_enf :: Float64
end

immutable KelvinKuttawFlap
    surf :: TwoDSurfwFlap
    field :: TwoDFlowField
end

immutable KelvinKuttaLLTldvm
    surf :: ThreeDSurfSimple
    field :: Vector{TwoDFlowField}
 end

function (kelv::KelvinKutta)(v_iter::Array{Float64})
    val = zeros(2)

    #Update the TEV and LEV strengths
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = v_iter[1]
    kelv.field.lev[nlev].s = v_iter[2]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf ,[kelv.field.u[1],kelv.field.w[1]])

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val[1] = val[1] + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val[1] = val[1] + kelv.field.lev[iv].s
    end

    if (kelv.surf.a0[1] > 0)
        lesp_cond = kelv.surf.lespcrit[1]
    else
        lesp_cond = -kelv.surf.lespcrit[1]
    end
    val[2] = kelv.surf.a0[1]-lesp_cond

    #Add kelv_enforced if necessary - merging will be better
    return val
end

function (kelv::KelvinKutta2DOF)(v_iter::Array{Float64})
    val = zeros(2)

    #Update the TEV and LEV strengths
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = v_iter[1]
    kelv.field.lev[nlev].s = v_iter[2]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val[1] = val[1] + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val[1] = val[1] + kelv.field.lev[iv].s
    end
    val[1] = val[1] + kelv.kelv_enf

    if (kelv.surf.a0[1] > 0)
        lesp_cond = kelv.surf.lespcrit[1]
    else
        lesp_cond = -kelv.surf.lespcrit[1]
    end
    val[2] = kelv.surf.a0[1]-lesp_cond

    #Add kelv_enforced if necessary - merging will be better
    return val
end

function (kelv::KelvinKutta2DFree)(v_iter::Array{Float64})
    val = zeros(2)

    #Update the TEV and LEV strengths
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = v_iter[1]
    kelv.field.lev[nlev].s = v_iter[2]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val[1] = val[1] + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val[1] = val[1] + kelv.field.lev[iv].s
    end
    val[1] = val[1] + kelv.kelv_enf

    if (kelv.surf.a0[1] > 0)
        lesp_cond = kelv.surf.lespcrit[1]
    else
        lesp_cond = -kelv.surf.lespcrit[1]
    end
    val[2] = kelv.surf.a0[1]-lesp_cond

    #Add kelv_enforced if necessary - merging will be better
    return val
end

function (kelv::KelvinKuttawFlap)(v_iter::Array{Float64})
    val = zeros(2)

    #Update the TEV and LEV strengths
    nlev = length(kelv.field.lev)
    ntev = length(kelv.field.tev)
    kelv.field.tev[ntev].s = v_iter[1]
    kelv.field.lev[nlev].s = v_iter[2]

    #Update incduced velocities on airfoil
    update_indbound(kelv.surf, kelv.field)

    #Calculate downwash
    update_downwash(kelv.surf)

    #Calculate first two fourier coefficients
    update_a0anda1(kelv.surf)

    val[1] = kelv.surf.uref*kelv.surf.c*pi*(kelv.surf.a0[1] + kelv.surf.aterm[1]/2.)

    for iv = 1:ntev
        val[1] = val[1] + kelv.field.tev[iv].s
    end
    for iv = 1:nlev
        val[1] = val[1] + kelv.field.lev[iv].s
    end

    if (kelv.surf.a0[1] > 0)
        lesp_cond = kelv.surf.lespcrit[1]
    else
        lesp_cond = -kelv.surf.lespcrit[1]
    end
    val[2] = kelv.surf.a0[1]-lesp_cond

    #Add kelv_enforced if necessary - merging will be better
    return val
end

function (kelv::KelvinKuttaLLTldvm)(tev_iter::Array{Float64})
    val = zeros(2*kelv.surf3d.nspan + kelv.nshed)
    bc = zeros(kelv.surf3d.nspan)
    AR = kelv.surf3d.bref/kelv.surf3d.cref
    lhs = zeros(kelv.surf3d.nspan,kelv.surf3d.nspan)
    rhs = zeros(kelv.surf3d.nspan)
    bcoeff = zeros(kelv.surf3d.nspan)
    a03d = zeros(kelv.surf3d.nspan)

    cntr = kelv.surf3d.nspan + 1

    #Assume symmetry condition for now
    for i = 1:kelv.surf3d.nspan
        nlev = length(kelv.field[i].lev)
        ntev = length(kelv.field[i].tev)
        kelv.field[i].tev[ntev].s = tev_iter[i]
        if kelv.surf[i].levflag[1] == 1
            kelv.field[i].lev[nlev].s = tev_iter[cntr]
            cntr += 1
        end

        #Update incduced velocities on airfoil
        update_indbound(kelv.surf[i], kelv.field[i])

        #Calculate downwash
        update_downwash(kelv.surf[i], [kelv.field[i].u[1],kelv.field[i].w[1]])

        #Calculate first two fourier coefficients
        update_a0anda1(kelv.surf[i])

        bc[i] = kelv.surf[i].a0[1] + 0.5*kelv.surf[i].aterm[1]

        for i = 1:kelv.surf3d.nspan
        a03d[i] = 0
            for n = 1:kelv.surf3d.nspan
                nn = 2*n - 1
                a03d[i] = a03d[i] - real(nn)*tev_iter[n+kelv.surf3d.nspan+kelv.nshed]*sin(nn*kelv.surf3d.psi[i])/sin(kelv.surf3d.psi[i])
            end
        end

        val[i] = kelv.surf[i].uref*kelv.surf[i].c*pi*(bc[i] + a03d[i])

        for iv = 1:ntev
            val[i] = val[i] + kelv.field[i].tev[iv].s
        end
        for iv = 1:nlev
            val[i] = val[i] + kelv.field[i].lev[iv].s
        end

        if kelv.surf[i].levflag[1] == 1
            if kelv.surf[i].a0[1] > 0
                lesp_cond = kelv.surf[i].lespcrit[1]
            else
                lesp_cond = -kelv.surf[i].lespcrit[1]
            end
            val[cntr-1] = kelv.surf[i].a0[1] + a03d[i] - lesp_cond
        end
    end

    for i = 1:kelv.surf3d.nspan
        for n = 1:kelv.surf3d.nspan
            nn = 2*n - 1
            lhs[i,n] = sin(nn*kelv.surf3d.psi[i])*(sin(kelv.surf3d.psi[i]) + (nn*pi/(2*AR)))
        end
        rhs[i] = pi*sin(kelv.surf3d.psi[i])*bc[i]/(2*AR)
    end

    val[kelv.surf3d.nspan+kelv.nshed+1:2*kelv.surf3d.nspan+kelv.nshed] = lhs*tev_iter[kelv.surf3d.nspan+kelv.nshed+1:2*kelv.surf3d.nspan+kelv.nshed] - rhs

    return val
end

# END KelvinKutta
# ---------------------------------------------------------------------------------------------


text:
using LinearAlgebra
using Transducers
n = 3
q = 4*rand(Float64)
r = 2*rand(Float64)
k = 10000
actions = 1:k

function empirical_var(diffs,n,k)
    sum(sum([d .^ 2 for d in diffs]))/(n*k)
end

lssm = LinearStateSpaceModel(Matrix{Float64}(I, n, n), q .* Matrix{Float64}(I, n, n))
states = collect(Simulator(lssm), actions)
lobs = LinearObservation(Matrix{Float64}(I, n, n), r .* Matrix{Float64}(I, n, n))
measures = states |> Map(s -> observe(s,lobs)) |> collect


text:
struct Dimensions
    q::Int         # configuration
    u::Int         # control
	w::Int         # disturbance
	c::Int         # contact points
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Glimpse
const Gl = Glimpse
import Glimpse: RGBAf0
using GeometryTypes
using LinearAlgebra
import LinearAlgebra: norm
using ColorSchemes

cscheme = ColorScheme([RGBAf0(1.0, 0.0, 0.0, 0.6), RGBAf0(0.0, 0.0, 0.0, 0.6),
                       RGBAf0(0.0, 0.0, 1.0, 0.6)], "custom", "red black blue")
sz = [1 0; 0 -1]
ψ = p -> [p[1] * p[2], (p[1] + 1im * p[2]) * p[3] / sqrt(2)] .* ℯ^-norm(p)
function densfunc(p, orb)
    ψ = orb(p)
    return convert(Float32, real(ψ' * ψ))
end

function colorfunc(p, orb, cscheme)
    ψ = orb(p)
    if norm(ψ) == 0
        return RGBAf0(0.0, 0.0, 0.0, 0.6)
    else
        c = get(cscheme, 0.5 + 0.5 * real(ψ' * sz * ψ / (ψ' * ψ)))
        return c::RGBAf0
    end
end
grid = Gl.Grid([Point3f0(a, b, c) for a in -30:0.3:30, b in -30:0.3:30, c in -30:0.3:30])

ccomp = Gl.FunctionColor(p -> colorfunc(p, ψ, cscheme))
dcomp = Gl.FunctionGeometry(p -> densfunc(p, ψ), 0.001f0)

dio = Diorama(; background = RGBAf0(0.0, 0.0, 0.0, 1.0))
Gl.new_entity!(dio;
               separate = [Gl.Spatial(), Gl.Material(), Gl.Shape(), ccomp, dcomp,
                           Gl.ProgramTag{Gl.PeelingProgram}()],
               shared = Gl.ComponentData[grid])
dio.loop = @async Gl.renderloop(dio)


text:
using Plots
using LaTeXStrings
#
using ThreeBodyDecay
#
# decay Λb ⟶ Jψ p K
ms = (Jψ = 3.09, p=0.938, K = 0.49367, Lb = 5.62) # masses of the particles
# create two-body system
tbs = ThreeBodySystem(ms.Jψ, ms.p, ms.K; m0=ms.Lb,   # masses m1,m2,m3,m0
            two_jps=([    1, 1//2,    0,  1//2] .|> x2,  # twice spin
                     [  '-',  '+',  '-',   '+'])) # parities
#
# chain-1, i.e. (2+3): Λs with the lowest ls, LS
Λ1520  = DecayChainLS(1, σ->BW(σ, 1.5195, 0.017); two_s = 3/2|>x2, tbs=tbs)
Λ1690  = DecayChainLS(1, σ->BW(σ, 1.685,  0.050); two_s = 1/2|>x2, tbs=tbs)
Λ1810  = DecayChainLS(1, σ->BW(σ, 1.80,   0.090); two_s = 5/2|>x2, tbs=tbs)
Λs = (Λ1520,Λ1690,Λ1810)
#
# chain-3, i.e. (1+2): Pentaquarks with the lowest ls, LS
Pc4312 = DecayChainLS(3, σ->BW(σ, 4.312, 0.015); two_s = 1/2|>x2, tbs=tbs)
Pc4440 = DecayChainLS(3, σ->BW(σ, 4.440, 0.010); two_s = 1/2|>x2, tbs=tbs)
Pc4457 = DecayChainLS(3, σ->BW(σ, 4.457, 0.020); two_s = 3/2|>x2, tbs=tbs)
Pcs = (Pc4312,Pc4440,Pc4457)
#
A(σs,two_λs,cs) = sum(c*amplitude(σs,two_λs,dc) for (c, dc) in zip(cs, (Λs...,Pcs...)))
I(σs,cs) = sum(abs2(A(σs,two_λs,cs)) for two_λs in itr(tbs.two_js))
#
Ini(σs,cs) = sum(abs2(sum(c*amplitude(σs,two_λs,dc) for (c, dc) in zip(cs, Λs ))) for two_λs in itr(tbs.two_js)) +
             sum(abs2(sum(c*amplitude(σs,two_λs,dc) for (c, dc) in zip(cs, Pcs))) for two_λs in itr(tbs.two_js))
#
I(σs)   = I(  σs,[1, 1.1, 3.9im, 2.2, 2.1im, -0.3im]) # set the couplings
Ini(σs) = Ini(σs,[1, 1.1, 3.9im, 2.2, 2.1im, -0.3im]) # set the couplings
#
dpp = randomPoint(tbs) # just a random point of the Dalitz Plot
@show I(dpp.σs)
@show Ini(dpp.σs)

let
    # generate data
    plot(size=(900,350), layout=grid(1,2),
        xlab=L"\cos\,\theta_{12}")
    #
    σ1v,σ2v,σ3v = flatDalitzPlotSample(tbs; Nev = 10000)
    cosθ12v = [cosθ12(σs,tbs.msq) for σs in zip(σ1v,σ2v,σ3v)]
    #
    weights    = [I(  σs) for σs in zip(σ1v,σ2v,σ3v)]
    weights_ni = [Ini(σs) for σs in zip(σ1v,σ2v,σ3v)]
    #
    plot!(sp=1, cosθ12v, bins=60, seriestype=:stephist, lab="ph.ps.")
    plot!(sp=2, cosθ12v, weights=weights, bins=60, seriestype=:stephist, lab="|L+P|^2")
    plot!(sp=2, cosθ12v, weights=weights_ni, bins=60, seriestype=:stephist, lab="|L|^2+|P|^2")
end

max_stars_repo_name: mmikhasenko/ThreeBodyDecay

max_issues_repo_name: mmikhasenko/ThreeBodyDecay

max_forks_repo_name: mmikhasenko/ThreeBodyDecay