diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 32641120..93668954 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,5 +1,5 @@
 # ========================================================================================
-#  (C) (or copyright) 2024. Triad National Security, LLC. All rights reserved.
+#  (C) (or copyright) 2024-2025. Triad National Security, LLC. All rights reserved.
 #
 #  This program was produced under U.S. Government contract 89233218CNA000001 for Los
 #  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -47,7 +47,32 @@ jobs:
           VERBOSE=1 CFM=clang-format ./style/format.sh
           git diff --exit-code --ignore-submodules
 
-  cpu:
+  unit-tests:
+    if: >
+      ${{ !contains(github.event.pull_request.title, 'Draft:') &&
+          !contains(github.event.pull_request.title, 'WIP:') }}
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v3
+        with:
+          submodules: recursive
+      - name: Install dependencies
+        run: |
+          sudo apt-get update -qq
+          sudo apt-get install -qq --no-install-recommends tzdata
+          sudo apt-get install -qq git
+          sudo apt-get install -qq make cmake g++
+          sudo apt-get install -qq libopenmpi-dev libhdf5-openmpi-dev
+      - name: Build and run unit tests
+        run: |
+          export MAKE_PROGRAM=make
+          mkdir -p build
+          cd build
+          cmake --preset=cpu-release ../
+          make -j 4
+          ctest -L unit --output-on-failure
+
+  regression-tests:
     if: >
       ${{ !contains(github.event.pull_request.title, 'Draft:') &&
           !contains(github.event.pull_request.title, 'WIP:') }}
@@ -65,7 +90,7 @@ jobs:
           sudo apt-get install -qq libopenmpi-dev libhdf5-openmpi-dev
           sudo apt-get install -qq openssh-client
           sudo apt-get install -qq python3 python3-numpy python3-h5py python3-matplotlib
-      - name: Run CPU tests
+      - name: Run regression tests
         run: |
           export MAKE_PROGRAM=make
           cd tst
diff --git a/.gitmodules b/.gitmodules
index fd85689b..4788484d 100644
--- a/.gitmodules
+++ b/.gitmodules
@@ -13,3 +13,6 @@
 [submodule "external/jaybenne"]
 	path = external/jaybenne
 	url = https://github.com/lanl/jaybenne.git
+[submodule "external/Catch2"]
+	path = external/Catch2
+	url = https://github.com/catchorg/Catch2.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 2e10a8bf..90bd8e1d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -244,3 +244,49 @@ add_subdirectory(external/jaybenne)
 # Add artemis
 message("\nConfiguring src")
 add_subdirectory(src)
