diff --git a/src/AtomicAndPhysicalConstants.jl b/src/AtomicAndPhysicalConstants.jl
index e846f47..83d50a3 100644
--- a/src/AtomicAndPhysicalConstants.jl
+++ b/src/AtomicAndPhysicalConstants.jl
@@ -7,6 +7,7 @@ using JSON
 using Reexport
 @reexport using Unitful
 
+include("units_definition.jl")
 include("physical_constants.jl")
 include("atomic_isotopes.jl")
 include("subatomic_species.jl")
@@ -17,13 +18,10 @@ include("update_iso_masses.jl")
 include("particle_functions.jl")
 include("set_units.jl")
 
-export setunits, printunits
-export PARTICLE_PHYSICS, MKS, CGS
-export mass, charge
-export c_light, m_electron, m_proton, m_neutron, m_muon, m_helion, m_deuteron, m_pion_0, m_pion_charged
-export r_e, e_charge, h_planck, mu_0_vac, eps_0_vac
-export kg_per_amu, eV_per_amu, N_avogadro, fine_structure, classical_radius_factor, r_p, h_bar_planck, kg_per_eV, eps_0_vac
-export mu_deuteron, mu_electron, mu_helion, mu_muon, mu_neutron, mu_proton, mu_triton
+export setunits
+export ACCELERATOR, MKS, CGS
+export massof, chargeof
+export C_LIGHT, H_PLANCK, H_BAR_PLANCK, R_E, R_P, E_CHARGE, MU_0_VAC, EPS_0_VAC, CLASSICAL_RADIUS_FACTOR, FINE_STRUCTURE, N_AVOGADRO
 export SubatomicSpecies
 export AtomicSpecies
 export SUBATOMIC_SPECIES
diff --git a/src/physical_constants.jl b/src/physical_constants.jl
index 344d5c0..6e5e273 100644
--- a/src/physical_constants.jl
+++ b/src/physical_constants.jl
@@ -2,24 +2,6 @@
 # Constants pulled from the NIST table of
 # the 2022 CODATA release
 
-const __b_e_charge = 1.602176634e-19 * u"C"
-# elementary charge [C]
-const __b_N_avogadro = 6.02214076e23
-# Avogadro's constant: Number / mole (exact)
-
-module NewUnits
-using Unitful
-AA = parentmodule(NewUnits)
-@unit e "e" elementary_charge AA.__b_e_charge false
-@unit amu "amu" Amu (1 / (AA.__b_N_avogadro)) * u"g" false
-end
-
-Unitful.register(NewUnits);
-using .NewUnits
-
-function __init__()
-    Unitful.register(NewUnits)
-end
 
 
 # Format is:
@@ -49,6 +31,10 @@ const __b_mu_0_vac = 1.25663706127e-6 * u"N/A^2"
 # Vacuum permeability in [N/A^2] (newtons per ampere squared)
 
 
+const __b_e_charge = 1.602176634e-19 * u"C"
+# elementary charge [C]
+const __b_N_avogadro = 6.02214076e23
+# Avogadro's constant: Number / mole (exact)
 
 
 
diff --git a/src/set_units.jl b/src/set_units.jl
index 85fc74d..47b1ffd 100644
--- a/src/set_units.jl
+++ b/src/set_units.jl
@@ -1,50 +1,26 @@
-"""
-    UnitSystem
-
-## Description:
-This defines an immutable struct for storing a specific
-system of units.
-
-## Fields:
-- `mass`    -- type:`Unitful.FreeUnits`, stores the unit for mass
-- `length`  -- type:`Unitful.FreeUnits`, stores the unit for length
-- `time`    -- type:`Unitful.FreeUnits`, stores the unit for time
-- `energy`  -- type:`Unitful.FreeUnits`, stores the unit for energy
-- `charge`  -- type:`Unitful.FreeUnits`, stores the unit for charge
-"""
-UnitSystem
-
-struct UnitSystem
-  mass::Unitful.FreeUnits
-  length::Unitful.FreeUnits
-  time::Unitful.FreeUnits
-  energy::Unitful.FreeUnits
-  charge::Unitful.FreeUnits
-end
-
 # Declare specific systems of units
 #   for particle physics
 """
-    PARTICLE_PHYSICS::UnitSystem
-## PARTICLE_PHYSICS units:
+    ACCELERATOR
+## ACCELERATOR units:
 - `mass`: eV/c^2
 - `length`: m
 - `time`: s
 - `energy`: eV
 - `charge`: elementary charge
 """
