text:
using Documenter, AlphaShapes

makedocs(sitename="AlphaShapes.jl",format = Documenter.HTML(prettyurls = true))

deploydocs(
    repo = "github.com/harveydevereux/AlphaShapes.jl.git"
)


text:
"""
`init_subspec!(m::Model1002)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::Model1002)
    if subspec(m) == "ss2"
        return
    elseif subspec(m) == "ss8"
        return ss8!(m)
    elseif subspec(m) == "ss9"
        return ss9!(m)
    elseif subspec(m) == "ss10"
        return ss10!(m)
    elseif subspec(m) == "ss11"
        # Normalize s.t. sum of variances of anticipated shocks equals the
        # variance of the contemporaneous shock
        # See measurement equation
        return ss10!(m)
    else
        error("This subspec is not defined.")
    end
end

"""
```
ss8!(m::Model1002)
```

Initializes subspec 8 of `Model1002`. In this subspecification, γ_gdi
and δ_gdi are fixed at 1 and 0 respectively.
"""
function ss8!(m::Model1002)
    m <= parameter(:γ_gdi, 1., (-10., 10.), (-10., -10.), Untransformed(), Normal(1., 2.), fixed = true,
                   tex_label = "\\gamma_{gdi}")

    m <= parameter(:δ_gdi, 0., (-10., 10.), (-10., -10.), Untransformed(), Normal(0.00, 2.), fixed = true,
                   tex_label = "\\delta_{gdi}")
end

"""
```
ss9!(m::Model1002)
```

Initializes subspec 9 of `Model1002`. This subspecification is ss8 + Iskander's
changes + the bounds for beta-distributed parameters have been changed from
(1e-5, 0.99) to (0.0, 1.0).
"""
function ss9!(m::Model1002)

    ss8!(m)

    m <= parameter(:ζ_p, 0.8940, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.1), fixed = false,
                   description = "ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label = "\\zeta_p")

    m <= parameter(:ι_p, 0.1865, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",
                   tex_label = "\\iota_p")

    m <= parameter(:h, 0.5347, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.7, 0.1), fixed = false,
                   description = "h: Consumption habit persistence.",
                   tex_label = "h")

    m <= parameter(:ppsi, 0.6862, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ppsi: Utilization costs.",
                   tex_label = "\\psi")

    m <= parameter(:ζ_w, 0.9291, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.1), fixed = false,
                   description = "ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label = "\\zeta_w")

    m <= parameter(:ι_w, 0.2992, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.15), fixed = false,
                   description = "ι_w: The weight attributed to last period's wage in wage indexation. (1-ι_w) is the weight attributed to steady-state wages.",
                   tex_label = "\\iota_w")

    m <= parameter(:ρ, 0.7126, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.10), fixed = false,
                   description = "ρ: The degree of inertia in the monetary policy rule.",
                   tex_label = "\\rho_R")

    # Financial frictions parameters
    m <= parameter(:Fω, 0.0300, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.03, 0.01), fixed = true,
                   scaling = x -> 1 - (1-x)^0.25,
                   description = "F(ω): The cumulative distribution function of ω (idiosyncratic iid shock that increases or decreases entrepreneurs' capital).",
                   tex_label = "F(\\bar{\\omega})")

    m <= parameter(:ζ_spb, 0.0559, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.05, 0.005), fixed = false,
                   description = "ζ_spb: The elasticity of the expected exess return on capital (or 'spread') with respect to leverage.",
                   tex_label = "\\zeta_{sp,b}")

    m <= parameter(:γ_star, 0.9900, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.99, 0.002),
                   fixed = true,
                   description = "γ_star: Fraction of entrepreneurs who survive and continue operating for another period.",
                   tex_label = "\\gamma_*")

    m <= parameter(:ρ_g, 0.9863, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label = "\\rho_g")

    m <= parameter(:ρ_b, 0.9410, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label = "\\rho_b")

    m <= parameter(:ρ_μ, 0.8735, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label = "\\rho_{\\mu}")

    m <= parameter(:ρ_z, 0.9446, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_z: AR(1) coefficient in the technology process.",
                   tex_label = "\\rho_z")

    m <= parameter(:ρ_λ_f, 0.8827, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w, 0.3884, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.",
                   tex_label = "\\rho_{\\lambda_w}")

    m <= parameter(:ρ_rm, 0.2135, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_rm: AR(1) coefficient in the monetary policy shock process.",
                   tex_label = "\\rho_{r^m}")

    m <= parameter(:ρ_σ_w, 0.9898, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = false,
                   description = "ρ_σ_w: The standard deviation of entrepreneurs' capital productivity follows an exogenous process with mean ρ_σ_w. Innovations to the process are called _spread shocks_.",
                   tex_label = "\\rho_{\\sigma_\\omega}")

    m <= parameter(:ρ_μ_e, 0.7500, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = true,
                   description = "ρ_μ_e: AR(1) coefficient in the exogenous bankruptcy cost process.",
                   tex_label = "\\rho_{\\mu_e}")

    m <= parameter(:ρ_γ, 0.7500, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.75, 0.15), fixed = true,
                   description = "ρ_γ: AR(1) coefficient in the process describing the fraction of entrepreneurs surviving period t.",
                   tex_label = "\\rho_{\\gamma}")

    m <= parameter(:ρ_π_star, 0.9900, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = true,
                   description = "ρ_π_star: AR(1) coefficient in the time-varying inflation target process.",
                   tex_label = "\\rho_{\\pi_*}")

    m <= parameter(:ρ_lr, 0.6936, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{10y}")

    m <= parameter(:ρ_z_p, 0.8910, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   description = "ρ_z_p: AR(1) coefficient in the process describing the permanent component of productivity.",
                   tex_label = "\\rho_{z^p}")

    m <= parameter(:ρ_tfp, 0.1953, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{tfp}")

    m <= parameter(:ρ_gdpdef, 0.5379, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{gdpdef}")

    m <= parameter(:ρ_corepce, 0.2320, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.5, 0.2), fixed = false,
                   tex_label = "\\rho_{pce}")

    m <= parameter(:η_gz, 0.8400, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   description = "η_gz: Correlate g and z shocks.",
                   tex_label = "\\eta_{gz}")

    m <= parameter(:η_λ_f, 0.7892, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   tex_label = "\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w, 0.4226, (0.0, 1.0), (0.0, 1.0), SquareRoot(), BetaAlt(0.50, 0.20), fixed = false,
                   description = "η_λ_w: AR(2) coefficient on wage markup shock process.",
                   tex_label = "\\eta_{\\lambda_w}")

    m <= parameter(:Iendoα, 0.0000, (0.0, 1.0), (0.0, 0.0), Untransformed(), BetaAlt(0.50, 0.20), fixed = true,
                   description = "Iendoα: Indicates whether to use the model's endogenous α in the capacity utilization adjustment of total factor productivity.",
                   tex_label = "I\\{\\alpha^{model}\\}")
end

"""
```
ss10!(m::Model1002)
```

