diff --git a/Project.toml b/Project.toml
index 8fe2296..1284997 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,12 +4,12 @@ authors = ["Zijin Zhang", "Luke Adams <luke@lukeclydeadams.com>"]
 version = "0.1.0"
 
 [deps]
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulEquivalences = "da9c4bc3-91c8-4f02-8a40-6b990d2a7e0c"
 
 [compat]
 Aqua = "0.8"
+LinearAlgebra = "1.10"
 Test = "1.10"
 Unitful = "1.0"
 UnitfulEquivalences = "0.2"
@@ -17,7 +17,8 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "Test"]
+test = ["Aqua", "LinearAlgebra", "Test"]
diff --git a/src/PlasmaFormulary.jl b/src/PlasmaFormulary.jl
index 08e3631..7fce756 100644
--- a/src/PlasmaFormulary.jl
+++ b/src/PlasmaFormulary.jl
@@ -6,7 +6,6 @@ using Unitful: μ0, ε0, c, q
 using Unitful: k, ħ
 using Unitful: me, mp, u
 using Unitful: Velocity, Mass, BField, Density, Charge
-using LinearAlgebra
 
 @derived_dimension NumberDensity Unitful.𝐋^-3
 
@@ -37,24 +36,20 @@ atmospheric_pressure = 1.01325e5 * u"Pa"
 standard_molar_volume = 2.24140e-2 * u"m^3 / mol"
 molar_weight_of_air = 2.89647e-2 * u"kg / mol"
 
-# Fundamental plasma parameters
-# These formulas have been converted to use SI units from the original Gaussian cgs units
-# that are used in the 2023 edition of the formulary.
-
 include("dimensionless.jl")
+export plasma_beta
+
 include("lengths.jl")
+export gyroradius, electron_gyroradius, electron_debroglie_length, classical_minimum_approach_distance, inertial_length, electron_inertial_length, debye_length
+
 include("speeds.jl")
+export alfven_velocity, thermal_speed, thermal_temperature, electron_thermal_velocity, ion_thermal_velocity
+
 include("frequencies.jl")
+export gyrofrequency, electron_gyrofrequency, ion_gyrofrequency, plasma_frequency
+
 include("misc.jl")
+export thermal_pressure, magnetic_pressure
 
-export plasma_beta
-export gyroradius, debye_length, inertial_length
-export Alfven_velocity, Alfven_speed
-export gyrofrequency, plasma_frequency
-
-# aliases
-const rc_ = gyroradius
-const di_ = ion_inertial_length
-export rc_, di_
 
 end
diff --git a/src/frequencies.jl b/src/frequencies.jl
index 18d2709..8ea905c 100644
--- a/src/frequencies.jl
+++ b/src/frequencies.jl
@@ -1,4 +1,4 @@
-function gyrofrequency(B::BField, mass, q)
+function gyrofrequency(B::BField, mass::Mass, q::Charge)
     upreferred(q * B / mass)
 end
 
@@ -6,7 +6,7 @@ function electron_gyrofrequency(B::BField)
     gyrofrequency(B, me, Unitful.q)
 end
 
-function ion_gyrofrequency(B::BField, Z, mass)
+function ion_gyrofrequency(B::BField, Z, mass::Mass)
     gyrofrequency(B, mass, Z * Unitful.q)
 end
 
@@ -29,8 +29,9 @@ plasma_frequency(n::NumberDensity) = plasma_frequency(n, Unitful.q, me)
 
 Ion plasma frequency.
 """
-plasma_frequency(n::NumberDensity, Z::Integer, mass_numb) =
-    plasma_frequency(n, Z * Unitful.q, mass_numb * Unitful.u)
+plasma_frequency(n::NumberDensity, Z::Integer, mass_number) =
+    plasma_frequency(n, Z * Unitful.q, mass_number * Unitful.u)
 
+# TODO: Do we need these?
 const electron_plasma_frequency = plasma_frequency
 const ion_plasma_frequency = plasma_frequency
diff --git a/src/lengths.jl b/src/lengths.jl
index 34741df..8a91dbe 100644
--- a/src/lengths.jl
+++ b/src/lengths.jl
@@ -9,27 +9,17 @@ julia> gyroradius(0.2u"T", 1e6u"K")
 0.006682528174870854 m
 ```
 """