-PARTICLE_PHYSICS
+ACCELERATOR
 
-const PARTICLE_PHYSICS = UnitSystem(
+const ACCELERATOR = [
   u"eV/c^2",
   u"m",
   u"s",
   u"eV",
-  u"e")
+  u"e"]
 
 #   MKS
 """
-    MKS::UnitSystem
+    MKS
 ## MKS units:
 - `mass`: kg
 - `length`: m
@@ -54,15 +30,15 @@ const PARTICLE_PHYSICS = UnitSystem(
 """
 MKS
 
-const MKS = UnitSystem(
+const MKS = [
   u"kg",
   u"m",
   u"s",
   u"J",
-  u"C")
+  u"C"]
 #   quasi-CGS
 """
-    CGS::UnitSystem
+    CGS
 ## CGS units:
 - `mass`: g
 - `length`: cm
@@ -72,76 +48,15 @@ const MKS = UnitSystem(
 """
 CGS
 
-const CGS = UnitSystem(
+const CGS = [
   u"g",
   u"cm",
   u"s",
   u"J",
-  u"C")
-
-
-"""
-    current_units :: UnitSystem
-
-## Description:
-This declares a `UnitSystem` that stores the units in current use.
-
-## Note:
-It is initialized when setunits() is called.
-"""
-# Note that the units don't have the proper dimension, this means set unit is not called
-current_units::UnitSystem = UnitSystem(u"m",
-  u"m",
-  u"m",
-  u"m",
-  u"m")
-
-
-
-"""
-    getunits(unit::Symbol)
-
-## Description:
-return the unit with the corresponding field `unit` in `current_units`
-
-## parameters:
-- `unit`  -- type:`Symbol`, the name of the field in current_units.
-
-"""
-getunit
-
-function getunit(unit::Symbol)
-  if dimension(current_units.mass) != dimension(u"kg")
-    throw(ErrorException("units are not set, call setunits() to initalize units and constants"))
-  else
-    return getfield(current_units, unit)
-  end
-end
-
-"""
-    printunits()
-
-## Description:
-This function returns nothing. It simply prints the set of units
-in current use.
-"""
-printunits
-
-function printunits()
-  if dimension(current_units.mass) != dimension(u"kg")
-    throw(ErrorException("units are not set, call setunits() to initalize units and constants"))
-  end
-  # prints the units for each dimensions
-  println("mass unit:\t", current_units.mass)
-  println("length unit:\t", current_units.length)
-  println("time unit:\t", current_units.time)
-  println("energy unit:\t", current_units.energy)
-  println("charge unit:\t", current_units.charge)
-  return
-end
+  u"C"]
 
 """
-    setunits(unitsystem::UnitSystem=PARTICLE_PHYSICS;
+    setunits(unitsystem::UnitSystem=ACCELERATOR;
       mass_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.mass,
       length_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.length,
       time_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.time,
@@ -184,30 +99,13 @@ Prints current units at the end (optional).
 """
 setunits
 