+
+# Configure Catch2 for unit testing
+if (NOT TARGET Catch2::Catch2WithMain)
+  message("Configuring Catch2")
+  set(CATCH_BUILD_TESTING OFF CACHE BOOL "" FORCE)
+  set(CATCH_INSTALL_DOCS OFF CACHE BOOL "" FORCE)
+  set(CATCH_INSTALL_EXTRAS OFF CACHE BOOL "" FORCE)
+  add_subdirectory(external/Catch2 Catch2)
+else()
+  message(STATUS "Catch2 already configured in this build, skipping.")
+endif()
+
+# Enable testing
+enable_testing()
+
+# Add unit tests
+message("\nConfiguring unit tests")
+add_subdirectory(tst/unit)
+
+# Add regression tests
+message("\nConfiguring regression tests")
+find_package(Python3 COMPONENTS Interpreter)
+if(Python3_FOUND)
+  # Determine which test suite to run based on GPU support
+  if(ARTEMIS_ENABLE_CUDA OR ARTEMIS_ENABLE_HIP)
+    set(REGRESSION_SUITE "gpu.suite")
+  else()
+    set(REGRESSION_SUITE "regression.suite")
+  endif()
+  
+  # Add the regression test
+  add_test(
+    NAME regression
+    COMMAND ${Python3_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/tst/run_tests.py
+            ${REGRESSION_SUITE} --cwd
+    WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+  )
+  
+  # Set label for regression test
+  set_tests_properties(regression PROPERTIES LABELS "regression")
+  
+  message(STATUS "Regression test suite: ${REGRESSION_SUITE}")
+  message(STATUS "Regression test command: ${Python3_EXECUTABLE} run_tests.py ${REGRESSION_SUITE} --cwd")
+else()
+  message(WARNING "Python3 not found. Regression tests will not be available via CTest.")
+endif()
diff --git a/external/Catch2 b/external/Catch2
new file mode 160000
index 00000000..97091636
--- /dev/null
+++ b/external/Catch2
@@ -0,0 +1 @@
+Subproject commit 97091636d012baa87711e7dbbd4cd778b96a211b
diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt
new file mode 100644
index 00000000..2704e440
--- /dev/null
+++ b/tst/unit/CMakeLists.txt
@@ -0,0 +1,47 @@
+# ========================================================================================
+#  (C) (or copyright) 2023-2025. Triad National Security, LLC. All rights reserved.
+#
+#  This program was produced under U.S. Government contract 89233218CNA000001 for Los
+#  Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+#  for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+#  in the program are reserved by Triad National Security, LLC, and the U.S. Department
+#  of Energy/National Nuclear Security Administration. The Government is granted for
+#  itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+#  license in this material to reproduce, prepare derivative works, distribute copies to
+#  the public, perform publicly and display publicly, and to permit others to do so.
+# ========================================================================================
+
+# This file was created in part or in whole by generative AI
+
+# Unit tests using Catch2
+
+# Helper function to add unit tests
+function(add_unit_test TEST_NAME TEST_SOURCE)
+  add_executable(${TEST_NAME} ${TEST_SOURCE})
+  
+  target_link_libraries(${TEST_NAME}
+    PRIVATE
+      Catch2::Catch2WithMain
+      parthenon
+      singularity-eos::singularity-eos
+      Kokkos::kokkos
+  )
+  
+  target_include_directories(${TEST_NAME}
+    PRIVATE
+      ${CMAKE_CURRENT_SOURCE_DIR}/../../src
+      ${ARTEMIS_SINGULARITY_INCLUDE_PATHS}
+  )
+  
+  # Discover and add tests to CTest
+  catch_discover_tests(${TEST_NAME}
+    PROPERTIES LABELS "unit")
+endfunction()
+
+# Include Catch2's module for test discovery
+include(Catch)
+
+# Add individual unit tests here
+add_unit_test(test_ideal_gas test_ideal_gas.cpp)
+add_unit_test(test_ideal_hhe test_ideal_hhe.cpp)
+add_unit_test(test_geometry test_geometry.cpp)
diff --git a/tst/unit/test_geometry.cpp b/tst/unit/test_geometry.cpp
new file mode 100644
index 00000000..b9276af8
--- /dev/null
+++ b/tst/unit/test_geometry.cpp
@@ -0,0 +1,606 @@
+//========================================================================================
+// (C) (or copyright) 2025. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+//  license in this material to reproduce, prepare derivative works, distribute copies to
+//  the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+// This file was created in part or in whole by generative AI
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <cmath>
+
+#include "artemis.hpp"
+#include "geometry/axisymmetric.hpp"
+#include "geometry/cylindrical.hpp"
+#include "geometry/geometry.hpp"
+#include "geometry/spherical.hpp"
+
+using Catch::Approx;
+using namespace geometry;
+
+// Helper to create a simple BBox for testing
+struct TestBBox {
+  Real x1[2];
+  Real x2[2];
+  Real x3[2];
+};
+
+TEST_CASE("Spherical 3D coordinate system", "[geometry][spherical]") {
+  // Create a test cell in spherical coordinates
+  // r: [1, 2], theta: [pi/4, pi/2], phi: [0, pi/2]
+
+  SECTION("Volume calculation") {
+    // Analytical volume: ∫∫∫ r² sin(θ) dr dθ dφ
+    // = [r³/3]₁² × [-cos(θ)]_{π/4}^{π/2} × [φ]₀^{π/2}
+    // = (8/3 - 1/3) × (0 - (-√2/2)) × π/2
+    // = 7/3 × √2/2 × π/2
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = M_PI / 4.0;
+    coords.bnds.x2[1] = M_PI / 2.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI / 2.0;
+
+    const Real expected_volume = (7.0 / 3.0) * (std::sqrt(2.0) / 2.0) * (M_PI / 2.0);
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Radial centroid (x1v)") {
+    // For spherical coordinates, the radial centroid is volume-weighted:
+    // <r> = ∫ r³ dr / ∫ r² dr
+    // For r ∈ [1, 2]: <r> = [r⁴/4]₁² / [r³/3]₁²
+    // = (16/4 - 1/4) / (8/3 - 1/3) = (15/4) / (7/3) = 45/28
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0 * M_PI;
+
+    const Real expected_x1v = 45.0 / 28.0;
+    const Real computed_x1v = coords.x1v();
+
+    REQUIRE(computed_x1v == Approx(expected_x1v).epsilon(1e-10));
+  }
+
+  SECTION("Theta centroid (x2v)") {
+    // Volume-weighted theta: ∫ θ sin(θ) dθ / ∫ sin(θ) dθ
+    // = [sin(θ) - θ cos(θ)] / [-cos(θ)]
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = M_PI / 4.0;
+    coords.bnds.x2[1] = 3.0 * M_PI / 4.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    const Real computed_x2v = coords.x2v();
+
+    // Should be close to pi/2 by symmetry
+    REQUIRE(computed_x2v == Approx(M_PI / 2.0).epsilon(1e-6));
+  }
+
+  SECTION("X1 face area") {
+    // Area of radial face: ∫∫ r² sin(θ) dθ dφ
+    // At r = 1.5, θ ∈ [0, π/2], φ ∈ [0, π]
+    // = 1.5² × [-cos(θ)]₀^{π/2} × π
+    // = 2.25 × 1 × π
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI / 2.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    const Real r = 1.5;
+    const Real expected_area = r * r * 1.0 * M_PI;
+    const Real computed_area = coords.AreaX1(r);
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("X2 face area") {
+    // Area of theta face: ∫∫ r sin(θ) dr dφ
+    // For r ∈ [1, 2], θ = π/4, φ ∈ [0, π]
+    // = [r²/2]₁² × sin(π/4) × π
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = M_PI / 6.0;
+    coords.bnds.x2[1] = M_PI / 3.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    const Real theta = M_PI / 4.0;
+    const Real expected_area = 1.5 * 1.0 * std::sin(theta) * M_PI;
+    const Real computed_area = coords.AreaX2(theta);
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("X3 face area") {
+    // Area of phi face: ∫∫ r dr dθ = [r²/2] × Δθ
+
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI / 2.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    const Real expected_area = 1.5 * 1.0 * (M_PI / 2.0);
+    const Real computed_area = coords.AreaX3(0.0); // phi value doesn't matter
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+}
+
+TEST_CASE("Spherical 2D coordinate system", "[geometry][spherical]") {
+  // 2D spherical (r, theta) with axisymmetry in phi
+
+  SECTION("Volume calculation") {
+    // Volume: ∫∫ r² sin(θ) dr dθ (no explicit phi integration in 2D)
+    // The 2D implementation doesn't include the 2π factor
+
+    Coords<Coordinates::spherical2D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = M_PI / 4.0;
+    coords.bnds.x2[1] = M_PI / 2.0;
+
+    // dx1 * rfac * |cos(x2[0]) - cos(x2[1])|
+    const Real dx1 = 1.0;
+    const Real rfac = (1.0 + 2.0 + 4.0) / 3.0; // (r₀² + r₀r₁ + r₁²) / 3
+    const Real dx2 = std::abs(std::cos(M_PI / 4.0) - std::cos(M_PI / 2.0));
+    const Real expected_volume = rfac * dx1 * dx2;
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Centroid calculations") {
+    Coords<Coordinates::spherical2D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+
+    const Real expected_x1v = 45.0 / 28.0;
+    REQUIRE(coords.x1v() == Approx(expected_x1v).epsilon(1e-10));
+  }
+}
+
+TEST_CASE("Cylindrical coordinate system", "[geometry][cylindrical]") {
+  // Cylindrical coordinates: (R, phi, z)
+
+  SECTION("Volume calculation") {
+    // Volume: ∫∫∫ R dR dφ dz = [R²/2] × Δφ × Δz
+    // For R ∈ [1, 2], φ ∈ [0, π], z ∈ [0, 3]
+
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 3.0;
+
+    // <R> = (R₀ + R₁) / 2 = 1.5
+    // Volume = <R> × ΔR × Δφ × Δz = 1.5 × 1.0 × π × 3.0
+    const Real expected_volume = 1.5 * 1.0 * M_PI * 3.0;
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Radial centroid (x1v)") {
+    // For cylindrical, <R> = ∫ R² dR / ∫ R dR
+    // = [R³/3]₁² / [R²/2]₁²
+    // = (8/3 - 1/3) / (4/2 - 1/2) = (7/3) / (3/2) = 14/9
+
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 2.0 * M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 1.0;
+
+    const Real expected_x1v = 14.0 / 9.0;
+    const Real computed_x1v = coords.x1v();
+
+    REQUIRE(computed_x1v == Approx(expected_x1v).epsilon(1e-10));
+  }
+
+  SECTION("X1 face area (radial face)") {
+    // Area: ∫∫ R dφ dz = R × Δφ × Δz
+
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI / 2.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0;
+
+    const Real R = 1.5;
+    const Real expected_area = R * (M_PI / 2.0) * 2.0;
+    const Real computed_area = coords.AreaX1(R);
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("X3 face area (z face)") {
+    // Area: ∫∫ R dR dφ = [R²/2] × Δφ = <R> × ΔR × Δφ
+
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 1.0;
+
+    const Real expected_area = 1.5 * 1.0 * M_PI;
+    const Real computed_area = coords.AreaX3(0.5); // z value doesn't matter
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("Scale factor hx2v") {
+    // In cylindrical coords, h_phi = R
+
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 2.0;
+    coords.bnds.x1[1] = 4.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 1.0;
+
+    const Real expected_hx2v = coords.x1v();
+    const Real computed_hx2v = coords.hx2v();
+
+    REQUIRE(computed_hx2v == Approx(expected_hx2v).epsilon(1e-10));
+  }
+}
+
+TEST_CASE("Axisymmetric coordinate system", "[geometry][axisymmetric]") {
+  // Axisymmetric coordinates: (R, z, phi) with phi as azimuthal direction
+
+  SECTION("Volume calculation") {
+    // Volume: ∫∫∫ R dR dz dφ = [R²/2] × Δz × Δφ
+    // Similar to cylindrical but with different ordering
+
+    Coords<Coordinates::axisymmetric> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 3.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    // <R> = (R₀ + R₁) / 2 = 1.5
+    // Volume = <R> × ΔR × Δz × Δφ
+    const Real expected_volume = 1.5 * 1.0 * 3.0 * M_PI;
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Radial centroid (x1v)") {
+    // Same as cylindrical
+
+    Coords<Coordinates::axisymmetric> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 1.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0 * M_PI;
+
+    const Real expected_x1v = 14.0 / 9.0;
+    const Real computed_x1v = coords.x1v();
+
+    REQUIRE(computed_x1v == Approx(expected_x1v).epsilon(1e-10));
+  }
+
+  SECTION("X1 face area (radial face)") {
+    // Area: ∫∫ R dz dφ = R × Δz × Δφ
+
+    Coords<Coordinates::axisymmetric> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 2.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI / 2.0;
+
+    const Real R = 1.5;
+    const Real expected_area = R * 2.0 * (M_PI / 2.0);
+    const Real computed_area = coords.AreaX1(R);
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("X2 face area (z face)") {
+    // Area: ∫∫ R dR dφ = [R²/2] × Δφ = <R> × ΔR × Δφ
+
+    Coords<Coordinates::axisymmetric> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 1.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = M_PI;
+
+    const Real expected_area = 1.5 * 1.0 * M_PI;
+    const Real computed_area = coords.AreaX2(0.5); // z value doesn't matter
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("Scale factor hx3v") {
+    // In axisymmetric coords, h_phi = R
+
+    Coords<Coordinates::axisymmetric> coords;
+    coords.bnds.x1[0] = 2.0;
+    coords.bnds.x1[1] = 4.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 1.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0 * M_PI;
+
+    const Real expected_hx3v = coords.x1v();
+    const Real computed_hx3v = coords.hx3v();
+
+    REQUIRE(computed_hx3v == Approx(expected_hx3v).epsilon(1e-10));
+  }
+}
+
+TEST_CASE("Cartesian coordinate system", "[geometry][cartesian]") {
+  // Cartesian coordinates: (x, y, z)
+
+  SECTION("Volume calculation") {
+    // Volume: ∫∫∫ dx dy dz = Δx × Δy × Δz
+
+    Coords<Coordinates::cartesian> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 3.0;
+    coords.bnds.x2[0] = 2.0;
+    coords.bnds.x2[1] = 5.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 4.0;
+
+    const Real expected_volume = 2.0 * 3.0 * 4.0;
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Centroids are cell centers") {
+    // In Cartesian coords, centroids are simple arithmetic means
+
+    Coords<Coordinates::cartesian> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 3.0;
+    coords.bnds.x2[0] = 2.0;
+    coords.bnds.x2[1] = 6.0;
+    coords.bnds.x3[0] = -1.0;
+    coords.bnds.x3[1] = 1.0;
+
+    REQUIRE(coords.x1v() == Approx(2.0).epsilon(1e-10));
+    REQUIRE(coords.x2v() == Approx(4.0).epsilon(1e-10));
+    REQUIRE(coords.x3v() == Approx(0.0).epsilon(1e-10));
+  }
+
+  SECTION("Face areas") {
+    Coords<Coordinates::cartesian> coords;
+    coords.bnds.x1[0] = 0.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 3.0;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 5.0;
+
+    // X1 face: dy × dz
+    const Real expected_area_x1 = 3.0 * 5.0;
+    REQUIRE(coords.AreaX1(1.0) == Approx(expected_area_x1).epsilon(1e-10));
+
+    // X2 face: dx × dz
+    const Real expected_area_x2 = 2.0 * 5.0;
+    REQUIRE(coords.AreaX2(1.5) == Approx(expected_area_x2).epsilon(1e-10));
+
+    // X3 face: dx × dy
+    const Real expected_area_x3 = 2.0 * 3.0;
+    REQUIRE(coords.AreaX3(2.5) == Approx(expected_area_x3).epsilon(1e-10));
+  }
+
+  SECTION("Scale factors are unity") {
+    // In Cartesian coordinates, all scale factors h_i = 1
+
+    Coords<Coordinates::cartesian> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 3.0;
+    coords.bnds.x2[1] = 4.0;
+    coords.bnds.x3[0] = 5.0;
+    coords.bnds.x3[1] = 6.0;
+
+    REQUIRE(coords.hx1v() == 1.0);
+    REQUIRE(coords.hx2v() == 1.0);
+    REQUIRE(coords.hx3v() == 1.0);
+
+    REQUIRE(coords.hx1(1.5, 3.5, 5.5) == 1.0);
+    REQUIRE(coords.hx2(1.5, 3.5, 5.5) == 1.0);
+    REQUIRE(coords.hx3(1.5, 3.5, 5.5) == 1.0);
+  }
+
+  SECTION("No coordinate dependence") {
+    Coords<Coordinates::cartesian> coords;
+
+    REQUIRE(coords.x1dep() == false);
+    REQUIRE(coords.x2dep() == false);
+    REQUIRE(coords.x3dep() == false);
+  }
+}
+
+TEST_CASE("Spherical 1D coordinate system", "[geometry][spherical]") {
+  // 1D spherical (r only) with spherical symmetry
+
+  SECTION("Volume calculation") {
+    // Volume: ∫ r² dr × 4π (full solid angle)
+    // = [r³/3] × 4π (but implementation doesn't include 4π explicitly)
+
+    Coords<Coordinates::spherical1D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+
+    // dx1 * rfac where rfac = (r₀² + r₀r₁ + r₁²) / 3
+    const Real dx1 = 1.0;
+    const Real rfac = (1.0 + 2.0 + 4.0) / 3.0;
+    const Real expected_volume = rfac * dx1;
+    const Real computed_volume = coords.Volume();
+
+    REQUIRE(computed_volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Radial centroid (x1v)") {
+    // Volume-weighted radial centroid
+    // Same formula as 3D spherical
+
+    Coords<Coordinates::spherical1D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+
+    const Real expected_x1v = 45.0 / 28.0;
+    const Real computed_x1v = coords.x1v();
+
+    REQUIRE(computed_x1v == Approx(expected_x1v).epsilon(1e-10));
+  }
+
+  SECTION("X1 face area (radial face)") {
+    // Area: 4π r² (but implementation returns just r²)
+
+    Coords<Coordinates::spherical1D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+
+    const Real r = 1.5;
+    const Real expected_area = r * r;
+    const Real computed_area = coords.AreaX1(r);
+
+    REQUIRE(computed_area == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("X2 and X3 face areas") {
+    // Theta and phi faces: [r²/2] × Δr
+
+    Coords<Coordinates::spherical1D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+
+    const Real dx1 = 1.0;
+    const Real expected_area = 0.5 * (1.0 + 2.0) * dx1;
+
+    REQUIRE(coords.AreaX2(0.0) == Approx(expected_area).epsilon(1e-10));
+    REQUIRE(coords.AreaX3(0.0) == Approx(expected_area).epsilon(1e-10));
+  }
+
+  SECTION("Coordinate dependence") {
+    Coords<Coordinates::spherical1D> coords;
+
+    // Only x1 (r) dependent
+    REQUIRE(coords.x1dep() == true);
+    REQUIRE(coords.x2dep() == false);
+    REQUIRE(coords.x3dep() == false);
+  }
+
+  SECTION("Thin shell approximation") {
+    // For a thin shell, volume ≈ 4πr² × Δr
+
+    Coords<Coordinates::spherical1D> coords;
+    coords.bnds.x1[0] = 100.0;
+    coords.bnds.x1[1] = 100.1;
+
+    const Real volume = coords.Volume();
+    const Real dr = 0.1;
+    const Real r_mid = 100.05;
+    const Real expected_volume =
+        (r_mid * r_mid + r_mid * r_mid / 3.0 * (dr * dr) / (r_mid * r_mid)) * dr;
+
+    // For thin shell: V ≈ r² × dr
+    REQUIRE(volume == Approx(r_mid * r_mid * dr).epsilon(0.01));
+  }
+}
+
+TEST_CASE("Coordinate system edge cases", "[geometry]") {
+
+  SECTION("Spherical coordinate at pole (theta = 0)") {
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 1.0;
+    coords.bnds.x1[1] = 2.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 1e-6; // Near pole
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0 * M_PI;
+
+    const Real volume = coords.Volume();
+
+    // Should be very small but positive
+    REQUIRE(volume > 0.0);
+    REQUIRE(volume < 1e-4);
+  }
+
+  SECTION("Cylindrical coordinate on axis (R = 0)") {
+    Coords<Coordinates::cylindrical> coords;
+    coords.bnds.x1[0] = 0.0;
+    coords.bnds.x1[1] = 1.0;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = 2.0 * M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 1.0;
+
+    const Real volume = coords.Volume();
+
+    // Volume should still be computed correctly
+    const Real expected_volume = 0.5 * 1.0 * 2.0 * M_PI * 1.0;
+    REQUIRE(volume == Approx(expected_volume).epsilon(1e-10));
+  }
+
+  SECTION("Very thin spherical shell") {
+    Coords<Coordinates::spherical3D> coords;
+    coords.bnds.x1[0] = 10.0;
+    coords.bnds.x1[1] = 10.01;
+    coords.bnds.x2[0] = 0.0;
+    coords.bnds.x2[1] = M_PI;
+    coords.bnds.x3[0] = 0.0;
+    coords.bnds.x3[1] = 2.0 * M_PI;
+
+    const Real volume = coords.Volume();
+    const Real surface_area = 4.0 * M_PI * 10.0 * 10.0;
+    const Real dr = 0.01;
+
+    // For thin shell: V ≈ surface_area × dr
+    REQUIRE(volume == Approx(surface_area * dr).epsilon(0.01));
+  }
+}
diff --git a/tst/unit/test_ideal_gas.cpp b/tst/unit/test_ideal_gas.cpp
new file mode 100644
index 00000000..a6aaef49
--- /dev/null
+++ b/tst/unit/test_ideal_gas.cpp
@@ -0,0 +1,165 @@
+//========================================================================================
+// (C) (or copyright) 2025. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+//  license in this material to reproduce, prepare derivative works, distribute copies to
+//  the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+// This file was created in part or in whole by generative AI
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+#include "artemis.hpp"
+#include "utils/eos/eos.hpp"
+
+using Catch::Approx;
+
+TEST_CASE("IdealGas EOS recovers input gamma", "[eos][ideal_gas]") {
+  // Create an ideal gas EOS with known gamma
+  constexpr Real gamma = 5.0 / 3.0; // Monoatomic ideal gas
+  constexpr Real cv = 1.0e7;        // Specific heat (arbitrary units)
+
+  // IdealGas constructor takes (Gamma - 1, Cv)
+  singularity::IdealGas eos(gamma - 1.0, cv);
+
+  SECTION("Gamma is correctly recovered") {
+    // The IdealGas EOS should return the gamma we set
+    Real recovered_gamma = eos.GruneisenParamFromDensityInternalEnergy(1.0, 1.0);
+
+    // Gruneisen parameter for ideal gas is (gamma - 1)
+    REQUIRE(recovered_gamma == Approx(gamma - 1.0).epsilon(1e-12));
+  }
+
+  SECTION("Pressure-density-temperature relation is consistent") {
+    // Test the ideal gas law: P = rho * (gamma - 1) * cv * T
+    Real rho = 1.0e-3;
+    Real temp = 1.0e4;
+
+    Real pressure = eos.PressureFromDensityTemperature(rho, temp);
+    Real expected_pressure = rho * (gamma - 1.0) * cv * temp;
+
+    REQUIRE(pressure == Approx(expected_pressure).epsilon(1e-12));
+  }
+
+  SECTION("Internal energy-density-temperature relation is consistent") {
+    // Test: sie = cv * T
+    Real rho = 1.0e-3;
+    Real temp = 1.0e4;
+
+    Real sie = eos.InternalEnergyFromDensityTemperature(rho, temp);
+    Real expected_sie = cv * temp;
+
+    REQUIRE(sie == Approx(expected_sie).epsilon(1e-12));
+  }
+
+  SECTION("Pressure from density and internal energy") {
+    // Test: P = rho * sie * (gamma - 1)
+    Real rho = 1.0e-3;
+    Real sie = 1.0e11;
+
+    Real pressure = eos.PressureFromDensityInternalEnergy(rho, sie);
+    Real expected_pressure = rho * sie * (gamma - 1.0);
+
+    REQUIRE(pressure == Approx(expected_pressure).epsilon(1e-12));
+  }
+
+  SECTION("Temperature from density and internal energy") {
+    // Test: T = sie / cv
+    Real rho = 1.0e-3;
+    Real sie = 1.0e11;
+
+    Real temp = eos.TemperatureFromDensityInternalEnergy(rho, sie);
+    Real expected_temp = sie / cv;
+
+    REQUIRE(temp == Approx(expected_temp).epsilon(1e-12));
+  }
+}
+
+TEST_CASE("IdealGas EOS with different gamma values", "[eos][ideal_gas]") {
+  constexpr Real cv = 1.0e7;
+
+  SECTION("Monoatomic gas (gamma = 5/3)") {
+    constexpr Real gamma = 5.0 / 3.0;
+    singularity::IdealGas eos(gamma - 1.0, cv);
+
+    Real gruneisen = eos.GruneisenParamFromDensityInternalEnergy(1.0, 1.0);
+    REQUIRE(gruneisen == Approx(gamma - 1.0).epsilon(1e-12));
+  }
+}
+
+TEST_CASE("IdealGas EOS with UnitSystem wrapper", "[eos][ideal_gas][units]") {
+  // Set up unit conversions (example values similar to gas.cpp)
+  constexpr Real gamma = 5.0 / 3.0;
+  constexpr Real cv_physical = 1.0e7; // Physical units (e.g., erg/g/K)
+  constexpr Real time_code_to_physical = 1.0;
+  constexpr Real mass_code_to_physical = 1.0;
+  constexpr Real length_code_to_physical = 1.0;
+  constexpr Real temperature_code_to_physical = 1.0;
+
+  // Create EOS with UnitSystem wrapper (as done in gas.cpp around line 98)
+  ArtemisUtils::EOS eos_host = singularity::UnitSystem<singularity::IdealGas>(
+      singularity::IdealGas(gamma - 1.0, cv_physical),
+      singularity::eos_units_init::LengthTimeUnitsInit(), time_code_to_physical,
+      mass_code_to_physical, length_code_to_physical, temperature_code_to_physical);
+
+  SECTION("Gamma is correctly recovered through UnitSystem") {
+    Real recovered_gamma = eos_host.GruneisenParamFromDensityInternalEnergy(1.0, 1.0);
+    REQUIRE(recovered_gamma == Approx(gamma - 1.0).epsilon(1e-12));
+  }
+
+  SECTION("Pressure calculation with UnitSystem") {
+    Real rho = 1.0e-3;
+    Real temp = 1.0e4;
+
+    Real pressure = eos_host.PressureFromDensityTemperature(rho, temp);
+    Real expected_pressure = rho * (gamma - 1.0) * cv_physical * temp;
+
+    REQUIRE(pressure == Approx(expected_pressure).epsilon(1e-12));
+  }
+
+  SECTION("Temperature recovery with UnitSystem") {
+    Real rho = 1.0e-3;
+    Real sie = 1.0e11;
+
+    Real temp = eos_host.TemperatureFromDensityInternalEnergy(rho, sie);
+    Real expected_temp = sie / cv_physical;
+
+    REQUIRE(temp == Approx(expected_temp).epsilon(1e-12));
+  }
+
+  SECTION("Unit conversions with different scale factors") {
+    // Test with non-trivial unit conversions
+    constexpr Real time_scale = 2.0;
+    constexpr Real mass_scale = 3.0;
+    constexpr Real length_scale = 4.0;
+    constexpr Real temp_scale = 5.0;
+
+    ArtemisUtils::EOS eos_scaled = singularity::UnitSystem<singularity::IdealGas>(
+        singularity::IdealGas(gamma - 1.0, cv_physical),
+        singularity::eos_units_init::LengthTimeUnitsInit(), time_scale, mass_scale,
+        length_scale, temp_scale);
+
+    // The EOS should still work with scaled units
+    Real rho_code = 1.0;
+    Real temp_code = 1.0;
+
+    Real pressure_code = eos_scaled.PressureFromDensityTemperature(rho_code, temp_code);
+
+    // Verify pressure is positive and physically reasonable
+    REQUIRE(pressure_code > 0.0);
+
+    // Verify round-trip consistency in code units
+    Real sie_code = eos_scaled.InternalEnergyFromDensityTemperature(rho_code, temp_code);
+    Real temp_recovered =
+        eos_scaled.TemperatureFromDensityInternalEnergy(rho_code, sie_code);
+
+    REQUIRE(temp_recovered == Approx(temp_code).epsilon(1e-10));
+  }
+}
diff --git a/tst/unit/test_ideal_hhe.cpp b/tst/unit/test_ideal_hhe.cpp
new file mode 100644
index 00000000..58afaed2
--- /dev/null
+++ b/tst/unit/test_ideal_hhe.cpp
@@ -0,0 +1,261 @@
+//========================================================================================
+// (C) (or copyright) 2025. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+//  license in this material to reproduce, prepare derivative works, distribute copies to
+//  the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+// This file was created in part or in whole by generative AI
+
+#include <catch2/catch_approx.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <cmath>
+
+#include "artemis.hpp"
+#include "utils/eos/eos.hpp"
+
+using ArtemisEOS::IdealHHe;
+using Catch::Approx;
+using singularity::IdealGas;
+
+#ifdef SPINER_USE_HDF
+TEST_CASE("IdealHHe EOS basic construction and properties", "[IdealHHe][EOS]") {
+  // Use solar composition: X = 0.7, Y = 0.28 (Z = 0.02)
+  const Real X = 0.7;
+  const Real Y = 0.28;
+
+  // Temperature range: 1 K to 1e6 K (ltmin=0, ltmax=6 in log10)
+  const Real ltmin = 0.0;
+  const Real ltmax = 6.0;
+  const int nt = 50;
+
+  // Density range: 1e-15 to 1e-3 g/cc (ldmin=-15, ldmax=-3 in log10)
+  const Real ldmin = -15.0;
+  const Real ldmax = -3.0;
+  const int nd = 50;
+
+  // Create EOS without table lookup for testing
+  const bool use_table = false;
+  const std::string save_to_file = "";
+
+  IdealHHe eos(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table);
+
+  SECTION("Temperature recovery from density and internal energy") {
+    const Real rho = 1.0e-10;    // g/cc
+    const Real T_input = 1000.0; // K
+
+    // Get internal energy at this state
+    const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T_input);
+
+    // Recover temperature from density and internal energy
+    const Real T_recovered = eos.TemperatureFromDensityInternalEnergy(rho, sie);
+
+    // Temperature should be recovered to within 1%
+    REQUIRE(T_recovered == Approx(T_input).epsilon(0.01));
+  }
+
+  SECTION("Pressure calculation from density and temperature") {
+    const Real rho = 1.0e-10; // g/cc
+    const Real T = 1000.0;    // K
+
+    const Real P = eos.PressureFromDensityTemperature(rho, T);
+
+    // Pressure should be positive
+    REQUIRE(P > 0.0);
+  }
+
+  SECTION("Round-trip consistency: rho, T -> sie -> T") {
+    const Real rho = 1.0e-12;       // g/cc
+    const Real T_original = 5000.0; // K
+
+    // Convert T to sie
+    const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T_original);
+
+    // Convert sie back to T
+    const Real T_recovered = eos.TemperatureFromDensityInternalEnergy(rho, sie);
+
+    // Should recover original temperature to within 1%
+    REQUIRE(T_recovered == Approx(T_original).epsilon(0.01));
+  }
+
+  SECTION("Internal energy increases with temperature") {
+    const Real rho = 1.0e-10; // g/cc
+    const Real T1 = 1000.0;   // K
+    const Real T2 = 2000.0;   // K
+
+    const Real sie1 = eos.InternalEnergyFromDensityTemperature(rho, T1);
+    const Real sie2 = eos.InternalEnergyFromDensityTemperature(rho, T2);
+
+    // Higher temperature should give higher internal energy
+    REQUIRE(sie2 > sie1);
+  }
+}
+
+TEST_CASE("IdealHHe EOS with different compositions", "[IdealHHe][EOS]") {
+  const Real ltmin = 0.0;
+  const Real ltmax = 6.0;
+  const int nt = 50;
+  const Real ldmin = -15.0;
+  const Real ldmax = -3.0;
+  const int nd = 50;
+  const bool use_table = false;
+  const std::string save_to_file = "";
+
+  SECTION("Hydrogen-rich composition (X=0.9, Y=0.1)") {
+    IdealHHe eos_H_rich(0.9, 0.1, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file,
+                        use_table);
+
+    const Real rho = 1.0e-10;
+    const Real T = 1000.0;
+
+    const Real sie = eos_H_rich.InternalEnergyFromDensityTemperature(rho, T);
+    const Real P = eos_H_rich.PressureFromDensityTemperature(rho, T);
+
+    REQUIRE(sie > 0.0);
+    REQUIRE(P > 0.0);
+  }
+
+  SECTION("Helium-rich composition (X=0.3, Y=0.7)") {
+    IdealHHe eos_He_rich(0.3, 0.7, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file,
+                         use_table);
+
+    const Real rho = 1.0e-10;
+    const Real T = 1000.0;
+
+    const Real sie = eos_He_rich.InternalEnergyFromDensityTemperature(rho, T);
+    const Real P = eos_He_rich.PressureFromDensityTemperature(rho, T);
+
+    REQUIRE(sie > 0.0);
+    REQUIRE(P > 0.0);
+  }
+}
+
+TEST_CASE("IdealHHe EOS with UnitSystem wrapper", "[IdealHHe][EOS][UnitSystem]") {
+  // Use solar composition
+  const Real X = 0.7;
+  const Real Y = 0.28;
+  const Real ltmin = 0.0;
+  const Real ltmax = 6.0;
+  const int nt = 50;
+  const Real ldmin = -15.0;
+  const Real ldmax = -3.0;
+  const int nd = 50;
+  const bool use_table = false;
+  const std::string save_to_file = "";
+
+  IdealHHe eos_base(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table);
+
+  SECTION("Unit conversion with CGS units") {
+    // Define unit system: time [s], mass [g], length [cm], temperature [K]
+    const Real time_cgs = 1.0;
+    const Real mass_cgs = 1.0;
+    const Real length_cgs = 1.0;
+    const Real temp_cgs = 1.0;
+
+    // Create IdealHHe inline as temporary for UnitSystem
+    auto eos = singularity::UnitSystem<IdealHHe>(
+        IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table),
+        singularity::eos_units_init::LengthTimeUnitsInit(), time_cgs, mass_cgs,
+        length_cgs, temp_cgs);
+
+    const Real rho = 1.0e-10; // g/cc
+    const Real T = 1000.0;    // K
+
+    const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T);
+    const Real P = eos.PressureFromDensityTemperature(rho, T);
+
+    REQUIRE(sie > 0.0);
+    REQUIRE(P > 0.0);
+  }
+
+  SECTION("Unit conversion with scaled units") {
+    // Define unit system with scale factors
+    // Time: 1 day = 86400 s
+    // Mass: 1 solar mass = 1.989e33 g
+    // Length: 1 AU = 1.496e13 cm
+    // Temperature: 1 K
+    const Real time_scale = 86400.0;
+    const Real mass_scale = 1.989e33;
+    const Real length_scale = 1.496e13;
+    const Real temp_scale = 1.0;
+
+    // Create IdealHHe inline as temporary for UnitSystem
+    auto eos = singularity::UnitSystem<IdealHHe>(
+        IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table),
+        singularity::eos_units_init::LengthTimeUnitsInit(), time_scale, mass_scale,
+        length_scale, temp_scale);
+
+    // Density in code units (Msun/AU^3)
+    // Convert 1e-10 g/cc to Msun/AU^3: 1e-10 * (AU^3/Msun)
+    const Real rho_cgs = 1.0e-10; // g/cc
+    const Real rho_code = rho_cgs * std::pow(length_scale, 3) / mass_scale;
+
+    // Temperature in code units (same as K in this case)
+    const Real T_code = 1000.0;
+
+    const Real sie_code = eos.InternalEnergyFromDensityTemperature(rho_code, T_code);
+    const Real P_code = eos.PressureFromDensityTemperature(rho_code, T_code);
+
+    // Convert back to CGS for comparison
+    // sie: [cm^2/s^2] -> code units: [length^2/time^2]
+    const Real sie_cgs = sie_code * std::pow(length_scale, 2) / std::pow(time_scale, 2);
+
+    // Verify physical values are positive
+    REQUIRE(sie_cgs > 0.0);
+    REQUIRE(P_code > 0.0);
+  }
+
+  SECTION("Consistency between wrapped and unwrapped EOS in CGS units") {
+    // Wrap with CGS units (scale factors = 1)
+    auto eos_wrapped = singularity::UnitSystem<IdealHHe>(
+        IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table),
+        singularity::eos_units_init::LengthTimeUnitsInit(), 1.0, 1.0, 1.0, 1.0);
+
+    const Real rho = 1.0e-10; // g/cc
+    const Real T = 1000.0;    // K
+
+    // Get results from both
+    const Real sie_unwrapped = eos_base.InternalEnergyFromDensityTemperature(rho, T);
+    const Real sie_wrapped = eos_wrapped.InternalEnergyFromDensityTemperature(rho, T);
+
+    const Real P_unwrapped = eos_base.PressureFromDensityTemperature(rho, T);
+    const Real P_wrapped = eos_wrapped.PressureFromDensityTemperature(rho, T);
+
+    // Should match exactly with unit scale factors = 1
+    REQUIRE(sie_wrapped == Approx(sie_unwrapped).epsilon(1e-10));
+    REQUIRE(P_wrapped == Approx(P_unwrapped).epsilon(1e-10));
+  }
+
+  SECTION("Round-trip with unit conversions") {
+    // Use code units
+    const Real time_scale = 1.0e3;
+    const Real mass_scale = 1.0e30;
+    const Real length_scale = 1.0e10;
+    const Real temp_scale = 1.0;
+
+    // Create IdealHHe inline as temporary for UnitSystem
+    auto eos = singularity::UnitSystem<IdealHHe>(
+        IdealHHe(X, Y, ltmin, ltmax, nt, ldmin, ldmax, nd, save_to_file, use_table),
+        singularity::eos_units_init::LengthTimeUnitsInit(), time_scale, mass_scale,
+        length_scale, temp_scale);
+
+    // Density in code units
+    const Real rho_cgs = 1.0e-10;
+    const Real rho_code = rho_cgs * std::pow(length_scale, 3) / mass_scale;
+    const Real T_code = 1000.0;
+
+    // Round-trip: T -> sie -> T
+    const Real sie = eos.InternalEnergyFromDensityTemperature(rho_code, T_code);
+    const Real T_recovered = eos.TemperatureFromDensityInternalEnergy(rho_code, sie);
+
+    // Should recover temperature to within 1%
+    REQUIRE(T_recovered == Approx(T_code).epsilon(0.01));
+  }
+}
+#endif // SPINER_USE_HDF