-function gyroradius(B, mass, q, Vperp::Velocity)
+function gyroradius(B::BField, mass::Mass, q::Charge, Vperp::Velocity)
     return upreferred(mass * Vperp / (q * B))
 end
 
-function gyroradius(B, Vperp::Velocity; mass_numb = 1, Z = 1)
-    return gyroradius(B, mass_numb * Unitful.u, Z * Unitful.q, Vperp)
+function gyroradius(B::BField, mass::Mass, q::Charge, T::EnergyOrTemp)
+    Vperp = thermal_speed(T, mass)
+    return gyroradius(B, mass, q, Vperp)
 end
 
-function gyroradius(B, T::EnergyOrTemp; mass_numb = 1, Z = 1)
-    Vperp = thermal_speed(T, mass_numb)
-    return gyroradius(B, Vperp; mass_numb, Z)
-end
-
-gyroradius(val::Union{Velocity,EnergyOrTemp}, B::BField; kw...) = gyroradius(B, val; kw...)
-
-electron_gyroradius(B, Vperp::Velocity) = gyroradius(B, me, Unitful.q, Vperp)
-electron_gyroradius(B, T::EnergyOrTemp) = electron_gyroradius(B, thermal_speed(T, me))
-
-# electron and ion trapping rates excluded
-
-#electron and ion collision rates excluded
+electron_gyroradius(B::BField, Vperp::Velocity) = gyroradius(B, me, Unitful.q, Vperp)
+electron_gyroradius(B::BField, T::EnergyOrTemp) = electron_gyroradius(B, thermal_speed(T, me))
 
 function electron_debroglie_length(eot::EnergyOrTemp)
     upreferred(sqrt(2 * pi * ħ^2 / me / energy(eot)))
