"Stores a set of coupling constants and meson masses that goes into the Lagrangian"
struct parameters
    # Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
    m_s::Float64 # MeV/c2
    m_ω::Float64 # MeV/c2
    m_ρ::Float64 # MeV/c2
    m_γ::Float64 # MeV/c2
    g2_s::Float64 # dimensionless
    g2_v::Float64 # dimensionless
    g2_ρ::Float64 # dimensionless
    g2_γ::Float64 # dimensionless
    κ_ss::Float64 # MeV
    λ::Float64 # dimensionless, aka LambdaSS
    ζ::Float64 # dimensionless, aka LambdaVV
    Λv::Float64 # dimensionless, aka LambdaVR

    "Dummy struct when parameters are not needed (for testing)"
    parameters() = new(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

    "Initialize parameters from a string with values provided in order of struct definition separated by commas"
    function parameters(s::String)
        values = parse.(Float64, strip.(split(s, ',')))
        @assert length(values) == 12 "String must contain exactly 12 values separated by commas"
        return new(values...)
    end
end

"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

    param::parameters
    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, parameters, r_max, divs) = new(Z, N, parameters, r_max, divs, unfilled_spectrum(), unfilled_spectrum(), [zeros(1 + divs) for _ in 1:8]...)

    "Dummy struct to define the mesh (for testing)"
    system(r_max, divs) = system(0, 0, parameters(), 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)