diff --git a/mesons.jl b/mesons.jl
index 96f90de..b3888e0 100644
--- a/mesons.jl
+++ b/mesons.jl
@@ -1,21 +1,6 @@
 include("common.jl")
 include("system.jl")
 
-# Values defined in C. J. Horowitz and J. Piekarewicz, Phys. Rev. Lett. 86, 5647 (2001)
-# Values taken from Hartree.f (FSUGarnet)
-const m_s = 496.939473213388 # MeV/c2
-const m_ω = 782.5 # MeV/c2
-const m_ρ = 763.0 # MeV/c2
-const m_γ = 0.000001000 # MeV/c2 # should be 0?
-const g2_s = 110.349189097820 # dimensionless
-const g2_v = 187.694676506801 # dimensionless
-const g2_ρ = 192.927428365698 # dimensionless
-const g2_γ = 0.091701236 # dimensionless # equal to 4πα
-const κ_ss = 3.260178893447 # MeV
-const λ = -0.003551486718 # dimensionless # LambdaSS
-const ζ = 0.023499504053 # dimensionless # LambdaVV
-const Λv = 0.043376933644 # dimensionless # LambdaVR
-
 "Green's function for Klein-Gordon equation in natural units"
 greensFunctionKG(m, r, rp) = -1 / (m * (r + r_reg) * (rp + r_reg)) * sinh(m * min(r, rp)) * exp(-m * max(r, rp))
 
@@ -42,6 +27,7 @@ end
     the inital solutions are read from s and the final solutions are saved in-place.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
 function solveMesonFields!(s::system, iterations=15, oscillation_control_parameter=0.3)
+    (m_s, m_ω, m_ρ, m_γ, g2_s, g2_v, g2_ρ, g2_γ, κ_ss, λ, ζ, Λv) = (s.param.m_s, s.param.m_ω, s.param.m_ρ, s.param.m_γ, s.param.g2_s, s.param.g2_v, s.param.g2_ρ, s.param.g2_γ, s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
     (Φ0, W0, B0, A0) = (s.Φ0, s.W0, s.B0, s.A0)
     (ρ_sp, ρ_vp, ρ_sn, ρ_vn) = (s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)
 
@@ -71,6 +57,8 @@ end
 
 "Calculate the total energy associated with meson fields"
 function meson_E(s::system)
+    (κ_ss, λ, ζ, Λv) = (s.param.κ_ss, s.param.λ, s.param.ζ, s.param.Λv)
+
     E_densities = map(zip(s.Φ0, s.W0, s.B0, s.A0, s.ρ_sp, s.ρ_vp, s.ρ_sn, s.ρ_vn)) do (Φ0, W0, B0, A0, ρ_sp, ρ_vp, ρ_sn, ρ_vn)
         E_σ = (1/2) * (Φ0/ħc) * (ρ_sp + ρ_sn) - ((κ_ss/ħc)/12 * (Φ0/ħc)^3 + (λ/24) * (Φ0/ħc)^4)
         E_ω = -(1/2) * (W0/ħc) * (ρ_vp + ρ_vn) + (ζ/24) * (W0/ħc)^4
diff --git a/system.jl b/system.jl
index a5066f7..0298438 100644
--- a/system.jl
+++ b/system.jl
@@ -1,3 +1,30 @@
+"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}
@@ -13,6 +40,7 @@ mutable struct system
     Z::Int
     N::Int
 
+    param::parameters
     r_max::Float64
     divs::Int
 
@@ -30,10 +58,10 @@ mutable struct system
     ρ_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]...)
+    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"
-    system(r_max, divs) = system(0, 0, r_max, divs)
+    "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"
diff --git a/test/Linear_nucleon_spectrum.jl b/test/Linear_nucleon_spectrum.jl
index 0c7d036..f9a27bb 100644
--- a/test/Linear_nucleon_spectrum.jl
+++ b/test/Linear_nucleon_spectrum.jl
@@ -3,7 +3,7 @@ include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
-test_data = readdlm("test/LinearFldsFSUGarnet.csv") 
+test_data = readdlm("test/LinearFldsFSUGarnet.csv")
 xs = test_data[:, 1]
 Ss = test_data[:, 2]
 Vs = test_data[:, 3]
