using FastGaussQuadrature

function gausslegendre_shifted(a, b, n)
    scale = (b - a) / 2
    shift = (a + b) / 2
    p, w = gausslegendre(n)
    p = p .* scale .+ shift
    w = w .* scale
    return (p, w)
end

function get_mesh(vertices, subdivisions)
    p = Vector{ComplexF64}()
    w = Vector{ComplexF64}()
    for (a, b) in zip(vertices, vertices[2:end])
        p_new, w_new = gausslegendre_shifted(a, b, subdivisions)
        append!(p, p_new)
        append!(w, w_new)
    end
    return (p, w)
end

get_V_matrix(V_pq, p, w) = V_pq.(p, transpose(p)) .* transpose(w)

get_K_matrix(p, μ) = collect(Diagonal(p.*p ./ (2*μ)))

get_H_matrix(V_pq, p, w, μ=0.5) = get_K_matrix(p, μ) + get_V_matrix(V_pq, p, w)