"Tabulates a nucleon spectrum (protons or neutrons) containing κ and occupancy"
struct spectrum
    κ::Vector{Int}
    E::Vector{Float64}
    occ::Vector{Int}
end

"Initializes an unfilled spectrum"
unfilled_spectrum() = spectrum(Int[], Float64[], Int[])

"Defines a nuclear system containing relevant parameters and meson/nucleon densities"
mutable struct system
    Z::Int
    N::Int

    r_max::Float64
    divs::Int

    p_spectrum::spectrum
    n_spectrum::spectrum

    Φ0::Vector{Float64}
    W0::Vector{Float64}
    B0::Vector{Float64}
    A0::Vector{Float64}

    ρ_sp::Vector{Float64}
    ρ_vp::Vector{Float64}
    ρ_sn::Vector{Float64}
    ρ_vn::Vector{Float64}

    "Initialize an unsolved system"
    system(Z, N, r_max, divs) = new(Z, N, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)

    "Dummy struct to define the mesh"
    system(r_max, divs) = system(0, 0, r_max, divs)
end

"Get mass number of nucleus"
A(s::system)::Int = s.Z + s.N

"Get r values in the mesh"
rs(s::system)::StepRangeLen = range(0, s.r_max, length=s.divs+1)

"Get Δr value for the mesh"
Δr(s::system)::Float64 = s.r_max / s.divs

"Get the number of protons or neutrons in the system"
Z_or_N(s::system, p::Bool)::Int = p ? s.Z : s.N

"Create an empty array for the size of the mesh"
zero_array(s::system) = zeros(1 + s.divs)