@@ -13,7 +13,7 @@ As = test_data[:, 5]
 p = false
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(8, 8, r_max, divs)
+s = system(8, 8, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
@@ -29,7 +29,7 @@ bench_data, _ = readdlm("test/LinearSpectrum.csv", ','; header=true)
 mask = findall(bench_data[:, 1] .== Int(p))
 bench_κs = bench_data[mask, 2]
 bench_Ebinds = bench_data[mask, 6]
-bench_filling = bench_data[mask, 3] 
+bench_filling = bench_data[mask, 3]
 bench_occ = Int.(bench_filling .* (2 .* j_κ.(Int.(bench_κs)) .+ 1))
 
 scatter(κs, Ebinds, label=(p ? "proton" : "neutron") * " spectrum")
diff --git a/test/Pb208.jl b/test/Pb208.jl
index 38fe28a..b79772d 100644
--- a/test/Pb208.jl
+++ b/test/Pb208.jl
@@ -2,9 +2,13 @@ include("../NuclearRMF.jl")
 
 Z = 82
 N = 126
+
+# Parameter values calibrated by FSUGarnet for Pb-208
+param = parameters("496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644")
+
 r_max = 20.0
 divs = 400
 
-s = system(Z, N, r_max, divs)
+s = system(Z, N, param, r_max, divs)
 
 solve_system!(s; live_plots=false)
diff --git a/test/Pb208ParamsFSUGarnet.txt b/test/Pb208ParamsFSUGarnet.txt
new file mode 100644
index 0000000..fd58700
--- /dev/null
+++ b/test/Pb208ParamsFSUGarnet.txt
@@ -0,0 +1 @@
+496.939473213388, 782.5, 763.0, 0.000001000, 110.349189097820, 187.694676506801, 192.927428365698, 0.091701236, 3.260178893447, -0.003551486718, 0.023499504053, 0.043376933644
\ No newline at end of file
diff --git a/test/Pb208_meson_flds.jl b/test/Pb208_meson_flds.jl
index ccab61b..9d1369d 100644
--- a/test/Pb208_meson_flds.jl
+++ b/test/Pb208_meson_flds.jl
@@ -3,7 +3,7 @@ include("../mesons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
-bench_data = readdlm("test/Pb208FldsFSUGarnet.csv") 
+bench_data = readdlm("test/Pb208FldsFSUGarnet.csv")
 xs_bench = bench_data[:, 1]
 Φ0_bench = bench_data[:, 2]
 W0_bench = bench_data[:, 3]
@@ -17,9 +17,12 @@ plot(xs_bench, hcat(Φ0_bench, W0_bench, B0_bench, A0_bench), layout=4, label=["
 test_data = readdlm("test/Pb208DensFSUGarnet.csv")
 xs = test_data[:, 1]
 
+N_p = 82
+N_n = 126
+param = parameters(read("test/Pb208ParamsFSUGarnet.txt", String))
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(r_max, divs)
+s = system(N_p, N_n, param, r_max, divs)
 
 s.ρ_sn = test_data[:, 2]
 s.ρ_vn = test_data[:, 3]
diff --git a/test/Pb208_nucleon_dens.jl b/test/Pb208_nucleon_dens.jl
index 1e11111..5ee0b13 100644
--- a/test/Pb208_nucleon_dens.jl
+++ b/test/Pb208_nucleon_dens.jl
@@ -3,7 +3,7 @@ include("../nucleons.jl")
 
 # test data generated from Hartree.f
 # format: x S(x) V(x) R(x) A(x)
-test_data = readdlm("test/Pb208FldsFSUGarnet.csv") 
+test_data = readdlm("test/Pb208FldsFSUGarnet.csv")
 xs = test_data[:, 1]
 Ss = test_data[:, 2]
 Vs = test_data[:, 3]
@@ -14,7 +14,7 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs
diff --git a/test/Pb208_nucleon_spectrum.jl b/test/Pb208_nucleon_spectrum.jl
index da55a59..eded7cc 100644
--- a/test/Pb208_nucleon_spectrum.jl
+++ b/test/Pb208_nucleon_spectrum.jl
@@ -15,7 +15,7 @@ N_p = 82
 N_n = 126
 r_max = maximum(xs)
 divs = length(xs) - 1
-s = system(N_p, N_n, r_max, divs)
+s = system(N_p, N_n, parameters(), r_max, divs)
 
 s.Φ0 = Ss
 s.W0 = Vs