diff --git a/.github/workflows/build_and_test_libshotscope.yml b/.github/workflows/build_and_test_libshotscope.yml
index 471288b..40bb45d 100644
--- a/.github/workflows/build_and_test_libshotscope.yml
+++ b/.github/workflows/build_and_test_libshotscope.yml
@@ -15,10 +15,10 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
 
     - name: Cache CMake build
-      uses: actions/cache@v4
+      uses: actions/cache@v5
       with:
         path: |
           build
diff --git a/.github/workflows/ci-multi-platform.yml b/.github/workflows/ci-multi-platform.yml
index 714e9ef..a714eeb 100644
--- a/.github/workflows/ci-multi-platform.yml
+++ b/.github/workflows/ci-multi-platform.yml
@@ -39,7 +39,7 @@ jobs:
             compiler: default
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
 
     - name: Set up compiler (Linux)
       if: runner.os == 'Linux' && matrix.compiler != 'default'
@@ -50,7 +50,7 @@ jobs:
         fi
 
     - name: Cache CMake build
-      uses: actions/cache@v4
+      uses: actions/cache@v5
       with:
         path: |
           build
@@ -79,7 +79,7 @@ jobs:
 
     - name: Upload build artifacts
       if: matrix.os == 'ubuntu-latest' && matrix.compiler == 'gcc'
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v7
       with:
         name: libgolf-${{ runner.os }}
         path: |
diff --git a/.github/workflows/code-coverage.yml b/.github/workflows/code-coverage.yml
index b0c967d..c3b4d90 100644
--- a/.github/workflows/code-coverage.yml
+++ b/.github/workflows/code-coverage.yml
@@ -11,7 +11,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
 
     - name: Install dependencies
       run: |
@@ -50,13 +50,13 @@ jobs:
         lcov --summary coverage.info 2>&1 | tee -a coverage-summary.txt
 
     - name: Upload coverage HTML
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v7
       with:
         name: coverage-report
         path: coverage_html/
 
     - name: Upload coverage summary
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v7
       with:
         name: coverage-summary
         path: coverage-summary.txt
diff --git a/.github/workflows/pages.yml b/.github/workflows/pages.yml
index cc2fd1d..67fed9e 100644
--- a/.github/workflows/pages.yml
+++ b/.github/workflows/pages.yml
@@ -29,13 +29,13 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
 
       - name: Configure Pages
-        uses: actions/configure-pages@v5
+        uses: actions/configure-pages@v6
 
       - name: Set up Emscripten
-        uses: mymindstorm/setup-emsdk@v14
+        uses: mymindstorm/setup-emsdk@v16
 
       - name: Configure CMake (wasm preset)
         run: cmake --preset wasm
@@ -63,7 +63,7 @@ jobs:
           destination: ./_site
 
       - name: Upload Pages artifact
-        uses: actions/upload-pages-artifact@v3
+        uses: actions/upload-pages-artifact@v5
         with:
           path: _site
 
@@ -76,4 +76,4 @@ jobs:
       url: ${{ steps.deployment.outputs.page_url }}
     steps:
       - id: deployment
-        uses: actions/deploy-pages@v4
+        uses: actions/deploy-pages@v5
diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml
index 3fb4798..5499557 100644
--- a/.github/workflows/release.yml
+++ b/.github/workflows/release.yml
@@ -30,7 +30,7 @@ jobs:
             platform: Windows
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
 
       - name: Get version from tag
         id: get_version
@@ -64,7 +64,7 @@ jobs:
           7z a ../libgolf-${{ matrix.platform }}-${{ steps.get_version.outputs.version }}.zip .
 
       - name: Upload release assets
-        uses: softprops/action-gh-release@v2
+        uses: softprops/action-gh-release@v3
         with:
           files: libgolf-${{ matrix.platform }}-${{ steps.get_version.outputs.version }}.${{ matrix.archive_ext }}
           body: |
@@ -81,14 +81,14 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v5
 
       - name: Get version from tag
         id: get_version
         run: echo "version=${GITHUB_REF#refs/tags/v}" >> $GITHUB_OUTPUT
 
       - name: Set up Emscripten
-        uses: mymindstorm/setup-emsdk@v14
+        uses: mymindstorm/setup-emsdk@v16
 
       - name: Configure CMake (wasm preset)
         run: cmake --preset wasm
@@ -105,6 +105,6 @@ jobs:
           zip -r ../libgolf-wasm-${{ steps.get_version.outputs.version }}.zip .
 
       - name: Upload release assets
-        uses: softprops/action-gh-release@v2
+        uses: softprops/action-gh-release@v3
         with:
           files: libgolf-wasm-${{ steps.get_version.outputs.version }}.zip
diff --git a/.github/workflows/static-analysis.yml b/.github/workflows/static-analysis.yml
index ffb95e1..d72794a 100644
--- a/.github/workflows/static-analysis.yml
+++ b/.github/workflows/static-analysis.yml
@@ -12,7 +12,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
 
     - name: Install dependencies
       run: |
@@ -31,7 +31,7 @@ jobs:
 
     - name: Upload clang-tidy report
       if: always()
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v7
       with:
         name: clang-tidy-report
         path: clang-tidy-report.txt
@@ -41,7 +41,7 @@ jobs:
     runs-on: ubuntu-latest
 
     steps:
-    - uses: actions/checkout@v4
+    - uses: actions/checkout@v5
 
     - name: Install cppcheck
       run: |
@@ -62,7 +62,7 @@ jobs:
 
     - name: Upload cppcheck report
       if: always()
-      uses: actions/upload-artifact@v4
+      uses: actions/upload-artifact@v7
       with:
         name: cppcheck-report
         path: cppcheck-report.txt
diff --git a/CMakeLists.txt b/CMakeLists.txt
index df346e6..c82b4a7 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.14)
-project(golf VERSION 4.5.2 LANGUAGES CXX)
+project(golf VERSION 4.6.0 LANGUAGES CXX)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -21,9 +21,6 @@ configure_file(
     @ONLY
 )
 
-include_directories(${PROJECT_INCLUDE_DIR})
-include_directories(${PROJECT_BINARY_DIR}/include)
-
 set(SOURCES
    ${PROJECT_SRC_DIR}/math_utils.cpp
    ${PROJECT_SRC_DIR}/ShotPhysicsContext.cpp
@@ -41,16 +38,28 @@ set(HEADERS
    ${PROJECT_INCLUDE_DIR}/FlightPhase.hpp
    ${PROJECT_INCLUDE_DIR}/FlightSimulator.hpp
    ${PROJECT_INCLUDE_DIR}/BallState.hpp
+   ${PROJECT_INCLUDE_DIR}/BallProperties.hpp
    ${PROJECT_INCLUDE_DIR}/ground_surface.hpp
    ${PROJECT_INCLUDE_DIR}/terrain_interface.hpp
    ${PROJECT_INCLUDE_DIR}/ground_physics.hpp
    ${PROJECT_INCLUDE_DIR}/physics_constants.hpp
    ${PROJECT_INCLUDE_DIR}/RollModel.hpp
    ${PROJECT_INCLUDE_DIR}/DefaultRollModel.hpp
+   ${PROJECT_INCLUDE_DIR}/Integrator.hpp
+   ${PROJECT_INCLUDE_DIR}/DefaultIntegrator.hpp
 )
 
 # Create the static library
 add_library(${PROJECT_NAME} STATIC ${SOURCES} ${HEADERS})
+add_library(golf::golf ALIAS ${PROJECT_NAME})
+
+# Publish the include paths as usage requirements so in-tree targets and
+# downstream consumers (find_package / add_subdirectory) both resolve headers.
+target_include_directories(${PROJECT_NAME} PUBLIC
+    $<BUILD_INTERFACE:${PROJECT_INCLUDE_DIR}>
+    $<BUILD_INTERFACE:${PROJECT_BINARY_DIR}/include>
+    $<INSTALL_INTERFACE:include>
+)
 
 # Set library version properties
 set_target_properties(${PROJECT_NAME} PROPERTIES
@@ -58,21 +67,12 @@ set_target_properties(${PROJECT_NAME} PROPERTIES
     SOVERSION ${PROJECT_VERSION_MAJOR}
 )
 
-# Set compilation flags based on compiler
+# No -ffast-math: -ffinite-math-only folds NaN comparisons to a constant and
+# breaks the convergence guards that rely on `NaN <= height` staying false.
 if(MSVC)
-    # MSVC-specific flags
     set(CMAKE_CXX_FLAGS "/W4")
-    set(CMAKE_CXX_FLAGS_RELEASE "/O2")
 else()
-    # GCC/Clang flags
     set(CMAKE_CXX_FLAGS "-Wall -Wextra")
-    # No -ffast-math: it implies -ffinite-math-only, which lets the compiler
-    # assume NaN/Inf never occur and fold NaN comparisons to a constant. The
-    # convergence guards (FlightSimulator run loop + *Phase::isPhaseComplete)
-    # rely on IEEE semantics — `NaN <= height` must stay false so a poisoned
-    # trajectory never falsely "completes". Apple Clang folded it the other way,
-    # breaking the hang guards on macOS.
-    set(CMAKE_CXX_FLAGS_RELEASE "-O3")
 endif()
 
 # Code coverage option
@@ -92,6 +92,7 @@ if(NOT EMSCRIPTEN)
       GIT_TAG v1.17.0
     )
     set(gtest_force_shared_crt ON CACHE BOOL "" FORCE)
+    set(INSTALL_GTEST OFF) # keep gtest/gmock out of our install prefix
     FetchContent_MakeAvailable(googletest)
 
     # Enable testing
@@ -163,6 +164,43 @@ if(EMSCRIPTEN)
     )
 endif()
 
