diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index be7fa7483..8dd814a29 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -88,13 +88,37 @@ public:
   //! clear data from a single coordinate level
   void reset();
 
-  Position r;  //!< particle position
-  Direction u; //!< particle direction
-  int cell {-1};
-  int universe {-1};
-  int lattice {-1};
-  array<int, 3> lattice_i {{-1, -1, -1}};
-  bool rotated {false}; //!< Is the level rotated?
+  // accessors
+  Position& r() { return r_; }
+  const Position& r() const { return r_; }
+
+  Direction& u() { return u_; }
+  const Direction& u() const { return u_; }
+
+  int& cell() { return cell_; }
+  const int& cell() const { return cell_; }
+
+  int& universe() { return universe_; }
+  const int& universe() const { return universe_; }
+
+  int& lattice() { return lattice_; }
+  int lattice() const { return lattice_; }
+
+  array<int, 3>& lattice_index() { return lattice_index_; }
+  const array<int, 3>& lattice_index() const { return lattice_index_; }
+
+  bool& rotated() { return rotated_; }
+  const bool& rotated() const { return rotated_; }
+
+private:
+  // Data members
+  Position r_;  //!< particle position
+  Direction u_; //!< particle direction
+  int cell_ {-1};
+  int universe_ {-1};
+  int lattice_ {-1};
+  array<int, 3> lattice_index_ {{-1, -1, -1}};
+  bool rotated_ {false}; //!< Is the level rotated?
 };
 
 //==============================================================================
@@ -301,20 +325,20 @@ public:
   const Position& u_last() const { return u_last_; }
 
   // Accessors for position in global coordinates
-  Position& r() { return coord_[0].r; }
-  const Position& r() const { return coord_[0].r; }
+  Position& r() { return coord_[0].r(); }
+  const Position& r() const { return coord_[0].r(); }
 
   // Accessors for position in local coordinates
-  Position& r_local() { return coord_[n_coord_ - 1].r; }
-  const Position& r_local() const { return coord_[n_coord_ - 1].r; }
+  Position& r_local() { return coord_[n_coord_ - 1].r(); }
+  const Position& r_local() const { return coord_[n_coord_ - 1].r(); }
 
   // Accessors for direction in global coordinates
-  Direction& u() { return coord_[0].u; }
-  const Direction& u() const { return coord_[0].u; }
+  Direction& u() { return coord_[0].u(); }
+  const Direction& u() const { return coord_[0].u(); }
 
   // Accessors for direction in local coordinates
-  Direction& u_local() { return coord_[n_coord_ - 1].u; }
-  const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
+  Direction& u_local() { return coord_[n_coord_ - 1].u(); }
+  const Direction& u_local() const { return coord_[n_coord_ - 1].u(); }
 
   // Surface token for the surface that the particle is currently on
   int& surface() { return surface_; }
diff --git a/include/openmc/plot.h b/include/openmc/plot.h
index 69e801c5f..7e27679ea 100644
--- a/include/openmc/plot.h
+++ b/include/openmc/plot.h
@@ -223,7 +223,7 @@ T SlicePlotBase::get_map() const
     GeometryState p;
     p.r() = xyz;
     p.u() = dir;
-    p.coord(0).universe = model::root_universe;
+    p.coord(0).universe() = model::root_universe;
     int level = slice_level_;
     int j {};
 
