text:
"""
### mmread(filename, infoonly::Bool=false, retcoord::Bool=false)
Read the contents of the Matrix Market file 'filename' into a matrix,
which will be either sparse or dense, depending on the Matrix Market format
indicated by 'coordinate' (coordinate sparse storage), or 'array' (dense
array storage).
If infoonly is true (default: false), only information on the size and
structure is returned from reading the header. The actual data for the
matrix elements are not parsed.
If retcoord is true (default: false), the rows, column and value vectors
are returned, if it is a sparse matrix, along with the header information.
"""
function mmread(mmstream::TranscodingStream, infoonly::Bool=false, retcoord::Bool=false)
    # Read first line
    firstline = chomp(readline(mmstream))
    tokens = split(firstline)
    if length(tokens) != 5
        throw(ParseError(string("Not enough words on first line: ", firstline)))
    end
    if tokens[1] != "%%MatrixMarket"
        throw(ParseError(string("Expected start of header `%%MatrixMarket`, got `$(tokens[1])`")))
    end
    (head1, rep, field, symm) = map(lowercase, tokens[2:5])
    if head1 != "matrix"
        throw(ParseError("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)"))
    end

    eltype = field == "real" ? Float64 :
             field == "complex" ? ComplexF64 :
             field == "integer" ? Int64 :
             field == "pattern" ? Bool :
             throw(ParseError("Unsupported field $field (only real and complex are supported)"))

    symlabel = symm == "general" ? identity :
               symm == "symmetric" ? symmetric! :
               symm == "hermitian" ? hermitian! :
               symm == "skew-symmetric" ? skewsymmetric! :
               throw(ParseError("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)"))

    # Skip all comments and empty lines
    ll   = readline(mmstream)
    while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
        ll = readline(mmstream)
    end
    # Read matrix dimensions (and number of entries) from first non-comment line
    dd = map(_parseint, split(ll))
    if length(dd) < (rep == "coordinate" ? 3 : 2)
        throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
    end
    rows = dd[1]
    cols = dd[2]
    entries = (rep == "coordinate") ? dd[3] : (rows * cols)
    infoonly && return (rows, cols, entries, rep, field, symm)

    rep == "coordinate" ||
        return symlabel(reshape([parse(Float64, readline(mmstream)) for i in 1:entries],
                                (rows,cols)))

    rr = Vector{Int}(undef, entries)
    cc = Vector{Int}(undef, entries)
    xx = Vector{eltype}(undef, entries)
    for i in 1:entries
        line = readline(mmstream)
        splits = find_splits(line, eltype == ComplexF64 ? 3 : (eltype == Bool ? 1 : 2))
        rr[i] = _parseint(line[1:splits[1]])
        cc[i] = _parseint(eltype == Bool
                          ? line[splits[1]:end]
                          : line[splits[1]:splits[2]])
        if eltype == ComplexF64
            real = parse(Float64, line[splits[2]:splits[3]])
            imag = parse(Float64, line[splits[3]:length(line)])
            xx[i] = ComplexF64(real, imag)
        elseif eltype == Bool
            xx[i] = true
        else
            xx[i] = parse(eltype, line[splits[2]:length(line)])
        end
    end
    (retcoord
     ? (rr, cc, xx, rows, cols, entries, rep, field, symm)
     : symlabel(sparse(rr, cc, xx, rows, cols)))
end


text:
using JuliaFormatter

# Test formatter on single script
format("./src/04_full-extent.jl", BlueStyle())
format("./src/04_full-extent.jl", BlueStyle(); pipe_to_function_call=false)

# Use formatter on main scripts
mainfiles = readdir("./src/"; join=true)
filter!(endswith(".jl"), mainfiles)
show(stdout, "text/plain", mainfiles)
problems = []
function formatfiles(files)
    for file in files
        @info "Formatting $(file)"
        try
            format(file, BlueStyle(); pipe_to_function_call=false)
        catch
            @warn "Error while $(file)"
            push!(problems, file)
            continue
        end
        return problems
    end
end
formatfiles(mainfiles)

# Investigate problematic files
format(problems, BlueStyle(); pipe_to_function_call=false)
format(mainfiles[5], BlueStyle(); pipe_to_function_call=false)

# Use formatter on lib scripts
libfiles = readdir("./src/lib/"; join=true)
filter!(endswith(".jl"), libfiles)
formatfiles(libfiles)

# Use formatter on other scripts
otherfiles = readdir("./src/others/"; join=true)
filter!(endswith(".jl"), otherfiles)
formatfiles(otherfiles)
unique(problems)

ext: jl

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
function compute_gap!(residuals::Residuals, pair::PrimalDual, a::AuxiliaryData, aff::AffineSets, p::Params)::Nothing

    # Inplace primal feasibility error
    if aff.p > 0
        residuals.equa_feasibility = 0.
        @simd for i in 1:aff.p
            @inbounds residuals.equa_feasibility = max(residuals.equa_feasibility, abs(a.Mx[i] - aff.b[i]))
        end
        residuals.equa_feasibility /= (1. + p.norm_b)
    end
    if aff.m > 0
        residuals.ineq_feasibility = 0.
        @simd for i in aff.p+1:aff.p+aff.m
            @inbounds residuals.ineq_feasibility = max(residuals.ineq_feasibility, a.Mx[i] - aff.h[i-aff.p])
        end
        residuals.ineq_feasibility /= (1. + p.norm_h)
    end
    residuals.feasibility[p.iter] = max(residuals.equa_feasibility, residuals.ineq_feasibility)

    # Primal-dual gap
    residuals.prim_obj[p.iter] = dot(aff.c, pair.x)
    residuals.dual_obj[p.iter] = 0.
    if aff.p > 0
        residuals.dual_obj[p.iter] -= dot(aff.b, @view pair.y[1:aff.p])
    end
    if aff.m > 0
        residuals.dual_obj[p.iter] -= dot(aff.h, @view pair.y[aff.p+1:end])
    end
    residuals.dual_gap[p.iter] =
        abs(residuals.prim_obj[p.iter] - residuals.dual_obj[p.iter]) /
            (1. + abs(residuals.prim_obj[p.iter]) + abs(residuals.dual_obj[p.iter]))

    return nothing
end

function compute_residual!(residuals::Residuals, pair::PrimalDual, a::AuxiliaryData, p::Params, aff::AffineSets)::Nothing

    # Primal residual
    # Px_old
    a.Mty_old .= pair.x_old .- p.primal_step .* a.Mty_old
    # Px
    pair.x_old .= pair.x .- p.primal_step .* a.Mty
    # Px - Px_old
    pair.x_old .-= a.Mty_old
    residuals.primal_residual[p.iter] = sqrt(aff.n) * norm(pair.x_old, Inf) / max(norm(a.Mty_old, Inf), p.norm_b, p.norm_h, 1.)

    # Dual residual
    # Py_old
    a.Mx_old .= pair.y_old .- p.dual_step .* a.Mx_old
    # Py
    pair.y_old .= pair.y .- p.dual_step .* a.Mx
    # Py - Py_old
    pair.y_old .-= a.Mx_old
    residuals.dual_residual[p.iter] = sqrt(aff.m + aff.p) * norm(pair.y_old, Inf) / max(norm(a.Mx_old, Inf), p.norm_c, 1.)

    # Compute combined residual
    residuals.comb_residual[p.iter] = max(residuals.primal_residual[p.iter], residuals.dual_residual[p.iter])

    # Keep track of previous iterates
    copyto!(pair.x_old, pair.x)
    copyto!(pair.y_old, pair.y)
    copyto!(a.Mty_old, a.Mty)
    copyto!(a.Mx_old, a.Mx)

    return nothing
end

function soc_convergence(a::AuxiliaryData, cones::ConicSets, pair::PrimalDual, opt::Options, p::Params)::Bool
    for (idx, soc) in enumerate(cones.socone)
        if soc_gap(a.soc_v[idx], a.soc_s[idx]) >= opt.tol_soc
            return false
        end
    end

    return true
end

function soc_gap(v::ViewVector, s::ViewScalar)

    return norm(v, 2) - s[]
end

function convergedrank(a, p::Params, cones::ConicSets, opt::Options)::Bool
    for (idx, sdp) in enumerate(cones.sdpcone)
        if !(
                sdp.sq_side < opt.min_size_krylov_eigs ||
                p.target_rank[idx] > opt.max_target_rank_krylov_eigs ||
                min_eig(a, idx, p) < opt.tol_psd
            )
            # @show min_eig(a, idx, p), -opt.tol_psd
            return false
        end
    end

    return true
end


text:
export OffPolicy

"""
    OffPolicy(π_target::P, π_behavior::B) -> OffPolicy{P,B}
"""
Base.@kwdef struct OffPolicy{P,B} <: AbstractPolicy
    π_target::P
    π_behavior::B
end

(π::OffPolicy)(obs) = π.π_behavior(obs)

function RLBase.update!(π::OffPolicy, t::AbstractTrajectory)
    experience = extract_experience(t, π)
    isnothing(experience) || update!(π, experience)
end

function RLBase.update!(π::OffPolicy{<:VBasedPolicy}, transitions::NamedTuple)
    # ??? define a `get_batch_prob` function for efficiency
    weights = [
        get_prob(π.π_target, (state = s,), a) / get_prob(π.π_behavior, (state = s,), a)
        for (s, a) in zip(transitions.states, transitions.actions)
    ]  # TODO: implement iterate interface for (SubArray of) CircularArrayBuffer
    experience = merge(transitions, (weights = weights,))
    update!(π.π_target, experience)
end

RLBase.extract_experience(t::AbstractTrajectory, π::OffPolicy) =
    extract_experience(t, π.π_target)


text:
import GLAbstraction: gluniform

@render_program DefaultProgram
@render_program InstancedDefaultProgram

@vao DefaultVao
@instanced_vao InstancedDefaultVao

