diff --git a/nucleons.jl b/nucleons.jl
index 0b0248d..271684c 100644
--- a/nucleons.jl
+++ b/nucleons.jl
@@ -16,7 +16,7 @@ const r_reg = 1E-8 # fm # regulator for the centrifugal term
     f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r) is a function of r in MeV (see optimized_dirac_potentials()),
     r is the radius in fm.
 Reference: P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, Frontiers in Physics 10, (2023)"
-function dirac!(du, u, (κ, E, f1, f2), r)
+function dirac!(du::Vector{Float64}, u::Vector{Float64}, (κ, E, f1, f2), r::Float64) # TODO: Static typing
     (g, f) = u
     @inbounds du[1] = -(κ/(r + r_reg)) * g + (E + f1(r)) * f / ħc
     @inbounds du[2] =  (κ/(r + r_reg)) * f - (E + f2(r)) * g / ħc
@@ -25,7 +25,7 @@ end
 "Get the potentials f1 and f2 that goes into the Dirac equation, defined as
     f1(r) =  M-Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r),
     f2(r) = -M+Φ0(r)-W0(r)-(p-0.5)*2B0(r)-p*A0(r)."
-function optimized_dirac_potentials(p, s::system)
+function optimized_dirac_potentials(p::Bool, s::system)
     M = p ? M_p : M_n
 
     f1s = zero_array(s)
@@ -44,7 +44,7 @@ end
     divs is the number of mesh divisions so solution would be returned as a 2×(1+divs) matrix,
     shooting method divides the interval into two partitions at r_max/2, ensuring convergence at both r=0 and r=r_max,
     the other parameters are the same from dirac!(...)."
-function solveNucleonWf(κ, p, E, s::system; shooting=true, normalize=true, algo=Tsit5())
+function solveNucleonWf(κ, p::Bool, E, s::system; shooting=true, normalize=true, algo=Tsit5())
     (f1, f2) = optimized_dirac_potentials(p, s)
 
     if shooting
@@ -81,7 +81,7 @@ end
 "Returns a function that solves the Dirac equation and returns g(r=r_max) where
     r_max is the outer boundary in fm,
     the other parameters are the same from dirac!(...)."
-function boundaryValueFunc(κ, p, s::system; dtype=Float64, algo=Tsit5())
+function boundaryValueFunc(κ, p::Bool, s::system; dtype=Float64, algo=Tsit5())
     prob = ODEProblem(dirac!, convert.(dtype, [0, 1]), (0, s.r_max))
     (f1, f2) = optimized_dirac_potentials(p, s)
     func(E) = solve(prob, algo, p=(κ, E, f1, f2), saveat=[s.r_max], save_idxs=[1])[1, 1]
@@ -91,7 +91,7 @@ end
 "Find all bound energies between E_min (=850.0) and E_max (=938.0) where
     tol_digits determines the precision for root finding and the threshold for identifying duplicates,
     the other parameters are the same from dirac!(...)."
-function findEs(κ, p, s::system, E_min=850.0, E_max=938.0, tol_digits=5)
+function findEs(κ, p::Bool, s::system, E_min=850.0, E_max=938.0, tol_digits=5)
     func = boundaryValueFunc(κ, p, s)
     Es = find_all_zeros(func, E_min, E_max; partitions=20, tol=1/10^tol_digits)
     return unique(E -> round(E; digits=tol_digits), Es)
@@ -99,7 +99,7 @@ end
 
 "Find all orbitals and return two lists containing κ values and corresponding energies for a single species where
     the other parameters are defined above"
-function findAllOrbitals(p, s::system, E_min=850.0, E_max=938.0)
+function findAllOrbitals(p::Bool, s::system, E_min=850.0, E_max=938.0)
     κs = Int[]
     Es = Float64[]
     # start from κ=-1 and go both up and down
@@ -140,7 +140,7 @@ l_κ(κ::Int) = abs(κ) - (κ < 0) # since true = 1 and false = 0
 "Calculate scalar and vector densities of a nucleon species on [0,r_max] divided into (divs+1) points and returns them as vectors (ρ_s, ρ_v) where
     the arrays κs, Es, occs tabulate the energies and occupation numbers corresponding to each κ,
     the other parameters are defined above"
-function calculateNucleonDensity(κs, Es, occs, p, s::system)
+function calculateNucleonDensity(κs, Es, occs, p::Bool, s::system)
     ρr2 = zeros(2, s.divs + 1) # ρ×r² values
     
     for (κ, E, occ) in zip(κs, Es, occs)
@@ -161,7 +161,7 @@ end
 
 "Solve the Dirac equation and calculate scalar and vector densities of a nucleon species where
     the other parameters are defined above"
-function solveNucleonDensity(p, s::system, E_min=850.0, E_max=938.0)
+function solveNucleonDensity(p::Bool, s::system, E_min=850.0, E_max=938.0)
     κs, Es = findAllOrbitals(p, s, E_min, E_max)
     occs = fillNucleons(Z_or_N(s, p), κs, Es)
     return calculateNucleonDensity(κs, Es, occs, p, s)