Initializes subspec 10 of `Model1002`. This subspecification is the same as ss9,
except `betabar` is defined with `m[:σ_c]` instead of `σ_ω_star`.
"""
function ss10!(m::Model1002)
    ss9!(m)
end

ext: jl


text:
# This file was generated, do not modify it. # hide
m = machine(tm, X, y)
fit!(m)


text:
import Base.+

function +(buffer::Array{String,1},content::String)
    push!(buffer,content)
end

max_stars_repo_name: varnerlab/CBModelTools

max_issues_repo_name: varnerlab/CBModelTools

max_forks_repo_name: varnerlab/CBModelTools


text:
let
	a = VecSE2()
	@test typeof(a) <: VecSE
	@test typeof(a) <: AbstractVec
	@test a.x == 0.0
	@test a.y == 0.0
	@test a.θ == 0.0

	b = VecSE2(0.0,0.0,0.0)
	@test b.x == 0.0
	@test b.y == 0.0
	@test b.θ == 0.0

	c = VecSE2(VecE2(0.0,0.0),0.0)
	@test c.x == 0.0
	@test c.y == 0.0
	@test c.θ == 0.0

	@test length(a) == 3
	@test a == b
	@test b == c
	@test a == c
	@test isequal(b, a)
	@test copy(a) == a
	@test convert(Vector{Float64}, a) == [0.0,0.0,0.0]
	@test convert(VecSE2, [1.0,2.0,3.0]) == VecSE2(1.0,2.0,3.0)
	@test convert(VecE2, VecSE2(1.0,2.0,3.0)) == VecE2(1.0,2.0)

	@test isapprox(polar( 1.0, 0.0, 0.5), VecSE2(1.0,0.0,0.5))
	@test isapprox(polar( 2.0, π/2,-0.5), VecSE2(0.0,2.0,-0.5))
	@test isapprox(polar( 3.0,-π/2,   1), VecSE2(0.0,-3.0,1.0))
	@test isapprox(polar(-0.5,  1π, -1), VecSE2(0.5,0.0,-1))

	a = VecSE2(0.0,1.0,π/2)
	b = VecSE2(0.5,2.0,-π)
	unit_e2 = VecE2(1.0, 1.0)

	@test a != b

	@test a + b == VecSE2(0.5, 3.0, -π/2)
	@test a + 1*unit_e2 == VecSE2(1.0, 2.0,  π/2)
	@test a + 0.5*unit_e2 == VecSE2(0.5, 1.5, π/2)
	@test 1*unit_e2 + a == VecSE2(1.0, 2.0,  π/2)
	@test 0.5*unit_e2 + a == VecSE2(0.5, 1.5, π/2)

	@test a - b == VecSE2(-0.5, -1.0, 3π/2)
	@test a - 1*unit_e2 == VecSE2(-1.0,  0.0,  π/2)
	@test a - 0.5*unit_e2 == VecSE2(-0.5, 0.5, π/2)
	@test 1*unit_e2 - a == VecSE2(1.0,  0.0, -π/2)
	@test 0.5*unit_e2 - a == VecSE2(0.5, -0.5, -π/2)

	@test scale_euclidean(a, 2) == VecSE2(0.0, 2.0, π/2)
	@test scale_euclidean(a, 0.5) == VecSE2(0.0, 0.5, π/2)

	@test scale_euclidean(a, 1/2) == VecSE2(0.0, 0.5, π/2)
	@test scale_euclidean(a, 1/0.5) == VecSE2(0.0, 2.0, π/2)

	@test clamp_euclidean(VecSE2(1.0, 10.0, 0.5), 0.0, 5.0) == VecSE2(1.0, 5.0, 0.5)
	@test clamp_euclidean(VecSE2(-1.0, 4.0, 5.5), 0.0, 5.0) == VecSE2(0.0, 4.0, 5.5)

	@test isfinite(a)
	@test !isfinite(VecSE2(Inf,0))
	@test !isfinite(VecSE2(0,Inf))
	@test !isfinite(VecSE2(0,-Inf))
	@test !isfinite(VecSE2(Inf,Inf))

	@test !isinf(a)
	@test isinf(VecSE2(Inf,0))
	@test isinf(VecSE2(0,Inf))
	@test isinf(VecSE2(0,-Inf))
	@test isinf(VecSE2(Inf,Inf))

	@test !isnan(a)
	@test isnan(VecSE2(NaN,0))
	@test isnan(VecSE2(0,NaN))
	@test isnan(VecSE2(NaN,NaN))

	c = VecSE2(3.0,4.0)

	@test isapprox(norm(VecE2(a)), 1.0)
	@test isapprox(norm(VecE2(b)), hypot(0.5, 2.0))
	@test isapprox(norm(VecE2(c)), hypot(3.0, 4.0))

	@test isapprox(normalize(VecE2(a)), VecE2(0.0,1.0))
	@test isapprox(normalize(VecE2(b)), VecE2(0.5/sqrt(4.25),2.0/sqrt(4.25)))
	@test isapprox(normalize(VecE2(c)), VecE2(3.0/5,4.0/5))
end


text:
module Arduino

import GetC.@getCFun

const ardLib = "libarduino"

@getCFun ardLib init serialport_init(serialport::Ptr{Uint8}, baud::Int32)::Int32
export init
@getCFun ardLib writebyte serialport_writebyte(fd::Int32, b::Uint8)::Int32
export writebyte
@getCFun ardLib write serialport_write(fd::Int32, str::Ptr{Uint8})::Int32
export write
@getCFun ardLib read_until serialport_read_until(fd::Int32, buf::Ptr{Uint8}, until::Char, buf_max::Int32, timeout::Int32)::Int32
export read_until
@getCFun ardLib flush serialport_flush(fd::Int32)::Int32
export flush
@getCFun ardLib close serialport_close(fd::Int32)::Int32
export close

end


text:
module ACEtb

using ACE, JuLIP, NeighbourLists, Reexport

include("bonds.jl")
@reexport using ACEtb.Bonds

include("SlaterKoster.jl")
@reexport using ACEtb.SlaterKoster

include("tbhelpers.jl")
@reexport using ACEtb.TBhelpers

include("utils.jl")
@reexport using ACEtb.Utils

include("predictions.jl")
@reexport using ACEtb.Predictions

include("bindings.jl")
@reexport using ACEtb.Bindings

include("fittools.jl")
@reexport using ACEtb.FitTools

include("analysistools.jl")
@reexport using ACEtb.AnalysisTools

end


text:
export Localization, LocElem, Loc

###############################################################################

#   Declaration types
#   Loc / LocElem
#
###############################################################################

import AbstractAlgebra: base_ring, parent, check_parent, isunit, inv, +, -, *, divexact, zero, Base.promote_rule,
       elem_type, parent_type, Base.one, Base.iszero, Base.isone, Base.==, Base.gcd, Base.deepcopy_internal, needs_parentheses, valuation, Base.show,
       displayed_with_minus_in_front, show_minus_one, Base.^, data, Base.numerator, Base.denominator, canonical_unit, Base.gcdx, Base.div, divides,
       Base.lcm, Base.rand

import Nemo.prime

#prime might be product of several primes if localized at several primes, those primes are in array primes
mutable struct Loc{T<: RingElem} <: AbstractAlgebra.Ring
   base_ring::AbstractAlgebra.Ring
   prime::T
   primes::Array{T,1}

   function Loc{T}(prime::T, primes::Array{T,1}, cached::Bool = true) where {T <: RingElem}
      length(primes) == 0 && error("No element to localize at since array of primes is empty")
      if cached && haskey(LocDict, (parent(prime), prime))
         return LocDict[parent(prime), prime]::Loc{T}
      else
         z = new(parent(prime), prime, primes)
         if cached
            LocDict[parent(prime), prime] = z
         end
         return z
      end
   end
end

const LocDict = Dict{Tuple{AbstractAlgebra.Ring, RingElement}, AbstractAlgebra.Ring}()

mutable struct LocElem{T<: RingElem} <: AbstractAlgebra.RingElem
   data::FracElem{T}
   parent::Loc{T}
   function LocElem{T}(data::FracElem{T}, par::Loc, checked::Bool = true) where {T <: RingElem}
      checked && !isone(gcd(denominator(data), prime(par))) && error("No valid element of localization since its denominator and the primes localized at are not coprime")
      return new{T}(data,par)
   end
end

###############################################################################
#
#   Unsafe operators and functions
#
###############################################################################

# TODO (easy): Do this properly
AbstractAlgebra.add!(c::LocElem, a::LocElem, b::LocElem) = a + b

AbstractAlgebra.mul!(c::LocElem, a::LocElem, b::LocElem) = a * b

AbstractAlgebra.addeq!(a::LocElem, b::LocElem) = a + b

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

elem_type(::Type{Loc{T}}) where {T <: RingElem} = LocElem{T}

parent_type(::Type{LocElem{T}}) where {T <: RingElem} = Loc{T}

base_ring(L::Loc{T})  where {T <: RingElem}  = L.base_ring::parent_type(T)

base_ring(a::LocElem{T}) where {T <: RingElem} = base_ring(parent(a))

parent(a::LocElem{T})  where {T <: RingElem} = a.parent

function check_parent(a::LocElem{T}, b::LocElem{T})  where {T <: RingElem}
    parent(a) !== parent(b) && error("Parent objects do not match")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

data(a::LocElem{T}) where {T <: RingElement} = a.data

numerator(a::LocElem{T}) where {T <: RingElement} = numerator(data(a))

denominator(a::LocElem{T}) where {T <: RingElement} = denominator(data(a))

prime(L::Loc{T}) where {T <: RingElement} = L.prime

primes(L::Loc{T}) where {T <: RingElement} = L.primes

zero(L::Loc{T}) where {T <: RingElement} = L(0)

one(L::Loc{T}) where {T <: RingElement} = L(1)

iszero(a::LocElem{T}) where {T <: RingElement} = iszero(data(a))

isone(a::LocElem{T}) where {T <: RingElement} = isone(data(a))

function isunit(a::LocElem{T})  where {T <: RingElem}
    return isone(gcd(numerator(data(a)),prime(parent(a))))
end

deepcopy_internal(a::LocElem{T}, dict::IdDict) where {T <: RingElem} = parent(a)(deepcopy(data(a)))

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::LocElem{T}) where {T <: RingElem}
   print(io, data(a))
end

function show(io::IO, L::Loc{T}) where {T <: RingElem}
   print(io, "Localization of ", base_ring(L), " at ", primes(L))
end

needs_parentheses(x::LocElem{T})  where {T <: RingElem} = needs_parentheses(data(x))

displayed_with_minus_in_front(x::LocElem{T})  where {T <: RingElem} = displayed_with_minus_in_front(data(x))

show_minus_one(::Type{LocElem{T}}) where {T <: RingElement} = true

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(a::LocElem{T})  where {T <: RingElem}
   parent(a)(-data(a))
end

###############################################################################
#
#   Binary operators
#
###############################################################################

function +(a::LocElem{T}, b::LocElem{T})  where {T <: RingElem}
   check_parent(a,b)
   return LocElem{T}(data(a) + data(b), parent(a), false)
end

function -(a::LocElem{T}, b::LocElem{T})  where {T <: RingElem}
   check_parent(a,b)
   return LocElem{T}(data(a) - data(b), parent(a), false)
end

function *(a::LocElem{T}, b::LocElem{T})  where {T <: RingElem}
   check_parent(a,b)
   return LocElem{T}(data(a) * data(b), parent(a), false)
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
   check_parent(a, b)
   return data(a) == data(b)
end

###############################################################################
#
#   Inversion
#
###############################################################################

@doc Markdown.doc"""
     inv(a::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
Returns the inverse element of $a$ if $a$ is a unit.
If 'checked = false' the invertibility of $a$ is not checked and the corresponding inverse element
of the Fraction Field is returned.
"""
function inv(a::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
   checked && !isunit(a) && error("$a not invertible in given localization")
   return LocElem{T}(inv(data(a)), parent(a), false)
end

###############################################################################
#
#   Exact division
#
###############################################################################

@doc Markdown.doc"""
     divides(a::LocElem{T}, b::LocElem{T}, checked::Bool = true) where {T <: RingElem}
Returns tuple (`true`,`c`) if $b$ divides $a$ where `c`*$b$ = $a$.
If 'checked = false' the corresponding element of the Fraction Field is returned and it is not
checked whether it is an element of the given localization.
"""
function divides(a::LocElem{T}, b::LocElem{T}, checked::Bool = true) where {T <: RingElem}
   check_parent(a,b)
   try
      true, parent(a)(divexact(data(a), data(b)), checked)
   catch
      false, parent(a)()
   end
end

@doc Markdown.doc"""
     divexact(a::LocElem{T}, b::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
Returns element 'c' of given localization s.th. `c`*$b$ = $a$ if such element exists.
If 'checked = false' the corresponding element of the Fraction Field is returned and it is not
checked whether it is an element of the given localization.
"""
function divexact(a::LocElem{T}, b::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
   d = divides(a, b, checked)
   d[1] ? d[2] : error("$a not divisible by $b in the given Localization")
end

@doc Markdown.doc"""
     div(a::LocElem{T}, b::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
Returns element `c` if $b$ divides $a$ where `c`* $b$ = $a$.
If $b$ does not divide $a$, `0`is returned.
If 'checked = false' the corresponding element of the Fraction Field is returned and it is not
checked whether it is an element of the given localization.
"""
function div(a::LocElem{T}, b::LocElem{T}, checked::Bool = true)  where {T <: RingElem}
   d = divides(a, b, checked)
   return d[2]
end

###############################################################################
#
#   GCD & LCM
#
###############################################################################

@doc Markdown.doc"""
    gcd(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
Returns gcd of $a$ and $b$ in canonical representation.
"""
function gcd(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
   check_parent(a,b)
   iszero(a) && return canonical_unit(b) * b
   iszero(b) && return canonical_unit(a) * a
   par = parent(a)
   elem = prod([pi^min(valuation(a,pi), valuation(b,pi)) for pi in primes(par)])
   return par(elem)
end

@doc Markdown.doc"""
    lcm(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
Returns lcm of $a$ and $b$ in canonical representation.
"""
function lcm(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
   check_parent(a,b)
   par = parent(a)
   (iszero(a) || iszero(b)) && return par()
   elem = prod([pi^max(valuation(a,pi), valuation(b,pi)) for pi in primes(par)])
   return par(elem)
end

###############################################################################
#
#   GCDX
#
###############################################################################

@doc Markdown.doc"""
    gcdx(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
Returns tuple `(g,u,v)` s.th. `g` = gcd($a$,$b$) and `g` = `u` * $a$ + `v` * $b$.
Requires method gcdx for ring that is localized.
"""
function gcdx(a::LocElem{T}, b::LocElem{T}) where {T <: RingElement}
   check_parent(a,b)
   par = parent(a)
   g = gcd(a,b)
   can_a = canonical_unit(a)
   can_b = canonical_unit(b)
   a_temp = numerator(divexact(can_a * a, g, false))
   b_temp = numerator(divexact(can_b * b, g, false))
   (_,u,v) = gcdx(a_temp, b_temp)
   return (g, par(u)*can_a, par(v)*can_b)
end

###############################################################################
#
#   Random Functions
#
###############################################################################

#mainly for testing
function rand(L::Loc{T}, num_scale = (1:1000), den_scale=(1:1000)) where {T <: fmpz}
   num = rand(num_scale)
   den = rand(den_scale)
   while gcd(den,prime(L)) != 1
      den = rand(den_scale)
   end
   return L(num//den)
end
###############################################################################
#
#   Powering
#
###############################################################################

function ^(a::LocElem{T}, b::Int) where {T <: RingElem}
   return LocElem{T}(data(a)^b, parent(a), false)
end

###############################################################################
#
#   Promotion rules
#
###############################################################################

promote_rule(::Type{LocElem{T}}, ::Type{LocElem{T}}) where {T <: RingElement} = LocElem{T}

###############################################################################
#
#   Parent object call overloading
#
###############################################################################

(L::Loc{T})() where {T <: RingElement} = L(zero(base_ring(L)))

(L::Loc{T})(a::Integer)  where {T <: RingElem} = L(base_ring(L)(a))

function (L::Loc{T})(data::FracElem{T}, checked::Bool = true) where {T <: RingElem}
   return LocElem{T}(data,L,checked)
end

function (L::Loc{T})(data::Rational{<: Integer}, checked::Bool = true) where {T <: RingElem}
   return LocElem{T}(base_ring(L)(numerator(data)) // base_ring(L)(denominator(data)),L,checked)
end

function (L::Loc{T})(data::T, checked::Bool = false) where {T <: RingElem}
   return LocElem{T}(data // parent(data)(1),L,checked)
end

function (L::Loc{T})(a::LocElem{T}) where {T <: RingElement}
   L != parent(a) && error("No element of $L")
   return a
end

################################################################################
#
#   Valuation
#
################################################################################

@doc Markdown.doc"""
    valuation(a::LocElem{T}) where {T <: RingElement}
Returns the valuation of $a$ at the prime localized at. Only implemented for localizations at one prime ideal.
"""
function valuation(a::LocElem{T}) where {T <: RingElement}
   length(primes(parent(a))) != 1 && error("Only implemented for localizations at one prime ideal")
   return valuation(data(a), prime(parent(a)))
end

@doc Markdown.doc"""
    valuation(a::LocElem{T}, p::T) where {T <: RingElement}
Returns the valuation `n` of $a$ at $p$, i.e. the integer `n` s.th $a$ = $p$^`n` * x, where x has valuation 0 at $p$
"""
valuation(a::LocElem{T}, p::T) where {T <: RingElement} = valuation(data(a), p)

###############################################################################
#
#   Canonicalisation
#
###############################################################################

@doc Markdown.doc"""
    canonical_unit(a::LocElem{T}) where {T <: RingElem}
Returns unit `b`::LocElem{T} s.th. $a$ * `b` only consists of powers of primes localized at.
"""
function canonical_unit(a::LocElem{T}) where {T <: RingElem}
   temp = data(a)
   for pi in primes(parent(a))
      _,temp = remove(temp, pi)
   end
   return parent(a)(inv(temp))
end

###############################################################################
#
#   Constructors
#
###############################################################################

@doc Markdown.doc"""
    Localization(R::AbstractAlgebra.Ring, prime::T; cached=true) where {T <: RingElement}
Returns the localization of the ring $R$ at the ideal generated by the ring element $prime$. Requires $R$ to
be an euclidean domain and $prime$ to be a prime element, both not checked.
If `cached == true` (the default) then the resulting
localization parent object is cached and returned for any subsequent calls
to the constructor with the same base ring $R$ and element $prime$.
"""
function Localization(R::AbstractAlgebra.Ring, prime::T; cached=true) where {T <: RingElement}
   primes = [R(prime)]
   return Loc{elem_type(R)}(R(prime), primes, cached)
end

@doc Markdown.doc"""
     Localization(R::AbstractAlgebra.Ring, primes::Array{T,1}; cached=true) where {T <: RingElement}
Returns the localization of the ring $R$ at the union of principal ideals that are generated by the ring elements in $primes$.
Requires $R$ to be an euclidean domain and $primes$ to be prime elements, both not checked.
If `cached == true` (the default) then the resulting
localization parent object is cached and returned for any subsequent calls
to the constructor with the same base ring $R$ and elements $primes$.
"""
function Localization(R::AbstractAlgebra.Ring, primes::Array{T,1}; cached=true) where {T <: RingElement}
   prime = R(prod(primes))
   return Loc{elem_type(R)}(prime, Array{elem_type(R),1}(primes), cached)
end


text:
using Mandrill
using Compat
using Base.Test

# write your own tests here
p = Mandrill.users_ping2()
@test p["PING"] == "PONG!"

p = Mandrill.users_info()
p["public_id"] == "WBdNPiLFZOoKiEZ-stgZOQ"

s = @compat Dict(
  "message" => Dict(
    "html" => "
 Example HTML Content
",
    "text" => "Example Test Content",
    "subject" => "Example Subject",
    "from_email" => "message.from@example.com",
    "from_name" => "Example Name",
    "to" => [
      Dict(
        "email" => "recipient.email@example.com",
        "name" => "Recipient Name",
        "type" => "to"
        )
    ],
    "headers" => Dict(
      "Reply-To" => "message.reply@example.com"
    ),
    "important" => false,
    "tags" => [
      "password-resets",
      "user-initiated"
    ],
    "metadata" => Dict(
      "website" => "www.example.com"
    )
  ),
  "async" => false,
  "ip_pool" => "Main Pool"
)

p=Mandrill.messages_send(s)
@test length(p) == 1
@test p[1]["status"] == "sent"

s = Dict{ASCIIString, Any}()
message = Dict{ASCIIString, Any}()
message["html"] = "
Example HTML Content
"
message["text"] = "Example text content"
message["subject"] = "message.from@example.com"
rcpt1 = Dict{ASCIIString, Any}()
rcpt1["email"] = "rcpt1@example.com"
rcpt1["name"] = "Recipient One"
rcpt1["type"] = "to"
rcpt2 = Dict{ASCIIString, Any}()
rcpt2["email"] = "rcpt2@example.com"
rcpt2["name"] = "Recipient Two"
rcpt2["type"] = "cc"
message["to"] = [rcpt1, rcpt2]
s["message"]=message
s["async"] = true

p=Mandrill.messages_send(s)
@test length(p) == 2
@test p[1]["status"] == "queued"
@test p[2]["status"] == "queued"


text:
#Function for estimating a problem's time step
function update_a2a3adot(surf::TwoDSurf,dt)
    for ia = 2:3
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:3
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

function update_atermdot(surf::TwoDSurf,dt)
    for ia = 2:naterm
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:naterm
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

function update_adot(surf::TwoDSurf,dt)
    surf.a0dot[1] = (surf.a0[1] - surf.a0prev[1])/dt
    for ia = 1:3
        surf.adot[ia] = (surf.aterm[ia]-surf.aprev[ia])/dt
    end
    return surf
end

# Function for updating the induced velocities
function update_indbound(surf::TwoDSurf, curfield::TwoDFlowField)
    surf.uind[1:surf.ndiv], surf.wind[1:surf.ndiv] = ind_vel([curfield.tev; curfield.lev; curfield.extv], surf.bnd_x, surf.bnd_z)
    return surf
end

function add_indbound_b(surf::TwoDSurf, surfj::TwoDSurf)
    uind, wind = ind_vel(surfj.bv, surf.bnd_x, surf.bnd_z)
    surf.uind[:] += uind
    surf.wind[:] += wind
    return surf

end

# Function for updating the downwash
function update_downwash(surf::TwoDSurf, vels::Vector{Float64})
    for ib = 1:surf.ndiv
        surf.downwash[ib] = -(surf.kinem.u + vels[1])*sin(surf.kinem.alpha) - surf.uind[ib]*sin(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*cos(surf.kinem.alpha) - surf.wind[ib]*cos(surf.kinem.alpha) - surf.kinem.alphadot*(surf.x[ib] - surf.pvt*surf.c) + surf.cam_slope[ib]*(surf.uind[ib]*cos(surf.kinem.alpha) + (surf.kinem.u + vels[1])*cos(surf.kinem.alpha) + (surf.kinem.hdot - vels[2])*sin(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))
    end
    return surf
end

# Function for a_0 and a_1 fourier coefficients
function update_a0anda1(surf::TwoDSurf)
    surf.a0[1] = simpleTrapz(surf.downwash,surf.theta)
    surf.aterm[1] = simpleTrapz(surf.downwash.*cos.(surf.theta),surf.theta)
    surf.a0[1] = -surf.a0[1]/(surf.uref*pi)
    surf.aterm[1] = 2. *surf.aterm[1]/(surf.uref*pi)
    return surf
end

# Function for calculating the fourier coefficients a_2 upwards to a_n
function update_a2toan(surf::TwoDSurf)
    for ia = 2:surf.naterm
        surf.aterm[ia] = simpleTrapz(surf.downwash.*cos.(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2. *surf.aterm[ia]/(surf.uref*pi)
    end
    return surf
end

#Function to update the external flowfield
function update_externalvel(curfield::TwoDFlowField, t)
    if (typeof(curfield.velX) == CosDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    elseif (typeof(curfield.velX) == SinDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    elseif (typeof(curfield.velX) == ConstDef)
        curfield.u[1] = curfield.velX(t)
        curfield.w[1] = curfield.velZ(t)
    end
    if (typeof(curfield.velX) == StepGustDef)
        curfield.u[1] = curfield.velX(t)
    end
    if (typeof(curfield.velZ) == StepGustDef)
        curfield.w[1] = curfield.velZ(t)
    end
end

# Function updating the dimensional kinematic parameters
function update_kinem(surf::TwoDSurf, t)

    # Pitch kinematics
    if (typeof(surf.kindef.alpha) == EldUpDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldUptstartDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldRampReturnDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == ConstDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = 0.
    elseif (typeof(surf.kindef.alpha) == SinDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == CosDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == VAWTalphaDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------

    # Plunge kinematics
    if (typeof(surf.kindef.h) == EldUpDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUptstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpIntDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpInttstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldRampReturnDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == ConstDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = 0.
    elseif (typeof(surf.kindef.h) == SinDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == CosDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot =  ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == VAWThDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    end
    # ---------------------------------------------------------------------------------------------

    # Forward velocity
    if (typeof(surf.kindef.u) == EldUpDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == EldRampReturnDef)
        surf.kinem.u, surf.kinem.udot = surf.kindef.u(t)
        surf.kinem.u = surf.kinem.u*surf.uref
        surf.kinem.udot = surf.kinem.udot*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == ConstDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = 0.
    elseif (typeof(surf.kindef.u) == SinDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == CosDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == LinearDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == VAWTuDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------
    return surf
end

# ---------------------------------------------------------------------------------------------
# Updates the bound vorticity distribution: eqn (2.1) in Ramesh et al. (2013)
# determines the strength of bound vortices
# determines the x, z components of the bound vortices
function update_bv(surf::TwoDSurf)
    gamma = zeros(surf.ndiv)
    for ib = 1:surf.ndiv
        gamma[ib] = (surf.a0[1]*(1 + cos(surf.theta[ib])))
        for ia = 1:surf.naterm
            gamma[ib] = gamma[ib] + surf.aterm[ia]*sin(ia*surf.theta[ib])*sin(surf.theta[ib])
        end
        gamma[ib] = gamma[ib]*surf.uref*surf.c
    end

    for ib = 2:surf.ndiv
        surf.bv[ib-1].s = (gamma[ib]+gamma[ib-1])*(surf.theta[2]-surf.theta[1])/2.
        surf.bv[ib-1].x = (surf.bnd_x[ib] + surf.bnd_x[ib-1])/2.
        surf.bv[ib-1].z = (surf.bnd_z[ib] + surf.bnd_z[ib-1])/2.
    end
end

# Function for calculating the wake rollup
function wakeroll(surf::TwoDSurf, curfield::TwoDFlowField, dt)

    nlev = length(curfield.lev)
    ntev = length(curfield.tev)
    nextv = length(curfield.extv)

    #Clean induced velocities
    for i = 1:ntev
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:nlev
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    for i = 1:nextv
        curfield.extv[i].vx = 0
        curfield.extv[i].vz = 0
    end

    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev; curfield.extv])

    #Add the influence of velocities induced by bound vortices
    utemp = zeros(ntev + nlev + nextv)
    wtemp = zeros(ntev + nlev + nextv)
    utemp, wtemp = ind_vel(surf.bv, [map(q -> q.x, curfield.tev); map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)], [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev); map(q -> q.z, curfield.extv) ])

    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Add the influence of freestream velocities
    for i = 1:ntev
        curfield.tev[i].vx += curfield.u[1]
        curfield.tev[i].vz += curfield.w[1]
    end
    for i = 1:nlev
        curfield.lev[i].vx += curfield.u[1]
        curfield.lev[i].vz += curfield.w[1]
    end
    for i = 1:nextv
        curfield.extv[i].vx += curfield.u[1]
        curfield.extv[i].vz += curfield.w[1]
    end

    #Convect free vortices with their induced velocities
    for i = 1:ntev
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:nlev
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:nextv
        curfield.extv[i].x += dt*curfield.extv[i].vx
        curfield.extv[i].z += dt*curfield.extv[i].vz
    end

    return curfield
end

function wakeroll(surf::Vector{TwoDSurf}, curfield::TwoDFlowField, dt)

    nlev = length(curfield.lev)
    ntev = length(curfield.tev)
    nextv = length(curfield.extv)
    nsurf = length(surf)

    #Clean induced velocities
    for i = 1:ntev
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:nlev
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    for i = 1:nextv
        curfield.extv[i].vx = 0
        curfield.extv[i].vz = 0
    end

    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev; curfield.extv])

    #Add the influence of velocities induced by bound vortices
    utemp = zeros(ntev + nlev + nextv)
    wtemp = zeros(ntev + nlev + nextv)
    for is = 1:nsurf
        ut, wt = ind_vel(surf[is].bv, [map(q -> q.x, curfield.tev); map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)], [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev); map(q -> q.z, curfield.extv) ])
        utemp += ut
        wtemp += wt
    end
    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Add the influence of freestream velocities
    for i = 1:ntev
        curfield.tev[i].vx += curfield.u[1]
        curfield.tev[i].vz += curfield.w[1]
    end
    for i = 1:nlev
        curfield.lev[i].vx += curfield.u[1]
        curfield.lev[i].vz += curfield.w[1]
    end
    for i = 1:nextv
        curfield.extv[i].vx += curfield.u[1]
        curfield.extv[i].vz += curfield.w[1]
    end

    #Convect free vortices with their induced velocities
    for i = 1:ntev
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:nlev
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:nextv
        curfield.extv[i].x += dt*curfield.extv[i].vx
        curfield.extv[i].z += dt*curfield.extv[i].vz
    end

    return curfield
end

# Places a trailing edge vortex
function place_tev(surf::TwoDSurf,field::TwoDFlowField,dt)
    ntev = length(field.tev)
    if ntev == 0
        xloc = surf.bnd_x[surf.ndiv] + 0.5*surf.kinem.u*dt
        zloc = surf.bnd_z[surf.ndiv]
    else
        xloc = surf.bnd_x[surf.ndiv]+(1. /3.)*(field.tev[ntev].x - surf.bnd_x[surf.ndiv])

        zloc = surf.bnd_z[surf.ndiv]+(1. /3.)*(field.tev[ntev].z - surf.bnd_z[surf.ndiv])
    end
    push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))
    return field
end

function place_tev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt)
    nsurf = length(surf)
    ntev = length(field.tev)
    if ntev < nsurf
        for i = 1:nsurf
            xloc = surf[i].bnd_x[end] + 0.5*surf[i].kinem.u*dt
            zloc = surf[i].bnd_z[end]
            push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
        end
    else
        for i = 1:nsurf
            xloc = surf[i].bnd_x[end]+(1. /3.)*(field.tev[ntev-nsurf+i].x - surf[i].bnd_x[end])
            zloc = surf[i].bnd_z[end]+(1. /3.)*(field.tev[ntev-nsurf+i].z - surf[i].bnd_z[end])
            push!(field.tev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
        end
    end
    return field
end

# Places a leading edge vortex
function place_lev(surf::TwoDSurf,field::TwoDFlowField,dt)
    nlev = length(field.lev)

    le_vel_x = surf.kinem.u - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1]
    le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c- surf.kinem.hdot + surf.wind[1]

    if (surf.levflag[1] == 0)
        xloc = surf.bnd_x[1] + 0.5*le_vel_x*dt
        zloc = surf.bnd_z[1] + 0.5*le_vel_z*dt
    else
        xloc = surf.bnd_x[1]+(1. /3.)*(field.lev[nlev].x - surf.bnd_x[1])
        zloc = surf.bnd_z[1]+(1. /3.)*(field.lev[nlev].z - surf.bnd_z[1])
    end

    push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf.c,0.,0.))

    return field
end

function place_lev(surf::Vector{TwoDSurf},field::TwoDFlowField,dt,shed_ind::Vector{Int})
    nsurf = length(surf)
    nlev = length(field.lev)

    for i = 1:nsurf
        if i in shed_ind
            if (surf[i].levflag[1] == 0)
                le_vel_x = surf[i].kinem.u - surf[i].kinem.alphadot*sin(surf[i].kinem.alpha)*surf[i].pvt*surf[i].c + surf[i].uind[1]
                le_vel_z = -surf[i].kinem.alphadot*cos(surf[i].kinem.alpha)*surf[i].pvt*surf[i].c- surf[i].kinem.hdot + surf[i].wind[1]
                xloc = surf[i].bnd_x[1] + 0.5*le_vel_x*dt
                zloc = surf[i].bnd_z[1] + 0.5*le_vel_z*dt
                push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
            else
                xloc = surf[i].bnd_x[1]+(1. /3.)*(field.lev[nlev-nsurf+i].x - surf[i].bnd_x[1])
                zloc = surf[i].bnd_z[1]+(1. /3.)*(field.lev[nlev-nsurf+i].z - surf[i].bnd_z[1])
                push!(field.lev,TwoDVort(xloc,zloc,0.,0.02*surf[i].c,0.,0.))
            end
        else
            push!(field.lev, TwoDVort(0., 0., 0., 0., 0., 0.))
        end
    end
    return field
end

# Function for updating the positions of the bound vortices
function update_boundpos(surf::TwoDSurf, dt::Float64)
    for i = 1:surf.ndiv
        surf.bnd_x[i] = surf.bnd_x[i] + dt*((surf.pvt*surf.c - surf.x[i])*sin(surf.kinem.alpha)*surf.kinem.alphadot - surf.kinem.u + surf.cam[i]*cos(surf.kinem.alpha)*surf.kinem.alphadot)
        surf.bnd_z[i] = surf.bnd_z[i] + dt*(surf.kinem.hdot + (surf.pvt*surf.c - surf.x[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot - surf.cam[i]*sin(surf.kinem.alpha)*surf.kinem.alphadot)
    end
    return surf
end

# Function to calculate induced velocities by a set of votices at a target location
function ind_vel(src::Vector{TwoDVort},t_x,t_z)
    #'s' stands for source and 't' for target
    uind = zeros(length(t_x))
    wind = zeros(length(t_x))

    for itr = 1:length(t_x)
        for isr = 1:length(src)
            xdist = src[isr].x - t_x[itr]
            zdist = src[isr].z - t_z[itr]
            distsq = xdist*xdist + zdist*zdist
            uind[itr] = uind[itr] - src[isr].s*zdist/(2*pi*sqrt(src[isr].vc*src[isr].vc*src[isr].vc*src[isr].vc + distsq*distsq))
            wind[itr] = wind[itr] + src[isr].s*xdist/(2*pi*sqrt(src[isr].vc*src[isr].vc*src[isr].vc*src[isr].vc + distsq*distsq))
        end
    end
    return uind, wind
end

# Function determining the effects of interacting vorticies - velocities induced on each other - classical n-body problem
function mutual_ind(vorts::Vector{TwoDVort})
    for i = 1:length(vorts)
        for j = i+1:length(vorts)
            dx = vorts[i].x - vorts[j].x
            dz = vorts[i].z - vorts[j].z
            #source- tar
            dsq = dx*dx + dz*dz

            magitr = 1. /(2*pi*sqrt(vorts[j].vc*vorts[j].vc*vorts[j].vc*vorts[j].vc + dsq*dsq))
            magjtr = 1. /(2*pi*sqrt(vorts[i].vc*vorts[i].vc*vorts[i].vc*vorts[i].vc + dsq*dsq))

            vorts[j].vx -= dz * vorts[i].s * magjtr
            vorts[j].vz += dx * vorts[i].s * magjtr

            vorts[i].vx += dz * vorts[j].s * magitr
            vorts[i].vz -= dx * vorts[j].s * magitr
        end
    end
    return vorts
end

"""
    controlVortCount(delvort, surf, curfield)

Performs merging or deletion operation on free vortices in order to
control computational cost according to algorithm specified.

Algorithms for parameter delvort

 - delNone

    Does nothing, no vortex count control.

 - delSpalart(limit=500, dist=10, tol=1e-6)

    Merges vortices according to algorithm given in Spalart,
    P. R. (1988). Vortex methods for separated flows.

     - limit: min number of vortices present for merging to occur

     - dist: small values encourage mergin near airfoil, large values in
         wake (see paper)

     - tol: tolerance for merging, merging is less likely to occur for low
        values (see paper)

There is no universal set of parameters that work for all problem. If
using vortex control, test and calibrate choice of parameters

"""
function controlVortCount(delvort :: delVortDef, surf_locx :: Float64, surf_locz :: Float64, curfield :: TwoDFlowField)

    if typeof(delvort) == delNone

    elseif typeof(delvort) == delSpalart
        D0 = delvort.dist
        V0 = delvort.tol
        if length(curfield.tev) > delvort.limit
            #Check possibility of merging the last vortex with the closest 20 vortices
            gamma_j = curfield.tev[1].s
            d_j = sqrt((curfield.tev[1].x- surf_locx)^2 + (curfield.tev[1].z - surf_locz)^2)
            z_j = sqrt(curfield.tev[1].x^2 + curfield.tev[1].z^2)

            for i = 2:20
                gamma_k = curfield.tev[i].s
                d_k = sqrt((curfield.tev[i].x - surf_locx)^2 + (curfield.tev[i].z - surf_locz)^2)
                z_k = sqrt(curfield.tev[i].x^2 + curfield.tev[i].z^2)

                fact = abs(gamma_j*gamma_k)*abs(z_j - z_k)/(abs(gamma_j + gamma_k)*(D0 + d_j)^1.5*(D0 + d_k)^1.5)

                if fact < V0
                    #Merge the 2 vortices into the one at k
                    curfield.tev[i].x = (abs(gamma_j)*curfield.tev[1].x + abs(gamma_k)*curfield.tev[i].x)/(abs(gamma_j + gamma_k))
                    curfield.tev[i].z = (abs(gamma_j)*curfield.tev[1].z + abs(gamma_k)*curfield.tev[i].z)/(abs(gamma_j + gamma_k))
                    curfield.tev[i].s += curfield.tev[1].s

                    popfirst!(curfield.tev)

                    break
                end
            end

            if length(curfield.lev) > delvort.limit
                #Check possibility of merging the last vortex with the closest 20 vortices

                gamma_j = curfield.lev[1].s
                d_j = sqrt((curfield.lev[1].x - surf_locx)^2 + (curfield.lev[1].z - surf_locz)^2)
                z_j = sqrt(curfield.lev[1].x^2 + curfield.lev[1].z^2)

                for i = 2:20
                    gamma_k = curfield.lev[i].s
                    d_k = sqrt((curfield.lev[i].x - surf_locx)^2 + (curfield.lev[i].z - surf_locz)^2)
                    z_k = sqrt(curfield.lev[i].x^2 + curfield.lev[i].z^2)

                    fact = abs(gamma_j*gamma_k)*abs(z_j - z_k)/(abs(gamma_j + gamma_k)*(D0 + d_j)^1.5*(D0 + d_k)^1.5)

                    if fact < V0
                        #Merge the 2 vortices into the one at k
                        curfield.lev[i].x = (abs(gamma_j)*curfield.lev[1].x + abs(gamma_k)*curfield.lev[i].x)/(abs(gamma_j + gamma_k))
                        curfield.lev[i].z = (abs(gamma_j)*curfield.lev[1].z + abs(gamma_k)*curfield.lev[i].z)/(abs(gamma_j + gamma_k))
                        curfield.lev[i].s += curfield.lev[1].s

                        popfirst!(curfield.lev)

                        break
                    end
                end
            end
        end
    end
end

function update_kinem2DOF(surf::TwoDSurf, strpar :: TwoDOFPar, kinem :: KinemPar2DOF, dt, cl, cm)
    #Update previous terms
    kinem.alpha_pr = kinem.alpha
    kinem.h_pr = kinem.h

    kinem.alphadot_pr3 = kinem.alphadot_pr2
    kinem.alphadot_pr2 = kinem.alphadot_pr
    kinem.alphadot_pr = kinem.alphadot

    kinem.hdot_pr3 = kinem.hdot_pr2
    kinem.hdot_pr2 = kinem.hdot_pr
    kinem.hdot_pr = kinem.hdot

    kinem.alphaddot_pr3 = kinem.alphaddot_pr2
    kinem.alphaddot_pr2 = kinem.alphaddot_pr

    kinem.hddot_pr3 = kinem.hddot_pr2
    kinem.hddot_pr2 = kinem.hddot_pr

    # Calculate hddot and alphaddot from forces based on 2DOF structural system
    m11 = 2. /surf.c
    m12 = -strpar.x_alpha*cos(kinem.alpha)
    m21 = -2. *strpar.x_alpha*cos(kinem.alpha)/surf.c
    m22 = strpar.r_alpha*strpar.r_alpha

    R1 = 4*strpar.kappa*surf.uref*surf.uref*cl/(pi*surf.c*surf.c) - 2*strpar.w_h*strpar.w_h*(strpar.cubic_h_1*kinem.h + strpar.cubic_h_3*kinem.h^3)/surf.c - strpar.x_alpha*sin(kinem.alpha)*kinem.alphadot*kinem.alphadot

    R2 = 8*strpar.kappa*surf.uref*surf.uref*cm/(pi*surf.c*surf.c) - strpar.w_alpha*strpar.w_alpha*strpar.r_alpha*strpar.r_alpha*(strpar.cubic_alpha_1*kinem.alpha + strpar.cubic_alpha_3*kinem.alpha^3)

    kinem.hddot_pr = (1/(m11*m22-m21*m12))*(m22*R1-m12*R2)
    kinem.alphaddot_pr = (1/(m11*m22-m21*m12))*(-m21*R1+m11*R2)

    #Kinematics are updated according to the 2DOF response
    kinem.alphadot = kinem.alphadot_pr + (dt/12.)*(23*kinem.alphaddot_pr - 16*kinem.alphaddot_pr2 + 5*kinem.alphaddot_pr3)
    kinem.hdot = kinem.hdot_pr + (dt/12)*(23*kinem.hddot_pr-16*kinem.hddot_pr2 + 5*kinem.hddot_pr3)

    kinem.alpha = kinem.alpha_pr + (dt/12)*(23*kinem.alphadot_pr-16*kinem.alphadot_pr2 + 5*kinem.alphadot_pr3)
    kinem.h = kinem.h_pr + (dt/12)*(23*kinem.hdot_pr - 16*kinem.hdot_pr2 + 5*kinem.hdot_pr3)

    #Copy over these values to the TwoDSurf
    surf.kinem.alphadot = kinem.alphadot
    surf.kinem.hdot = kinem.hdot
    surf.kinem.alpha = kinem.alpha
    surf.kinem.h = kinem.h

    return surf, kinem
end

function update_indbound(surf::TwoDSurfThick, curfield::TwoDFlowField)
    surf.uind_u[1:surf.ndiv], surf.wind_u[1:surf.ndiv] = ind_vel([curfield.tev;
                                                                  curfield.lev; curfield.extv], surf.bnd_x_u, surf.bnd_z_u)

    surf.uind_l[1:surf.ndiv], surf.wind_l[1:surf.ndiv] = ind_vel([curfield.tev;
                                                                  curfield.lev; curfield.extv], surf.bnd_x_l, surf.bnd_z_l)

    return surf
end

function update_kinem(surf::TwoDSurfThick, t)

    # Pitch kinematics
    if (typeof(surf.kindef.alpha) == EldUpDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldUptstartDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == EldRampReturnDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == ConstDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = 0.
    elseif (typeof(surf.kindef.alpha) == SinDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == CosDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.alpha) == VAWTalphaDef)
        surf.kinem.alpha = surf.kindef.alpha(t)
        surf.kinem.alphadot = ForwardDiff.derivative(surf.kindef.alpha,t)*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------

    # Plunge kinematics
    if (typeof(surf.kindef.h) == EldUpDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUptstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpIntDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldUpInttstartDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == EldRampReturnDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == ConstDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = 0.
    elseif (typeof(surf.kindef.h) == SinDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == CosDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot =  ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    elseif (typeof(surf.kindef.h) == VAWThDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = ForwardDiff.derivative(surf.kindef.h,t)*surf.uref
    end
    # ---------------------------------------------------------------------------------------------

    # Forward velocity
    if (typeof(surf.kindef.u) == EldUpDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == EldRampReturnDef)
        surf.kinem.u, surf.kinem.udot = surf.kindef.u(t)
        surf.kinem.u = surf.kinem.u*surf.uref
        surf.kinem.udot = surf.kinem.udot*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == ConstDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = 0.
    elseif (typeof(surf.kindef.u) == SinDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == CosDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == LinearDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    elseif (typeof(surf.kindef.u) == VAWTuDef)
        surf.kinem.u = surf.kindef.u(t)*surf.uref
        surf.kinem.udot = ForwardDiff.derivative(surf.kindef.u,t)*surf.uref*surf.uref/surf.c
    end
    # ---------------------------------------------------------------------------------------------
    return surf
end

function wakeroll(surf::TwoDSurfThick, curfield::TwoDFlowField, dt)

    nlev = length(curfield.lev)
    ntev = length(curfield.tev)
    nextv = length(curfield.extv)

    #Clean induced velocities
    for i = 1:ntev
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:nlev
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    for i = 1:nextv
        curfield.extv[i].vx = 0
        curfield.extv[i].vz = 0
    end

    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev; curfield.extv])

    #Add the influence of velocities induced by bound vortices
    utemp = zeros(ntev + nlev + nextv)
    wtemp = zeros(ntev + nlev + nextv)
    utemp, wtemp = ind_vel(surf.bv, [map(q -> q.x, curfield.tev); map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)], [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev); map(q -> q.z, curfield.extv) ])

    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Add effect of sources
    utemp, wtemp = ind_vel_src(surf.src, [map(q -> q.x, curfield.tev);
                                          map(q -> q.x, curfield.lev); map(q -> q.x, curfield.extv)],
                               [map(q -> q.z, curfield.tev); map(q -> q.z, curfield.lev);
                                map(q -> q.z, curfield.extv)])

    for i = 1:ntev
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = ntev+1:ntev+nlev
        curfield.lev[i-ntev].vx += utemp[i]
        curfield.lev[i-ntev].vz += wtemp[i]
    end
    for i = ntev+nlev+1:ntev+nlev+nextv
        curfield.extv[i-ntev-nlev].vx += utemp[i]
        curfield.extv[i-ntev-nlev].vz += wtemp[i]
    end

    #Convect free vortices with their induced velocities
    for i = 1:ntev
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:nlev
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:nextv
        curfield.extv[i].x += dt*curfield.extv[i].vx
        curfield.extv[i].z += dt*curfield.extv[i].vz
    end

    return curfield
end

function update_boundpos(surf::TwoDSurfThick, dt::Float64)
    for i = 1:surf.ndiv
        surf.bnd_x_u[i] = surf.bnd_x_u[i] + dt*((surf.pvt*surf.c - surf.x[i])*
                                                sin(surf.kinem.alpha)*surf.kinem.alphadot - surf.kinem.u + (surf.cam[i] +
                                                                                surf.thick[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot)

        surf.bnd_z_u[i] = surf.bnd_z_u[i] + dt*(surf.kinem.hdot + (surf.pvt*surf.c -
                                                                   surf.x[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot - (surf.cam[i] +
                                                                                surf.thick[i])*sin(surf.kinem.alpha)*surf.kinem.alphadot)

        surf.bnd_x_l[i] = surf.bnd_x_l[i] + dt*((surf.pvt*surf.c - surf.x[i])*
                                                sin(surf.kinem.alpha)*surf.kinem.alphadot - surf.kinem.u + (surf.cam[i] -
                                                                                surf.thick[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot)

        surf.bnd_z_l[i] = surf.bnd_z_l[i] + dt*(surf.kinem.hdot + (surf.pvt*
                                                                   surf.c - surf.x[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot - (surf.cam[i] -
                                                                                surf.thick[i])*sin(surf.kinem.alpha)*surf.kinem.alphadot)

        surf.bnd_x_chord[i] = surf.bnd_x_chord[i] + dt*((surf.pvt*surf.c - surf.x[i])*
                                                        sin(surf.kinem.alpha)*surf.kinem.alphadot - surf.kinem.u)
        surf.bnd_z_chord[i] = surf.bnd_z_chord[i] + dt*(surf.kinem.hdot + (surf.pvt*surf.c -
                                                                           surf.x[i])*cos(surf.kinem.alpha)*surf.kinem.alphadot)
    end

    return surf
end

function update_thickLHS(surf::TwoDSurfThick, curfield::TwoDFlowField, dt::Float64, vcore::Float64)
    #Calculate the missing column in LHS that depends on last shed vortex location

    ntev = length(curfield.tev)

    if ntev == 0
        xloc_tev = surf.bnd_x_chord[surf.ndiv] + 0.5*surf.kinem.u*dt
        zloc_tev = surf.bnd_z_chord[surf.ndiv]
    else
        xloc_tev = surf.bnd_x_chord[surf.ndiv] + (1. /3.)*(curfield.tev[ntev].x - surf.bnd_x_chord[surf.ndiv])
        zloc_tev = surf.bnd_z_chord[surf.ndiv] + (1. /3.)*(curfield.tev[ntev].z - surf.bnd_z_chord[surf.ndiv])
    end

    dummyvort = TwoDVort(xloc_tev, zloc_tev, 1., vcore, 0., 0.)

    uu, wu = ind_vel([dummyvort], surf.bnd_x_u, surf.bnd_z_u)
    ul, wl = ind_vel([dummyvort], surf.bnd_x_l, surf.bnd_z_l)

    for i = 2:surf.ndiv-1
        wlz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) +
                   wl[i]*cos(surf.kinem.alpha) + ul[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) -
                   wl[i]*cos(surf.kinem.alpha) - ul[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) +
                   ul[i]*cos(surf.kinem.alpha) - wl[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) -
                   ul[i]*cos(surf.kinem.alpha) + wl[i]*sin(surf.kinem.alpha))

        surf.LHS[i-1,2+2*surf.naterm] = -surf.cam_slope[i]*wtx -
            surf.thick_slope[i]*wlx + wlz

        surf.LHS[surf.ndiv+i-3,2+2*surf.naterm] = -surf.cam_slope[i]*wlx -
            surf.thick_slope[i]*wtx + wtz

    end

    return surf, xloc_tev, zloc_tev
end

function update_thickLHSLEV(surf::TwoDSurfThick, curfield::TwoDFlowField, dt::Float64, vcore::Float64)
    #Calculate the missing column in LHS that depends on last shed vortex location

    ntev = length(curfield.tev)
    nlev = length(curfield.lev)

    if ntev == 0
        xloc_tev = surf.bnd_x_chord[surf.ndiv] + 0.5*surf.kinem.u*dt
        zloc_tev = surf.bnd_z_chord[surf.ndiv]
    else
        xloc_tev = surf.bnd_x_chord[surf.ndiv] + (1. /3.)*(curfield.tev[ntev].x - surf.bnd_x_chord[surf.ndiv])
        zloc_tev = surf.bnd_z_chord[surf.ndiv] + (1. /3.)*(curfield.tev[ntev].z - surf.bnd_z_chord[surf.ndiv])
    end

    if surf.levflag[1] == 0

        le_vel_x = (surf.kinem.u + curfield.u[1]) - surf.kinem.alphadot*sin(surf.kinem.alpha)*surf.pvt*surf.c + surf.uind[1]
        le_vel_z = -surf.kinem.alphadot*cos(surf.kinem.alpha)*surf.pvt*surf.c - (surf.kinem.hdot - curfield.w[1]) + surf.wind[1]

        xloc = surf.bnd_x[1] + 0.5*le_vel_x*dt
        zloc = surf.bnd_z[1] + 0.5*le_vel_z*dt
    else
        xloc = surf.bnd_x[1]+(1. /3.)*(curfield.lev[nlev].x - surf.bnd_x[1])
        zloc = surf.bnd_z[1]+(1. /3.)*(curfield.lev[nlev].z - surf.bnd_z[1])
    end

    dummyvort = TwoDVort(xloc_tev, zloc_tev, 1., vcore, 0., 0.)

    uu, wu = ind_vel([dummyvort], surf.bnd_x_u, surf.bnd_z_u)
    ul, wl = ind_vel([dummyvort], surf.bnd_x_l, surf.bnd_z_l)

    for i = 2:surf.ndiv-1
        wlz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) +
                   wl[i]*cos(surf.kinem.alpha) + ul[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) -
                   wl[i]*cos(surf.kinem.alpha) - ul[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) +
                   ul[i]*cos(surf.kinem.alpha) - wl[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) -
                   ul[i]*cos(surf.kinem.alpha) + wl[i]*sin(surf.kinem.alpha))

        surf.LHS[i-1,2+2*surf.naterm] = -surf.cam_slope[i]*wtx -
            surf.thick_slope[i]*wlx + wlz

        surf.LHS[surf.ndiv+i-3,2+2*surf.naterm] = -surf.cam_slope[i]*wlx -
            surf.thick_slope[i]*wtx + wtz

    end

    dummyvort = TwoDVort(xloc_lev, zloc_lev, 1., vcore, 0., 0.)

    for i = 2:surf.ndiv-1
        wlz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) +
                   wl[i]*cos(surf.kinem.alpha) + ul[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(wu[i]*cos(surf.kinem.alpha) + uu[i]*sin(surf.kinem.alpha) -
                   wl[i]*cos(surf.kinem.alpha) - ul[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) +
                   ul[i]*cos(surf.kinem.alpha) - wl[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(uu[i]*cos(surf.kinem.alpha) - wu[i]*sin(surf.kinem.alpha) -
                   ul[i]*cos(surf.kinem.alpha) + wl[i]*sin(surf.kinem.alpha))

        surf.LHS[i-1,3+2*surf.naterm] = -surf.cam_slope[i]*wtx -
            surf.thick_slope[i]*wlx + wlz

        surf.LHS[surf.ndiv+i-3,3+2*surf.naterm] = -surf.cam_slope[i]*wlx -
            surf.thick_slope[i]*wtx + wtz
    end

    LHS[2*ndiv-3,2*naterm+3] = 1.

    LHS[2*ndiv-2,1] = 1.

    return surf, xloc_tev, zloc_tev, xloc_lev, zloc_lev
end

function update_thickRHS(surf::TwoDSurfThick, curfield::TwoDFlowField)
    for i = 2:surf.ndiv-1
        #RHS for lifting equation

        wlz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) +
                   surf.wind_l[i]*cos(surf.kinem.alpha) + surf.uind_l[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) -
                   surf.wind_l[i]*cos(surf.kinem.alpha) - surf.uind_l[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) +
                   surf.uind_l[i]*cos(surf.kinem.alpha) - surf.wind_l[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) -
                   surf.uind_l[i]*cos(surf.kinem.alpha) + surf.wind_l[i]*sin(surf.kinem.alpha))

        surf.RHS[i-1] = -(surf.kinem.u + curfield.u[1])*sin(surf.kinem.alpha)/surf.uref -
            surf.kinem.alphadot*(surf.x[i] - surf.pvt*surf.c)/surf.uref +
            (surf.kinem.hdot - curfield.w[1])*cos(surf.kinem.alpha)/surf.uref - wlz/surf.uref +
            surf.cam_slope[i]*((surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) +
                               (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha) + wtx -
                               surf.kinem.alphadot*surf.cam[i]) +
                               surf.thick_slope[i]*(wlx - surf.kinem.alphadot*surf.thick[i])

        surf.RHS[surf.ndiv+i-3] = surf.cam_slope[i]*(wlx - surf.kinem.alphadot*surf.thick[i]) +
            surf.thick_slope[i]*((surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) +
                                 (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha) + wtx -
                                 surf.kinem.alphadot*surf.cam[i]) - wtz
    end

    #RHS for Kelvin condition (negative strength of all previously shed vortices)
    #surf.RHS[2*surf.ndiv-3] = -(sum(map(q->q.s, curfield.tev)) + sum(map(q->q.s, curfield.lev)))/(surf.uref*surf.c)
    surf.RHS[2*surf.ndiv-3] = pi*(surf.a0prev[1] + 0.5*surf.aprev[1])

    surf.RHS[2*surf.ndiv-2] = -(surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) - (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha)

    surf.RHS[2*surf.ndiv-1] = -(surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) - (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha)

end

function update_thickRHSLEV(surf::TwoDSurfThick, curfield::TwoDFlowField)
    for i = 2:surf.ndiv-1
        #RHS for lifting equation

        wlz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) +
                   surf.wind_l[i]*cos(surf.kinem.alpha) + surf.uind_l[i]*sin(surf.kinem.alpha))

        wtz = 0.5*(surf.wind_u[i]*cos(surf.kinem.alpha) + surf.uind_u[i]*sin(surf.kinem.alpha) -
                   surf.wind_l[i]*cos(surf.kinem.alpha) - surf.uind_l[i]*sin(surf.kinem.alpha))

        wtx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) +
                   surf.uind_l[i]*cos(surf.kinem.alpha) - surf.wind_l[i]*sin(surf.kinem.alpha))

        wlx = 0.5*(surf.uind_u[i]*cos(surf.kinem.alpha) - surf.wind_u[i]*sin(surf.kinem.alpha) -
                   surf.uind_l[i]*cos(surf.kinem.alpha) + surf.wind_l[i]*sin(surf.kinem.alpha))

        surf.RHS[i-1] = -(surf.kinem.u + curfield.u[1])*sin(surf.kinem.alpha)/surf.uref -
            surf.kinem.alphadot*(surf.x[i] - surf.pvt*surf.c)/surf.uref +
            (surf.kinem.hdot - curfield.w[1])*cos(surf.kinem.alpha)/surf.uref - wlz/surf.uref +
            surf.cam_slope[i]*((surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) +
                               (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha) + wtx -
                               surf.kinem.alphadot*surf.cam[i]) +
                               surf.thick_slope[i]*(wlx - surf.kinem.alphadot*surf.thick[i])

        surf.RHS[surf.ndiv+i-3] = surf.cam_slope[i]*(wlx - surf.kinem.alphadot*surf.thick[i]) +
            surf.thick_slope[i]*((surf.kinem.u + curfield.u[1])*cos(surf.kinem.alpha) +
                                 (surf.kinem.hdot - curfield.w[1])*sin(surf.kinem.alpha) + wtx -
                                 surf.kinem.alphadot*surf.cam[i]) - wtz
    end

    #RHS for Kelvin condition (negative strength of all previously shed vortices)
    surf.RHS[2*surf.ndiv-3] = -(sum(map(q->q.s, curfield.tev)) + sum(map(q->q.s, curfield.lev)))/(surf.uref*surf.c)

    if surf.a0[1] > 0.
        surf.RHS[2*surf.ndiv-2] = surf.lespcrit[1]
    else
        surf.RHS[2*surf.ndiv-2] = -surf.lespcrit[1]
    end
end

#Update position and strengths of bound vortices and sources
function update_bv_src(surf::TwoDSurfThick)
    gamma = zeros(surf.ndiv)
    src = zeros(surf.ndiv)
    for ib = 1:surf.ndiv
        gamma[ib] = surf.a0[1]*(1 + cos(surf.theta[ib]))
        for ia = 1:surf.naterm
            gamma[ib] += surf.aterm[ia]*sin(ia*surf.theta[ib])*sin(surf.theta[ib])
            src[ib] += surf.bterm[ia]*sin(ia*surf.theta[ib])*sin(surf.theta[ib])
        end
        gamma[ib] = gamma[ib]*surf.uref*surf.c
        src[ib] = src[ib]*surf.uref*surf.c
    end

    for ib = 2:surf.ndiv
        surf.bv[ib-1].s = (gamma[ib]+gamma[ib-1])*(surf.theta[2]-surf.theta[1])/2.
        surf.bv[ib-1].x = (surf.bnd_x_chord[ib] + surf.bnd_x_chord[ib-1])/2.
        surf.bv[ib-1].z = (surf.bnd_z_chord[ib] + surf.bnd_z_chord[ib-1])/2.
        surf.src[ib-1].s = (src[ib] + src[ib-1])*(surf.theta[2]-surf.theta[1])/2.
        surf.src[ib-1].x = surf.bv[ib-1].x
        surf.src[ib-1].z = surf.bv[ib-1].z
    end
end

# Function to calculate induced velocities by a set of sources at a target location
function ind_vel_src(src::Vector{TwoDSource},t_x,t_z)
    #'s' stands for source and 't' for target
    uind = zeros(length(t_x))
    wind = zeros(length(t_x))

    for itr = 1:length(t_x)
	for isr = 1:length(src)
            xdist = src[isr].x - t_x[itr]
            zdist = src[isr].z - t_z[itr]
            distsq = xdist*xdist + zdist*zdist
            uind[itr] = uind[itr] - src[isr].s*xdist/(2*pi*distsq)
            wind[itr] = wind[itr] - src[isr].s*zdist/(2*pi*distsq)
        end
    end
    return uind, wind
end


text:
# ========== Perform 1 step of discrete wavelet transform ==========
# ----- 1 step of dwt for 1D signals -----
"""
    dwt_step(v, h, g)

Perform one level of the discrete wavelet transform (DWT) on the vector `v`, which is the
`d`-th level scaling coefficients (Note the 0th level scaling coefficients is the raw
signal). The vectors `h` and `g` are the detail and scaling filters.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of DWT
DWT.dwt_step(v, h, g)
```

**See also:** [`dwt_step!`](@ref), [`idwt_step`](@ref)
"""
function dwt_step(v::AbstractVector{T}, h::Array{S,1}, g::Array{S,1}) where
                 {T<:Number, S<:Number}
    n = length(v)
    w₁ = zeros(T, n÷2)
    w₂ = zeros(T, n÷2)

    dwt_step!(w₁, w₂, v, h, g)
    return w₁, w₂
end

"""
    dwt_step!(w₂, w₂, v, h, g)

Same as `dwt_step` but without array allocation.

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.
- `v::AbstractVector{T} where T<:Number`: Vector of coefficients from a node at level `d`.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `w₁::Vector{T}`: Output from the low pass filter.
- `w₂::Vector{T}`: Output from the high pass filter.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)
w₁ = zeros(8)
w₂ = zeros(8)

# One step of DWT
DWT.dwt_step!(w₁, w₂, v, 0, h, g)
```

**See also:** [`dwt_step`](@ref), [`idwt_step`](@ref)
"""
function dwt_step!(w₁::AbstractVector{T},
                   w₂::AbstractVector{T},
                   v::AbstractVector{T},
                   h::Array{S,1},
                   g::Array{S,1}) where {T<:Number, S<:Number}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)÷2
    @assert length(h) == length(g)

    # Setup
    n = length(v)           # Parent length
    n₁ = length(w₁)         # Child length
    filtlen = length(h)     # Filter length

    # One step of discrete transform
    for i in 1:n₁
        k₁ = 2*i-1          # Start index for low pass filtering
        k₂ = 2*i            # Start index for high pass filtering
        @inbounds w₁[i] = g[end] * v[k₁]
        @inbounds w₂[i] = h[1] * v[k₂]
        for j in 2:filtlen
            k₁ = k₁+1 |> k₁ -> k₁>n ? mod1(k₁,n) : k₁
            k₂ = k₂-1 |> k₂ -> k₂≤0 ? mod1(k₂,n) : k₂
            @inbounds w₁[i] += g[end-j+1] * v[k₁]
            @inbounds w₂[i] += h[j] * v[k₂]
        end
    end
    return w₁, w₂
end

"""
    idwt_step(w₁, w₂, h, g)

Perform one level of the inverse discrete wavelet transform (IDWT) on the vectors `w₁` and
`w₂`, which are the scaling and detail coefficients. The vectors `h` and `g` are the detail
and scaling filters.

# Arguments
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `v::Vector{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of SDWT
w₁, w₂ = DWT.dwt_step(v, h, g)

# One step of ISDWT
v̂ = DWT.idwt_step(w₁, w₂, h, g)
```

**See also:** [`idwt_step!`](@ref)
"""
function idwt_step(w₁::AbstractVector{T},
                   w₂::AbstractVector{T},
                   h::Array{S,1},
                   g::Array{S,1}) where {T<:Number, S<:Number}
    n = length(w₁)
    v = Vector{T}(undef, 2*n)
    idwt_step!(v, w₁, w₂, h, g)
    return v
end

"""
    idwt_step!(v, w₁, w₂, h, g)

Same as `idwt_step` but without array allocation.

# Arguments
- `v::AbstractVector{T} where T<:Number`: Vector allocation for reconstructed coefficients.
- `w₁::AbstractVector{T} where T<:Number`: Vector allocation for output from low pass
  filter.
- `w₂::AbstractVector{T} where T<:Number`: Vector allocation for output from high pass
  filter.
- `h::Vector{S} where S<:Number`: High pass filter.
- `g::Vector{S} where S<:Number`: Low pass filter.

# Returns
- `v::Vector{T}`: Reconstructed coefficients.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8)
v̂ = similar(v)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of SDWT
w₁, w₂ = DWT.dwt_step(v, h, g)

# One step of ISDWT
DWT.idwt_step!(v̂, w₁, w₂, h, g)
```

**See also:** [`idwt_step`](@ref)
"""
function idwt_step!(v::AbstractVector{T},
                    w₁::AbstractVector{T},
                    w₂::AbstractVector{T},
                    h::Array{S,1},
                    g::Array{S,1}) where {T<:Number, S<:Number}
    # Sanity check
    @assert length(w₁) == length(w₂) == length(v)÷2
    @assert length(h) == length(g)

    # Setup
    n = length(v)           # Parent length
    n₁ = length(w₁)         # Child length
    filtlen = length(h)     # Filter length

    # One step of inverse discrete transform
    for i in 1:n
        j₀ = mod1(i,2)      # Pivot point to determine start index for filter
        j₁ = filtlen-j₀+1   # Index for low pass filter g
        j₂ = mod1(i+1,2)    # Index for high pass filter h
        k₁ = (i+1)>>1       # Index for approx coefs w₁
        k₂ = (i+1)>>1       # Index for detail coefs w₂
        @inbounds v[i] = g[j₁] * w₁[k₁] + h[j₂] * w₂[k₂]
        for j in (j₀+2):2:filtlen
            j₁ = filtlen-j+1
            j₂ = j + isodd(j) - iseven(j)
            k₁ = k₁-1 |> k₁ -> k₁≤0 ? mod1(k₁,n₁) : k₁
            k₂ = k₂+1 |> k₂ -> k₂>n₁ ? mod1(k₂,n₁) : k₂
            @inbounds v[i] += g[j₁] * w₁[k₁] + h[j₂] * w₂[k₂]
        end
    end
    return v
end

# ----- 1 step of dwt for 2D signals -----
"""
    dwt_step(v, h, g)

Compute 1 step of 2D discrete wavelet transform (DWT).

# Arguments
- `v::AbstractArray{T,2} where T<:Number`: Array of coefficients of size ``(n,m)``.
- `h::Array{S,1} where S<:Number`: High pass filter.
- `g::Array{S,1} where S<:Number`: Low pass filter.

# Returns
- `w₁::Array{T,2}`: Top left output. Result of low pass filter on columns + low pass filter
  on rows.
- `w₂::Array{T,2}`: Top left output. Result of low pass filter on columns + high pass filter
  on rows.
- `w₃::Array{T,2}`: Top left output. Result of high pass filter on columns + low pass filter
  on rows.
- `w₄::Array{T,2}`: Top left output. Result of high pass filter on columns + high pass filter
  on rows.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8,8)
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of DWT
DWT.dwt_step(v, h, g)
```

**See also:** [`dwt_step`](@ref), [`dwt_step!`](@ref), [`idwt_step`](@ref)
"""
function dwt_step(v::AbstractArray{T,2}, h::Array{S,1}, g::Array{S,1}) where
                 {T<:Number, S<:Number}
    n,m = size(v)
    n₁ = n÷2
    m₁ = m÷2
    w₁ = Array{T,2}(undef, (n₁,m₁))
    w₂ = Array{T,2}(undef, (n₁,m₁))
    w₃ = Array{T,2}(undef, (n₁,m₁))
    w₄ = Array{T,2}(undef, (n₁,m₁))
    temp = Array{T,2}(undef, (n,m))
    dwt_step!(w₁, w₂, w₃, w₄, v, h, g, temp)
    return w₁, w₂, w₃, w₄
end

"""
    dwt_step!(v, w₁, w₂, w₃, w₄, h, g, temp)

Same as 2D version of `dwt_step` but without array allocation.

# Arguments
- `w₁::AbstractArray{T,2} where T<:Number`: Array allocation for top left output.
- `w₂::AbstractArray{T,2} where T<:Number`: Array allocation for top right output.
- `w₃::AbstractArray{T,2} where T<:Number`: Array allocation for bottom left output.
- `w₄::AbstractArray{T,2} where T<:Number`: Array allocation for bottom right output.
- `v::AbstractArray{T,2} where T<:Number`: Array of coefficients to be transformed.
- `h::Array{S,1} where S<:Number`: High pass filter.
- `g::Array{S,1} where S<:Number`: Low pass filter.
- `temp::AbstractArray{T,2} where T<:Number`: Array allocation for intermediate
  computations.

# Returns
- `w₁::Array{T,2}`: Top left output. Result of low pass filter on columns + low pass filter
  on rows.
- `w₂::Array{T,2}`: Top left output. Result of low pass filter on columns + high pass filter
  on rows.
- `w₃::Array{T,2}`: Top left output. Result of high pass filter on columns + low pass filter
  on rows.
- `w₄::Array{T,2}`: Top left output. Result of high pass filter on columns + high pass filter
  on rows.

# Examples
```julia
using Wavelets, WaveletsExt

# Setup
v = randn(8,8)
temp = similar(v)
w₁ = Array{Float64,2}(undef, (4,4))
w₂ = Array{Float64,2}(undef, (4,4))
w₃ = Array{Float64,2}(undef, (4,4))
w₄ = Array{Float64,2}(undef, (4,4))
wt = wavelet(WT.haar)
g, h = WT.makereverseqmfpair(wt, true)

# One step of DWT
DWT.dwt_step!(w₁, w₂, w₃, w₄, v, h, g, temp)
```
"""
function dwt_step!(w₁::AbstractArray{T,2}, w₂::AbstractArray{T,2},
                   w₃::AbstractArray{T,2}, w₄::AbstractArray{T,2},
                   v::AbstractArray{T,2},
                   h::Array{S,1}, g::Array{S,1},
                   temp::AbstractArray{T,2}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(temp)
    @assert size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert size(w₁,1)*2 == size(v,1)
    @assert size(w₁,2)*2 == size(v,2)

    # Setup
    n,m = size(w₁)

    # --- Transform ---
    # Compute dwt for all columns
    for j in 1:(2*m)
        @inbounds temp₁ⱼ = @view temp[1:n,j]
        @inbounds temp₂ⱼ = @view temp[(n+1):end,j]
        @inbounds vⱼ = @view v[:,j]
        @inbounds dwt_step!(temp₁ⱼ, temp₂ⱼ, vⱼ, h, g)
    end
    # Compute dwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds dwt_step!(w₁ᵢ, w₂ᵢ, temp₁ᵢ, h, g)

        @inbounds temp₂ᵢ = @view temp[n+i,:]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds dwt_step!(w₃ᵢ, w₄ᵢ, temp₂ᵢ, h, g)
    end
    return w₁, w₂, w₃, w₄
end

"""
    idwt_step(w₁, w₂, w₃, w₄, h, g)

Computes one step of inverse discrete wavelet transform on 2D-signals.

# Arguments
- `w₁::AbstractArray{T,2} where T<:Number`: Top left child coefficients.
- `w₂::AbstractArray{T,2} where T<:Number`: Top right child coefficients.
- `w₃::AbstractArray{T,2} where T<:Number`: Bottom left child coefficients.
- `w₄::AbstractArray{T,2} where T<:Number`: Bottom right child coefficients.
- `h::Array{S,1} where S<:Number`: High pass filter.
- `g::Array{S,1} where S<:Number`: Low pass filter.

# Returns
- `::Array{T,2}`: Reconstructed coefficients.
"""
function idwt_step(w₁::AbstractArray{T,2}, w₂::AbstractArray{T,2},
                   w₃::AbstractArray{T,2}, w₄::AbstractArray{T,2},
                   h::Array{S,1}, g::Array{S,1}) where {T<:Number, S<:Number}
    n,m = size(w₁)
    v = Array{T,2}(undef, (2*n,2*m))
    temp = Array{T,2}(undef, (2*n,2*m))
    idwt_step!(v, w₁, w₂, w₃, w₄, h, g, temp)
    return v
end

"""
    idwt_step!(v, w₁, w₂, w₃, w₄, h, g, temp)

Same as 2D version of `idwt_step` but without array allocation.

# Arguments
- `v::AbstractArray{T,2} where T<:Number`: Array allocation for inverse transformed output.
- `w₁::AbstractArray{T,2} where T<:Number`: Array allocation for top left coefficients.
- `w₂::AbstractArray{T,2} where T<:Number`: Array allocation for top right coefficients.
- `w₃::AbstractArray{T,2} where T<:Number`: Array allocation for bottom left coefficients.
- `w₄::AbstractArray{T,2} where T<:Number`: Array allocation for bottom right coefficients.
- `h::Array{S,1} where S<:Number`: High pass filter.
- `g::Array{S,1} where S<:Number`: Low pass filter.
- `temp::AbstractArray{T,2} where T<:Number`: Array allocation for intermediate
  computations.

# Returns
- `v::Array{T,2}`: Reconstructed coefficients.
"""
function idwt_step!(v::AbstractArray{T,2},
                    w₁::AbstractArray{T,2}, w₂::AbstractArray{T,2},
                    w₃::AbstractArray{T,2}, w₄::AbstractArray{T,2},
                    h::Array{S,1}, g::Array{S,1},
                    temp::AbstractArray{T,2}) where {T<:Number, S<:Number}
    # Sanity check
    @assert size(v) == size(temp)
    @assert size(w₁) == size(w₂) == size(w₃) == size(w₄)
    @assert size(w₁,1)*2 == size(v,1)
    @assert size(w₁,2)*2 == size(v,2)

    # Setup
    n,m = size(w₁)

    # --- Inverse transform ---
    # Compute idwt for all rows
    for i in 1:n
        @inbounds temp₁ᵢ = @view temp[i,:]
        @inbounds w₁ᵢ = @view w₁[i,:]
        @inbounds w₂ᵢ = @view w₂[i,:]
        @inbounds idwt_step!(temp₁ᵢ, w₁ᵢ, w₂ᵢ, h, g)

        @inbounds temp₂ᵢ = @view temp[n+i,:]
        @inbounds w₃ᵢ = @view w₃[i,:]
        @inbounds w₄ᵢ = @view w₄[i,:]
        @inbounds idwt_step!(temp₂ᵢ, w₃ᵢ, w₄ᵢ, h, g)
    end
    # Compute idwt for all columns
    for j in 1:(2*m)
        @inbounds vⱼ = @view v[:,j]
        @inbounds temp₁ⱼ = @view temp[1:n,j]
        @inbounds temp₂ⱼ = @view temp[(n+1):end,j]
        @inbounds idwt_step!(vⱼ, temp₁ⱼ, temp₂ⱼ, h, g)
    end
    return v
end

# ----- 1 step of dwt for nD signals -----
# TODO: Implement this function