struct DefaultRenderer <: AbstractRenderSystem end

function Overseer.requested_components(::DefaultRenderer)
    return (DefaultVao, DefaultProgram, InstancedDefaultVao, InstancedDefaultProgram,
            ModelMat, Material, PointLight, UniformColor, BufferColor, Spatial, Camera3D,
            IOTarget)
end

function Overseer.update(::DefaultRenderer, m::AbstractLedger)
    fbo   = singleton(m, IOTarget)
    prog  = singleton(m, DefaultProgram)
    iprog = singleton(m, InstancedDefaultProgram)

    bind(fbo)
    draw(fbo)
    glDisable(GL_BLEND)
    glEnable(GL_DEPTH_TEST)
    glDepthFunc(GL_LEQUAL)
    glEnableCullFace(:back)
    #Render instanced-renderables

    bind(prog)

    light, ucolor, spat, cam, modelmat, material, vao, vis = m[PointLight], m[UniformColor],
                                                             m[Spatial], m[Camera3D],
                                                             m[ModelMat], m[Material],
                                                             m[DefaultVao], m[Visible]

    set_light_camera_uniforms = (prog) -> begin
        for e in @entities_in(m, PointLight && UniformColor && Spatial)
            gluniform(prog, e[PointLight], e[UniformColor], e[Spatial])
        end

        for e in @entities_in(spat && cam)
            gluniform(prog, spat[e], cam[e])
        end
    end
    set_light_camera_uniforms(prog)

    set_model_material = (e_modelmat, e_material) -> begin
        gluniform(prog, :material, Vec2(e_material.specpow, e_material.specint))
        gluniform(prog, :modelmat, e_modelmat.modelmat)
    end

    idcolor = m[Selectable]
    #Colors inside Vao
    for e in @entities_in(vao && modelmat && material && !ucolor)
        evao = vao[e]
        if vis[e].visible
            set_model_material(modelmat[e], material[e])
            if e in idcolor
                gluniform(prog, :object_id_color, idcolor[e].color)
            end
            GLA.bind(evao)
            GLA.draw(evao)
        end
    end

    bind(iprog)
    set_light_camera_uniforms(iprog)
    ivao = m[InstancedDefaultVao]
    for evao in ivao.data
        GLA.bind(evao)
        GLA.draw(evao)
    end
end


text:
using Documenter, UnROOT

makedocs(;
    modules=[UnROOT],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets=String[],
    ),
    pages=[
        "Introduction" => "index.md",
        "Example Usage" => "exampleusage.md",
        "Performance Tips" => "performancetips.md",
        "Advanced Usage" => [
            "Parse Custom Branch" => "advanced/custom_branch.md",
            "Reduce startup latency" => "advanced/reduce_latency.md",
        ],
        "For Contributors" => "devdocs.md",
        "APIs" => "internalapis.md",
    ],
    repo="https://github.com/JuliaHEP/UnROOT.jl/blob/{commit}{path}#L{line}",
    sitename="UnROOT.jl",
    authors="Tamas Gal and contributors",
)

deploydocs(;
    repo="github.com/JulieHEP/UnROOT.jl",
)


text:
#Function for estimating a problem's time step

#Simple linear interpolation function
function interp(x1 ::Float64, x2 :: Float64, y1 :: Float64, y2 :: Float64, x::Float64)
    y = y1 + (y2 - y1)*(x - x1)/(x2 - x1)
    return y
end

function find_tstep(kin:: Array{CosDef})
    dtstar = 0.015
    for i = 1:length(kin)
        dt_tmp = 0.015*0.2/(kin[i].k*kin[i].amp)
        dtstar = minimum([dt_tmp dtstar])
    end
    return dtstar
end

function find_tstep(kin:: EldUpDef)
    dtstar = 0.015
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

function find_tstep(kin:: EldUptstartDef)
    dtstar = 0.015
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

function find_tstep(kin:: EldRampReturnDef)
    dtstar = 0.015
    dtstar = minimum([0.015*0.2/kin.K 0.015])
    return dtstar
end

#Function for calculating a simple 3D correction based on quasi-steady LLT

function simple_LLT(mat::Array{Float64,2}, surf::TwoDSurf, s::Float64, n_bterm::Int64, n_span::Int64)
    # s is the semi-span
    nsteps = length(mat[:,1])

    bnd_circ2d = mat[:,9]
    t = mat[:,1]
    alpha = mat[:,2]
    h = mat[:,3]
    hdot = zeros(nsteps)
    hdot[1] = 0
    hdot[2:end] = diff(h)./diff(t)
    u = mat[:,4]

    lhs = zeros(n_span,n_bterm)
    rhs = zeros(n_span)
    b_coeff = zeros(n_bterm)
    bcoeff_prev = zeros(n_bterm)
    bdot = zeros(n_bterm)
    sp_gam = zeros(nsteps,n_span)
    psi = zeros(n_span)
    dt = mat[2,1]-mat[1,1]
    cnc_f = zeros(nsteps)
    cnnc_f = zeros(nsteps)

    for i = 1:n_span
        psi[i] = (real(i)/(n_span+1))*pi
    end

    for i = 1:nsteps
        for j = 1:n_span
            for n = 1:n_bterm
                lhs[j,n] = sin(n*psi[j])*(sin(psi[j]) + (n*surf.c*pi/(4*s)))
            end
            rhs[j] = surf.c*pi*sin(psi[j])*bnd_circ2d[i]/(4*s)
        end
        bcoeff_prev = b_coeff
        b_coeff = \(lhs, rhs)
        bdot = (b_coeff - bcoeff_prev)/dt

        # for j = 1:n_span
        #     sp_gam[i,j] = 0
        #     for n = 1:n_bterm
        #         sp_gam[i,j] = sp_gam[i,j] + 4*s*u*b_coeff[n]*sin(n*psi[j])
        #     end
        # end

        sum_bcoeff = 0
        for n = 1:n_bterm
            if rem(n,2) != 0
                sum_bcoeff = sum_bcoeff + b_coeff[n]
            end
        end
        cnc_f[i] = -2*pi*(u[i]*cos(alpha[i])/surf.uref + hdot[i]*sin(alpha[i])/surf.uref)*(sum_bcoeff)
        sum_bdot = 0
        for n = 1:n_bterm
            if rem(n,2) != 0
                sum_bdot = sum_bdot + bdot[n]
            end
        end
        cnnc_f[i] = -(2*pi*surf.c/(surf.uref))*(3*sum_bdot/4)
    end
    return cnc_f, cnnc_f
end