-install(TARGETS golf DESTINATION lib)
-install(DIRECTORY include/ DESTINATION include)
-install(FILES "${PROJECT_BINARY_DIR}/include/version.hpp" DESTINATION include)
\ No newline at end of file
+# ----------------------------------------------------------------------------
+# Install + export: make find_package(golf) and the golf::golf target work for
+# downstream projects.
+# ----------------------------------------------------------------------------
+include(GNUInstallDirs)
+include(CMakePackageConfigHelpers)
+
+install(TARGETS ${PROJECT_NAME}
+    EXPORT golfTargets
+    ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+    INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+)
+install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+install(FILES "${PROJECT_BINARY_DIR}/include/version.hpp"
+        DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+
+set(GOLF_CMAKE_CONFIG_DIR ${CMAKE_INSTALL_LIBDIR}/cmake/golf)
+
+install(EXPORT golfTargets
+    FILE golfTargets.cmake
+    NAMESPACE golf::
+    DESTINATION ${GOLF_CMAKE_CONFIG_DIR}
+)
+
+configure_package_config_file(
+    "${PROJECT_SOURCE_DIR}/cmake/golfConfig.cmake.in"
+    "${PROJECT_BINARY_DIR}/golfConfig.cmake"
+    INSTALL_DESTINATION ${GOLF_CMAKE_CONFIG_DIR}
+)
+write_basic_package_version_file(
+    "${PROJECT_BINARY_DIR}/golfConfigVersion.cmake"
+    VERSION ${PROJECT_VERSION}
+    COMPATIBILITY SameMajorVersion
+)
+install(FILES
+    "${PROJECT_BINARY_DIR}/golfConfig.cmake"
+    "${PROJECT_BINARY_DIR}/golfConfigVersion.cmake"
+    DESTINATION ${GOLF_CMAKE_CONFIG_DIR}
+)
diff --git a/README.md b/README.md
index 9fb108e..37267bb 100644
--- a/README.md
+++ b/README.md
@@ -19,6 +19,22 @@ chmod +x build.sh
 ./build.sh
 ```
 
+## Using libgolf in your project
+
+After installing (`cmake --install build`), consume it from another CMake
+project with `find_package`:
+
+```cmake
+find_package(golf REQUIRED)
+
+add_executable(my_app main.cpp)
+target_link_libraries(my_app PRIVATE golf::golf)
+```
+
+`golf::golf` carries its include paths, so `#include <libgolf.hpp>` works with
+no extra configuration. The same target name is available via
+`add_subdirectory(libgolf)` for in-tree builds.
+
 ## Features
 
 - Full trajectory simulation with automatic phase transitions (aerial → bounce → roll)
diff --git a/build.sh b/build.sh
index ece80b8..b965457 100755
--- a/build.sh
+++ b/build.sh
@@ -1,5 +1,12 @@
 #!/bin/bash
 
+RUN_TESTS=0
+for arg in "$@"; do
+    case "$arg" in
+        -t|--test) RUN_TESTS=1 ;;
+    esac
+done
+
 # Use Clang if available and not overridden, otherwise use system default
 if [ -z "$CC" ] && [ -z "$CXX" ]; then
     if command -v clang &> /dev/null; then
@@ -23,6 +30,11 @@ cmake ..
 # Build the project
 cmake --build .
 
+if [ "$RUN_TESTS" -eq 1 ]; then
+    echo "Running tests"
+    ./libgolf_tests
+fi
+
 # Optionally install (uncomment if you want to install system-wide)
 # echo "Installing to system directories (password may be required)"
 # sudo cmake --install .