@@ -44,7 +34,7 @@ The inertial length is the characteristic length scale for a particle to be acce
 
 References: [PlasmaPy API Documentation](https://docs.plasmapy.org/en/latest/api/plasmapy.formulary.lengths.inertial_length.html)
 """
-function inertial_length(n::NumberDensity, q, mass::Mass)
+function inertial_length(n::NumberDensity, q::Charge, mass::Mass)
     upreferred(c / plasma_frequency(n, q, mass))
 end
 
@@ -54,8 +44,6 @@ function ion_inertial_length(n::NumberDensity, Z, mass::Mass)
     inertial_length(n, Z * Unitful.q, mass)
 end
 
-ion_inertial_length(n::NumberDensity; Z = 1, mass = mp) = ion_inertial_length(n, Z, mass)
-
 function debye_length(density::NumberDensity, eot::EnergyOrTemp)
     upreferred(sqrt(ε0 * energy(eot) / density / q^2))
-end
\ No newline at end of file
+end
diff --git a/src/misc.jl b/src/misc.jl
index 3048c1e..8d6ec24 100644
--- a/src/misc.jl
+++ b/src/misc.jl
@@ -11,8 +11,6 @@ function thermal_pressure(T::EnergyOrTemp, n::NumberDensity)
     return n * k * temperature(T) |> upreferred
 end
 
-const p_th = thermal_pressure
-
 """
     magnetic_pressure(B)
 
@@ -20,4 +18,4 @@ Calculate the magnetic pressure.
 """
 function magnetic_pressure(B::BField)
     return (B^2) / (2 * μ0)
-end
\ No newline at end of file
+end
diff --git a/src/speeds.jl b/src/speeds.jl
index f15fb5a..e7f300d 100644
--- a/src/speeds.jl
+++ b/src/speeds.jl
@@ -1,93 +1,73 @@
 # TODO: Sound speed
 
-function Alfven_velocity(B::BField, ρ::Density)
+"""
+The typical propagation speed of magnetic disturbances in a quasineutral plasma.
+
+References: [PlasmaPy API Documentation](https://docs.plasmapy.org/en/stable/api/plasmapy.formulary.speeds.Alfven_speed.html)
+"""
+function alfven_velocity(B::BField, ρ::Density)
     return B / sqrt(μ0 * ρ) |> upreferred
 end
 
-function Alfven_velocity(B::BField, n::NumberDensity; mass_numb = 1, Z = 1)
-    ρ = n * (mass_numb * mp + Z * me)
-    return Alfven_velocity(B, ρ)
+function alfven_velocity(B::BField, n::NumberDensity, mass_number)
+    ρ = n * mass_number * mp
+    return alfven_velocity(B, ρ)
 end
 
-Alfven_velocity(B::AbstractVector, args...; kwargs...) =
-    [Alfven_velocity(Bi, args...; kwargs...) for Bi in B]
+# TODO: Add docstrings
+abstract type ThermalSpeedMethod end
+struct MostProbable <: ThermalSpeedMethod end
+struct RMS <: ThermalSpeedMethod end
+struct MeanMagnitude <: ThermalSpeedMethod end
+struct NRL <: ThermalSpeedMethod end
 
 """
-    Alfven_speed(args...; kwargs...)
+    thermal_speed_coefficients(method::ThermalSpeedMethod, ndim::Int)
 
-The typical propagation speed of magnetic disturbances in a quasineutral plasma.
+Get the thermal speed coefficient corresponding to the desired thermal speed definition.
 
-References: [PlasmaPy API Documentation](https://docs.plasmapy.org/en/stable/api/plasmapy.formulary.speeds.Alfven_speed.html)
+# Arguments
+- `method::ThermalSpeedMethod`: Method to be used for calculating the thermal speed. Valid values are `MostProbable()`, `RMS()`, `MeanMagnitude()`, and `NRL()`.
+- `ndim::Val{Int}`: Dimensionality (1D, 2D, 3D) of space in which to calculate thermal speed. Valid values are `Val(1)`, `Val(2)`, or `Val{3}`.
 """
-Alfven_speed(args...; kwargs...) = norm(Alfven_velocity(args...; kwargs...))
-
+thermal_speed_coefficients(::MostProbable, ::Val{1}) = 0.0
+thermal_speed_coefficients(::MostProbable, ::Val{2}) = 1.0
+thermal_speed_coefficients(::MostProbable, ::Val{3}) = sqrt(2)
+thermal_speed_coefficients(::RMS, ::Val{1}) = 1.0
+thermal_speed_coefficients(::RMS, ::Val{2}) = sqrt(2)
+thermal_speed_coefficients(::RMS, ::Val{3}) = sqrt(3)
+thermal_speed_coefficients(::MeanMagnitude, ::Val{1}) = sqrt(2 / π)
+thermal_speed_coefficients(::MeanMagnitude, ::Val{2}) = sqrt(π / 2)
+thermal_speed_coefficients(::MeanMagnitude, ::Val{3}) = sqrt(8 / π)
+thermal_speed_coefficients(::NRL, ::Val{1}) = 1.0
+thermal_speed_coefficients(::NRL, ::Val{2}) = 1.0
+thermal_speed_coefficients(::NRL, ::Val{3}) = 1.0
 
 function thermal_speed(
     T::EnergyOrTemp,
     mass::Unitful.Mass,
-    method = "most_probable",
+    method::ThermalSpeedMethod = MostProbable(),
     ndim = 3,
 )
-    coeff = thermal_speed_coefficients(method, ndim)
+    coeff = thermal_speed_coefficients(method, Val(ndim))
     return coeff * sqrt(k * temperature(T) / mass)
 end
 
-thermal_speed(T::EnergyOrTemp, mass_numb, args...) =
-    thermal_speed(T, mass_numb * u, args...)
-
 function thermal_temperature(
     V::Unitful.Velocity,
     mass::Unitful.Mass,
-    method = "most_probable",
+    method::ThermalSpeedMethod = MostProbable(),
     ndim = 3,
 )
-    coeff = thermal_speed_coefficients(method, ndim)
+    coeff = thermal_speed_coefficients(method, Val(ndim))
     return mass * V^2 / (k * coeff^2) |> upreferred
 end
 
-"""
-    thermal_speed_coefficients(method::String, ndim::Int)
-
-Get the thermal speed coefficient corresponding to the desired thermal speed definition.
-
-# Arguments
-- `method::String`: Method to be used for calculating the thermal speed. Valid values are `"most_probable"`, `"rms"`, `"mean_magnitude"`, and `"nrl"`.
-- `ndim::Int`: Dimensionality (1D, 2D, 3D) of space in which to calculate thermal speed. Valid values are `1`, `2`, or `3`.
-"""
-function thermal_speed_coefficients(method::String, ndim::Int)
-    # Define the coefficient dictionary
-    _coefficients = Dict(
-        (1, "most_probable") => 0.0,
-        (2, "most_probable") => 1.0,
-        (3, "most_probable") => sqrt(2),
-        (1, "rms") => 1.0,
-        (2, "rms") => sqrt(2),
-        (3, "rms") => sqrt(3),
-        (1, "mean_magnitude") => sqrt(2 / π),
-        (2, "mean_magnitude") => sqrt(π / 2),
-        (3, "mean_magnitude") => sqrt(8 / π),
-        (1, "nrl") => 1.0,
-        (2, "nrl") => 1.0,
-        (3, "nrl") => 1.0,
-    )
-
-    # Attempt to retrieve the coefficient
-    try
-        return _coefficients[(ndim, method)]
-    catch
-        throw(
-            ArgumentError(
-                "Value for (ndim, method) pair not valid, got '($ndim, $method)'.",
-            ),
-        )
-    end
-end
-
-
 function electron_thermal_velocity(eot::EnergyOrTemp)
-    upreferred(sqrt(k * temperature(eot) / me))
+    thermal_speed(eot, me, NRL())
 end
 
+# TODO: Do we need this, or is the thermal_speed function sufficient?
 function ion_thermal_velocity(eot::EnergyOrTemp, ion_mass::Unitful.Mass)
-    upreferred(sqrt(k * temperature(eot) / ion_mass))
-end
\ No newline at end of file
+    thermal_speed(eot, ion_mass, NRL())
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 76536fc..f0ab2ae 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,13 +2,30 @@ using PlasmaFormulary
 using Test
 using Aqua
 using Unitful
+using LinearAlgebra
 
 @testset "PlasmaFormulary.jl" begin
     @testset "Code quality (Aqua.jl)" begin
         Aqua.test_all(PlasmaFormulary)
     end
 
-    @testset "Conversions" begin
-        @test isapprox(uconvert(u"erg/K", PlasmaFormulary.boltzmann_constant), 1.38065e-16 * u"erg/K")
+    @testset "Misc" begin
+        @test gyroradius(0.2u"T", Unitful.mp, Unitful.q, 1e6u"K") ≈ 0.0067u"m" rtol=0.01
+        @test_broken electron_gyroradius(1e6u"K", 0.2u"T", Unitful.q, Unitful.mp) ≈ 0.0067u"m" rtol=0.01
+
+        B = [-0.014u"T", 0.0u"T", 0.0u"T"]
+        Bmag = norm(B)
+        n = 5e19 * u"m^-3"
+        mass_number = 1
+        rho = n * mass_number * Unitful.mp
+        @test norm(alfven_velocity.(B, rho)) ==
+              alfven_velocity(Bmag, rho) ==
+              alfven_velocity(Bmag, n, mass_number) ≈
+              43185.625u"m/s"
+
+        @test alfven_velocity.(B, rho) ≈ [-43185.625u"m/s", 0.0u"m/s", 0.0u"m/s"]
+        @test plasma_frequency(1e19u"m^-3", Unitful.q, Unitful.mp) ≈ 4163294534.0u"s^-1"
+        @test plasma_frequency(1e19u"m^-3", 1, Unitful.mp / Unitful.u) ≈ 4163294534.0u"s^-1"
+        @test plasma_frequency(1e19u"m^-3") ≈ 1.783986365e11u"s^-1"
     end
 end