# ---------------------------------------------------------------------------------------------
# Function for calculating a_2 and a_3 fourier coefficients
function update_a2a3adot(surf::TwoDSurf,dt)
    for ia = 2:3
        surf.aterm[ia] = trapz(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_a2a3adot(surf::TwoDSurf_2DOF,dt)
    for ia = 2:3
        surf.aterm[ia] = trapz(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_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 update_adot(surf::TwoDSurf_2DOF,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 update_a2a3adot(surf::TwoDFreeSurf,dt)
    for ia = 2:3
        surf.aterm[ia] = trapz(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_a2a3adot(surf::TwoDSurfwFlap,dt)
    for ia = 2:3
        surf.aterm[ia] = trapz(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 for updating the induced velocities u, w: eqns (2.12), (2.13) in Ramesh et al.
function update_indbound(surf::TwoDSurf, curfield::TwoDFlowField)
    surf.uind[1:surf.ndiv], surf.wind[1:surf.ndiv] = ind_vel([curfield.tev; curfield.lev], surf.bnd_x, surf.bnd_z)
    return surf
end

function update_indbound(surf::TwoDSurfwFlap, curfield::TwoDFlowField)
    surf.uind[1:surf.ndiv], surf.wind[1:surf.ndiv] = ind_vel([curfield.tev; curfield.lev], surf.bnd_x, surf.bnd_z)
    return surf
end

function update_indbound(surf::TwoDSurf_2DOF, curfield::TwoDFlowField)
    surf.uind[1:surf.ndiv], surf.wind[1:surf.ndiv] = ind_vel([curfield.tev; curfield.lev], surf.bnd_x, surf.bnd_z)
    return surf
end

function update_indbound(surf::TwoDFreeSurf, 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 for updating the downwash W: eqn (2.5) in Ramesh et al.
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 update_downwash(surf::TwoDSurf_2DOF)
    for ib = 1:surf.ndiv
        surf.downwash[ib] = -surf.kinem.u*sin(surf.kinem.alpha) - surf.uind[ib]*sin(surf.kinem.alpha) + surf.kinem.hdot*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*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))
    end
    return surf
end

function update_downwash(surf::TwoDFreeSurf)
    for ib = 1:surf.ndiv
        surf.downwash[ib] = -surf.kinem.u*sin(surf.kinem.alpha) - surf.uind[ib]*sin(surf.kinem.alpha) + surf.kinem.hdot*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*cos(surf.kinem.alpha) + surf.kinem.hdot*sin(surf.kinem.alpha) - surf.wind[ib]*sin(surf.kinem.alpha))
    end
    return surf
end
# ---------------------------------------------------------------------------------------------

# Function for updating the downwash W: eqn (2.5) in Ramesh et al.
# Additional terms: alphadot*eta and time derivative of eta to account for deformations
# Added by Laura Merchant 2016
function update_downwash(surf::TwoDSurfwFlap)
    for ib = 1:surf.ndiv
        surf.downwash[ib] = -surf.kinem.u*sin(surf.kinem.alpha) - surf.uind[ib]*sin(surf.kinem.alpha) + surf.kinem.hdot*cos(surf.kinem.alpha) - surf.wind[ib]*cos(surf.kinem.alpha) - surf.kinem.alphadot*(surf.x[ib] - surf.pvt*surf.c) + surf.cam_tder[ib] + surf.cam_slope[ib]*(surf.uind[ib]*cos(surf.kinem.alpha) + surf.kinem.u*cos(surf.kinem.alpha) - surf.kinem.alphadot*surf.cam[ib] + surf.kinem.hdot*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] = trapz(surf.downwash,surf.theta)
    surf.aterm[1] = trapz(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 update_a0anda1(surf::TwoDSurf_2DOF)
    surf.a0[1] = trapz(surf.downwash,surf.theta)
    surf.aterm[1] = trapz(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 update_a0anda1(surf::TwoDFreeSurf)
    surf.a0[1] = trapz(surf.downwash,surf.theta)
    surf.aterm[1] = trapz(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 update_a0anda1(surf::TwoDSurfwFlap)
    surf.a0[1] = trapz(surf.downwash,surf.theta)
    surf.aterm[1] = trapz(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] = trapz(surf.downwash.*cos(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2.*surf.aterm[ia]/(surf.uref*pi)
    end
    return surf
end

function update_a2toan(surf::TwoDSurf_2DOF)
    for ia = 2:surf.naterm
        surf.aterm[ia] = trapz(surf.downwash.*cos(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2.*surf.aterm[ia]/(surf.uref*pi)
    end
    return surf
end

function update_a2toan(surf::TwoDFreeSurf)
    for ia = 2:surf.naterm
        surf.aterm[ia] = trapz(surf.downwash.*cos(ia*surf.theta),surf.theta)
        surf.aterm[ia] = 2.*surf.aterm[ia]/(surf.uref*pi)
    end
    return surf
end

function update_a2toan(surf::TwoDSurfwFlap)
    for ia = 2:surf.naterm
        surf.aterm[ia] = trapz(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
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
    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
    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.
    end
    # ---------------------------------------------------------------------------------------------
    return surf
end
# END kinem function
# ---------------------------------------------------------------------------------------------

function update_kinem(surf::TwoDSurf_2DOF, dt)
    #Kinematics are updated according to the 2DOF response
    surf.kinem.alphadot = surf.kinem.alphadot_pr + (dt/12.)*(23*surf.kinem.alphaddot_pr - 16*surf.kinem.alphaddot_pr2 + 5*surf.kinem.alphaddot_pr3)
    surf.kinem.hdot = surf.kinem.hdot_pr + (dt/12)*(23*surf.kinem.hddot_pr-16*surf.kinem.hddot_pr2 + 5*surf.kinem.hddot_pr3)
    surf.kinem.alpha = surf.kinem.alpha_pr + (dt/12)*(23*surf.kinem.alphadot_pr-16*surf.kinem.alphadot_pr2 + 5*surf.kinem.alphadot_pr3)
    surf.kinem.h = surf.kinem.h_pr + (dt/12)*(23*surf.kinem.hdot_pr - 16*surf.kinem.hdot_pr2 + 5*surf.kinem.hdot_pr3)
end

function update_kinem(surf::TwoDFreeSurf, dt)
    #Kinematics are updated according to the 2DOF response
    surf.kinem.alphadot = surf.kinem.alphadot_pr + (dt/12.)*(23*surf.kinem.alphaddot_pr - 16*surf.kinem.alphaddot_pr2 + 5*surf.kinem.alphaddot_pr3)
    surf.kinem.hdot = surf.kinem.hdot_pr + (dt/12)*(23*surf.kinem.hddot_pr-16*surf.kinem.hddot_pr2 + 5*surf.kinem.hddot_pr3)
    surf.kinem.alpha = surf.kinem.alpha_pr + (dt/12)*(23*surf.kinem.alphadot_pr-16*surf.kinem.alphadot_pr2 + 5*surf.kinem.alphadot_pr3)
    surf.kinem.h = surf.kinem.h_pr + (dt/12)*(23*surf.kinem.hdot_pr - 16*surf.kinem.hdot_pr2 + 5*surf.kinem.hdot_pr3)
    surf.kinem.u = surf.kinem.u_pr + (dt/12)*(23*surf.kinem.udot_pr - 16*surf.kinem.udot_pr2 + 5*surf.kinem.udot_pr3)
end

# ---------------------------------------------------------------------------------------------
# Function updating the kinematic parameters
function update_kinem(surf::TwoDSurfwFlap, t)
# Added by Laura Merchant 2016
    # ---------------------------------------------------------------------------------------------
    # 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) == 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
    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) == EldUpIntDef)
        surf.kinem.h = surf.kindef.h(t)*surf.c
        surf.kinem.hdot = EldUpDef(30,0.2,0.8)(t)*surf.kindef.h.amp*surf.uref
        #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
    end

    # ---------------------------------------------------------------------------------------------
    # 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.
    end

    # ---------------------------------------------------------------------------------------------
    # Deformation kinematics
    if (typeof(surf.kindef.n) == CosDef)
        surf.kinem.n = surf.kindef.n(t)
        surf.kinem.ndot = ForwardDiff.derivative(surf.kindef.n,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.n) == EldRampReturnDef)
        surf.kinem.n = surf.kindef.n(t)
        surf.kinem.ndot = ForwardDiff.derivative(surf.kindef.n,t)*surf.uref/surf.c
    elseif (typeof(surf.kindef.n) == ConstDef)
        surf.kinem.n = surf.kindef.n(t)
        surf.kinem.ndot = 0.
    elseif (typeof(surf.kindef.n) == SinDef)
        surf.kinem.n = surf.kindef.n(t)
        surf.kinem.ndot = ForwardDiff.derivative(surf.kindef.n,t)*surf.uref/surf.c
    end

    return surf
end
# END kinemwFlap function
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# 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 update_bv(surf::TwoDSurf_2DOF)
    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 update_bv(surf::TwoDFreeSurf)
    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

# ---------------------------------------------------------------------------------------------
# Update KinemPar2DOF for 2DOF simulations
function update_kinem2DOF(surf::TwoDSurf_2DOF)
    surf.kinem.alpha_pr = surf.kinem.alpha
    surf.kinem.alpha_pr2 = surf.kinem.alpha_pr
    surf.kinem.alpha_pr3 = surf.kinem.alpha_pr2
    surf.kinem.h_pr = surf.kinem.h
    surf.kinem.h_pr2 = surf.kinem.h_pr
    surf.kinem.h_pr3 = surf.kinem.h_pr2
    surf.kinem.alphadot_pr = surf.kinem.alphadot
    surf.kinem.alphadot_pr2 = surf.kinem.alphadot_pr
    surf.kinem.alphadot_pr3 = surf.kinem.alphadot_pr2
    surf.kinem.hdot_pr = surf.kinem.hdot
    surf.kinem.hdot_pr2 = surf.kinem.hdot_pr
    surf.kinem.hdot_pr3 = surf.kinem.hdot_pr2
    surf.kinem.alphaddot_pr = surf.kinem.alphaddot
    surf.kinem.alphaddot_pr2 = surf.kinem.alphaddot_pr
    surf.kinem.alphaddot_pr3 = surf.kinem.alphaddot_pr2
    surf.kinem.hddot_pr = surf.kinem.hddot
    surf.kinem.hddot_pr2 = surf.kinem.hddot_pr
    surf.kinem.hddot_pr3 = surf.kinem.hddot_pr2
    return surf
end
# ---------------------------------------------------------------------------------------------
# Update KinemPar2DOF for Free simulations
function update_kinem2DOF(surf::TwoDFreeSurf)
    surf.kinem.alpha_pr = surf.kinem.alpha
    surf.kinem.h_pr = surf.kinem.h
    surf.kinem.u_pr = surf.kinem.u
    surf.kinem.alphadot_pr = surf.kinem.alphadot
    surf.kinem.alphadot_pr2 = surf.kinem.alphadot_pr
    surf.kinem.alphadot_pr3 = surf.kinem.alphadot_pr2
    surf.kinem.udot_pr = surf.kinem.udot
    surf.kinem.udot_pr2 = surf.kinem.udot_pr
    surf.kinem.udot_pr3 = surf.kinem.udot_pr2
    surf.kinem.hdot_pr = surf.kinem.hdot
    surf.kinem.hdot_pr2 = surf.kinem.hdot_pr
    surf.kinem.hdot_pr3 = surf.kinem.hdot_pr2
    surf.kinem.alphaddot_pr = surf.kinem.alphaddot
    surf.kinem.alphaddot_pr2 = surf.kinem.alphaddot_pr
    surf.kinem.alphaddot_pr3 = surf.kinem.alphaddot_pr2
    surf.kinem.hddot_pr = surf.kinem.hddot
    surf.kinem.hddot_pr2 = surf.kinem.hddot_pr
    surf.kinem.hddot_pr3 = surf.kinem.hddot_pr2
    return surf
end
# ---------------------------------------------------------------------------------------------

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

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

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

    surf.kinem.hddot = (1/(m11*m22-m21*m12))*(m22*R1-m12*R2)
    surf.kinem.alphaddot = (1/(m11*m22-m21*m12))*(-m21*R1+m11*R2)
    return surf
end

function calc_moveFree(surf::TwoDFreeSurf, cl::Float64, cd :: Float64,
cm :: Float64, cf :: Float64)
    accl_g = 9.8

    surf.kinem.udot = -2*surf.strpar.kappa*surf.uref*surf.uref*(cd+cf*sin(surf.kinem.alpha))/(pi*surf.c)
    surf.kinem.hddot = 2*surf.strpar.kappa*surf.uref*surf.uref*(cl+cf*cos(surf.kinem.alpha))/(pi*surf.c) - accl_g
    surf.kinem.alphaddot = 8*surf.strpar.kappa*surf.uref*surf.uref*cm/(pi*surf.c*surf.c*surf.strpar.r_g*surf.strpar.r_g)
end

function update_bv(surf::TwoDSurfwFlap)
    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 = 1:surf.ndiv-1
       surf.bv_prev[ib].s = surf.bv[ib].s
    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
    return surf
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# Numerical integration method: Trapezoidal
function trapz{T<:Real}(y::Vector{T}, x::Vector{T})
    local len = length(y)
    if (len != length(x))
        error("Vectors must be of same length")
    end
    r = 0.0
    for i in 2:len
        r += (x[i] - x[i-1]) * (y[i] + y[i-1])
    end
    r/2.0
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# Aerofoil camber calculation from coordinate file
function camber_calc(x::Vector,airfoil::String)
    #Determine camber and camber slope on airfoil from airfoil input file

    ndiv = length(x);
    c = x[ndiv];

    cam = Array(Float64,ndiv)
    cam_slope = Array(Float64,ndiv)

    in_air = readdlm(airfoil);
    xcoord = in_air[:,1];
    ycoord = in_air[:,2];
    ncoord = length(xcoord);
    xcoord_sum = Array(Float64,ncoord);
    xcoord_sum[1] = 0;
    for i = 1:ncoord-1
        xcoord_sum[i+1] = xcoord_sum[i] + abs(xcoord[i+1]-xcoord[i]);
    end
    y_spl = Spline1D(xcoord_sum,ycoord);
    y_ans = Array(Float64,2*ndiv);

    for i=1:ndiv
        y_ans[i] = evaluate(y_spl,x[i]/c);
    end
    for i=ndiv+1:2*ndiv
        y_ans[i] = evaluate(y_spl,(x[ndiv]/c) + (x[i-ndiv]/c));
    end
    cam[1:ndiv] = [(y_ans[i] + y_ans[(2*ndiv) + 1 - i])*c/2 for i = ndiv:-1:1];
    cam[1] = 0;
    cam_spl = Spline1D(x,cam);
    cam_slope[1:ndiv] = Dierckx.derivative(cam_spl,x);
    return cam, cam_slope
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# Updating the deformation of the aerofoil: camber, spatial and time derivatives (from flap)
# Added by Laura Merchant 2016
function update_deform(surf::TwoDSurfwFlap, t)
    ndiv = length(surf.x);

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

    cam_spl = Spline1D(surf.x,surf.cam);
    surf.cam_slope[1:surf.ndiv] = Dierckx.derivative(cam_spl,surf.x);
    return surf
end
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# Function for calculating the wake rollup
function wakeroll(surf::TwoDSurf, curfield::TwoDFlowField, dt)
    #Clean induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:length(curfield.lev)
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

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

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

    for i = 1:length(curfield.tev)
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = length(curfield.tev)+1:length(utemp)
        curfield.lev[i-length(curfield.tev)].vx += utemp[i]
        curfield.lev[i-length(curfield.tev)].vz += wtemp[i]
    end

    #Convect free vortices with their induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:length(curfield.lev)
            curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    return curfield
end

function wakeroll(surf::TwoDSurf_2DOF, curfield::TwoDFlowField, dt)
    #Clean induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:length(curfield.lev)
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

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

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

    for i = 1:length(curfield.tev)
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = length(curfield.tev)+1:length(utemp)
        curfield.lev[i-length(curfield.tev)].vx += utemp[i]
        curfield.lev[i-length(curfield.tev)].vz += wtemp[i]
    end

    #Convect free vortices with their induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:length(curfield.lev)
            curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    return curfield
end

function wakeroll(surf::TwoDFreeSurf, curfield::TwoDFlowField, dt)
    #Clean induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:length(curfield.lev)
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

    #Assume external vortex is not affected by flow vortices
    #Velocities induced by free vortices on each other
    mutual_ind([curfield.tev; curfield.lev])

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

    for i = 1:length(curfield.tev)
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = length(curfield.tev)+1:length(utemp)
        curfield.lev[i-length(curfield.tev)].vx += utemp[i]
        curfield.lev[i-length(curfield.tev)].vz += wtemp[i]
    end

    #Convect free vortices with their induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:length(curfield.lev)
        curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[i].vz
    end
    for i = 1:length(curfield.extv)
        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::TwoDSurfwFlap, curfield::TwoDFlowField, dt)
    #Clean induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].vx = 0
        curfield.tev[i].vz = 0
    end

    for i = 1:length(curfield.lev)
        curfield.lev[i].vx = 0
        curfield.lev[i].vz = 0
    end

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

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

    for i = 1:length(curfield.tev)
        curfield.tev[i].vx += utemp[i]
        curfield.tev[i].vz += wtemp[i]
    end
    for i = length(curfield.tev)+1:length(utemp)
        curfield.lev[i-length(curfield.tev)].vx += utemp[i]
        curfield.lev[i-length(curfield.tev)].vz += wtemp[i]
    end

    #Convect free vortices with their induced velocities
    for i = 1:length(curfield.tev)
        curfield.tev[i].x += dt*curfield.tev[i].vx
        curfield.tev[i].z += dt*curfield.tev[i].vz
    end
    for i = 1:length(curfield.lev)
            curfield.lev[i].x += dt*curfield.lev[i].vx
        curfield.lev[i].z += dt*curfield.lev[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::TwoDSurf_2DOF,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::TwoDFreeSurf,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::TwoDSurfwFlap,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
# ---------------------------------------------------------------------------------------------

# ---------------------------------------------------------------------------------------------
# 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::TwoDSurf_2DOF,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::TwoDFreeSurf,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) then
        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::TwoDSurfwFlap,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) then
        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 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

function fgfromrhoGauss(rho::Float64)
    f = sqrt(2./pi)*exp(-rho^2/2.)
    g = erf(rho/sqrt(2.)) - rho*f
    return f, g
end

function mutual_ind(vorts::Vector{ThreeDVort})
      for i = 1:length(vorts)
          for j = i+1:length(vorts)
              dx = vorts[i].x[1] - vorts[j].x[1]
              dy = vorts[i].x[2] - vorts[j].x[2]
              dz = vorts[i].x[3] - vorts[j].x[3]
              #source - tar

              r = sqrt(dx*dx + dz*dz + dy*dy)
              rhoi = r/vorts[i].vc
              rhoj = r/vorts[i].vc

              f, g = fgfromrhoGauss(rhoi)

              vorts[j].v += -g*cross([dx; dy; dz],vorts[i].s)/(4.*pi*r^3)

              vorts[j].ds += -(-g*cross(vorts[j].s, vorts[i].s)/rhoi^3 +
                ((3*g/rhoi^3 - f)*dot(vorts[j].s, [dx;dy;dz])*cross([dx;dy;dz],vorts[i].s))/r^2)/(4*pi*vorts[i].vc^3)

            f, g = fgfromrhoGauss(rhoj)
            vorts[i].v += -g*cross([dx; dy; dz],vorts[j].s)/(4.*pi*r^3)

            vorts[i].ds += -(-g*cross(vorts[i].s, vorts[j].s)/rhoj^3 +
            ((3*g/rhoj^3 - f)*dot(vorts[i].s, [dx;dy;dz])*cross([dx;dy;dz],vorts[j].s))/r^2)/(4*pi*vorts[j].vc^3)

          end
      end
      return vorts
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 update_boundpos(surf::TwoDSurf_2DOF, 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 update_boundpos(surf::TwoDFreeSurf, 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 update_boundpos(surf::TwoDSurfwFlap, 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
# ---------------------------------------------------------------------------------------------
#Routines for integral BL calculation (Matsushita et al.)

function FfromE(E::Float64)
      F = 4.8274*E^4 - 5.9816*E^3 + 4.0274*E^2 + 0.23247*E + 0.15174
end

function BfromE(E::Float64)
    if E < -0.0616
        B = -225.86*E^3 - 3016.6*E^2 - 208.68*E - 17.915
    elseif E > -0.0395
        B = 131.9*E^3 - 167.32*E^2 + 76.642*E - 11.068
    else
        B = 0.5*(-225.86*E^3 - 3016.6*E^2 - 208.68*E - 17.915 + 131.9*E^3 - 167.32*E^2 + 76.642*E - 11.068)
    end
end

function SfromE(E::Float64)
    if E < -0.0582
        S = 451.55*E^3 + 2010.*E^2 + 138.96*E + 11.296
    elseif E > -0.042
        S = -96.739*E^3 + 117.74*E^2 - 46.432*E + 6.8074
    else
        S = 0.5*(451.55*E^3 + 2010.*E^2 + 138.96*E + 11.296 - 96.739*E^3 + 117.74*E^2 - 46.432*E + 6.8074)
    end
end

function dfdefromE(E::Float64)
    dfde = 4*4.8274*E^3 - 3*5.9816*E^2 + 2*4.0274*E + 0.23247
end

function calcdt(dx::Float64, cfl::Float64, lamb::Array{Float64,2})
    dt = 10000
    for ic = 2:size(lamb,2)-1
        dti = cfl*(dx/(abs(lamb[1,ic]) + abs(lamb[2,ic])))
        if dti < dt
            dt=dti
        end
    end
    return dt
end

function calcEigen(ue::Vector{Float64}, E::Vector{Float64}, F::Vector{Float64}, dfde::Vector{Float64})
    ncell = length(E) - 2
    lamb = zeros(2,ncell+2)
    for ic = 1:ncell+2
        aq = 1.
        bq = -ue[ic]*(dfde[ic] - 1.)
        cq = ue[ic]*ue[ic]*(E[ic]*dfde[ic] - F[ic])
        lamb[1,ic] = (-bq + sqrt(bq*bq - 4*aq*cq))/(2*aq)
        lamb[2,ic] = (-bq - sqrt(bq*bq - 4*aq*cq))/(2*aq)

        #Always have lamb1 > lamb2
        if lamb[2,ic] > lamb[1,ic]
            temp  = lamb[2,ic]
            lamb[2,ic] = lamb[1,ic]
            lamb[1,ic] = temp
        end
    end
    return lamb
end

function calc_flux(lamb::Array{Float64,2}, ue::Vector{Float64}, E::Vector{Float64}, del::Vector{Float64}, F::Vector{Float64})
    ncell = length(ue) - 2
    flux = zeros(2,2,ncell+2)
    Apos = zeros(2,2)
    Aneg = zeros(2,2)

    for ic = 1:ncell+2
        if lamb[1,ic] >= 0. && lamb[2,ic] >= 0.
            flux[1,1,ic] = ue[ic]*E[ic]*del[ic]
            flux[1,2,ic] = ue[ic]*F[ic]*del[ic]
        elseif lamb[1,ic] < 0. && lamb[2,ic] < 0.
            flux[1,:,ic] = 0.
            else
            Apos[1,1] = (ue[ic]*lamb[1,ic]/(lamb[1,ic] - lamb[2,ic]))*(-1. - lamb[2,ic]/ue[ic])
            Apos[1,2] = ue[ic]*lamb[1,ic]/(lamb[1,ic] - lamb[2,ic])
            Apos[2,1] = -(ue[ic]*lamb[1,ic]/(lamb[1,ic] - lamb[2,ic]))*(1 + lamb[1,ic]/ue[ic])*(1 + lamb[2,ic]/ue[ic])
            Apos[2,2] = (ue[ic]*lamb[1,ic]/(lamb[1,ic] - lamb[2,ic]))*(1 + lamb[1,ic]/ue[ic])
            flux[1,1,ic] =  Apos[1,1]*del[ic] + Apos[1,2]*(E[ic] + 1.)*del[ic]
            flux[1,2,ic] = Apos[2,1]*del[ic] + Apos[2,2]*(E[ic] + 1.)*del[ic]
        end
    end
    for ic = 1:ncell+2
        if lamb[1,ic] >= 0. && lamb[2,ic] >= 0.
            flux[2,:,ic] = 0.
         elseif lamb[1,ic] < 0. && lamb[2,ic] < 0.
            flux[2,1,ic] = ue[ic]*E[ic]*del[ic]
            flux[2,2,ic] = ue[ic]*F[ic]*del[ic]
         else
            Aneg[1,1] = (ue[ic]*lamb[2,ic]/(lamb[1,ic] - lamb[2,ic]))*(1. + lamb[1,ic]/ue[ic])
            Aneg[1,2] = -ue[ic]*lamb[2,ic]/(lamb[1,ic] - lamb[2,ic])
            Aneg[2,1] = (ue[ic]*lamb[2,ic]/(lamb[1,ic] - lamb[2,ic]))*(1. + lamb[1,ic]/ue[ic])*(1. + lamb[2,ic]/ue[ic])
            Aneg[2,2] = (ue[ic]*lamb[2,ic]/(lamb[1,ic] - lamb[2,ic]))*(-1. - lamb[2,ic]/ue[ic])
            flux[2,1,ic] = Aneg[1,1]*del[ic] + Aneg[1,2]*(E[ic] + 1.)*del[ic]
            flux[2,2,ic] = Aneg[2,1]*del[ic] + Aneg[2,2]*(E[ic] + 1.)*del[ic]
         end
    end
    return flux
end

function IBLm_cyl(ntime = 1000, cfl = 0.5, ncell = 64, ttime = 1.4)

    sepflag = Int(0)
    dx = pi/(real(ncell) + 1.)
    t = 0.005

    x = zeros(ncell+2)
    ue = zeros(ncell+2)
    uex = zeros(ncell+2)
    uet = zeros(ncell+2)

    #Set sources
    for i = 1:ncell+2
        x[i] = real(i-1)*dx
        ue[i] = 2*sin(x[i])
        uex[i] = 2*cos(x[i])
    end

    uet[1:ncell+2] = 0

    E = zeros(ncell+2)
    B = zeros(ncell+2)
    del = zeros(ncell+2)
    F = zeros(ncell+2)
    dfde = zeros(ncell+2)
    S = zeros(ncell+2)
    unk = zeros(2,ncell+2)

    lamb = zeros(2,ncell+2)
    unkh = zeros(2,ncell+2)
    rhs = zeros(2,ncell+2)
    flux = zeros(2,2,ncell+2)
    crit = zeros(ncell+2)
    Apos = zeros(2,2)
    Aneg = zeros(2,2)

    #Set initial conditions
    E[:] = 0.4142
    for ic = 1:ncell+2
        B[ic] = BfromE(E[ic])
        F[ic] = FfromE(E[ic])
    end
    del = sqrt(B*t)

    unk[1,:] = del
    unk[2,:] = del.*(E + 1.)

    #Main loop over time steps
    for i = 1:ntime
        #i = 1
        unk[:,1] = 2*unk[:,2] - unk[:,3]
        unk[:,ncell+2] = 2*unk[:,ncell+1] - unk[:,ncell]
        #Calculate derived quantities
        for ic = 1:ncell+2
            del[ic] = unk[1,ic]
            E[ic] = (unk[2,ic]./del[ic]) - 1.
            F[ic] = FfromE(E[ic])
            B[ic] = BfromE(E[ic])
            S[ic] = SfromE(E[ic])
            dfde[ic] = dfdefromE(E[ic])
        end

        #Compute eigenvalues
        lamb = calcEigen(ue, E, F, dfde)

        #Compute timestep
        dt = calcdt(dx, cfl, lamb)

        #Compute fluxes
        flux = calc_flux(lamb, ue, E, del, F)

        #compute rhs
        for ic = 2:ncell+1
        rhs[1,ic] = B[ic]/(2*del[ic]) - del[ic]*uet[ic]/ue[ic] - (E[ic] + 1.)*del[ic]*uex[ic]
            rhs[2,ic] = S[ic]/del[ic] - 2*E[ic]*del[ic]*uet[ic]/ue[ic] - 2*F[ic]*del[ic]*uex[ic]
        end

        for ic = 2:ncell+1
            unkh[1,ic] = unk[1,ic] - (dt/dx)*(flux[1,1,ic] - flux[1,1,ic-1] + flux[2,1,ic+1]
                                              - flux[2,1,ic]) + dt*rhs[1,ic]
            unkh[2,ic] = unk[2,ic] - (dt/dx)*(flux[1,2,ic] - flux[1,2,ic-1] + flux[2,2,ic+1]
                                              - flux[2,2,ic]) + dt*rhs[2,ic]
        end
        ic = ncell+1
        #unkh[1,ic] = unk[1,ic] - (dt/dx)*(flux[1,1,ic] - flux[1,1,ic-1]) + dt*rhs[1,ic]
        #unkh[2,ic] = unk[2,ic] - (dt/dx)*(flux[1,2,ic] - flux[1,2,ic-1]) + dt*rhs[2,ic]

        #Update ghost cells
        unkh[:,1] = 2*unkh[:,2] - unkh[:,3]
        unkh[:,ncell+2] = 2*unkh[:,ncell+1] - unkh[:,ncell]

        #Calculate derived quantities
        for ic = 1:ncell+2
            del[ic] = unkh[1,ic]
            E[ic] = (unkh[2,ic]./del[ic]) - 1.
            F[ic] = FfromE(E[ic])
            B[ic] = BfromE(E[ic])
            S[ic] = SfromE(E[ic])
            dfde[ic] = dfdefromE(E[ic])
        end

        #Compute eigenvalues
        lamb = calcEigen(ue, E, F, dfde)

        #Compute fluxes
        fluxhalf = calc_flux(lamb, ue, E, del, F)

        #compute rhs
        for ic = 2:ncell+1
            rhs[1,ic] = B[ic]/(2*del[ic]) - del[ic]*uet[ic]/ue[ic] - (E[ic] + 1.)*del[ic]*uex[ic]
            rhs[2,ic] = S[ic]/del[ic] - 2*E[ic]*del[ic]*uet[ic]/ue[ic] - 2*F[ic]*del[ic]*uex[ic]
        end

        ic = 2

        unk[1,ic] = 0.5*(unk[1,ic] + unkh[1,ic]) -
        (0.5*dt/dx)*(flux[1,1,ic] - flux[1,1,ic-1] - flux[2,1,ic] +
        2*flux[2,1,ic+1] - flux[2,1,ic+2] + fluxhalf[1,1,ic] -
        fluxhalf[1,1,ic-1] + fluxhalf[2,1,ic+1] - fluxhalf[2,1,ic]) +
        0.5*dt*rhs[1,ic]

        unk[2,ic] = 0.5*(unk[2,ic] + unkh[2,ic]) -
        (0.5*dt/dx)*(flux[1,2,ic] - flux[1,2,ic-1] - flux[2,2,ic] +
        2*flux[2,2,ic+1] - flux[2,2,ic+2] + fluxhalf[1,2,ic] -
        fluxhalf[1,2,ic-1] + fluxhalf[2,2,ic+1] - fluxhalf[2,2,ic]) +
        0.5*dt*rhs[2,ic]

        for ic = 3:ncell

            unk[1,ic] = 0.5*(unk[1,ic] + unkh[1,ic]) -
            (0.5*dt/dx)*(flux[1,1,ic] - 2*flux[1,1,ic-1] +
            flux[1,1,ic-2] - flux[2,1,ic] + 2*flux[2,1,ic+1] -
            flux[2,1,ic+2] + fluxhalf[1,1,ic] -fluxhalf[1,1,ic-1] +
            fluxhalf[2,1,ic+1] - fluxhalf[2,1,ic]) + 0.5*dt*rhs[1,ic]

            unk[2,ic] = 0.5*(unk[2,ic] + unkh[2,ic]) -
            (0.5*dt/dx)*(flux[1,2,ic] - 2*flux[1,2,ic-1] +
            flux[1,2,ic-2] - flux[2,2,ic] + 2*flux[2,2,ic+1] -
            flux[2,2,ic+2] + fluxhalf[1,2,ic] - fluxhalf[1,2,ic-1] +
            fluxhalf[2,2,ic+1] - fluxhalf[2,2,ic]) + 0.5*dt*rhs[2,ic]

        end
        ic = ncell+1

        unk[1,ic] = 0.5*(unk[1,ic] + unkh[1,ic]) -
        (0.5*dt/dx)*(flux[1,1,ic] - 2*flux[1,1,ic-1] + flux[1,1,ic-2]
        +fluxhalf[1,1,ic] - fluxhalf[1,1,ic-1]) + 0.5*dt*rhs[1,ic]

        unk[2,ic] = 0.5*(unk[2,ic] + unkh[2,ic]) -
        (0.5*dt/dx)*(flux[1,2,ic] - 2*flux[1,2,ic-1] + flux[1,2,ic-2]
        +fluxhalf[1,2,ic] - fluxhalf[1,2,ic-1]) + 0.5*dt*rhs[2,ic]

        t = t + dt
        if t > ttime
            break
        end

        for ic = 2:ncell+1
            crit[ic] = abs((del[ic+1] - del[ic])/(del[ic] - del[ic-1]))
            if abs(crit[ic]) > 10.
                sepflag = 1
                println(t, " ",x[ic]," ", ue[ic]," ", crit[ic])
                break
            end
        end
        if t > ttime || sepflag == 1
            break
        end
    end
t, x, del, E
end


text:
#####
##### High level training procedure
#####

"""
    Env{Game, Network, Board}

Type for an AlphZero environment.

The environment features the current neural network, the best neural network
seen so far that is used for data generation, a memory buffer
and an iteration counter.

# Constructor

    Env{Game}(params, curnn, bestnn=copy(curnn), experience=[], itc=0)

Construct a new AlphaZero environment.
- `Game` is the type of the game being played
- `params` has type [`Params`](@ref)
- `curnn` is the current neural network and has type [`AbstractNetwork`](@ref)
- `bestnn` is the best neural network so far, which is used for data generation
- `experience` is the initial content of the memory buffer
   as a vector of [`TrainingSample`](@ref)
- `itc` is the value of the iteration counter (0 at the start of training)
"""
mutable struct Env{Game, Network, Board}
  params :: Params
  curnn  :: Network
  bestnn :: Network
  memory :: MemoryBuffer{Board}
  itc    :: Int
  function Env{Game}(
      params, curnn, bestnn=copy(curnn), experience=[], itc=0) where Game
    Board = GI.State(Game)
    msize = max(params.mem_buffer_size[itc], length(experience))
    memory = MemoryBuffer{Board}(msize, experience)
    return new{Game, typeof(curnn), Board}(
      params, curnn, bestnn, memory, itc)
  end
end

GameType(env::Env{Game}) where Game = Game

#####
##### Training handlers
#####

"""
    Handlers

Namespace for the callback functions that are used during training.
This enables logging, saving and plotting to be implemented separately.
An example handler object is `Session`.

All callback functions take a handler object `h` as their first argument
and sometimes a second argment `r` that consists in a report.

| Callback                    | Comment                                        |
|:----------------------------|:-----------------------------------------------|
| `iteration_started(h)`      | called at the beggining of an iteration        |
| `self_play_started(h)`      | called once per iter before self play starts   |
| `game_played(h)`            | called after each game of self play            |
| `self_play_finished(h, r)`  | sends report: [`Report.SelfPlay`](@ref)        |
| `memory_analyzed(h, r)`     | sends report: [`Report.Memory`](@ref)          |
| `learning_started(h, r)`    | sends report: [`Report.LearningStatus`](@ref)  |
| `updates_started(h)`        | called before each series of batch updates     |
| `updates_finished(h, r)`    | sends report: [`Report.LearningStatus`](@ref)  |
| `checkpoint_started(h)`     | called before a checkpoint evaluation starts   |
| `checkpoint_game_played(h)` | called after each arena game                   |
| `checkpoint_finished(h, r)` | sends report: [`Report.Checkpoint`](@ref)      |
| `learning_finished(h, r)`   | sends report: [`Report.Learning`](@ref)        |
| `iteration_finished(h, r)`  | sends report: [`Report.Iteration`](@ref)       |
| `training_finished(h)`      | called once at the end of training             |
"""
module Handlers

  import ..Report

  function iteration_started(h)      return end
  function self_play_started(h)      return end
  function game_played(h)            return end
  function self_play_finished(h, r)  return end
  function memory_analyzed(h, r)     return end
  function learning_started(h, r)    return end
  function updates_started(h)        return end
  function updates_finished(h, r)    return end
  function checkpoint_started(h)     return end
  function checkpoint_game_played(h) return end
  function checkpoint_finished(h, r) return end
  function learning_finished(h, r)   return end
  function iteration_finished(h, r)  return end
  function training_finished(h)      return end

end

import .Handlers

#####
##### Public utilities
#####

"""
    get_experience(env::Env)

Return the content of the agent's memory as a
vector of [`TrainingSample`](@ref).
"""
get_experience(env::Env) = get_experience(env.memory)

"""
    initial_report(env::Env)

Return a report summarizing the configuration of agent before training starts,
as an object of type [`Report.Initial`](@ref).
"""
function initial_report(env::Env)
  num_network_parameters = Network.num_parameters(env.curnn)
  num_reg_params = Network.num_regularized_parameters(env.curnn)
  player = MctsPlayer(env.curnn, env.params.self_play.mcts)
  mcts_footprint_per_node = MCTS.memory_footprint_per_node(player.mcts)
  return Report.Initial(
    num_network_parameters, num_reg_params, mcts_footprint_per_node)
end

#####
##### Training loop
#####

function resize_memory!(env::Env{G,N,B}, n) where {G,N,B}
  exp = get_experience(env.memory)
  env.memory = MemoryBuffer{B}(n, exp)
  return
end

function evaluate_network(contender, baseline, params, handler)
  use_gpu = params.arena.use_gpu
  contender = Network.copy(contender, on_gpu=use_gpu, test_mode=true)
  baseline = Network.copy(baseline, on_gpu=use_gpu, test_mode=true)
  simulator = Simulator((contender, baseline), record_trace) do oracles
    white = MctsPlayer(oracles[1], params.arena.mcts)
    black = MctsPlayer(oracles[2], params.arena.mcts)
    return TwoPlayers(white, black)
  end
  samples = simulate(
    simulator,
    num_games=params.arena.num_games,
    num_workers=params.arena.num_workers,
    handler=(trace -> Handlers.checkpoint_game_played(handler)),
    reset_every=params.arena.reset_mcts_every,
    flip_probability=params.arena.flip_probability,
    color_policy=ALTERNATE_COLORS)
  gamma = params.self_play.mcts.gamma
  rewards, redundancy = rewards_and_redundancy(samples, gamma=gamma)
  return mean(rewards), redundancy
end

function learning_step!(env::Env, handler)
  ap = env.params.arena
  lp = env.params.learning
  checkpoints = Report.Checkpoint[]
  losses = Float32[]
  tloss, teval, ttrain = 0., 0., 0.
  experience = get_experience(env.memory)
  if env.params.use_symmetries
    experience = augment_with_symmetries(GameType(env), experience)
  end
  trainer, tconvert = @timed Trainer(env.curnn, experience, lp)
  init_status = learning_status(trainer)
  Handlers.learning_started(handler, init_status)
  # Compute the number of batches between each checkpoint
  nbatches = lp.max_batches_per_checkpoint
  if !iszero(lp.min_checkpoints_per_epoch)
    ntotal = num_batches_total(trainer)
    nbatches = min(nbatches, ntotal ÷ lp.min_checkpoints_per_epoch)
  end
  # Loop state variables
  best_evalz = ap.update_threshold
  nn_replaced = false

  for k in 1:lp.num_checkpoints
    # Execute a series of batch updates
    Handlers.updates_started(handler)
    dlosses, dttrain = @timed batch_updates!(trainer, nbatches)
    status, dtloss = @timed learning_status(trainer)
    Handlers.updates_finished(handler, status)
    tloss += dtloss
    ttrain += dttrain
    append!(losses, dlosses)
    # Run a checkpoint evaluation
    Handlers.checkpoint_started(handler)
    env.curnn = get_trained_network(trainer)
    (evalz, redundancy), dteval =
      @timed evaluate_network(env.curnn, env.bestnn, env.params, handler)
    teval += dteval
    # If eval is good enough, replace network
    success = (evalz >= best_evalz)
    if success
      nn_replaced = true
      env.bestnn = copy(env.curnn)
      best_evalz = evalz
    end
    checkpoint_report = Report.Checkpoint(
      k * nbatches, status, evalz, redundancy, success)
    push!(checkpoints, checkpoint_report)
    Handlers.checkpoint_finished(handler, checkpoint_report)
  end

  report = Report.Learning(
    tconvert, tloss, ttrain, teval,
    init_status, losses, checkpoints, nn_replaced)
  Handlers.learning_finished(handler, report)
  return report
end

function simple_memory_stats(env)
  mem = get_experience(env)
  nsamples = length(mem)
  ndistinct = length(merge_by_state(mem))
  return nsamples, ndistinct
end

# To be given as an argument to `Simulator`
function self_play_measurements(trace, _, player)
  mem = MCTS.approximate_memory_footprint(player.mcts)
  edepth = MCTS.average_exploration_depth(player.mcts)
  return (trace=trace, mem=mem, edepth=edepth)
end

function self_play_step!(env::Env{G}, handler) where G
  params = env.params.self_play
  Handlers.self_play_started(handler)
  network = Network.copy(env.bestnn, on_gpu=params.use_gpu, test_mode=true)
  simulator = Simulator(network, self_play_measurements) do oracle
    return MctsPlayer(oracle, params.mcts)
  end
  # Run the simulations
  results, elapsed = @timed simulate(
    simulator,
    num_games=params.num_games,
    num_workers=params.num_workers,
    handler=_->Handlers.game_played(handler),
    reset_every=params.reset_mcts_every,
    fill_batches=true,
    flip_probability=0.,
    color_policy=nothing)
  # Add the collected samples in memory
  new_batch!(env.memory)
  for x in results
    push_game!(env.memory, x.trace, params.mcts.gamma)
  end
  speed = cur_batch_size(env.memory) / elapsed
  edepth = mean([x.edepth for x in results])
  mem_footprint = maximum([x.mem for x in results])
  memsize, memdistinct = simple_memory_stats(env)
  report = Report.SelfPlay(
    speed, edepth, mem_footprint, memsize, memdistinct)
  Handlers.self_play_finished(handler, report)
  return report
end

function memory_report(env::Env, handler)
  if isnothing(env.params.memory_analysis)
    return nothing
  else
    report = memory_report(
      env.memory, env.curnn, env.params.learning, env.params.memory_analysis)
    Handlers.memory_analyzed(handler, report)
    return report
  end
end

"""
    train!(env::Env, handler=nothing)

Start or resume the training of an AlphaZero agent.

A `handler` object can be passed that implements a subset of the callback
functions defined in [`Handlers`](@ref).
"""
function train!(env::Env, handler=nothing)
  while env.itc < env.params.num_iters
    Handlers.iteration_started(handler)
    resize_memory!(env, env.params.mem_buffer_size[env.itc])
    sprep, spperfs = Report.@timed self_play_step!(env, handler)
    mrep, mperfs = Report.@timed memory_report(env, handler)
    lrep, lperfs = Report.@timed learning_step!(env, handler)
    rep = Report.Iteration(spperfs, mperfs, lperfs, sprep, mrep, lrep)
    env.itc += 1
    Handlers.iteration_finished(handler, rep)
  end
  Handlers.training_finished(handler)
end


text:
"""
Script for generating local documentation in docs/build.

To run, enter `;` in the Julia REPL and `cd` into `docs`.
Then run
```shell
shell> julia --project=@. makeLocal.jl
```
"""

using RandomStrainDistributions
using Documenter

DocMeta.setdocmeta!(RandomStrainDistributions, :DocTestSetup, :(using RandomStrainDistributions); recursive=true)

makedocs(;
    modules=[RandomStrainDistributions],
    authors="W. Joe Meese <meese022@umn.edu> and contributors",
    repo="https://github.com/meese-wj/RandomStrainDistributions.jl/blob/{commit}{path}#{line}",
    sitename="RandomStrainDistributions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
abstract type Thermostat
end

struct NullThermostat <: Thermostat
end

struct AndersenThermostat{tType <: Real, νType <: Real} <: Thermostat
    T::tType
    ν::νType
end

struct BerendsenThermostat{τType <: Real,tType <: Real} <: Thermostat
    T::tType
    τ::τType
    γ::τType
end

function BerendsenThermostat(T::Real, τ::Real)
    BerendsenThermostat(T, τ, 0.5 / τ)
end

function berendsen_acceleration!(dv, v, ms, kb, N, Nc, p::BerendsenThermostat)
    T = md_temperature(v, ms, kb, N, Nc)
    if inv(T) == Inf
        @. dv += p.γ * v
    else
        @. dv += p.γ * (p.T / T - 1) * v
    end
end

# N - number of particles
# Nc - number of constraints
function md_temperature(vs, ms, kb, N, Nc)
    e_kin = sum(dot(ms, vec(sum(vs.^2, dims=1))))
    temperature = e_kin / (kb * (3 * N - Nc))
    return temperature
end

struct NoseHooverThermostat{tType <: Real,τType <: Real} <: Thermostat
    T::tType
    τ::τType
end

function nosehoover_acceleration!(dv, u, v, ms, kb, N, Nc, ζind, p::NoseHooverThermostat)
    @. dv -= u[ζind] * v
    @. dv[:,end] = 0
    T = md_temperature(v[:, 1:N], ms, kb, N, Nc)
    ndf = 3 * N - Nc
    v[ζind] = inv(p.τ)^2 * (T/p.T - (ndf+1)/ndf)
    #v[ζind] = inv(p.Q) * ( sum(dot(ms, vec(sum(v[:,1:N].^2, 1)))) - (3 * N - Nc) * kb * p.T)
end

struct LangevinThermostat{tType <: Real, gType <: Real} <: Thermostat
    T::tType
    γ::gType
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "FastHadamardStructuredTransforms"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://gitlab.inria.fr/SketchedLearning/fasttransformsc.git", "5e7c023fbb89dc3c4d1bcaa41a97fc4a7bd37078"),
    GitSource("https://github.com/FALCONN-LIB/FFHT.git", "fba727a3ff72c862260460631ae1cc8e3e44e861")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/FFHT/
CFLAGS="-O3 -std=c99 -pedantic -Wall -Wshadow -Wpointer-arith -Wcast-qual -Wstrict-prototypes -Wmissing-prototypes"
if [[ ${target} == x86_64-* ]] || [[ ${target} == i686-* ]]; then
    CFLAGS="$CFLAGS -msse -msse2 -mno-avx"
fi

$CC fht.c -o fht.o -c $CFLAGS -fPIC
$CC fast_copy.c -o fast_copy.o -c $CFLAGS -fPIC
cd $WORKSPACE/srcdir/fasttransformsc
$CC fast_transforms.c -o fast_transforms.o -I$WORKSPACE/srcdir -c $CFLAGS -fPIC -lm
mkdir -p "${libdir}"
$CC -shared fast_transforms.o ../FFHT/fht.o ../FFHT/fast_copy.o -o "${libdir}/libfasttransforms.${dlext}" $CFLAGS -lm
install_license $WORKSPACE/srcdir/fasttransformsc/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "freebsd"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libfasttransforms", :libfasttransforms)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module PowerSystemsUnits

import Unitful
import Unitful: J, W, hr, 𝐋, 𝐌, 𝐓
using Unitful: @unit, @derived_dimension, @dimension, @refunit, @u_str, uconvert, Quantity

export asqtype, fustrip, UnitfulMissing

# Power Units
@derived_dimension PowerHour 𝐋^2*𝐌*𝐓^-2
@unit Wh "Wh" WattHour 3600J true

@derived_dimension ReactivePowerHour 𝐋^2*𝐌*𝐓^-2
@unit VARh "VARh" VARHour 3600J true

# Monetary Units
@dimension Money "Money" Currency
@refunit USD "USD" Currency Money false

# Monetary and Power Units
@derived_dimension MoneyPerPowerHour Money*𝐋^-2*𝐌^-1*𝐓^2
@unit USDPerMWh "USDPerMWh" DollarPerMegaWattHour USD/(1000000*Wh) false

include("utils.jl")

# Some gymnastics needed to get this to work at run-time.
# Sourced from https://github.com/ajkeller34/UnitfulUS.jl
const localunits = Unitful.basefactors
const localpromotion = Unitful.promotion
function __init__()
    merge!(Unitful.basefactors, localunits)
    merge!(Unitful.promotion, localpromotion)  # only if you've used @dimension
    Unitful.register(PowerSystemsUnits)
end

end

max_stars_repo_name: invenia/PowerSystemsUnits

max_issues_repo_name: invenia/PowerSystemsUnits

max_forks_repo_name: invenia/PowerSystemsUnits


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

function generate_deprecated(varargs...; kwargs...)
    Base.depwarn("Pkg.generate is deprecated. Please use PkgTemplates.jl instead.", Symbol("Pkg.generate"))
    return generate(varargs...; kwargs...)
end

generate(path::String; kwargs...) = generate(Context(), path; kwargs...)
function generate(ctx::Context, path::String; kwargs...)
    Context!(ctx; kwargs...)
    dir, pkg = dirname(path), basename(path)
    Base.isidentifier(pkg) || pkgerror("$(repr(pkg)) is not a valid package name")
    isdir(path) && pkgerror("$(abspath(path)) already exists")
    printpkgstyle(ctx, :Generating, " project $pkg:")
    uuid = project(ctx, pkg, dir)
    entrypoint(ctx, pkg, dir)
    return Dict(pkg => uuid)
end

function genfile(f::Function, ctx::Context, pkg::String, dir::String, file::String)
    path = joinpath(dir, pkg, file)
    println(ctx.io, "    $(Base.contractuser(path))")
    mkpath(dirname(path))
    open(f, path, "w")
    return
end

function project(ctx::Context, pkg::String, dir::String)
    name = email = nothing
    gitname = LibGit2.getconfig("user.name", "")
    isempty(gitname) || (name = gitname)
    gitmail = LibGit2.getconfig("user.email", "")
    isempty(gitmail) || (email = gitmail)

    if name === nothing
        for env in ["GIT_AUTHOR_NAME", "GIT_COMMITTER_NAME", "USER", "USERNAME", "NAME"]
            name = get(ENV, env, nothing)
            name !== nothing && break
        end
    end

    name === nothing && (name = "Unknown")

    if email === nothing
        for env in ["GIT_AUTHOR_EMAIL", "GIT_COMMITTER_EMAIL", "EMAIL"];
            email = get(ENV, env, nothing)
            email !== nothing && break
        end
    end

    authors = ["$name " * (email === nothing ? "" : "<$email>")]

    uuid = UUIDs.uuid4()
    genfile(ctx, pkg, dir, "Project.toml") do io
        toml = Dict{String,Any}("authors" => authors,
                                "name" => pkg,
                                "uuid" => string(uuid),
                                "version" => "0.1.0",
                                )
        TOML.print(io, toml, sorted=true, by=key -> (Types.project_key_order(key), key))
    end
    return uuid
end

function entrypoint(ctx::Context, pkg::String, dir)
    genfile(ctx, pkg, dir, "src/$pkg.jl") do io
        print(io,
           """
            module $pkg

            greet() = print("Hello World!")

            end # module
            """
        )
    end
end


text:
module tank_examples
using FinEtools
using FinEtoolsAcoustics
using FinEtoolsAcoustics.AlgoAcoustModule

function tank_piston_platten()
    println("""
    Axially symmetric geometry  of the tank plus piston plus platten
    Version: 08/26/2017
    """)

    t0 = time()

    rho = 1001*phun("kg/m^3");# mass density
    c  = 1500.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  200*phun("rev/s");      # frequency of the piston
    a_piston = 1.0*phun("m/s")

    axisymmetric = true
    Platten_radius = 12*phun("mm")
    Piston_radius = 33*phun("mm")
    @assert Piston_radius > Platten_radius
    Tank_radius = 50*phun("mm")
    @assert Tank_radius > Piston_radius
    Platten_height = 100*phun("mm")
    Piston_height = 33*phun("mm")
    Tank_height = Platten_height + Piston_height + 75*phun("mm")
    Gap_height = Tank_height-(Platten_height + Piston_height)
    Tolerance = Platten_radius/1.0e6

    xs = vec([0.0 Platten_radius Piston_radius Tank_radius])
    ys = vec(vcat(vec(linearspace(0.0, Platten_height, 4)),
    vec(linearspace(Platten_height+Gap_height/4, Platten_height+Gap_height, 3)),
    vec(linearspace(Platten_height+Gap_height+Piston_height/4, Tank_height, 3))))

    fens, fes = Q4blockx(xs, ys)

    l1 = selectelem(fens, fes, box=[0.0 Piston_radius-10*Tolerance Tank_height-Piston_height+10*Tolerance Tank_height],
    inflate=Tolerance, allin=false);
    l2 = selectelem(fens, fes, box=[0.0 Platten_radius-10*Tolerance 0.0 Platten_height-10*Tolerance],
    inflate=Tolerance, allin=false);
    fes = subset(fes, setdiff(1:count(fes), vcat(l1, l2)))
    for ixxxx in 1:3
        fens, fes = Q4refine(fens, fes)
    end
    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    # Postprocessing
    # File = "tank_mesh.vtk"
    # vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
    # @async run(`"paraview.exe" $File`)

    edge_fes = meshboundary(fes);

    material = MatAcoustFluid(bulk, rho)
    # The pressure boundary condition
    l1 = selectelem(fens, edge_fes, box=[0.0 Piston_radius Tank_height-Piston_height Tank_height-Piston_height]);
    flux1  =  FDataDict("femm"=>FEMMAcoustSurf(IntegDomain(subset(edge_fes, l1), GaussRule(1, 2), axisymmetric),
    material),  "normal_flux"=> -rho*a_piston+0.0im);
    l2 = selectelem(fens, edge_fes, box=[Piston_radius Tank_radius Tank_height Tank_height]);
    ebc2 = FDataDict("node_list"=>connectednodes(subset(edge_fes, l2)),
    "pressure"=>x -> 0.0) # entering the domain

    femm = FEMMAcoust(IntegDomain(fes,  GaussRule(2, 2), axisymmetric),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,"omega"=>omega,
    "regions"=>[region1], "flux_bcs"=>[flux1], "essential_bcs"=>[ebc2]);

    # Call the solver
    modeldata = AlgoAcoustModule.steadystate(modeldata)
    geom = modeldata["geom"]
    P = modeldata["P"]
    println("Minimum/maximum pressure, real= $(minimum(real(P.values)))/$(maximum(real(P.values))))")
    println("Minimum/maximum pressure, imag= $(minimum(imag(P.values)))/$(maximum(imag(P.values))))")

    println("Total time elapsed = ",time() - t0,"s")

    # Postprocessing
    File = "tank_piston_platten.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.Q4; scalars=[("Pre", real(P.values)), ("Pim", imag(P.values))])
    @async run(`"paraview.exe" $File`)

end # tank_piston_platten

function tank_piston_platten_pressure()
   println("""
    Axially symmetric geometry  of the tank plus piston plus platten
    Version: 08/26/2017
    """)

    t0 = time()

    rho = 1001*phun("kg/m^3");# mass density
    c  = 1500.0*phun("m/s");# sound speed
    bulk =  c^2*rho;
    omega =  200*phun("rev/s");      # frequency of the piston

    axisymmetric = true
    Platten_radius = 15*phun("mm")
    Piston_radius = 35*phun("mm")
    @assert Piston_radius > Platten_radius
    Tank_radius = 50*phun("mm")
    @assert Tank_radius > Piston_radius
    Platten_height = 75*phun("mm")
    Piston_height = 75*phun("mm")
    Tank_height = Platten_height + Piston_height + 75*phun("mm")
    Gap_height = Tank_height-(Platten_height + Piston_height)
    Tolerance = Platten_radius/1.0e6

    xs = vec([0.0 Platten_radius Piston_radius Tank_radius])
    ys = vec(vcat(vec(linearspace(0.0, Platten_height, 4)),
    vec(linearspace(Platten_height+Gap_height/4, Platten_height+Gap_height, 3)),
    vec(linearspace(Platten_height+Gap_height+Piston_height/4, Tank_height, 3))))

    fens, fes = Q4blockx(xs, ys)

    l1 = selectelem(fens, fes, box=[0.0 Piston_radius-10*Tolerance Tank_height-Piston_height+10*Tolerance Tank_height],
    inflate=Tolerance, allin=false);
    l2 = selectelem(fens, fes, box=[0.0 Platten_radius-10*Tolerance 0.0 Platten_height-10*Tolerance],
    inflate=Tolerance, allin=false);
    fes = subset(fes, setdiff(1:count(fes), vcat(l1, l2)))
    for ixxxx in 1:3
        fens, fes = Q4refine(fens, fes)
    end
    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    # Postprocessing
    # File = "tank_mesh.vtk"
    # vtkexportmesh(File, fes.conn, fens.xyz, FinEtools.MeshExportModule.Q4)
    # @async run(`"paraview.exe" $File`)

    edge_fes = meshboundary(fes);

    # The pressure boundary condition
    l1 = selectelem(fens, edge_fes, box=[0.0 Piston_radius Tank_height-Piston_height Tank_height-Piston_height]);
    ebc1 = FDataDict("node_list"=>connectednodes(subset(edge_fes, l1)),
    "pressure"=>x -> 1.0) # entering the domain
    l2 = selectelem(fens, edge_fes, box=[Piston_radius Tank_radius Tank_height Tank_height]);
    ebc2 = FDataDict("node_list"=>connectednodes(subset(edge_fes, l2)),
    "pressure"=>x -> 0.0) # entering the domain
    material = MatAcoustFluid(bulk, rho)
    femm = FEMMAcoust(IntegDomain(fes,  GaussRule(2, 2), axisymmetric),  material)
    region1 = FDataDict("femm"=>femm)

    # Make model data
    modeldata = FDataDict("fens"=>fens,"omega"=>omega,
    "regions"=>[region1], "essential_bcs"=>[ebc1 ebc2]);

    # Call the solver
    modeldata = AlgoAcoustModule.steadystate(modeldata)
    geom = modeldata["geom"]
    P = modeldata["P"]
    println("Minimum/maximum pressure, real= $(minimum(real(P.values)))/$(maximum(real(P.values)))")
    println("Minimum/maximum pressure, imag= $(minimum(imag(P.values)))/$(maximum(imag(P.values)))")

    println("Total time elapsed = ",time() - t0,"s")

    # Postprocessing
    File = "tank_piston_platten_pressure.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.Q4; scalars=[("Pre", real(P.values)), ("Pim", imag(P.values))])
    @async run(`"paraview.exe" $File`)

end # tank_piston_platten_pressure

function allrun()
    println("#####################################################")
    println("# tank_piston_platten ")
    tank_piston_platten()
    println("#####################################################")
    println("# tank_piston_platten_pressure ")
    tank_piston_platten_pressure()
    return true
end # function allrun

end # module tank_examples


text:
push!(LOAD_PATH, ".")
using Statistics
using Points

function mean_var(ps :: Array{P}) where {P<:Point}
    m = mean(ps)
    ps2 = map(p->abs2(p-m),ps)
    v = sum(ps2)/(length(ps2)-1)
    return abs(m),v
end

function randomwalk(::Type{P}, num::Integer, nsteps::Integer) where {P<:Point}
    means, vars = zeros(nsteps+1), zeros(nsteps+1)
    ps = zeros(P,num)
    means[1], vars[1] = mean_var(ps)
    for t in 0:nsteps
        map!(update, ps, ps)
        means[t+1], vars[t+1] = mean_var(ps)
    end
    means, vars
end

function main(dim, num, nsteps, filename)
    if dim == 1
        P = Point1D
    elseif dim == 2
        P = Point2D
    else
        error("dimension error!")
    end

    if isempty(filename)
        io = stdout
    else
        io = open(filename, "w")
    end

    means, vars = randomwalk(P, num, nsteps)

    for (i, (m, v)) in enumerate(zip(means, vars))
        println(io, i-1, " ", m, " ", v)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    dim = length(ARGS)>0 ? parse(Int, ARGS[1]) : 1
    num = length(ARGS)>1 ? parse(Int, ARGS[2]) : 1000
    nsteps = length(ARGS)>2 ? parse(Int, ARGS[3]) : 1000
    filename = length(ARGS)>3 ? ARGS[4] : ""

    main(dim, num, nsteps, filename)
end

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
"""
Abstract interactor. The top node of the Interactor type hierarchy.
"""
abstract type AbstractInteractor end

"""
Abstract speaker: an interactor that is an individual.
"""
abstract type AbstractSpeaker <: AbstractInteractor end

"""
Abstract language: an interactor that is a language (community).
"""
abstract type AbstractLanguage <: AbstractInteractor end

"""
Abstract probabilistic speaker: a speaker containing a probabilistic replicator.
"""
abstract type AbstractProbabilisticSpeaker <: AbstractSpeaker end

"""
Abstract utterance selector: an interactor of the Utterance Selection Model type.
"""
abstract type AbstractUtteranceSelector <: AbstractProbabilisticSpeaker end