\ No newline at end of file
diff --git a/cmake/golfConfig.cmake.in b/cmake/golfConfig.cmake.in
new file mode 100644
index 0000000..12adb8c
--- /dev/null
+++ b/cmake/golfConfig.cmake.in
@@ -0,0 +1,5 @@
+@PACKAGE_INIT@
+
+include("${CMAKE_CURRENT_LIST_DIR}/golfTargets.cmake")
+
+check_required_components(golf)
diff --git a/docs/how.md b/docs/how.md
index 8f3c440..f02ca8d 100644
--- a/docs/how.md
+++ b/docs/how.md
@@ -68,6 +68,15 @@ const AtmosphericData atmos{
 };
 ```
 
+Every field defaults to a sea-level standard day (59°F, 29.92 inHg, no wind,
+dry air), so `AtmosphericData{}` is a valid baseline and you only need to set
+the fields that differ from standard:
+
+```c++
+AtmosphericData atmos{}; // standard day
+atmos.elevation = 5280.0f; // mile-high course
+```
+
 Field definitions are documented in `include/atmospheric_data.hpp`.
 
 ### 3. Ground Surface Properties
@@ -248,14 +257,62 @@ FlightSimulator sim(ball, atmos, ground, /*aero*/ nullptr, /*bounce*/ nullptr, r
 
 See [Roll Models](roll_model.md) for the interface and a worked example.
 
+### Custom Ball Properties
+
+The simulator models a standard golf ball by default. Pass a `BallProperties`
+to simulate a different ball; its mass and circumference feed the aerodynamic
+coefficients and the radius the force models receive:
+
+```c++
+BallProperties ball{.massOz = 1.80f, .circumferenceIn = 5.30f};
+FlightSimulator sim(launch, atmos, ground,
+                    /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr, ball);
+```
+
+A default-constructed `BallProperties{}` reproduces the standard ball exactly,
+so omitting the argument leaves results unchanged.
+
+### Custom Gravity
+
+Gravity defaults to Earth (`32.174 ft/s²`) and can be set through the
+`FlightSimulator` constructor. It is applied to the aerial and between-bounce
+flight integration:
+
+```c++
+constexpr float kMoonGravity = 5.31f; // ft/s²
+FlightSimulator sim(launch, atmos, ground,
+                    /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr,
+                    /*ball*/ BallProperties{}, kMoonGravity);
+```
+
+The built-in roll model decelerates under Earth gravity; a custom `RollModel`
+can use any value.
+
+### Custom Integrator
+
+The aerial and between-bounce flight integration uses a semi-implicit Euler
+scheme by default. Implement `Integrator` to substitute your own (e.g. RK4 or
+an adaptive step) and pass it to `FlightSimulator`:
+
+```c++
+auto integrator = std::make_shared<MyRK4Integrator>();
+FlightSimulator sim(launch, atmos, ground,
+                    /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr,
+                    /*ball*/ BallProperties{},
+                    physics_constants::GRAVITY_FT_PER_S2, integrator);
+```
+
+The flight phase owns spin decay, wind, and the acceleration model; the
+integrator owns only how position and velocity advance. It receives an
+acceleration field it can sample at trial states. The roll phase runs its own
+integrator inside `RollModel`.
+
 ### What Isn't Pluggable
 
 You can replace the three per-phase physics models (aerodynamics, bounce, roll) and the terrain. Everything else is fixed in the current release:
 
-- **Gravity** — a fixed constant, not a constructor parameter.
-- **Ball properties** — mass, circumference, and radius use standard golf-ball constants. A custom `AerodynamicModel` *reads* `ballRadius` and `c0` through `AerodynamicState`, but cannot change the constants the simulator itself bakes into `c0` and the spin-rate scaling.
 - **Air model** — the air-density, viscosity, and saturation-vapor-pressure formulas are fixed. You supply `AtmosphericData` inputs; you cannot swap the model that converts them into density.
-- **Integrator and phase machine** — the aerial time integration and the aerial → bounce → roll transition logic are internal. You can replace what each phase *computes*, not how it is stepped or sequenced.
+- **Phase machine** — the aerial → bounce → roll transition logic and the criteria for when each transition fires are internal. You can replace what each phase *computes* and how the flight phases step (see Custom Integrator), but not how the phases are sequenced.
 - **Launch transform** — the mapping from `LaunchData` (launch-monitor inputs) to the initial state vector is fixed.
 
 ### Example Programs
diff --git a/include/AerodynamicModel.hpp b/include/AerodynamicModel.hpp
index 92896cd..79e783f 100644
--- a/include/AerodynamicModel.hpp
+++ b/include/AerodynamicModel.hpp
@@ -79,8 +79,9 @@ struct AerodynamicState
  *
  * This model sets aerodynamic forces, not the integration. AerialPhase runs the
  * step loop and picks the timestep; you supply acceleration and a spin-decay
- * constant for it to integrate. You cannot swap the aerial integrator itself.
- * RollModel works the other way: it gets dt and runs its own integrator.
+ * constant for it to integrate. The integration scheme itself is swappable
+ * separately through the Integrator interface. RollModel works differently
+ * again: it gets dt and runs its own integrator internally.
  *
  * computeSpinDecayTau only affects spin in the air. Bounce spin comes from
  * BounceModel (BounceResult::newSpinVector) and roll spin from RollModel
diff --git a/include/BallProperties.hpp b/include/BallProperties.hpp
new file mode 100644
index 0000000..a3d20ab
--- /dev/null
+++ b/include/BallProperties.hpp
@@ -0,0 +1,39 @@
+#ifndef BALL_PROPERTIES_HPP
+#define BALL_PROPERTIES_HPP
+
+#include "physics_constants.hpp"
+
+/**
+ * @brief Physical properties of the ball being simulated.
+ *
+ * Defaults reproduce a standard golf ball, so a default-constructed
+ * BallProperties{} leaves the simulation identical to the historical
+ * hardcoded behaviour. Supply custom values to simulate a different ball:
+ *
+ * @code
+ * BallProperties heavy{.massOz = 1.80f, .circumferenceIn = 5.30f};
+ * FlightSimulator sim(launch, atmos, ground,
+ *                     nullptr, nullptr, nullptr, heavy);
+ * @endcode
+ *
+ * Mass and circumference feed the lumped aerodynamic coefficients (c0, the
+ * Reynolds reference, and surface spin speed); the radius the phases hand to
+ * the force models is derived from circumference.
+ */
+struct BallProperties
+{
+	/// Ball mass (ounces).
+	float massOz = physics_constants::STD_BALL_MASS_OZ;
+
+	/// Ball circumference (inches).
+	float circumferenceIn = physics_constants::STD_BALL_CIRCUMFERENCE_IN;
+
+	/// Ball radius (feet), derived from circumference.
+	[[nodiscard]] float radiusFt() const
+	{
+		return circumferenceIn / (2.0F * physics_constants::PI) /
+		       physics_constants::INCHES_PER_FOOT;
+	}
+};
+
+#endif // BALL_PROPERTIES_HPP
diff --git a/include/DefaultAerodynamicModel.hpp b/include/DefaultAerodynamicModel.hpp
index af26663..950530f 100644
--- a/include/DefaultAerodynamicModel.hpp
+++ b/include/DefaultAerodynamicModel.hpp
@@ -17,8 +17,8 @@
  *
  * where vw = |v - v_wind| and omega is the current spin vector.
  *
- * Drag model (piecewise-linear through the drag crisis):
- *   Re <= RE_THRESHOLD_LOW:                  Cd = CD_LOW
+ * Drag model (piecewise-linear through the drag crisis, continuous in Re):
+ *   Re <= RE_THRESHOLD_LOW:                  Cd = CD_LOW  + CD_SPIN * S
  *   RE_THRESHOLD_LOW < Re < RE_THRESHOLD_HIGH: linear + CD_SPIN * S
  *   Re >= RE_THRESHOLD_HIGH:                 Cd = CD_HIGH + CD_SPIN * S
  *
@@ -46,18 +46,23 @@ class DefaultAerodynamicModel : public AerodynamicModel
 {
 public:
 	// ========================================================================
-	// DRAG FORCE LUMPED-SCALAR DERIVATION (consumed by ShotPhysicsContext)
+	// DRAG FORCE LUMPED-SCALAR DERIVATION
 	// ========================================================================
-	// c0 = DRAG_FORCE_CONST * rho * (REF_BALL_MASS_OZ / mass) * (circ / REF_BALL_CIRC_IN)^2
+	// The lumped coefficient c0 is derived by ShotPhysicsContext from the
+	// reference-ball constants in physics_constants:
+	//   c0 = DRAG_FORCE_CONST * rho * (REF_BALL_MASS_OZ / mass) * (circ / REF_BALL_CIRC_IN)^2
+	// These aliases are retained for inspection by derived models; the single
+	// source of truth lives in physics_constants so the context layer need not
+	// depend on this concrete model.
 
 	/// Drag force constant. c0 = 0.07182 * rho * (5.125/mass) * (circ/9.125)²
-	static constexpr float DRAG_FORCE_CONST = 0.07182F;
+	static constexpr float DRAG_FORCE_CONST = physics_constants::DRAG_FORCE_CONST;
 
 	/// Reference golf ball mass for drag calculation (oz)
-	static constexpr float REF_BALL_MASS_OZ = 5.125F;
+	static constexpr float REF_BALL_MASS_OZ = physics_constants::REF_BALL_MASS_OZ;
 
 	/// Reference golf ball circumference for drag calculation (inches)
-	static constexpr float REF_BALL_CIRC_IN = 9.125F;
+	static constexpr float REF_BALL_CIRC_IN = physics_constants::REF_BALL_CIRC_IN;
 
 	// ========================================================================
 	// SPIN DECAY
@@ -158,7 +163,7 @@ class DefaultAerodynamicModel : public AerodynamicModel
 		const float vRelZ = state.velocity[2] - state.windVelocity[2];
 		const double vw = std::sqrt(static_cast<double>(vRelX * vRelX + vRelY * vRelY + vRelZ * vRelZ));
 
-		if (vw < static_cast<double>(physics_constants::MIN_VELOCITY_THRESHOLD))
+		if (vw < static_cast<double>(physics_constants::MIN_SPEED))
 		{
 			return {0.0F, 0.0F, 0.0F};
 		}
@@ -186,7 +191,7 @@ class DefaultAerodynamicModel : public AerodynamicModel
 
 		// Magnus: C0 * (Cl / omega) * vw * (spinVector × vRel)
 		float magnusX = 0.0F, magnusY = 0.0F, magnusZ = 0.0F;
-		if (omegaMag > static_cast<double>(physics_constants::MIN_VELOCITY_THRESHOLD))
+		if (omegaMag > static_cast<double>(physics_constants::MIN_SPIN))
 		{
 			const double magnusScale = static_cast<double>(state.c0) * (Cl / omegaMag) * vw;
 			magnusX = static_cast<float>(magnusScale * (omegaY * vRelZ - omegaZ * vRelY));
@@ -213,7 +218,9 @@ class DefaultAerodynamicModel : public AerodynamicModel
 
 		if (Re_x_e5 <= reLow)
 		{
-			return cdLow;
+			// Carry the spin term through the low-Re branch so Cd is continuous
+			// at reLow: the linear branch evaluates to cdLow + cdSpin*S there.
+			return cdLow + cdSpin * spinFactor;
 		}
 		else if (Re_x_e5 < reHigh)
 		{
diff --git a/include/DefaultBounceModel.hpp b/include/DefaultBounceModel.hpp
index 5d8d1cc..6f39e67 100644
--- a/include/DefaultBounceModel.hpp
+++ b/include/DefaultBounceModel.hpp
@@ -118,7 +118,7 @@ class DefaultBounceModel : public BounceModel
 		const Vector3D vNormalAfter = vNormal * -effectiveCor;
 
 		float impactAngle = 0.0F;
-		if (impactSpeed > physics_constants::MIN_VELOCITY_THRESHOLD)
+		if (impactSpeed > physics_constants::MIN_SPEED)
 		{
 			const float sinAngle = std::clamp(-vDotN / impactSpeed, -1.0F, 1.0F);
 			impactAngle = std::asin(sinAngle);
@@ -131,7 +131,7 @@ class DefaultBounceModel : public BounceModel
 		Vector3D vTangentAfter{};
 
 		if (steepImpact && energeticImpact &&
-		    tangentMag > physics_constants::MIN_VELOCITY_THRESHOLD)
+		    tangentMag > physics_constants::MIN_SPEED)
 		{
 			// Backspin scalar: positive when spin axis aligns with t̂ × n̂.
 			// For ball moving +Y, normal +Z: t̂ × n̂ = +X — matches golf
diff --git a/include/DefaultIntegrator.hpp b/include/DefaultIntegrator.hpp
new file mode 100644
index 0000000..7f156ac
--- /dev/null
+++ b/include/DefaultIntegrator.hpp
@@ -0,0 +1,35 @@
+#ifndef DEFAULT_INTEGRATOR_HPP
+#define DEFAULT_INTEGRATOR_HPP
+
+#include "Integrator.hpp"
+#include "physics_constants.hpp"
+
+/**
+ * @brief Built-in semi-implicit (symplectic) Euler integrator.
+ *
+ * Advances on the acceleration from the start of the step:
+ *   position += velocity * dt + 0.5 * a * dt²
+ *   velocity += a * dt
+ *
+ * Position carries a 2nd-order term; velocity is forward-Euler (1st-order
+ * accurate). This scheme needs only the start-of-step acceleration already in
+ * `state.acceleration`, so it does not sample the acceleration field.
+ */
+class DefaultIntegrator : public Integrator
+{
+public:
+	void step(BallState &state, float dt, const AccelerationField &accel) const override
+	{
+		(void)accel; // start-of-step acceleration is sufficient for this scheme
+
+		const Vector3D a = state.acceleration;
+		for (int i = 0; i < 3; ++i)
+		{
+			state.position[i] += state.velocity[i] * dt +
+			                     physics_constants::HALF * a[i] * dt * dt;
+			state.velocity[i] += a[i] * dt;
+		}
+	}
+};
+
+#endif // DEFAULT_INTEGRATOR_HPP
diff --git a/include/DefaultRollModel.hpp b/include/DefaultRollModel.hpp
index 910bc43..d57ee24 100644
--- a/include/DefaultRollModel.hpp
+++ b/include/DefaultRollModel.hpp
@@ -97,7 +97,7 @@ class DefaultRollModel : public RollModel
 
 		// Stationary: friction direction undefined, gravity along slope is
 		// balanced by static friction at rest.
-		if (vHorizontal < physics_constants::MIN_VELOCITY_THRESHOLD)
+		if (vHorizontal < physics_constants::MIN_SPEED)
 		{
 			return acceleration;
 		}
diff --git a/include/FlightPhase.hpp b/include/FlightPhase.hpp
index 40a179e..2038ef3 100644
--- a/include/FlightPhase.hpp
+++ b/include/FlightPhase.hpp
@@ -2,8 +2,10 @@
 #define FLIGHTPHASE_HPP
 
 #include "AerodynamicModel.hpp"
+#include "BallProperties.hpp"
 #include "BallState.hpp"
 #include "BounceModel.hpp"
+#include "Integrator.hpp"
 #include "RollModel.hpp"
 #include "ShotPhysicsContext.hpp"
 #include "atmospheric_data.hpp"
@@ -68,11 +70,13 @@ class AerialPhase : public FlightPhase
 	            const LaunchData &launch,
 	            const AtmosphericData &atmos,
 	            std::shared_ptr<TerrainInterface> terrain,
-	            std::shared_ptr<AerodynamicModel> model = nullptr);
+	            std::shared_ptr<AerodynamicModel> model = nullptr,
+	            const BallProperties &ball = {},
+	            float gravity = physics_constants::GRAVITY_FT_PER_S2,
+	            std::shared_ptr<Integrator> integrator = nullptr);
 
 	void initialize(BallState &state);
 	void calculateStep(BallState &state, float dt) override;
-	void calculateAccelerations(BallState &state);
 	bool isPhaseComplete(const BallState &state) const override;
 
 	// Getters for observable flight quantities (useful for testing and diagnostics)
@@ -88,6 +92,14 @@ class AerialPhase : public FlightPhase
 	AtmosphericData atmos;
 	std::shared_ptr<TerrainInterface> terrain;
 	std::shared_ptr<AerodynamicModel> model;
+	std::shared_ptr<Integrator> integrator;
+	float ballRadius;
+	float gravity;
+
+	// Acceleration field handed to the integrator. Built once and reused every
+	// step; stepSpin holds the spin (fixed across a step) the field samples with.
+	std::shared_ptr<Vector3D> stepSpin;
+	AccelerationField accelField;
 
 	// Cached scalar quantities derived from BallState each step
 	float v;
@@ -99,14 +111,10 @@ class AerialPhase : public FlightPhase
 	Vector3D velocity3D_w;
 
 	// Private calculation methods
-	void calculatePosition(BallState &state, float dt);
-	void calculateV(BallState &state, float dt);
 	void calculateVelocityw(const BallState &state);
 	void calculateTau(const BallState &state);
 	void calculateRw(const BallState &state);
 	void calculateAccel(BallState &state);
-
-	[[nodiscard]] AerodynamicState buildAerodynamicState(const BallState &state) const;
 };
 
 /**
@@ -125,15 +133,28 @@ class BouncePhase : public FlightPhase
 	            const AtmosphericData &atmos,
 	            std::shared_ptr<TerrainInterface> terrain,
 	            std::shared_ptr<AerodynamicModel> aeroModel = nullptr,
-	            std::shared_ptr<BounceModel> bounceModel = nullptr);
+	            std::shared_ptr<BounceModel> bounceModel = nullptr,
+	            const BallProperties &ball = {},
+	            float gravity = physics_constants::GRAVITY_FT_PER_S2,
+	            std::shared_ptr<Integrator> integrator = nullptr);
 
 	void calculateStep(BallState &state, float dt) override;
 	bool isPhaseComplete(const BallState &state) const override;
 
 private:
+	ShotPhysicsContext &physicsVars;
+	AtmosphericData atmos;
 	std::shared_ptr<TerrainInterface> terrain;
+	std::shared_ptr<AerodynamicModel> model;
 	std::shared_ptr<BounceModel> bounceModel;
-	AerialPhase aerialPhase; // Used for aerodynamic calculations between bounces
+	std::shared_ptr<Integrator> integrator;
+	float ballRadius;
+	float gravity;
+
+	// Acceleration field handed to the integrator. Built once and reused every
+	// step; stepSpin holds the spin (fixed across a step) the field samples with.
+	std::shared_ptr<Vector3D> stepSpin;
+	AccelerationField accelField;
 };
 
 /**
@@ -147,7 +168,8 @@ class RollPhase : public FlightPhase
 {
 public:
 	explicit RollPhase(std::shared_ptr<TerrainInterface> terrain,
-	                   std::shared_ptr<RollModel> model = nullptr);
+	                   std::shared_ptr<RollModel> model = nullptr,
+	                   const BallProperties &ball = {});
 
 	void calculateStep(BallState &state, float dt) override;
 	bool isPhaseComplete(const BallState &state) const override;
@@ -155,6 +177,7 @@ class RollPhase : public FlightPhase
 private:
 	std::shared_ptr<TerrainInterface> terrain;
 	std::shared_ptr<RollModel> model;
+	float ballRadius;
 	bool atRest = false;
 };
 
diff --git a/include/FlightSimulator.hpp b/include/FlightSimulator.hpp
index 603c375..6e2f833 100644
--- a/include/FlightSimulator.hpp
+++ b/include/FlightSimulator.hpp
@@ -1,9 +1,11 @@
 #ifndef FLIGHT_SIMULATOR_HPP
 #define FLIGHT_SIMULATOR_HPP
 
+#include "BallProperties.hpp"
 #include "BallState.hpp"
 #include "BounceModel.hpp"
 #include "FlightPhase.hpp"
+#include "Integrator.hpp"
 #include "RollModel.hpp"
 #include "ShotPhysicsContext.hpp"
 #include "atmospheric_data.hpp"
@@ -51,13 +53,20 @@ class FlightSimulator
 	 * @param aeroModel Aerodynamic coefficient model (nullptr uses DefaultAerodynamicModel)
 	 * @param bounceModel Bounce model (nullptr uses DefaultBounceModel)
 	 * @param rollModel Roll model (nullptr uses DefaultRollModel)
+	 * @param ball Ball properties (defaults to a standard golf ball)
+	 * @param gravity Gravitational acceleration in ft/s² (defaults to Earth);
+	 *        applied to the aerial and between-bounce flight integration
+	 * @param integrator Time integrator for the flight phases (nullptr uses DefaultIntegrator)
 	 */
 	FlightSimulator(const LaunchData &launch,
 	                const AtmosphericData &atmos,
 	                const GroundSurface &ground,
 	                std::shared_ptr<AerodynamicModel> aeroModel = nullptr,
 	                std::shared_ptr<BounceModel> bounceModel = nullptr,
-	                std::shared_ptr<RollModel> rollModel = nullptr);
+	                std::shared_ptr<RollModel> rollModel = nullptr,
+	                const BallProperties &ball = {},
+	                float gravity = physics_constants::GRAVITY_FT_PER_S2,
+	                std::shared_ptr<Integrator> integrator = nullptr);
 
 	/**
 	 * @brief Constructs a flight simulator with a custom terrain.
@@ -72,13 +81,20 @@ class FlightSimulator
 	 * @param aeroModel Aerodynamic coefficient model (nullptr uses DefaultAerodynamicModel)
 	 * @param bounceModel Bounce model (nullptr uses DefaultBounceModel)
 	 * @param rollModel Roll model (nullptr uses DefaultRollModel)
+	 * @param ball Ball properties (defaults to a standard golf ball)
+	 * @param gravity Gravitational acceleration in ft/s² (defaults to Earth);
+	 *        applied to the aerial and between-bounce flight integration
+	 * @param integrator Time integrator for the flight phases (nullptr uses DefaultIntegrator)
 	 */
 	FlightSimulator(const LaunchData &launch,
 	                const AtmosphericData &atmos,
 	                std::shared_ptr<TerrainInterface> terrain,
 	                std::shared_ptr<AerodynamicModel> aeroModel = nullptr,
 	                std::shared_ptr<BounceModel> bounceModel = nullptr,
-	                std::shared_ptr<RollModel> rollModel = nullptr);
+	                std::shared_ptr<RollModel> rollModel = nullptr,
+	                const BallProperties &ball = {},
+	                float gravity = physics_constants::GRAVITY_FT_PER_S2,
+	                std::shared_ptr<Integrator> integrator = nullptr);
 
 	/**
 	 * @brief Runs the simulation to completion.
@@ -149,6 +165,7 @@ class FlightSimulator
 	BallState state;
 
 	Vector3D startPosition_{0.0F, 0.0F, 0.0F};
+	float gravity_;
 
 	// Must be declared before phases since phases hold a reference to it
 	ShotPhysicsContext physicsVars_;
diff --git a/include/Integrator.hpp b/include/Integrator.hpp
new file mode 100644
index 0000000..ba5ec49
--- /dev/null
+++ b/include/Integrator.hpp
@@ -0,0 +1,81 @@
+#ifndef INTEGRATOR_HPP
+#define INTEGRATOR_HPP
+
+#include "BallState.hpp"
+#include "math_utils.hpp"
+
+#include <functional>
+
+/**
+ * @brief Kinematic sample the integrator evaluates the acceleration field at.
+ *
+ * A trial position/velocity/time that a higher-order or adaptive scheme may
+ * probe between the start and end of a step.
+ */
+struct IntegratorState
+{
+	Vector3D position; ///< Trial position (ft)
+	Vector3D velocity; ///< Trial velocity (ft/s)
+	float    time;     ///< Trial simulation time (s)
+};
+
+/**
+ * @brief Total acceleration field (aerodynamic + gravity, ft/s²) for the flight
+ *        phases, evaluated at an arbitrary trial state.
+ */
+using AccelerationField = std::function<Vector3D(const IntegratorState &)>;
+
+/**
+ * @brief Pluggable time integrator for the aerial and between-bounce flight.
+ *
+ * Implement this to replace the built-in semi-implicit Euler scheme with, e.g.,
+ * RK4 or an adaptive step. The flight phase owns spin decay, wind, and the
+ * acceleration model; the integrator owns only how position and velocity
+ * advance over a step.
+ *
+ * `step` receives the current ball state — with `state.acceleration` already
+ * holding the acceleration at the start of the step — and an @ref
+ * AccelerationField the implementation may sample at trial states. It must
+ * update `state.position` and `state.velocity` in place; the phase recomputes
+ * the start-of-step acceleration for the next step afterwards.
+ *
+ * @code
+ * class ForwardEuler : public Integrator {
+ * public:
+ *     void step(BallState& s, float dt, const AccelerationField&) const override {
+ *         const Vector3D a = s.acceleration;
+ *         for (int i = 0; i < 3; ++i) {
+ *             s.position[i] += s.velocity[i] * dt;
+ *             s.velocity[i] += a[i] * dt;
+ *         }
+ *     }
+ * };
+ * FlightSimulator sim(launch, atmos, ground, nullptr, nullptr, nullptr,
+ *                     BallProperties{}, physics_constants::GRAVITY_FT_PER_S2,
+ *                     std::make_shared<ForwardEuler>());
+ * @endcode
+ */
+class Integrator
+{
+public:
+	virtual ~Integrator() = default;
+
+	Integrator(const Integrator &) = delete;
+	Integrator &operator=(const Integrator &) = delete;
+	Integrator(Integrator &&) = delete;
+	Integrator &operator=(Integrator &&) = delete;
+
+	/**
+	 * @brief Advances position and velocity by one step of size dt.
+	 *
+	 * @param state Ball state to update; `state.acceleration` is the start-of-step value
+	 * @param dt Time step (s)
+	 * @param accel Acceleration field sampleable at trial states
+	 */
+	virtual void step(BallState &state, float dt, const AccelerationField &accel) const = 0;
+
+protected:
+	Integrator() = default;
+};
+
+#endif // INTEGRATOR_HPP
diff --git a/include/ShotPhysicsContext.hpp b/include/ShotPhysicsContext.hpp
index 15c459c..c99e292 100644
--- a/include/ShotPhysicsContext.hpp
+++ b/include/ShotPhysicsContext.hpp
@@ -1,6 +1,7 @@
 #ifndef SHOT_PHYSICS_CONTEXT_HPP
 #define SHOT_PHYSICS_CONTEXT_HPP
 
+#include "BallProperties.hpp"
 #include "atmospheric_data.hpp"
 #include "launch_data.hpp"
 #include "math_utils.hpp"
@@ -9,7 +10,8 @@ class ShotPhysicsContext
 {
 public:
 	ShotPhysicsContext(const LaunchData &launch,
-							 const AtmosphericData &atmos);
+							 const AtmosphericData &atmos,
+							 const BallProperties &ball = {});
 
 	// Getters
 	// Mixed unit convention, matching AerodynamicState: imperial for kinematics
@@ -59,6 +61,7 @@ class ShotPhysicsContext
 	void calculateAllVariables();
 	LaunchData launch;
 	AtmosphericData atmos;
+	BallProperties ball;
 
 	// Member variables
 	float rhoImperial = 0.0F;
diff --git a/include/atmospheric_data.hpp b/include/atmospheric_data.hpp
index 9724a7c..6bc6496 100644
--- a/include/atmospheric_data.hpp
+++ b/include/atmospheric_data.hpp
@@ -2,7 +2,11 @@
 #define ATMOSPHERIC_DATA_HPP
 
 /**
- * @brief Atmospheric conditions for flight simulation
+ * @brief Atmospheric conditions for flight simulation.
+ *
+ * Defaults describe a sea-level standard day (59°F, 29.92 inHg, no wind, dry
+ * air), so a default-constructed AtmosphericData{} is a usable baseline and
+ * designated initializers need only override the fields that differ.
  */
 struct AtmosphericData
 {
@@ -12,7 +16,7 @@ struct AtmosphericData
 	 * Affects air density and ball flight.
 	 * Typical range: 0-120°F (Earth conditions)
 	 */
-	float temp;
+	float temp = 59.0F;
 
 	/**
 	 * @brief Elevation above sea level (in feet).
@@ -20,7 +24,7 @@ struct AtmosphericData
 	 * Higher elevations have lower air density.
 	 * Typical range: -500 to 15000 ft (most golf courses 0-8000 ft)
 	 */
-	float elevation;
+	float elevation = 0.0F;
 
 	/**
 	 * @brief Wind speed (in mph).
@@ -28,7 +32,7 @@ struct AtmosphericData
 	 * Magnitude of wind velocity.
 	 * Typical range: 0-40 mph
 	 */
-	float vWind;
+	float vWind = 0.0F;
 
 	/**
 	 * @brief Wind direction (in degrees).
@@ -46,14 +50,14 @@ struct AtmosphericData
 	 *
 	 * Range: -180 to 180 deg
 	 */
-	float phiWind;
+	float phiWind = 0.0F;
 
 	/**
 	 * @brief Height at which wind acts (in feet).
 	 *
 	 * Wind affects ball above this altitude.
 	 */
-	float hWind;
+	float hWind = 0.0F;
 
 	/**
 	 * @brief Relative humidity (in percent).
@@ -61,7 +65,7 @@ struct AtmosphericData
 	 * Affects air density slightly.
 	 * Range: 0-100%
 	 */
-	float relHumidity;
+	float relHumidity = 0.0F;
 
 	/**
 	 * @brief Barometric pressure (in inches of mercury).
@@ -69,7 +73,7 @@ struct AtmosphericData
 	 * Standard sea level pressure is 29.92 inHg.
 	 * Typical range: 28-31 inHg
 	 */
-	float pressure;
+	float pressure = 29.92F;
 };
 
-#endif // ATMOSPHERIC_DATA_HPP
\ No newline at end of file
+#endif // ATMOSPHERIC_DATA_HPP
diff --git a/include/libgolf.hpp b/include/libgolf.hpp
index d2f1168..bb48d3e 100644
--- a/include/libgolf.hpp
+++ b/include/libgolf.hpp
@@ -2,21 +2,24 @@
 #define LIBGOLF_HPP
 
 #include "AerodynamicModel.hpp"
+#include "BallProperties.hpp"
+#include "BallState.hpp"
 #include "BounceModel.hpp"
 #include "DefaultAerodynamicModel.hpp"
 #include "DefaultBounceModel.hpp"
+#include "DefaultIntegrator.hpp"
 #include "DefaultRollModel.hpp"
-#include "RollModel.hpp"
-#include "FlightSimulator.hpp"
 #include "FlightPhase.hpp"
+#include "FlightSimulator.hpp"
+#include "Integrator.hpp"
+#include "RollModel.hpp"
 #include "ShotPhysicsContext.hpp"
-#include "BallState.hpp"
 #include "atmospheric_data.hpp"
-#include "launch_data.hpp"
-#include "ground_surface.hpp"
-#include "terrain_interface.hpp"
 #include "ground_physics.hpp"
+#include "ground_surface.hpp"
+#include "launch_data.hpp"
 #include "math_utils.hpp"
 #include "physics_constants.hpp"
+#include "terrain_interface.hpp"
 
 #endif // LIBGOLF_HPP
\ No newline at end of file
diff --git a/include/physics_constants.hpp b/include/physics_constants.hpp
index 01c1572..93c4ec2 100644
--- a/include/physics_constants.hpp
+++ b/include/physics_constants.hpp
@@ -30,6 +30,23 @@ namespace physics_constants
     /// Standard golf ball mass (ounces)
     constexpr float STD_BALL_MASS_OZ = 1.62F;
 
+    // ------------------------------------------------------------------------
+    // Lumped drag-coefficient (c0) derivation inputs
+    // ------------------------------------------------------------------------
+    // c0 = DRAG_FORCE_CONST * rho * (REF_BALL_MASS_OZ / mass) * (circ / REF_BALL_CIRC_IN)^2
+    // These describe the reference ball the empirical drag fit was taken
+    // against; they are an input to the derivation, not a tunable of any one
+    // model, so they live here rather than on a concrete model class.
+
+    /// Drag force constant for the c0 derivation
+    constexpr float DRAG_FORCE_CONST = 0.07182F;
+
+    /// Reference golf ball mass for the c0 derivation (oz)
+    constexpr float REF_BALL_MASS_OZ = 5.125F;
+
+    /// Reference golf ball circumference for the c0 derivation (inches)
+    constexpr float REF_BALL_CIRC_IN = 9.125F;
+
     // ========================================================================
     // FUNDAMENTAL PHYSICAL CONSTANTS
     // ========================================================================
@@ -188,10 +205,18 @@ namespace physics_constants
     // ========================================================================
     // NUMERICAL STABILITY THRESHOLDS
     // ========================================================================
+    // Dimensionally distinct floors guarding against division by zero and
+    // degenerate cases. Equal in value today, but split by dimension so each
+    // can be retuned independently.
+
+    /// Minimum speed magnitude before a velocity is treated as zero (ft/s)
+    constexpr float MIN_SPEED = 0.01F;
+
+    /// Minimum spin magnitude before spin is treated as zero (rad/s)
+    constexpr float MIN_SPIN = 0.01F;
 
-    /// Minimum velocity magnitude to avoid division by zero in calculations (ft/s)
-    /// Used in spin factor and friction calculations
-    constexpr float MIN_VELOCITY_THRESHOLD = 0.01F;
+    /// Minimum vector length before it is treated as zero (ft)
+    constexpr float MIN_LENGTH = 0.01F;
 
     // ========================================================================
     // PHASE TRANSITION THRESHOLDS
diff --git a/src/FlightPhase.cpp b/src/FlightPhase.cpp
index 2c0ac69..a4a4804 100644
--- a/src/FlightPhase.cpp
+++ b/src/FlightPhase.cpp
@@ -12,6 +12,7 @@
 #include "FlightPhase.hpp"
 #include "DefaultAerodynamicModel.hpp"
 #include "DefaultBounceModel.hpp"
+#include "DefaultIntegrator.hpp"
 #include "DefaultRollModel.hpp"
 #include "atmospheric_data.hpp"
 #include "math_utils.hpp"
@@ -30,6 +31,72 @@ namespace
 	{
 		return p ? std::move(p) : std::make_shared<Default>();
 	}
+
+	// Effective wind at the ball: the launch wind above hWind, otherwise none.
+	Vector3D effectiveWind(const ShotPhysicsContext &physicsVars,
+	                       const AtmosphericData &atmos, float heightZ)
+	{
+		if (heightZ >= atmos.hWind)
+		{
+			return {physicsVars.getVw()[0], physicsVars.getVw()[1], 0.0F};
+		}
+		return {0.0F, 0.0F, 0.0F};
+	}
+
+	// Snapshot the ball + launch atmosphere into the model-facing state.
+	AerodynamicState buildAeroState(const BallState &state,
+	                                const ShotPhysicsContext &physicsVars,
+	                                const AtmosphericData &atmos, float ballRadius)
+	{
+		const Vector3D wind = effectiveWind(physicsVars, atmos, state.position[2]);
+		return AerodynamicState{
+		    .velocity          = state.velocity,
+		    .windVelocity      = {wind[0], wind[1], 0.0F}, // vertical wind not modelled
+		    .spinVector        = state.spinVector,
+		    .position          = state.position,
+		    .currentTime       = state.currentTime,
+		    .ballRadius        = ballRadius,
+		    .airDensityKgPerM3 = physicsVars.getRhoMetric(),
+		    .airViscosity      = physicsVars.getAirViscosity(),
+		    .tempKelvin        = physicsVars.getTempKelvin(),
+		    .pressureMmHg      = physicsVars.getBarometricPressure(),
+		    .relHumidity       = physicsVars.getRelHumidity(),
+		    .c0                = physicsVars.getC0(),
+		    .re100             = physicsVars.getRe100()
+		};
+	}
+
+	// Aerodynamic acceleration plus gravity, shared by the aerial and
+	// between-bounce flight steps.
+	Vector3D flightAcceleration(const AerodynamicModel &model,
+	                            const AerodynamicState &aeroState, float gravity)
+	{
+		const Vector3D aero = model.computeAcceleration(aeroState);
+		return {aero[0], aero[1], aero[2] - gravity};
+	}
+
+	// Acceleration field an Integrator can sample at trial states, built once per
+	// phase and reused every step. The per-step spin (decayed once, before the
+	// step, then held fixed) is read through stepSpin so the field need not be
+	// rebuilt. Captures are move-safe: the model is held by shared_ptr, atmos by
+	// value, and physicsVars outlives the phase, so a moved-from phase does not
+	// dangle the field.
+	AccelerationField makeAccelField(std::shared_ptr<AerodynamicModel> model,
+	                                 const ShotPhysicsContext &physicsVars,
+	                                 const AtmosphericData &atmos, float ballRadius,
+	                                 float gravity, std::shared_ptr<Vector3D> stepSpin)
+	{
+		return [model = std::move(model), &physicsVars, atmos, ballRadius, gravity,
+		        stepSpin = std::move(stepSpin)](const IntegratorState &s) -> Vector3D {
+			BallState trial;
+			trial.position    = s.position;
+			trial.velocity    = s.velocity;
+			trial.spinVector  = *stepSpin;
+			trial.currentTime = s.time;
+			return flightAcceleration(
+			    *model, buildAeroState(trial, physicsVars, atmos, ballRadius), gravity);
+		};
+	}
 }
 
 // ============================================================================
@@ -39,9 +106,14 @@ namespace
 AerialPhase::AerialPhase(
 	ShotPhysicsContext &physicsVars, [[maybe_unused]] const LaunchData &launch,
 	const AtmosphericData &atmos, std::shared_ptr<TerrainInterface> terrain,
-	std::shared_ptr<AerodynamicModel> model)
+	std::shared_ptr<AerodynamicModel> model, const BallProperties &ball,
+	float gravity, std::shared_ptr<Integrator> integrator)
 	: physicsVars(physicsVars), atmos(atmos), terrain(terrain),
-	  model(orDefault<DefaultAerodynamicModel>(std::move(model)))
+	  model(orDefault<DefaultAerodynamicModel>(std::move(model))),
+	  integrator(orDefault<DefaultIntegrator>(std::move(integrator))),
+	  ballRadius(ball.radiusFt()), gravity(gravity),
+	  stepSpin(std::make_shared<Vector3D>()),
+	  accelField(makeAccelField(this->model, physicsVars, atmos, ballRadius, gravity, stepSpin))
 {
 	if (!terrain)
 	{
@@ -59,7 +131,7 @@ AerialPhase::AerialPhase(
 
 void AerialPhase::initialize(BallState &state)
 {
-	if (math_utils::magnitude(state.spinVector) < physics_constants::MIN_VELOCITY_THRESHOLD)
+	if (math_utils::magnitude(state.spinVector) < physics_constants::MIN_SPIN)
 	{
 		state.spinVector = physicsVars.getW();
 	}
@@ -76,19 +148,6 @@ void AerialPhase::initialize(BallState &state)
 	calculateAccel(state);
 }
 
-void AerialPhase::calculateAccelerations(BallState &state)
-{
-	v    = std::sqrt(state.velocity[0] * state.velocity[0] +
-	                 state.velocity[1] * state.velocity[1] +
-	                 state.velocity[2] * state.velocity[2]);
-	vMph = v / physics_constants::MPH_TO_FT_PER_S;
-
-	calculateVelocityw(state);
-	calculateTau(state);
-	calculateRw(state);
-	calculateAccel(state);
-}
-
 void AerialPhase::calculateStep(BallState &state, float dt)
 {
 	state.currentTime += dt;
@@ -102,11 +161,17 @@ void AerialPhase::calculateStep(BallState &state, float dt)
 	state.spinVector[1] *= decay;
 	state.spinVector[2] *= decay;
 
-	calculatePosition(state, dt);
-	calculateV(state, dt);       // updates v, vMph, state.velocity
+	// state.acceleration holds the start-of-step acceleration; advance position
+	// and velocity through the (pluggable) integrator. The field samples the
+	// just-decayed spin held in stepSpin.
+	*stepSpin = state.spinVector;
+	integrator->step(state, dt, accelField);
+
+	v    = math_utils::magnitude(state.velocity);
+	vMph = v / physics_constants::MPH_TO_FT_PER_S;
 	calculateVelocityw(state);   // updates velocity3D_w, vw, vwMph
 	calculateRw(state);
-	calculateAccel(state);       // calls model, then adds gravity
+	calculateAccel(state);       // start-of-step acceleration for the next step
 }
 
 bool AerialPhase::isPhaseComplete(const BallState &state) const
@@ -115,28 +180,6 @@ bool AerialPhase::isPhaseComplete(const BallState &state) const
 	return state.position[2] <= terrainHeight;
 }
 
-void AerialPhase::calculatePosition(BallState &state, float dt)
-{
-	state.position[0] = state.position[0] + state.velocity[0] * dt +
-						physics_constants::HALF * state.acceleration[0] * dt * dt;
-	state.position[1] = state.position[1] + state.velocity[1] * dt +
-						physics_constants::HALF * state.acceleration[1] * dt * dt;
-	state.position[2] = state.position[2] + state.velocity[2] * dt +
-						physics_constants::HALF * state.acceleration[2] * dt * dt;
-}
-
-void AerialPhase::calculateV(BallState &state, float dt)
-{
-	float vx = state.velocity[0] + state.acceleration[0] * dt;
-	float vy = state.velocity[1] + state.acceleration[1] * dt;
-	float vz = state.velocity[2] + state.acceleration[2] * dt;
-
-	state.velocity = {vx, vy, vz};
-
-	v    = std::sqrt(vx * vx + vy * vy + vz * vz);
-	vMph = v / physics_constants::MPH_TO_FT_PER_S;
-}
-
 void AerialPhase::calculateVelocityw(const BallState &state)
 {
 	if (state.position[2] >= atmos.hWind)
@@ -159,13 +202,13 @@ void AerialPhase::calculateVelocityw(const BallState &state)
 
 void AerialPhase::calculateTau(const BallState &state)
 {
-	if (v < physics_constants::MIN_VELOCITY_THRESHOLD)
+	if (v < physics_constants::MIN_SPEED)
 	{
 		tau = 1e6F;
 		return;
 	}
 
-	tau = model->computeSpinDecayTau(buildAerodynamicState(state));
+	tau = model->computeSpinDecayTau(buildAeroState(state, physicsVars, atmos, ballRadius));
 }
 
 void AerialPhase::calculateRw(const BallState &state)
@@ -175,29 +218,9 @@ void AerialPhase::calculateRw(const BallState &state)
 
 void AerialPhase::calculateAccel(BallState &state)
 {
-	Vector3D aeroAccel = model->computeAcceleration(buildAerodynamicState(state));
-	state.acceleration[0] = aeroAccel[0];
-	state.acceleration[1] = aeroAccel[1];
-	state.acceleration[2] = aeroAccel[2] - physics_constants::GRAVITY_FT_PER_S2;
-}
-
-AerodynamicState AerialPhase::buildAerodynamicState(const BallState &state) const
-{
-	return AerodynamicState{
-	    .velocity          = state.velocity,
-	    .windVelocity      = {velocity3D_w[0], velocity3D_w[1], 0.0F}, // vertical wind not modelled
-	    .spinVector        = state.spinVector,
-	    .position          = state.position,
-	    .currentTime       = state.currentTime,
-	    .ballRadius        = physics_constants::STD_BALL_RADIUS_FT,
-	    .airDensityKgPerM3 = physicsVars.getRhoMetric(),
-	    .airViscosity      = physicsVars.getAirViscosity(),
-	    .tempKelvin        = physicsVars.getTempKelvin(),
-	    .pressureMmHg      = physicsVars.getBarometricPressure(),
-	    .relHumidity       = physicsVars.getRelHumidity(),
-	    .c0                = physicsVars.getC0(),
-	    .re100             = physicsVars.getRe100()
-	};
+	*stepSpin = state.spinVector;
+	state.acceleration =
+	    accelField(IntegratorState{state.position, state.velocity, state.currentTime});
 }
 
 // ============================================================================
@@ -205,13 +228,19 @@ AerodynamicState AerialPhase::buildAerodynamicState(const BallState &state) cons
 // ============================================================================
 
 BouncePhase::BouncePhase(
-	ShotPhysicsContext &physicsVars, const LaunchData &launch,
+	ShotPhysicsContext &physicsVars, [[maybe_unused]] const LaunchData &launch,
 	const AtmosphericData &atmos, std::shared_ptr<TerrainInterface> terrain,
 	std::shared_ptr<AerodynamicModel> aeroModel,
-	std::shared_ptr<BounceModel> bounceModel)
-	: terrain(terrain),
+	std::shared_ptr<BounceModel> bounceModel,
+	const BallProperties &ball,
+	float gravity, std::shared_ptr<Integrator> integrator)
+	: physicsVars(physicsVars), atmos(atmos), terrain(terrain),
+	  model(orDefault<DefaultAerodynamicModel>(std::move(aeroModel))),
 	  bounceModel(orDefault<DefaultBounceModel>(std::move(bounceModel))),
-	  aerialPhase(physicsVars, launch, atmos, terrain, std::move(aeroModel))
+	  integrator(orDefault<DefaultIntegrator>(std::move(integrator))),
+	  ballRadius(ball.radiusFt()), gravity(gravity),
+	  stepSpin(std::make_shared<Vector3D>()),
+	  accelField(makeAccelField(this->model, physicsVars, atmos, ballRadius, gravity, stepSpin))
 {
 	if (!terrain)
 	{
@@ -235,7 +264,7 @@ void BouncePhase::calculateStep(BallState &state, float dt)
 			state.velocity,
 			surfaceNormal,
 			state.spinVector,
-			physics_constants::STD_BALL_RADIUS_FT
+			ballRadius
 		};
 		BounceResult result = bounceModel->resolveBounce(bounceState, surface);
 		state.velocity = result.newVelocity;
@@ -243,15 +272,15 @@ void BouncePhase::calculateStep(BallState &state, float dt)
 		state.position[2] = terrainHeight;
 	}
 
-	aerialPhase.calculateAccelerations(state);
-
-	state.position[0] += state.velocity[0] * dt + 0.5F * state.acceleration[0] * dt * dt;
-	state.position[1] += state.velocity[1] * dt + 0.5F * state.acceleration[1] * dt * dt;
-	state.position[2] += state.velocity[2] * dt + 0.5F * state.acceleration[2] * dt * dt;
+	// Aerodynamic acceleration for the free-flight arc between bounces, advanced
+	// through the same (pluggable) integrator the aerial phase uses. Both the
+	// start-of-step acceleration and the integrator's trial samples come from the
+	// one cached field, so they cannot drift apart. Spin is fixed across the arc.
+	*stepSpin = state.spinVector;
+	state.acceleration =
+	    accelField(IntegratorState{state.position, state.velocity, state.currentTime});
 
-	state.velocity[0] += state.acceleration[0] * dt;
-	state.velocity[1] += state.acceleration[1] * dt;
-	state.velocity[2] += state.acceleration[2] * dt;
+	integrator->step(state, dt, accelField);
 
 	state.currentTime += dt;
 
@@ -277,9 +306,11 @@ bool BouncePhase::isPhaseComplete(const BallState &state) const
 
 RollPhase::RollPhase(
 	std::shared_ptr<TerrainInterface> terrain,
-	std::shared_ptr<RollModel> model)
+	std::shared_ptr<RollModel> model,
+	const BallProperties &ball)
 	: terrain(terrain),
-	  model(orDefault<DefaultRollModel>(std::move(model)))
+	  model(orDefault<DefaultRollModel>(std::move(model))),
+	  ballRadius(ball.radiusFt())
 {
 	if (!terrain)
 	{
@@ -297,7 +328,7 @@ void RollPhase::calculateStep(BallState &state, float dt)
 		state.velocity,
 		state.spinVector,
 		surfaceNormal,
-		physics_constants::STD_BALL_RADIUS_FT,
+		ballRadius,
 		dt,
 		terrain.get()
 	};
diff --git a/src/FlightSimulator.cpp b/src/FlightSimulator.cpp
index b230d30..78fd9ef 100644
--- a/src/FlightSimulator.cpp
+++ b/src/FlightSimulator.cpp
@@ -24,13 +24,17 @@ FlightSimulator::FlightSimulator(
 	const GroundSurface &ground,
 	std::shared_ptr<AerodynamicModel> aeroModel,
 	std::shared_ptr<BounceModel> bounceModel,
-	std::shared_ptr<RollModel> rollModel)
+	std::shared_ptr<RollModel> rollModel,
+	const BallProperties &ball,
+	float gravity,
+	std::shared_ptr<Integrator> integrator)
 	: currentPhase(Phase::Aerial),
-	  physicsVars_(launch, atmos),
+	  gravity_(gravity),
+	  physicsVars_(launch, atmos, ball),
 	  terrainStorage_(std::make_shared<FlatTerrain>(ground)),
-	  aerialPhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel),
-	  bouncePhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, bounceModel),
-	  rollPhase(terrainStorage_, rollModel)
+	  aerialPhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, ball, gravity, integrator),
+	  bouncePhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, bounceModel, ball, gravity, integrator),
+	  rollPhase(terrainStorage_, rollModel, ball)
 {
 	initializeFromLaunch(launch);
 }
@@ -41,13 +45,17 @@ FlightSimulator::FlightSimulator(
 	std::shared_ptr<TerrainInterface> terrain,
 	std::shared_ptr<AerodynamicModel> aeroModel,
 	std::shared_ptr<BounceModel> bounceModel,
-	std::shared_ptr<RollModel> rollModel)
+	std::shared_ptr<RollModel> rollModel,
+	const BallProperties &ball,
+	float gravity,
+	std::shared_ptr<Integrator> integrator)
 	: currentPhase(Phase::Aerial),
-	  physicsVars_(launch, atmos),
+	  gravity_(gravity),
+	  physicsVars_(launch, atmos, ball),
 	  terrainStorage_(terrain),
-	  aerialPhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel),
-	  bouncePhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, bounceModel),
-	  rollPhase(terrainStorage_, rollModel)
+	  aerialPhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, ball, gravity, integrator),
+	  bouncePhase(physicsVars_, launch, atmos, terrainStorage_, aeroModel, bounceModel, ball, gravity, integrator),
+	  rollPhase(terrainStorage_, rollModel, ball)
 {
 	initializeFromLaunch(launch);
 }
@@ -63,7 +71,7 @@ void FlightSimulator::initializeFromLaunch(const LaunchData &launch)
 		launch.launchAngleDeg,
 		launch.directionDeg,
 		startPos,
-		physics_constants::GRAVITY_FT_PER_S2,
+		gravity_,
 		physicsVars_.getW());
 
 	aerialPhase.initialize(state);
diff --git a/src/ShotPhysicsContext.cpp b/src/ShotPhysicsContext.cpp
index c10f03f..d6a30f8 100644
--- a/src/ShotPhysicsContext.cpp
+++ b/src/ShotPhysicsContext.cpp
@@ -15,7 +15,6 @@
  */
 
 #include "ShotPhysicsContext.hpp"
-#include "DefaultAerodynamicModel.hpp"
 #include "launch_data.hpp"
 #include "atmospheric_data.hpp"
 #include "math_utils.hpp"
@@ -34,8 +33,8 @@
  *       physically reasonable ranges. Passing invalid or out-of-range data may lead to
  *       unexpected behavior or incorrect calculations.
  */
-ShotPhysicsContext::ShotPhysicsContext(const LaunchData &launch, const AtmosphericData &atmos)
-    : launch(launch), atmos(atmos)
+ShotPhysicsContext::ShotPhysicsContext(const LaunchData &launch, const AtmosphericData &atmos, const BallProperties &ball)
+    : launch(launch), atmos(atmos), ball(ball)
 {
     tempC = math_utils::convertFahrenheitToCelsius(atmos.temp);
     elevationM = math_utils::convertFeetToMeters(atmos.elevation);
@@ -77,8 +76,8 @@ void ShotPhysicsContext::calculateRhoImperial()
 
 void ShotPhysicsContext::calculateC0()
 {
-    c0 = DefaultAerodynamicModel::DRAG_FORCE_CONST * rhoImperial * (DefaultAerodynamicModel::REF_BALL_MASS_OZ / physics_constants::STD_BALL_MASS_OZ) *
-         std::pow(physics_constants::STD_BALL_CIRCUMFERENCE_IN / DefaultAerodynamicModel::REF_BALL_CIRC_IN, 2);
+    c0 = physics_constants::DRAG_FORCE_CONST * rhoImperial * (physics_constants::REF_BALL_MASS_OZ / ball.massOz) *
+         std::pow(ball.circumferenceIn / physics_constants::REF_BALL_CIRC_IN, 2);
 }
 
 void ShotPhysicsContext::calculateV0()
@@ -109,7 +108,7 @@ void ShotPhysicsContext::calculateOmega()
 
 void ShotPhysicsContext::calculateROmega()
 {
-    rOmega = (physics_constants::STD_BALL_CIRCUMFERENCE_IN / (2 * physics_constants::PI)) * (omega / physics_constants::INCHES_PER_FOOT);
+    rOmega = (ball.circumferenceIn / (2 * physics_constants::PI)) * (omega / physics_constants::INCHES_PER_FOOT);
 }
 
 void ShotPhysicsContext::calculateVw()
@@ -140,7 +139,7 @@ void ShotPhysicsContext::calculateAirViscosity()
 void ShotPhysicsContext::calculateRe100()
 {
     // Re = ρ · v · D / μ, evaluated at v = RE100_VELOCITY_M_PER_S.
-    const float diameterM = physics_constants::STD_BALL_CIRCUMFERENCE_IN /
+    const float diameterM = ball.circumferenceIn /
                             (physics_constants::PI * physics_constants::INCHES_PER_METER);
     Re100 = rhoMetric * physics_constants::RE100_VELOCITY_M_PER_S * diameterM / airViscosity;
 }
diff --git a/src/math_utils.cpp b/src/math_utils.cpp
index 9653dff..e0675f2 100644
--- a/src/math_utils.cpp
+++ b/src/math_utils.cpp
@@ -135,7 +135,7 @@ Vector3D math_utils::normalize(const Vector3D& v)
 {
     float mag = magnitude(v);
 
-    if (mag < physics_constants::MIN_VELOCITY_THRESHOLD)
+    if (mag < physics_constants::MIN_LENGTH)
     {
         throw std::invalid_argument(std::string(__func__) + ": Cannot normalize zero-length vector");
     }
@@ -155,7 +155,7 @@ Vector3D math_utils::project(const Vector3D& v, const Vector3D& onto)
 {
     float ontoMagSquared = dot(onto, onto);
 
-    if (ontoMagSquared < physics_constants::MIN_VELOCITY_THRESHOLD * physics_constants::MIN_VELOCITY_THRESHOLD)
+    if (ontoMagSquared < physics_constants::MIN_LENGTH * physics_constants::MIN_LENGTH)
     {
         throw std::invalid_argument(std::string(__func__) + ": Cannot project onto zero-length vector");
     }
diff --git a/test/test_aerodynamic_model.cpp b/test/test_aerodynamic_model.cpp
index 1a63239..5fbc734 100644
--- a/test/test_aerodynamic_model.cpp
+++ b/test/test_aerodynamic_model.cpp
@@ -30,16 +30,36 @@ class DefaultModelTest : public ::testing::Test
 
 TEST_F(DefaultModelTest, CdBelowLowReThreshold)
 {
-	// Re_x_e5 < RE_THRESHOLD_LOW (0.5) → CD_LOW regardless of spin
+	// Re_x_e5 < RE_THRESHOLD_LOW (0.5), no spin → CD_LOW
 	EXPECT_NEAR(model.computeCd(0.25, 0.0), DefaultAerodynamicModel::CD_LOW, 1e-6);
 }
 
 TEST_F(DefaultModelTest, CdAtLowReThresholdIsInclusive)
 {
-	// Re_x_e5 == RE_THRESHOLD_LOW uses <= branch → still CD_LOW
+	// Re_x_e5 == RE_THRESHOLD_LOW uses <= branch; no spin → CD_LOW
 	EXPECT_NEAR(model.computeCd(0.5, 0.0), DefaultAerodynamicModel::CD_LOW, 1e-6);
 }
 
+TEST_F(DefaultModelTest, CdLowReCarriesSpinTerm)
+{
+	// The low-Re branch includes CD_SPIN * S, matching the linear branch at
+	// the boundary so Cd is continuous there.
+	const double S = 0.2;
+	const double expected = static_cast<double>(DefaultAerodynamicModel::CD_LOW) +
+	                        static_cast<double>(DefaultAerodynamicModel::CD_SPIN) * S;
+	EXPECT_NEAR(model.computeCd(0.25, S), expected, 1e-6);
+}
+
+TEST_F(DefaultModelTest, CdIsContinuousAcrossLowReThreshold)
+{
+	// No step at RE_THRESHOLD_LOW: approaching 0.5 from below and from above
+	// converges to the same value even with spin.
+	const double S = 0.2;
+	const double below = model.computeCd(0.5 - 1e-6, S);
+	const double above = model.computeCd(0.5 + 1e-6, S);
+	EXPECT_NEAR(below, above, 1e-4);
+}
+
 TEST_F(DefaultModelTest, CdMidRangeNoSpin)
 {
 	// Re = 0.75, S = 0:
diff --git a/test/test_ball_physics_vars.cpp b/test/test_ball_physics_vars.cpp
index 336f096..bccee16 100644
--- a/test/test_ball_physics_vars.cpp
+++ b/test/test_ball_physics_vars.cpp
@@ -87,3 +87,69 @@ TEST(GolfTest, initVarsNotDefault)
   EXPECT_NEAR(vars.getBarometricPressure(), 759.97, 0.1);
   EXPECT_NEAR(vars.getRe100(), 123200, 100);
 }
+
+TEST(GolfTest, defaultAtmosphereIsStandardDay)
+{
+  const AtmosphericData atmos{};
+
+  // Sea-level standard day: 59°F, 29.92 inHg, no wind, dry air.
+  EXPECT_FLOAT_EQ(atmos.temp, 59.0f);
+  EXPECT_FLOAT_EQ(atmos.elevation, 0.0f);
+  EXPECT_FLOAT_EQ(atmos.vWind, 0.0f);
+  EXPECT_FLOAT_EQ(atmos.phiWind, 0.0f);
+  EXPECT_FLOAT_EQ(atmos.hWind, 0.0f);
+  EXPECT_FLOAT_EQ(atmos.relHumidity, 0.0f);
+  EXPECT_FLOAT_EQ(atmos.pressure, 29.92f);
+
+  // Default-constructed atmosphere is usable and yields a sane air density.
+  const LaunchData launch{
+      .ballSpeedMph = 160.0f,
+      .launchAngleDeg = 11.0f,
+      .directionDeg = 0.0f,
+      .backspinRpm = 3000.0f,
+      .sidespinRpm = 0.0f,
+  };
+  const ShotPhysicsContext vars(launch, atmos);
+  EXPECT_NEAR(vars.getRhoMetric(), 1.225, 0.02); // ~ISA sea-level density
+  EXPECT_GT(vars.getC0(), 0.0f);
+}
+
+TEST(GolfTest, ballPropertiesThreadIntoDerivation)
+{
+  const LaunchData launch{
+      .ballSpeedMph = 160.0f,
+      .launchAngleDeg = 11.0f,
+      .directionDeg = 0.0f,
+      .backspinRpm = 3000.0f,
+      .sidespinRpm = 0.0f,
+  };
+  const AtmosphericData atmos{
+      .temp = 70.0f,
+      .elevation = 0.0f,
+      .vWind = 0.0f,
+      .phiWind = 0.0f,
+      .hWind = 0.0f,
+      .relHumidity = 50.0f,
+      .pressure = 29.92f,
+  };
+
+  // A default-constructed BallProperties must reproduce the standard ball.
+  const ShotPhysicsContext implicitDefault(launch, atmos);
+  const ShotPhysicsContext explicitDefault(launch, atmos, BallProperties{});
+  EXPECT_FLOAT_EQ(implicitDefault.getC0(), explicitDefault.getC0());
+  EXPECT_FLOAT_EQ(implicitDefault.getROmega(), explicitDefault.getROmega());
+  EXPECT_FLOAT_EQ(implicitDefault.getRe100(), explicitDefault.getRe100());
+
+  // c0 scales as 1/mass: a heavier ball of the same size drags less per unit speed.
+  const BallProperties heavy{.massOz = 2.0f};
+  const ShotPhysicsContext heavyVars(launch, atmos, heavy);
+  EXPECT_LT(heavyVars.getC0(), implicitDefault.getC0());
+  EXPECT_FLOAT_EQ(heavyVars.getROmega(), implicitDefault.getROmega());
+
+  // A larger circumference raises c0 (∝ area), surface speed, and the Reynolds reference.
+  const BallProperties big{.circumferenceIn = 6.0f};
+  const ShotPhysicsContext bigVars(launch, atmos, big);
+  EXPECT_GT(bigVars.getC0(), implicitDefault.getC0());
+  EXPECT_GT(bigVars.getROmega(), implicitDefault.getROmega());
+  EXPECT_GT(bigVars.getRe100(), implicitDefault.getRe100());
+}
diff --git a/test/test_flight_simulator.cpp b/test/test_flight_simulator.cpp
index c899989..36bcf4e 100644
--- a/test/test_flight_simulator.cpp
+++ b/test/test_flight_simulator.cpp
@@ -9,9 +9,12 @@
 #include <stdexcept>
 #include "FlightSimulator.hpp"
 #include "AerodynamicModel.hpp"
+#include "Integrator.hpp"
+#include "DefaultIntegrator.hpp"
 #include "terrain_interface.hpp"
 #include "math_utils.hpp"
 #include "physics_constants.hpp"
+#include <memory>
 
 namespace
 {
@@ -70,6 +73,23 @@ namespace
 			return 1e6F;
 		}
 	};
+
+	// Default-equivalent integrator that counts how many steps it drove, used to
+	// prove a custom integrator is actually wired into the flight phases.
+	class CountingIntegrator : public Integrator
+	{
+	public:
+		explicit CountingIntegrator(std::shared_ptr<int> calls) : calls_(std::move(calls)) {}
+
+		void step(BallState &state, float dt, const AccelerationField &accel) const override
+		{
+			++(*calls_);
+			DefaultIntegrator{}.step(state, dt, accel);
+		}
+
+	private:
+		std::shared_ptr<int> calls_;
+	};
 }
 
 class FlightSimulatorTest : public ::testing::Test
@@ -117,6 +137,69 @@ TEST_F(FlightSimulatorTest, RunsToCompletion)
 	EXPECT_STREQ(sim.getCurrentPhaseName(), "complete");
 }
 
+TEST_F(FlightSimulatorTest, ExplicitDefaultIntegratorMatchesImplicit)
+{
+	FlightSimulator implicitSim(ball, atmos, ground);
+	implicitSim.run(0.01F);
+
+	FlightSimulator explicitSim(ball, atmos, ground,
+	                            /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr,
+	                            /*ball*/ BallProperties{},
+	                            physics_constants::GRAVITY_FT_PER_S2,
+	                            std::make_shared<DefaultIntegrator>());
+	explicitSim.run(0.01F);
+
+	// Supplying the default integrator explicitly must change nothing.
+	EXPECT_FLOAT_EQ(explicitSim.getLandingResult().distance,
+	                implicitSim.getLandingResult().distance);
+	EXPECT_FLOAT_EQ(explicitSim.getLandingResult().timeOfFlight,
+	                implicitSim.getLandingResult().timeOfFlight);
+}
+
+TEST_F(FlightSimulatorTest, CustomIntegratorDrivesTheFlightSteps)
+{
+	auto calls = std::make_shared<int>(0);
+	FlightSimulator sim(ball, atmos, ground,
+	                    /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr,
+	                    /*ball*/ BallProperties{},
+	                    physics_constants::GRAVITY_FT_PER_S2,
+	                    std::make_shared<CountingIntegrator>(calls));
+	sim.run(0.01F);
+
+	// The injected integrator must have advanced the aerial/bounce steps.
+	EXPECT_GT(*calls, 0);
+	EXPECT_STREQ(sim.getCurrentPhaseName(), "complete");
+}
+
+TEST_F(FlightSimulatorTest, GravityParameterAffectsFlight)
+{
+	// Earth-gravity reference shot.
+	FlightSimulator earth(ball, atmos, ground);
+	earth.run(0.01F);
+	const float earthApex = [&] {
+		FlightSimulator s(ball, atmos, ground);
+		float apex = 0.0F;
+		for (const auto &st : s.runAndGetTrajectory(0.01F))
+			apex = std::max(apex, st.position[2]);
+		return apex;
+	}();
+	const float earthDistance = earth.getLandingResult().distance;
+
+	// Same shot under reduced gravity must fly higher and farther — proving the
+	// knob reaches the per-step integration, not just the unused initial accel.
+	const float reducedGravity = physics_constants::GRAVITY_FT_PER_S2 / 6.0F;
+	FlightSimulator low(ball, atmos, ground,
+	                    /*aero*/ nullptr, /*bounce*/ nullptr, /*roll*/ nullptr,
+	                    /*ball*/ BallProperties{}, reducedGravity);
+	float lowApex = 0.0F;
+	for (const auto &st : low.runAndGetTrajectory(0.01F))
+		lowApex = std::max(lowApex, st.position[2]);
+	low.run(0.01F);
+
+	EXPECT_GT(lowApex, earthApex);
+	EXPECT_GT(low.getLandingResult().distance, earthDistance);
+}
+
 TEST_F(FlightSimulatorTest, TransitionsThroughAllPhases)
 {
 	FlightSimulator sim(ball, atmos, ground);
diff --git a/test/test_roll_phase.cpp b/test/test_roll_phase.cpp
index 9cc2491..7f01ffc 100644
--- a/test/test_roll_phase.cpp
+++ b/test/test_roll_phase.cpp
@@ -336,12 +336,12 @@ TEST_F(RollPhaseTest, BallCanStartRollingFromNearZeroVelocityOnSlope)
 
 	RollPhase roll(slopedTerrain);
 
-	// Ball starting with very small velocity (just above MIN_VELOCITY_THRESHOLD)
+	// Ball starting with very small velocity (just above MIN_SPEED)
 	// This tests that the velocity reversal fix allows the ball to accelerate
 	// from near-zero velocity without getting stuck
 	BallState state;
 	state.position = {0.0F, 0.0F, 0.0F};
-	state.velocity = {0.0F, 0.02F, 0.0F};  // Small velocity above MIN_VELOCITY_THRESHOLD
+	state.velocity = {0.0F, 0.02F, 0.0F};  // Small velocity above MIN_SPEED
 	state.acceleration = {0.0F, 0.0F, 0.0F};
 	state.currentTime = 0.0F;