diff --git a/src/cell.cpp b/src/cell.cpp
index 4b2699229..eb5c228bc 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -1314,10 +1314,10 @@ vector<ParentCell> Cell::find_parent_cells(
   bool cell_found = false;
   for (auto it = coords.begin(); it != coords.end(); it++) {
     const auto& coord = *it;
-    const auto& cell = model::cells[coord.cell];
+    const auto& cell = model::cells[coord.cell()];
     // if the cell at this level matches the current cell, stop adding to the
     // stack
-    if (coord.cell == model::cell_map[this->id_]) {
+    if (coord.cell() == model::cell_map[this->id_]) {
       cell_found = true;
       break;
     }
@@ -1327,10 +1327,10 @@ vector<ParentCell> Cell::find_parent_cells(
     int lattice_idx = C_NONE;
     if (cell->type_ == Fill::LATTICE) {
       const auto& next_coord = *(it + 1);
-      lattice_idx = model::lattices[next_coord.lattice]->get_flat_index(
-        next_coord.lattice_i);
+      lattice_idx = model::lattices[next_coord.lattice()]->get_flat_index(
+        next_coord.lattice_index());
     }
-    stack.push(coord.universe, {coord.cell, lattice_idx});
+    stack.push(coord.universe(), {coord.cell(), lattice_idx});
   }
 
   // if this loop finished because the cell was found and
diff --git a/src/dagmc.cpp b/src/dagmc.cpp
index 134360886..1d92d9d57 100644
--- a/src/dagmc.cpp
+++ b/src/dagmc.cpp
@@ -430,7 +430,7 @@ bool DAGUniverse::find_cell(GeometryState& p) const
   // cells, place it in the implicit complement
   bool found = Universe::find_cell(p);
   if (!found && model::universe_map[this->id_] != model::root_universe) {
-    p.lowest_coord().cell = implicit_complement_idx();
+    p.lowest_coord().cell() = implicit_complement_idx();
     found = true;
   }
   return found;
@@ -676,7 +676,7 @@ std::pair<double, int32_t> DAGCell::distance(
     p->history().reset();
   }
 
-  const auto& univ = model::universes[p->lowest_coord().universe];
+  const auto& univ = model::universes[p->lowest_coord().universe()];
 
   DAGUniverse* dag_univ = static_cast<DAGUniverse*>(univ.get());
   if (!dag_univ)
diff --git a/src/geometry.cpp b/src/geometry.cpp
index 9e975c00d..cefe7d364 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -38,17 +38,17 @@ bool check_cell_overlap(GeometryState& p, bool error)
 
   // Loop through each coordinate level
   for (int j = 0; j < n_coord; j++) {
-    Universe& univ = *model::universes[p.coord(j).universe];
+    Universe& univ = *model::universes[p.coord(j).universe()];
 
     // Loop through each cell on this level
     for (auto index_cell : univ.cells_) {
       Cell& c = *model::cells[index_cell];
-      if (c.contains(p.coord(j).r, p.coord(j).u, p.surface())) {
-        if (index_cell != p.coord(j).cell) {
+      if (c.contains(p.coord(j).r(), p.coord(j).u(), p.surface())) {
+        if (index_cell != p.coord(j).cell()) {
           if (error) {
             fatal_error(
               fmt::format("Overlapping cells detected: {}, {} on universe {}",
-                c.id_, model::cells[p.coord(j).cell]->id_, univ.id_));
+                c.id_, model::cells[p.coord(j).cell()]->id_, univ.id_));
           }
           return true;
         }
@@ -73,7 +73,7 @@ int cell_instance_at_level(const GeometryState& p, int level)
   }
 
   // determine the cell instance
-  Cell& c {*model::cells[p.coord(level).cell]};
+  Cell& c {*model::cells[p.coord(level).cell()]};
 
   // quick exit if this cell doesn't have distribcell instances
   if (c.distribcell_index_ == C_NONE)
@@ -82,13 +82,13 @@ int cell_instance_at_level(const GeometryState& p, int level)
   // compute the cell's instance
   int instance = 0;
   for (int i = 0; i < level; i++) {
-    const auto& c_i {*model::cells[p.coord(i).cell]};
+    const auto& c_i {*model::cells[p.coord(i).cell()]};
     if (c_i.type_ == Fill::UNIVERSE) {
       instance += c_i.offset_[c.distribcell_index_];
     } else if (c_i.type_ == Fill::LATTICE) {
       instance += c_i.offset_[c.distribcell_index_];
-      auto& lat {*model::lattices[p.coord(i + 1).lattice]};
-      const auto& i_xyz {p.coord(i + 1).lattice_i};
+      auto& lat {*model::lattices[p.coord(i + 1).lattice()]};
+      const auto& i_xyz {p.coord(i + 1).lattice_index()};
       if (lat.are_valid_indices(i_xyz)) {
         instance += lat.offset(c.distribcell_index_, i_xyz);
       }
@@ -111,7 +111,7 @@ bool find_cell_inner(
       i_cell = *it;
 
       // Make sure the search cell is in the same universe.
-      int i_universe = p.lowest_coord().universe;
+      int i_universe = p.lowest_coord().universe();
       if (model::cells[i_cell]->universe_ != i_universe)
         continue;
 
@@ -120,7 +120,7 @@ bool find_cell_inner(
       Direction u {p.u_local()};
       auto surf = p.surface();
       if (model::cells[i_cell]->contains(r, u, surf)) {
-        p.lowest_coord().cell = i_cell;
+        p.lowest_coord().cell() = i_cell;
         found = true;
         break;
       }
@@ -146,7 +146,7 @@ bool find_cell_inner(
     // code below this conditional, we set i_cell back to C_NONE to indicate
     // that.
     if (i_cell == C_NONE) {
-      int i_universe = p.lowest_coord().universe;
+      int i_universe = p.lowest_coord().universe();
       const auto& univ {model::universes[i_universe]};
       found = univ->find_cell(p);
     }
@@ -154,7 +154,7 @@ bool find_cell_inner(
     if (!found) {
       return found;
     }
-    i_cell = p.lowest_coord().cell;
+    i_cell = p.lowest_coord().cell();
 
     // Announce the cell that the particle is entering.
     if (found && verbose) {
@@ -186,14 +186,14 @@ bool find_cell_inner(
 
       // Set the lower coordinate level universe.
       auto& coord {p.coord(p.n_coord())};
-      coord.universe = c.fill_;
+      coord.universe() = c.fill_;
 
       // Set the position and direction.
-      coord.r = p.r_local();
-      coord.u = p.u_local();
+      coord.r() = p.r_local();
+      coord.u() = p.u_local();
 
       // Apply translation.
-      coord.r -= c.translation_;
+      coord.r() -= c.translation_;
 
       // Apply rotation.
       if (!c.rotation_.empty()) {
@@ -208,11 +208,11 @@ bool find_cell_inner(
 
       // Set the position and direction.
       auto& coord {p.coord(p.n_coord())};
-      coord.r = p.r_local();
-      coord.u = p.u_local();
+      coord.r() = p.r_local();
+      coord.u() = p.u_local();
 
       // Apply translation.
-      coord.r -= c.translation_;
+      coord.r() -= c.translation_;
 
       // Apply rotation.
       if (!c.rotation_.empty()) {
@@ -220,21 +220,21 @@ bool find_cell_inner(
       }
 
       // Determine lattice indices.
-      auto& i_xyz {coord.lattice_i};
-      lat.get_indices(coord.r, coord.u, i_xyz);
+      auto& i_xyz {coord.lattice_index()};
+      lat.get_indices(coord.r(), coord.u(), i_xyz);
 
       // Get local position in appropriate lattice cell
-      coord.r = lat.get_local_position(coord.r, i_xyz);
+      coord.r() = lat.get_local_position(coord.r(), i_xyz);
 
       // Set lattice indices.
-      coord.lattice = c.fill_;
+      coord.lattice() = c.fill_;
 
       // Set the lower coordinate level universe.
       if (lat.are_valid_indices(i_xyz)) {
-        coord.universe = lat[i_xyz];
+        coord.universe() = lat[i_xyz];
       } else {
         if (lat.outer_ != NO_OUTER_UNIVERSE) {
-          coord.universe = lat.outer_;
+          coord.universe() = lat.outer_;
         } else {
           p.mark_as_lost(fmt::format(
             "Particle {} left lattice {}, but it has no outer definition.",
@@ -261,7 +261,7 @@ bool neighbor_list_find_cell(GeometryState& p, bool verbose)
 
   // Get the cell this particle was in previously.
   auto coord_lvl = p.n_coord() - 1;
-  auto i_cell = p.coord(coord_lvl).cell;
+  auto i_cell = p.coord(coord_lvl).cell();
   Cell& c {*model::cells[i_cell]};
 
   // Search for the particle in that cell's neighbor list.  Return if we
@@ -275,15 +275,15 @@ bool neighbor_list_find_cell(GeometryState& p, bool verbose)
   // neighboring cell.
   found = find_cell_inner(p, nullptr, verbose);
   if (found)
-    c.neighbors_.push_back(p.coord(coord_lvl).cell);
+    c.neighbors_.push_back(p.coord(coord_lvl).cell());
   return found;
 }
 
 bool exhaustive_find_cell(GeometryState& p, bool verbose)
 {
-  int i_universe = p.lowest_coord().universe;
+  int i_universe = p.lowest_coord().universe();
   if (i_universe == C_NONE) {
-    p.coord(0).universe = model::root_universe;
+    p.coord(0).universe() = model::root_universe;
     p.n_coord() = 1;
     i_universe = model::root_universe;
   }
@@ -299,32 +299,32 @@ bool exhaustive_find_cell(GeometryState& p, bool verbose)
 void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
 {
   auto& coord {p.lowest_coord()};
-  auto& lat {*model::lattices[coord.lattice]};
+  auto& lat {*model::lattices[coord.lattice()]};
 
   if (verbose) {
     write_message(
       fmt::format("    Crossing lattice {}. Current position ({},{},{}). r={}",
-        lat.id_, coord.lattice_i[0], coord.lattice_i[1], coord.lattice_i[2],
-        p.r()),
+        lat.id_, coord.lattice_index()[0], coord.lattice_index()[1],
+        coord.lattice_index()[2], p.r()),
       1);
   }
 
   // Set the lattice indices.
-  coord.lattice_i[0] += boundary.lattice_translation[0];
-  coord.lattice_i[1] += boundary.lattice_translation[1];
-  coord.lattice_i[2] += boundary.lattice_translation[2];
+  coord.lattice_index()[0] += boundary.lattice_translation[0];
+  coord.lattice_index()[1] += boundary.lattice_translation[1];
+  coord.lattice_index()[2] += boundary.lattice_translation[2];
 
   // Set the new coordinate position.
   const auto& upper_coord {p.coord(p.n_coord() - 2)};
-  const auto& cell {model::cells[upper_coord.cell]};
-  Position r = upper_coord.r;
+  const auto& cell {model::cells[upper_coord.cell()]};
+  Position r = upper_coord.r();
   r -= cell->translation_;
   if (!cell->rotation_.empty()) {
     r = r.rotate(cell->rotation_);
   }
-  p.r_local() = lat.get_local_position(r, coord.lattice_i);
+  p.r_local() = lat.get_local_position(r, coord.lattice_index());
 
-  if (!lat.are_valid_indices(coord.lattice_i)) {
+  if (!lat.are_valid_indices(coord.lattice_index())) {
     // The particle is outside the lattice.  Search for it from the base coords.
     p.n_coord() = 1;
     bool found = exhaustive_find_cell(p);
@@ -337,7 +337,7 @@ void cross_lattice(GeometryState& p, const BoundaryInfo& boundary, bool verbose)
 
   } else {
     // Find cell in next lattice element.
-    p.lowest_coord().universe = lat[coord.lattice_i];
+    p.lowest_coord().universe() = lat[coord.lattice_index()];
     bool found = exhaustive_find_cell(p);
 
     if (!found) {
@@ -367,9 +367,9 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
   // Loop over each coordinate level.
   for (int i = 0; i < p.n_coord(); i++) {
     const auto& coord {p.coord(i)};
-    const Position& r {coord.r};
-    const Direction& u {coord.u};
-    Cell& c {*model::cells[coord.cell]};
+    const Position& r {coord.r()};
+    const Direction& u {coord.u()};
+    Cell& c {*model::cells[coord.cell()]};
 
     // Find the oncoming surface in this cell and the distance to it.
     auto surface_distance = c.distance(r, u, p.surface(), &p);
@@ -377,24 +377,24 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
     level_surf_cross = surface_distance.second;
 
     // Find the distance to the next lattice tile crossing.
-    if (coord.lattice != C_NONE) {
-      auto& lat {*model::lattices[coord.lattice]};
+    if (coord.lattice() != C_NONE) {
+      auto& lat {*model::lattices[coord.lattice()]};
       // TODO: refactor so both lattice use the same position argument (which
       // also means the lat.type attribute can be removed)
       std::pair<double, array<int, 3>> lattice_distance;
       switch (lat.type_) {
       case LatticeType::rect:
-        lattice_distance = lat.distance(r, u, coord.lattice_i);
+        lattice_distance = lat.distance(r, u, coord.lattice_index());
         break;
       case LatticeType::hex:
-        auto& cell_above {model::cells[p.coord(i - 1).cell]};
-        Position r_hex {p.coord(i - 1).r};
+        auto& cell_above {model::cells[p.coord(i - 1).cell()]};
+        Position r_hex {p.coord(i - 1).r()};
         r_hex -= cell_above->translation_;
-        if (coord.rotated) {
+        if (coord.rotated()) {
           r_hex = r_hex.rotate(cell_above->rotation_);
         }
-        r_hex.z = coord.r.z;
-        lattice_distance = lat.distance(r_hex, u, coord.lattice_i);
+        r_hex.z = coord.r().z;
+        lattice_distance = lat.distance(r_hex, u, coord.lattice_index());
         break;
       }
       d_lat = lattice_distance.first;
@@ -468,7 +468,7 @@ extern "C" int openmc_find_cell(
     return OPENMC_E_GEOMETRY;
   }
 
-  *index = geom_state.lowest_coord().cell;
+  *index = geom_state.lowest_coord().cell();
   *instance = geom_state.cell_instance();
   return 0;
 }
diff --git a/src/mesh.cpp b/src/mesh.cpp
index becbbd3fd..1b2c71eb2 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -368,11 +368,11 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 
           // Set birth cell attribute
           if (p.cell_born() == C_NONE)
-            p.cell_born() = p.lowest_coord().cell;
+            p.cell_born() = p.lowest_coord().cell();
 
           // Initialize last cells from current cell
           for (int j = 0; j < p.n_coord(); ++j) {
-            p.cell_last(j) = p.coord(j).cell;
+            p.cell_last(j) = p.coord(j).cell();
           }
           p.n_coord_last() = p.n_coord();
 
@@ -411,7 +411,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size,
 
             // cross next geometric surface
             for (int j = 0; j < p.n_coord(); ++j) {
-              p.cell_last(j) = p.coord(j).cell;
+              p.cell_last(j) = p.coord(j).cell();
             }
             p.n_coord_last() = p.n_coord();
 
diff --git a/src/output.cpp b/src/output.cpp
index baaf682a1..0316e313e 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -177,25 +177,26 @@ void print_particle(Particle& p)
   for (auto i = 0; i < p.n_coord(); i++) {
     fmt::print("  Level {}\n", i);
 
-    if (p.coord(i).cell != C_NONE) {
-      const Cell& c {*model::cells[p.coord(i).cell]};
+    if (p.coord(i).cell() != C_NONE) {
+      const Cell& c {*model::cells[p.coord(i).cell()]};
       fmt::print("    Cell             = {}\n", c.id_);
     }
 
-    if (p.coord(i).universe != C_NONE) {
-      const Universe& u {*model::universes[p.coord(i).universe]};
+    if (p.coord(i).universe() != C_NONE) {
+      const Universe& u {*model::universes[p.coord(i).universe()]};
       fmt::print("    Universe         = {}\n", u.id_);
     }
 
-    if (p.coord(i).lattice != C_NONE) {
-      const Lattice& lat {*model::lattices[p.coord(i).lattice]};
+    if (p.coord(i).lattice() != C_NONE) {
+      const Lattice& lat {*model::lattices[p.coord(i).lattice()]};
       fmt::print("    Lattice          = {}\n", lat.id_);
-      fmt::print("    Lattice position = ({},{},{})\n", p.coord(i).lattice_i[0],
-        p.coord(i).lattice_i[1], p.coord(i).lattice_i[2]);
+      fmt::print("    Lattice position = ({},{},{})\n",
+        p.coord(i).lattice_index()[0], p.coord(i).lattice_index()[1],
+        p.coord(i).lattice_index()[2]);
     }
 
-    fmt::print("    r = {}\n", p.coord(i).r);
-    fmt::print("    u = {}\n", p.coord(i).u);
+    fmt::print("    r = {}\n", p.coord(i).r());
+    fmt::print("    u = {}\n", p.coord(i).u());
   }
 
   // Display miscellaneous info.
diff --git a/src/particle.cpp b/src/particle.cpp
index 6cfdf2175..4758cded5 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -176,7 +176,7 @@ void Particle::event_calculate_xs()
   // If the cell hasn't been determined based on the particle's location,
   // initiate a search for the current cell. This generally happens at the
   // beginning of the history and again for any secondary particles
-  if (lowest_coord().cell == C_NONE) {
+  if (lowest_coord().cell() == C_NONE) {
     if (!exhaustive_find_cell(*this)) {
       mark_as_lost(
         "Could not find the cell containing particle " + std::to_string(id()));
@@ -185,11 +185,11 @@ void Particle::event_calculate_xs()
 
     // Set birth cell attribute
     if (cell_born() == C_NONE)
-      cell_born() = lowest_coord().cell;
+      cell_born() = lowest_coord().cell();
 
     // Initialize last cells from current cell
     for (int j = 0; j < n_coord(); ++j) {
-      cell_last(j) = coord(j).cell;
+      cell_last(j) = coord(j).cell();
     }
     n_coord_last() = n_coord();
   }
@@ -248,7 +248,7 @@ void Particle::event_advance()
   // Short-term solution until the surface source is revised and we can use
   // this->move_distance(distance)
   for (int j = 0; j < n_coord(); ++j) {
-    coord(j).r += distance * coord(j).u;
+    coord(j).r() += distance * coord(j).u();
   }
   double dt = distance / this->speed();
   this->time() += dt;
@@ -293,7 +293,7 @@ void Particle::event_cross_surface()
 {
   // Saving previous cell data
   for (int j = 0; j < n_coord(); ++j) {
-    cell_last(j) = coord(j).cell;
+    cell_last(j) = coord(j).cell();
   }
   n_coord_last() = n_coord();
 
@@ -393,14 +393,14 @@ void Particle::event_collide()
   // Set all directions to base level -- right now, after a collision, only
   // the base level directions are changed
   for (int j = 0; j < n_coord() - 1; ++j) {
-    if (coord(j + 1).rotated) {
+    if (coord(j + 1).rotated()) {
       // If next level is rotated, apply rotation matrix
-      const auto& m {model::cells[coord(j).cell]->rotation_};
-      const auto& u {coord(j).u};
-      coord(j + 1).u = u.rotate(m);
+      const auto& m {model::cells[coord(j).cell()]->rotation_};
+      const auto& u {coord(j).u()};
+      coord(j + 1).u() = u.rotate(m);
     } else {
       // Otherwise, copy this level's direction
-      coord(j + 1).u = coord(j).u;
+      coord(j + 1).u() = coord(j).u();
     }
   }
 
@@ -445,7 +445,7 @@ void Particle::event_revive_from_secondary()
       // Since the birth cell of the particle has not been set we
       // have to determine it before the energy of the secondary particle can be
       // removed from the pulse-height of this cell.
-      if (lowest_coord().cell == C_NONE) {
+      if (lowest_coord().cell() == C_NONE) {
         bool verbose = settings::verbosity >= 10 || trace();
         if (!exhaustive_find_cell(*this, verbose)) {
           mark_as_lost("Could not find the cell containing particle " +
@@ -454,11 +454,11 @@ void Particle::event_revive_from_secondary()
         }
         // Set birth cell attribute
         if (cell_born() == C_NONE)
-          cell_born() = lowest_coord().cell;
+          cell_born() = lowest_coord().cell();
 
         // Initialize last cells from current cell
         for (int j = 0; j < n_coord(); ++j) {
-          cell_last(j) = coord(j).cell;
+          cell_last(j) = coord(j).cell();
         }
         n_coord_last() = n_coord();
       }
@@ -516,7 +516,7 @@ void Particle::pht_collision_energy()
 
   // determine index of cell in pulse_height_cells
   auto it = std::find(model::pulse_height_cells.begin(),
-    model::pulse_height_cells.end(), lowest_coord().cell);
+    model::pulse_height_cells.end(), lowest_coord().cell());
 
   if (it != model::pulse_height_cells.end()) {
     int index = std::distance(model::pulse_height_cells.begin(), it);
@@ -571,13 +571,13 @@ void Particle::cross_surface(const Surface& surf)
   // in DAGMC, we know what the next cell should be
   if (surf.geom_type() == GeometryType::DAG) {
     int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
-                       lowest_coord().universe) -
+                       lowest_coord().universe()) -
                      1;
     // save material and temp
     material_last() = material();
     sqrtkT_last() = sqrtkT();
     // set new cell value
-    lowest_coord().cell = i_cell;
+    lowest_coord().cell() = i_cell;
     auto& cell = model::cells[i_cell];
 
     cell_instance() = 0;
@@ -681,7 +681,7 @@ void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
   u() = new_u;
 
   // Reassign particle's cell and surface
-  coord(0).cell = cell_last(0);
+  coord(0).cell() = cell_last(0);
   surface() = -surface();
 
   // If a reflective surface is coincident with a lattice or universe
@@ -942,7 +942,7 @@ void add_surf_source_to_bank(Particle& p, const Surface& surf)
     // Check if the cell of interest has been entered
     bool entered = false;
     for (int i = 0; i < p.n_coord(); ++i) {
-      if (p.coord(i).cell == cell_idx) {
+      if (p.coord(i).cell() == cell_idx) {
         entered = true;
       }
     }
diff --git a/src/particle_data.cpp b/src/particle_data.cpp
index fef8359a3..aff2b68ba 100644
--- a/src/particle_data.cpp
+++ b/src/particle_data.cpp
@@ -32,20 +32,20 @@ void GeometryState::mark_as_lost(const std::stringstream& message)
 
 void LocalCoord::rotate(const vector<double>& rotation)
 {
-  r = r.rotate(rotation);
-  u = u.rotate(rotation);
-  rotated = true;
+  r_ = r_.rotate(rotation);
+  u_ = u_.rotate(rotation);
+  rotated_ = true;
 }
 
 void LocalCoord::reset()
 {
-  cell = C_NONE;
-  universe = C_NONE;
-  lattice = C_NONE;
-  lattice_i[0] = 0;
-  lattice_i[1] = 0;
-  lattice_i[2] = 0;
-  rotated = false;
+  cell_ = C_NONE;
+  universe_ = C_NONE;
+  lattice_ = C_NONE;
+  lattice_index_[0] = 0;
+  lattice_index_[1] = 0;
+  lattice_index_[2] = 0;
+  rotated_ = false;
 }
 
 GeometryState::GeometryState()
@@ -64,7 +64,7 @@ void GeometryState::advance_to_boundary_from_void()
 
   for (auto c_i : root_universe->cells_) {
     auto dist =
-      model::cells.at(c_i)->distance(root_coord.r, root_coord.u, 0, this);
+      model::cells.at(c_i)->distance(root_coord.r(), root_coord.u(), 0, this);
     if (dist.first < boundary().distance) {
       boundary().distance = dist.first;
       boundary().surface = dist.second;
@@ -86,7 +86,7 @@ void GeometryState::advance_to_boundary_from_void()
 void GeometryState::move_distance(double length)
 {
   for (int j = 0; j < n_coord(); ++j) {
-    coord(j).r += length * coord(j).u;
+    coord(j).r() += length * coord(j).u();
   }
 }
 
@@ -123,7 +123,7 @@ TrackState ParticleData::get_track_state() const
   state.E = this->E();
   state.time = this->time();
   state.wgt = this->wgt();
-  state.cell_id = model::cells[this->lowest_coord().cell]->id_;
+  state.cell_id = model::cells[this->lowest_coord().cell()]->id_;
   state.cell_instance = this->cell_instance();
   if (this->material() != MATERIAL_VOID) {
     state.material_id = model::materials[material()]->id_;
diff --git a/src/plot.cpp b/src/plot.cpp
index dbc25b21e..db10c7600 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -54,14 +54,14 @@ void IdData::set_value(size_t y, size_t x, const GeometryState& p, int level)
     data_(y, x, 0) = NOT_FOUND;
     data_(y, x, 1) = NOT_FOUND;
   } else {
-    data_(y, x, 0) = model::cells.at(p.coord(level).cell)->id_;
+    data_(y, x, 0) = model::cells.at(p.coord(level).cell())->id_;
     data_(y, x, 1) = level == p.n_coord() - 1
                        ? p.cell_instance()
                        : cell_instance_at_level(p, level);
   }
 
   // set material data
-  Cell* c = model::cells.at(p.lowest_coord().cell).get();
+  Cell* c = model::cells.at(p.lowest_coord().cell()).get();
   if (p.material() == MATERIAL_VOID) {
     data_(y, x, 2) = MATERIAL_VOID;
     return;
@@ -83,7 +83,7 @@ PropertyData::PropertyData(size_t h_res, size_t v_res)
 void PropertyData::set_value(
   size_t y, size_t x, const GeometryState& p, int level)
 {
-  Cell* c = model::cells.at(p.lowest_coord().cell).get();
+  Cell* c = model::cells.at(p.lowest_coord().cell()).get();
   data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
   if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
     Material* m = model::materials.at(p.material()).get();
@@ -1692,7 +1692,7 @@ void Ray::trace()
 
     // Advance particle, prepare for next intersection
     for (int lev = 0; lev < n_coord(); ++lev) {
-      coord(lev).r += boundary().distance * coord(lev).u;
+      coord(lev).r() += boundary().distance * coord(lev).u();
     }
     surface() = boundary().surface;
     n_coord_last() = n_coord();
@@ -1743,7 +1743,7 @@ void ProjectionRay::on_intersection()
   line_segments_.emplace_back(
     plot_.color_by_ == PlottableInterface::PlotColorBy::mats
       ? material()
-      : lowest_coord().cell,
+      : lowest_coord().cell(),
     traversal_distance_, boundary().surface_index());
 }
 
@@ -1752,7 +1752,7 @@ void PhongRay::on_intersection()
   // Check if we hit an opaque material or cell
   int hit_id = plot_.color_by_ == PlottableInterface::PlotColorBy::mats
                  ? material()
-                 : lowest_coord().cell;
+                 : lowest_coord().cell();
 
   // If we are reflected and have advanced beyond the camera,
   // the ray is done. This is checked here because we should
@@ -1798,8 +1798,8 @@ void PhongRay::on_intersection()
     // Need to apply translations to find the normal vector in
     // the base level universe's coordinate system.
     for (int lev = n_coord() - 2; lev >= 0; --lev) {
-      if (coord(lev + 1).rotated) {
-        const Cell& c {*model::cells[coord(lev).cell]};
+      if (coord(lev + 1).rotated()) {
+        const Cell& c {*model::cells[coord(lev).cell()]};
         normal = normal.inverse_rotate(c.rotation_);
       }
     }
diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp
index 4a5b43036..b04b72981 100644
--- a/src/random_ray/flat_source_domain.cpp
+++ b/src/random_ray/flat_source_domain.cpp
@@ -794,7 +794,7 @@ void FlatSourceDomain::output_to_vtk() const
             continue;
           }
 
-          int i_cell = p.lowest_coord().cell;
+          int i_cell = p.lowest_coord().cell();
           int64_t sr = source_region_offsets_[i_cell] + p.cell_instance();
           if (RandomRay::mesh_subdivision_enabled_) {
             int mesh_idx = base_source_regions_.mesh(sr);
@@ -1070,7 +1070,7 @@ void FlatSourceDomain::convert_external_sources()
                                 "point source at {}",
           sp->r()));
       }
-      int i_cell = gs.lowest_coord().cell;
+      int i_cell = gs.lowest_coord().cell();
       int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance();
 
       if (RandomRay::mesh_subdivision_enabled_) {
@@ -1475,7 +1475,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle(
   gs.r() = r + TINY_BIT * u;
   gs.u() = {1.0, 0.0, 0.0};
   exhaustive_find_cell(gs);
-  int gs_i_cell = gs.lowest_coord().cell;
+  int gs_i_cell = gs.lowest_coord().cell();
   int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance();
   if (sr_found != sr) {
     discovered_source_regions_.unlock(sr_key);
diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp
index 36eab0628..2eca5b0ab 100644
--- a/src/random_ray/random_ray.cpp
+++ b/src/random_ray/random_ray.cpp
@@ -330,14 +330,14 @@ void RandomRay::event_advance_ray()
 
   // Advance particle
   for (int j = 0; j < n_coord(); ++j) {
-    coord(j).r += distance * coord(j).u;
+    coord(j).r() += distance * coord(j).u();
   }
 }
 
 void RandomRay::attenuate_flux(double distance, bool is_active, double offset)
 {
   // Determine source region index etc.
-  int i_cell = lowest_coord().cell;
+  int i_cell = lowest_coord().cell();
 
   // The base source region is the spatial region index
   int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
@@ -799,7 +799,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
   this->from_source(&site);
 
   // Locate ray
-  if (lowest_coord().cell == C_NONE) {
+  if (lowest_coord().cell() == C_NONE) {
     if (!exhaustive_find_cell(*this)) {
       this->mark_as_lost(
         "Could not find the cell containing particle " + std::to_string(id()));
@@ -807,12 +807,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain)
 
     // Set birth cell attribute
     if (cell_born() == C_NONE)
-      cell_born() = lowest_coord().cell;
+      cell_born() = lowest_coord().cell();
   }
 
   // Initialize ray's starting angular flux to starting location's isotropic
   // source
-  int i_cell = lowest_coord().cell;
+  int i_cell = lowest_coord().cell();
   int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance();
 
   SourceRegionHandle srh;
diff --git a/src/source.cpp b/src/source.cpp
index 65873f5a7..ae3daa404 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -246,9 +246,10 @@ bool Source::satisfies_spatial_constraints(Position r) const
       }
     } else {
       for (int i = 0; i < geom_state.n_coord(); i++) {
-        auto id = (domain_type_ == DomainType::CELL)
-                    ? model::cells[geom_state.coord(i).cell]->id_
-                    : model::universes[geom_state.coord(i).universe]->id_;
+        auto id =
+          (domain_type_ == DomainType::CELL)
+            ? model::cells[geom_state.coord(i).cell()].get()->id_
+            : model::universes[geom_state.coord(i).universe()].get()->id_;
         if ((accepted = contains(domain_ids_, id)))
           break;
       }
diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp
index b545801ea..7a6394956 100644
--- a/src/tallies/filter_cell.cpp
+++ b/src/tallies/filter_cell.cpp
@@ -49,7 +49,7 @@ void CellFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
   for (int i = 0; i < p.n_coord(); i++) {
-    auto search = map_.find(p.coord(i).cell);
+    auto search = map_.find(p.coord(i).cell());
     if (search != map_.end()) {
       match.bins_.push_back(search->second);
       match.weights_.push_back(1.0);
diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp
index 0634175dd..316a758d1 100644
--- a/src/tallies/filter_cell_instance.cpp
+++ b/src/tallies/filter_cell_instance.cpp
@@ -73,7 +73,7 @@ void CellInstanceFilter::set_cell_instances(span<CellInstance> instances)
 void CellInstanceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  int64_t index_cell = p.lowest_coord().cell;
+  int64_t index_cell = p.lowest_coord().cell();
   int64_t instance = p.cell_instance();
 
   if (cells_.count(index_cell) > 0) {
@@ -89,7 +89,7 @@ void CellInstanceFilter::get_all_bins(
     return;
 
   for (int i = 0; i < p.n_coord() - 1; i++) {
-    int64_t index_cell = p.coord(i).cell;
+    int64_t index_cell = p.coord(i).cell();
     // if this cell isn't used on the filter, move on
     if (cells_.count(index_cell) == 0)
       continue;
diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp
index 821e843da..02e88dde2 100644
--- a/src/tallies/filter_distribcell.cpp
+++ b/src/tallies/filter_distribcell.cpp
@@ -43,18 +43,18 @@ void DistribcellFilter::get_all_bins(
   int offset = 0;
   auto distribcell_index = model::cells[cell_]->distribcell_index_;
   for (int i = 0; i < p.n_coord(); i++) {
-    auto& c {*model::cells[p.coord(i).cell]};
+    auto& c {*model::cells[p.coord(i).cell()]};
     if (c.type_ == Fill::UNIVERSE) {
       offset += c.offset_[distribcell_index];
     } else if (c.type_ == Fill::LATTICE) {
-      auto& lat {*model::lattices[p.coord(i + 1).lattice]};
-      const auto& i_xyz {p.coord(i + 1).lattice_i};
+      auto& lat {*model::lattices[p.coord(i + 1).lattice()]};
+      const auto& i_xyz {p.coord(i + 1).lattice_index()};
       if (lat.are_valid_indices(i_xyz)) {
         offset +=
           lat.offset(distribcell_index, i_xyz) + c.offset_[distribcell_index];
       }
     }
-    if (cell_ == p.coord(i).cell) {
+    if (cell_ == p.coord(i).cell()) {
       match.bins_.push_back(offset);
       match.weights_.push_back(1.0);
       return;
diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp
index 48114cfad..f4b22decd 100644
--- a/src/tallies/filter_universe.cpp
+++ b/src/tallies/filter_universe.cpp
@@ -48,7 +48,7 @@ void UniverseFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
   for (int i = 0; i < p.n_coord(); i++) {
-    auto search = map_.find(p.coord(i).universe);
+    auto search = map_.find(p.coord(i).universe());
     if (search != map_.end()) {
       match.bins_.push_back(search->second);
       match.weights_.push_back(1.0);
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 71b103cfb..e73fb90f3 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2615,7 +2615,7 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
 
   // Save original cell/energy information
   int orig_n_coord = p.n_coord();
-  int orig_cell = p.coord(0).cell;
+  int orig_cell = p.coord(0).cell();
   double orig_E_last = p.E_last();
 
   for (auto i_tally : tallies) {
@@ -2634,7 +2634,7 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
 
           // Temporarily change cell of particle
           p.n_coord() = 1;
-          p.coord(0).cell = cell_id;
+          p.coord(0).cell() = cell_id;
 
           // Determine index of cell in model::pulse_height_cells
           auto it = std::find(model::pulse_height_cells.begin(),
@@ -2674,7 +2674,7 @@ void score_pulse_height_tally(Particle& p, const vector<int>& tallies)
     }
     // Restore cell/energy
     p.n_coord() = orig_n_coord;
-    p.coord(0).cell = orig_cell;
+    p.coord(0).cell() = orig_cell;
     p.E_last() = orig_E_last;
   }
 }
diff --git a/src/universe.cpp b/src/universe.cpp
index b4ef6b4b2..78ddd54d4 100644
--- a/src/universe.cpp
+++ b/src/universe.cpp
@@ -45,14 +45,14 @@ bool Universe::find_cell(GeometryState& p) const
   Position r {p.r_local()};
   Position u {p.u_local()};
   auto surf = p.surface();
-  int32_t i_univ = p.lowest_coord().universe;
+  int32_t i_univ = p.lowest_coord().universe();
 
   for (auto i_cell : cells) {
     if (model::cells[i_cell]->universe_ != i_univ)
       continue;
     // Check if this cell contains the particle
     if (model::cells[i_cell]->contains(r, u, surf)) {
-      p.lowest_coord().cell = i_cell;
+      p.lowest_coord().cell() = i_cell;
       return true;
     }
   }
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index 8b5c27f14..1deffb804 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -179,7 +179,7 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         } else if (domain_type_ == TallyDomain::CELL) {
           for (int level = 0; level < p.n_coord(); ++level) {
             for (int i_domain = 0; i_domain < n; i_domain++) {
-              if (model::cells[p.coord(level).cell]->id_ ==
+              if (model::cells[p.coord(level).cell()]->id_ ==
                   domain_ids_[i_domain]) {
                 this->check_hit(
                   p.material(), indices[i_domain], hits[i_domain]);
@@ -190,7 +190,7 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         } else if (domain_type_ == TallyDomain::UNIVERSE) {
           for (int level = 0; level < p.n_coord(); ++level) {
             for (int i_domain = 0; i_domain < n; ++i_domain) {
-              if (model::universes[p.coord(level).universe]->id_ ==
+              if (model::universes[p.coord(level).universe()]->id_ ==
                   domain_ids_[i_domain]) {
                 check_hit(p.material(), indices[i_domain], hits[i_domain]);
                 break;