-function setunits(unitsystem::UnitSystem=PARTICLE_PHYSICS;
-  mass_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.mass,
-  length_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.length,
-  time_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.time,
-  energy_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.energy,
-  charge_unit::Union{Unitful.FreeUnits,AbstractString}=unitsystem.charge,
-  print_units::Bool=true
+function setunits(unitsystem=ACCELERATOR;
+  mass_unit::Unitful.FreeUnits=unitsystem[1],
+  length_unit::Unitful.FreeUnits=unitsystem[2],
+  time_unit::Unitful.FreeUnits=unitsystem[3],
+  energy_unit::Unitful.FreeUnits=unitsystem[4],
+  charge_unit::Unitful.FreeUnits=unitsystem[5],
 )
-  # convert types to Unitful.FreeUnits
-  if mass_unit isa AbstractString
-    mass_unit = uparse(mass_unit)
-  end
-  if length_unit isa AbstractString
-    length_unit = uparse(length_unit)
-  end
-  if time_unit isa AbstractString
-    time_unit = uparse(time_unit)
-  end
-  if energy_unit isa AbstractString
-    energy_unit = uparse(energy_unit)
-  end
-  if charge_unit isa AbstractString
-    charge_unit = uparse(charge_unit)
-  end
   # check dimensions of units
   if dimension(mass_unit) != dimension(u"kg")
     throw(ErrorException("unit for mass does not have proper dimension"))
@@ -224,70 +122,21 @@ function setunits(unitsystem::UnitSystem=PARTICLE_PHYSICS;
   if dimension(charge_unit) != dimension(u"C")
     throw(ErrorException("unit for charge does not have proper dimension"))
   end
-  # record what units is currently being used
-  global current_units = UnitSystem(mass_unit, length_unit, time_unit, energy_unit, charge_unit)
-
-  # convert all the variables with dimension mass
-  AtomicAndPhysicalConstants.m_electron = (__b_m_electron |> mass_unit).val         # Electron Mass 
-  AtomicAndPhysicalConstants.m_proton = (__b_m_proton |> mass_unit).val             # Proton Mass 
-  AtomicAndPhysicalConstants.m_neutron = (__b_m_neutron |> mass_unit).val           # Neutron Mass 
-  AtomicAndPhysicalConstants.m_muon = (__b_m_muon |> mass_unit).val                 # Muon Mass 
-  AtomicAndPhysicalConstants.m_helion = (__b_m_helion |> mass_unit).val             # Helion Mass He3 nucleus 
-  AtomicAndPhysicalConstants.m_deuteron = (__b_m_deuteron |> mass_unit).val         # Deuteron Mass 
-  AtomicAndPhysicalConstants.m_pion_0 = (__b_m_pion_0 |> mass_unit).val             # Pion 0 mass
-  AtomicAndPhysicalConstants.m_pion_charged = (__b_m_pion_charged |> mass_unit).val # Pion+- Mass 
-
-  # convert all the variables with dimension length
-  AtomicAndPhysicalConstants.r_e = (__b_r_e |> length_unit).val                     # classical electron radius
-
-  # convert all the variables with dimension time
-
-  # convert all the variables with dimension speed
-  AtomicAndPhysicalConstants.c_light = (__b_c_light |> length_unit / time_unit).val          # speed of light
-
-  # convert all the variables with dimension energy
-
-  # convert all the variables with dimension charge
-  AtomicAndPhysicalConstants.e_charge = (__b_e_charge |> charge_unit).val                                # elementary charge
-
-  # constants with special dimenisions
-  # convert Planck's constant with dimension energy * time
-  AtomicAndPhysicalConstants.h_planck = (__b_h_planck |> energy_unit * time_unit).val        # Planck's constant 
-  # convert Vacuum permeability with dimension force / (current)^2
-  AtomicAndPhysicalConstants.mu_0_vac = __b_mu_0_vac.val
-  # convert Vacuum permeability with dimension capacitance / distance
-  AtomicAndPhysicalConstants.eps_0_vac = __b_eps_0_vac.val
-
-  # convert magnet moments dimension: energy / magnetic field strength
-  AtomicAndPhysicalConstants.mu_deuteron = (__b_mu_deuteron |> energy_unit / u"T").val    # deuteron magnetic moment
-  AtomicAndPhysicalConstants.mu_electron = (__b_mu_electron |> energy_unit / u"T").val    # electron magnetic moment
-  AtomicAndPhysicalConstants.mu_helion = (__b_mu_helion |> energy_unit / u"T").val        # helion magnetic moment
-  AtomicAndPhysicalConstants.mu_muon = (__b_mu_muon |> energy_unit / u"T").val            # muon magnetic moment
-  AtomicAndPhysicalConstants.mu_neutron = (__b_mu_neutron |> energy_unit / u"T").val      # neutron magnetic moment
-  AtomicAndPhysicalConstants.mu_proton = (__b_mu_proton |> energy_unit / u"T").val        # proton magnetic moment
-  AtomicAndPhysicalConstants.mu_triton = (__b_mu_triton |> energy_unit / u"T").val        # triton magnetic moment
-
-  # convert unitless variables
-  AtomicAndPhysicalConstants.kg_per_amu = __b_kg_per_amu.val               # kg per standard atomic mass unit (dalton)
-  AtomicAndPhysicalConstants.eV_per_amu = __b_eV_per_amu.val                  # eV per standard atomic mass unit (dalton)
-  AtomicAndPhysicalConstants.N_avogadro = __b_N_avogadro                # Number / mole  (exact)
-  AtomicAndPhysicalConstants.fine_structure = __b_fine_structure                # fine structure constant
-
-  # values calculated from other constants
-  AtomicAndPhysicalConstants.classical_radius_factor = r_e * m_electron                 # e^2 / (4 pi eps_0) = classical_radius * mass * c^2. Is same for all particles of charge +/- 1.
-  AtomicAndPhysicalConstants.r_p = r_e * m_electron / m_proton      # proton radius
-  AtomicAndPhysicalConstants.h_bar_planck = h_planck / 2pi                   # h_planck/twopi
-  AtomicAndPhysicalConstants.kg_per_eV = kg_per_amu / eV_per_amu
-  AtomicAndPhysicalConstants.eps_0_vac = 1 / (c_light^2 * mu_0_vac)       # Permeability of free space
-  if print_units
-    printunits()
-  end
 
-  return
+  eval(:(c_light() = uconvert($length_unit / $time_unit, __b_c_light)))
+  eval(:(h_planck() = uconvert($energy_unit * $time_unit, __b_h_planck)))
+  eval(:(h_bar_planck() = uconvert($energy_unit * $time_unit, __b_h_bar_planck)))
+  eval(:(r_e() = uconvert($length_unit, __b_r_e)))
+  eval(:(r_p() = uconvert($length_unit, __b_r_p)))
+  eval(:(e_charge() = uconvert($charge_unit, __b_e_charge)))
+  eval(:(massof(species::Species) = uconvert($mass_unit, species.mass_in_eV * u"eV/c^2")))
+  eval(:(chargeof(species::Species) = uconvert($charge_unit, species.charge * u"e")))
+  return [mass_unit, length_unit, time_unit, energy_unit, charge_unit]
+
 end
 
 """
-    mass(
+    massof(
       species::Species,
       unit::Union{Unitful.FreeUnits,AbstractString}=current_units.mass
     )
@@ -300,20 +149,21 @@ return mass of 'species' in current unit or unit of the user's choice
 - `unit`        -- type:`Union{Unitful.FreeUnits,AbstractString}`, default to the unit set from setunits(), the unit of the mass variable
 
 """
-mass
+massof
 
-function mass(species::Species, unit::Union{Unitful.FreeUnits,AbstractString}=getunit(:mass))
+function massof(species::Species, unit::Union{Unitful.FreeUnits,AbstractString})
   if unit isa AbstractString
     unit = uparse(unit)
   end
   if dimension(unit) != dimension(u"kg")
-    throw(ErrorException("mass unit doesn't have proper dimension"))
+    error("mass unit doesn't have proper dimension")
   end
   return (species.mass_in_eV * u"eV/c^2" |> unit).val
 end
 
+
 """
-    charge(
+    chargeof(
       species::Species,
       unit::Union{Unitful.FreeUnits,AbstractString}=current_units.charge
     )
@@ -326,9 +176,9 @@ return charge of 'species' in current unit or unit of the user's choice
 - `unit`        -- type:`Union{Unitful.FreeUnits,AbstractString}`, default to the unit set from setunits(), the unit of the charge variable
 
 """
-charge
+chargeof
 
-function charge(species::Species, unit::Union{Unitful.FreeUnits,AbstractString}=getunit(:charge))
+function chargeof(species::Species, unit::Union{Unitful.FreeUnits,AbstractString})
   if unit isa AbstractString
     unit = uparse(unit)
   end
@@ -340,41 +190,59 @@ end
 
 
 
-# Define global constants
-m_electron::Float64 = NaN
-m_proton::Float64 = NaN
-m_neutron::Float64 = NaN
-m_muon::Float64 = NaN
-m_helion::Float64 = NaN
-m_deuteron::Float64 = NaN
-m_pion_0::Float64 = NaN
-m_pion_charged::Float64 = NaN
 
-r_e::Float64 = NaN
+function c_light(unit::Unitful.FreeUnits)
+  return __b_c_light |> unit
+end
+
+function h_planck(unit::Unitful.FreeUnits)
+  return __b_h_planck |> unit
+end
+
+
+function h_bar_planck(unit::Unitful.FreeUnits)
+  return __b_h_bar_planck |> unit
+end
+
 
-c_light::Float64 = NaN
+function r_e(unit::Unitful.FreeUnits)
+  return __b_r_e |> unit
+end
 
-e_charge::Float64 = NaN
+function r_p(unit::Unitful.FreeUnits)
+  return __b_r_p |> unit
+end
 
-h_planck::Float64 = NaN
-mu_0_vac::Float64 = NaN
-eps_0_vac::Float64 = NaN
+function e_charge(unit::Unitful.FreeUnits)
+  return __b_e_charge |> unit
+end
 
-mu_deuteron::Float64 = NaN
-mu_electron::Float64 = NaN
-mu_helion::Float64 = NaN
-mu_muon::Float64 = NaN
-mu_neutron::Float64 = NaN
-mu_proton::Float64 = NaN
-mu_triton::Float64 = NaN
 
-kg_per_amu::Float64 = NaN
-eV_per_amu::Float64 = NaN
-N_avogadro::Float64 = NaN
-fine_structure::Float64 = NaN
+function mu_0_vac()
+  return __b_mu_0_vac
+end
+
+function mu_0_vac(unit::Unitful.FreeUnits)
+  return __b_mu_0_vac |> unit
+end
+
+
+function eps_0_vac()
+  return __b_eps_0_vac
+end
+
+function eps_0_vac(unit::Unitful.FreeUnits)
+  return __b_eps_0_vac |> unit
+end
+
+function classical_radius_factor()
+  return __b_classical_radius_factor
+end
+
+function fine_structure()
+  return __b_classical_radius_factor
+end
 
-classical_radius_factor::Float64 = NaN
-r_p::Float64 = NaN
-h_bar_planck::Float64 = NaN
-kg_per_eV::Float64 = NaN
-eps_0_vac::Float64 = NaN
\ No newline at end of file
+function N_avogadro()
+  return __b_N_avogadro
+end
\ No newline at end of file
diff --git a/src/units_definition.jl b/src/units_definition.jl
new file mode 100644
index 0000000..ce4365f
--- /dev/null
+++ b/src/units_definition.jl
@@ -0,0 +1,19 @@
+
+
+# Units pulled from the NIST table of
+# the 2022 CODATA release
+
+module NewUnits
+using Unitful
+AA = parentmodule(NewUnits)
+@unit amu "amu" Amu 1.66053906892 * 10^(-27) * u"kg" false
+@unit e "e" Elementary_charge 1.602176634e-19 * u"C" false
+end
+
+using Unitful
+Unitful.register(NewUnits);
+using .NewUnits
+
+function __init__()
+    Unitful.register(NewUnits)
+end
\ No newline at end of file