diff --git a/docs/sphinx/spelling-wordlist.txt b/docs/sphinx/spelling-wordlist.txt
index b42ba2decd..c6000fab5c 100644
--- a/docs/sphinx/spelling-wordlist.txt
+++ b/docs/sphinx/spelling-wordlist.txt
@@ -117,6 +117,9 @@ Prandtl
 pre
 priori
 proj
+quadcopter
+quaternion
+quaternions
 reStructuredText
 regrid
 regridding
@@ -140,6 +143,8 @@ teardown
 timestep
 timestepping
 timesteps
+timetable
+timetables
 tke
 vel
 vof
@@ -154,4 +159,4 @@ ylo
 yhi
 yt
 zlo
-zhi
\ No newline at end of file
+zhi
diff --git a/docs/sphinx/user/features.rst b/docs/sphinx/user/features.rst
index 8ae7c69fc1..241d6c9c77 100644
--- a/docs/sphinx/user/features.rst
+++ b/docs/sphinx/user/features.rst
@@ -76,6 +76,8 @@ Flow physics
 
    * Actuator turbine representations: Joukowsky disks, uniform disks, actuator line [:ref:`doc <turbine_models>`, :ref:`inp <inputs_actuator>`]
 
+   * Prescribed-motion rotor and multi-rotor Drone actuator representations [:ref:`inp <inputs_actuator>`]
+
    * Coupling with OpenFAST
 
    * Coupling with Nalu-Wind for blade resolved simulations
@@ -164,4 +166,3 @@ High performance computing
    * native AMReX solvers such as MLMG [:ref:`inp <inputs_mlmg>`]
 
    * hypre
-
diff --git a/docs/sphinx/user/inputs_Actuator.rst b/docs/sphinx/user/inputs_Actuator.rst
index 8abd285ba1..4f9b0b2538 100644
--- a/docs/sphinx/user/inputs_Actuator.rst
+++ b/docs/sphinx/user/inputs_Actuator.rst
@@ -26,8 +26,8 @@ turbines as actuator disks and actuator line models.
 
    This string identifies the type of actuator to use. The ones currently
    supported are: ``UniformCtDisk``, ``JoukowskyDisk``, ``TurbineFastLine``,
-   ``TurbineFastDisk``, ``TurbineKynemaLine``, ``TurbineKynemaFMBLine``,
-   ``TurbineKynemaFMBDisk``, ``FixedWingLine``, and ``ActuatorSector``.
+   ``TurbineFastDisk``, ``TurbineKynemaFMBLine``, ``TurbineKynemaFMBDisk``, ``FixedWingLine``,
+   ``ActuatorSector``, and ``Drone``.
 
    The Kynema-FMB coupled actuator turbine types are documented in the
    dedicated reference page :doc:`inputs_Actuator_KynemaFMB`.
@@ -275,6 +275,111 @@ Example for ``FixedWingLine``::
    This input allows actuator force coordinate directions to be deactivated by specifying a 0.0 in
    for the x, y, or z component of this vector.
 
+Prescribed actuator motion
+""""""""""""""""""""""""""
+
+``ActuatorSector`` and ``Drone`` support constant or time-dependent rigid-body
+motion. Timetable files are whitespace-delimited text files. The first line is
+a header and each subsequent line begins with time in seconds. Times must be
+strictly increasing. Values are linearly interpolated; velocity histories are
+integrated with the corresponding piecewise-linear, trapezoidal rule.
+
+Translation must be specified using no more than one of
+``translation_velocity``, ``velocity_timetable``, and ``position_timetable``.
+When a velocity source is used, ``center`` is the position at time zero.
+``center`` must not be combined with ``position_timetable`` because the
+position history supplies the complete position.
+
+Similarly, orientation must be specified using no more than one of a constant
+angular velocity, ``angular_velocity_timetable``, and
+``orientation_timetable``. Angular velocities are expressed in the global
+frame; this is currently the only supported angular-velocity frame.
+
+For example, a position history has three coordinates in meters::
+
+   Time X Y Z
+   0.0  0.0 0.0 1.0
+   1.0  0.0 0.0 2.0
+
+The recommended orientation format is roll, pitch, and yaw in degrees::
+
+   Time Roll Pitch Yaw
+   0.0  0.0 0.0  0.0
+   1.0  0.0 5.0 10.0
+
+These angles define the body-to-global rotation
+``Rz(yaw) Ry(pitch) Rx(roll)``. A quaternion history is also supported using
+scalar-first ``W X Y Z`` unit quaternions. Input quaternions are normalized
+before they are interpolated. Both orientation formats are interpolated with
+quaternion spherical linear interpolation.
+
+.. input_param:: Actuator.ActuatorSector.position_timetable
+
+   **type:** String, optional
+
+   File containing ``Time X Y Z``, where position is in meters in the global
+   frame.
+
+.. input_param:: Actuator.ActuatorSector.velocity_timetable
+
+   **type:** String, optional
+
+   File containing ``Time U V W``, where translational velocity is in m/s in
+   the global frame. ``center`` is required and supplies the initial position.
+
+.. input_param:: Actuator.ActuatorSector.orientation_timetable
+
+   **type:** String, optional
+
+   File containing either ``Time Roll Pitch Yaw`` in degrees or
+   ``Time W X Y Z``, as selected by ``orientation_format``.
+
+.. input_param:: Actuator.ActuatorSector.orientation_format
+
+   **type:** String, optional, default = ``roll_pitch_yaw``
+
+   Format of ``orientation_timetable``. Valid values are
+   ``roll_pitch_yaw`` and ``quaternion``.
+
+.. input_param:: Actuator.ActuatorSector.angular_velocity_timetable
+
+   **type:** String, optional
+
+   File containing ``Time OmegaX OmegaY OmegaZ``, where angular velocity is in
+   rad/s in the global frame.
+
+.. input_param:: Actuator.ActuatorSector.angular_velocity_frame
+
+   **type:** String, optional, default = ``global``
+
+   Coordinate frame for constant and tabulated body angular velocity. Only
+   ``global`` is currently supported.
+
+.. input_param:: Actuator.ActuatorSector.rotor_speed_timetable
+
+   **type:** String, optional
+
+   Rotor-speed history in rad/s. For ``ActuatorSector`` the file contains
+   ``Time Omega`` and replaces ``omega``. For ``Drone`` it contains one
+   angular-speed column per rotor and replaces ``rotor_omegas``. Rotor azimuth
+   is obtained by integrating this history.
+
+.. input_param:: Actuator.ActuatorSector.initial_azimuth_degrees
+
+   **type:** Real number, optional, default = 0.0
+
+   Initial blade azimuth in degrees. A Drone accepts either one value shared by
+   every rotor or ``num_rotors`` values.
+
+.. input_param:: Actuator.ActuatorSector.timetable_extrapolation
+
+   **type:** String, optional, default = ``hold``
+
+   Behavior outside a timetable's time range. ``hold`` holds the nearest
+   endpoint value and prints a warning the first time each bound is crossed.
+   ``error`` terminates the simulation instead.
+
+
 ActuatorSector
 """"""""""""""
 
@@ -341,13 +446,13 @@ Example for ``ActuatorSector``::
 
    **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0
 
-   Initial rotor hub center in meters in the CFD frame.
+   Initial rotor hub center in meters in the global frame.
 
 .. input_param:: Actuator.ActuatorSector.translation_velocity
 
    **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0
 
-   Prescribed drone/body translational velocity in m/s in the CFD frame.
+   Prescribed drone/body translational velocity in m/s in the global frame.
 
 .. input_param:: Actuator.ActuatorSector.rotor_normal
 
@@ -371,7 +476,7 @@ Example for ``ActuatorSector``::
 
    **type:** List of 3 real numbers, optional
 
-   Rotor-frame/body angular velocity in rad/s in the CFD frame. If specified,
+   Rotor-frame/body angular velocity in rad/s in the global frame. If specified,
    this overrides ``rotor_rotation_degrees_per_revolution``.
 
 .. input_param:: Actuator.ActuatorSector.epsilon
@@ -519,6 +624,187 @@ Example for ``ActuatorSector``::
    airfoil coefficients.
 
 
+Drone
+"""""
+
+The ``Drone`` actuator is a rigid-body composite of ``ActuatorSector`` rotors.
+All rotors share the drone position and orientation, while each rotor retains
+its own signed speed, azimuth, hub location, and blade-mirroring choice. Rotor
+aerodynamic and force-projection inputs use the ``Actuator.Drone`` namespace
+and are passed directly to every child sector.
+
+The drone center is the origin of its body frame. Arms lie in the body x-y
+plane, and each rotor axis initially points along body +z. Viewed from body +z,
+positive arm angles proceed from body +x toward body +y::
+
+                         body +y
+                            ^
+                            |
+                       R2   |   R1
+                          \ | /
+             body -x <-----+-----> body +x
+                          / | \
+                       R3   |   R4
+                            |
+
+                         body +z: out of page
+
+The sketch shows a four-rotor X layout with ``arm_phase_degrees = 45``.
+``R1`` is placed at the phase angle and subsequent rotors are placed at
+increasing, uniformly spaced angles. A positive arm angle therefore denotes a
+counterclockwise geometric placement in this view; it does not specify rotor
+spin direction.
+
+Example for a stationary X-layout quadcopter::
+
+   incflo.physics = FreeStream Actuator
+   ICNS.source_terms = ActuatorForcing
+   Actuator.labels = Q1
+   Actuator.Q1.type = Drone
+
+   Actuator.Drone.num_rotors = 4
+   Actuator.Drone.arm_length = 0.075
+   Actuator.Drone.arm_phase_degrees = 45.0
+   Actuator.Drone.rotor_omegas = 2500.0 -2500.0 2500.0 -2500.0
+   Actuator.Drone.initial_azimuth_degrees = 0.0 90.0 180.0 270.0
+   Actuator.Drone.mirror_blades = false true false true
+
+   Actuator.Q1.center = 0.0 0.0 1.0
+   Actuator.Q1.body_orientation_degrees = 0.0 0.0 0.0
+
+   Actuator.Drone.rotor_diameter = 0.10
+   Actuator.Drone.root_radius_fraction = 0.18
+   Actuator.Drone.num_blades = 2
+   Actuator.Drone.epsilon_chord = 0.5
+   Actuator.Drone.epsilon_dr = 1.5
+   Actuator.Drone.epsilon_dl = 1.5
+   Actuator.Drone.min_chord_dr = 2.0
+   Actuator.Drone.span_locs = 0.0 1.0
+   Actuator.Drone.chord = 0.01 0.006
+   Actuator.Drone.twist = 12.0 4.0
+   Actuator.Drone.airfoil_table = airfoil.txt
+   Actuator.Drone.airfoil_type = openfast
+
+For unequal arm lengths, replace the scalar with one value per rotor::
+
+   Actuator.Drone.arm_length = 0.075 0.080 0.075 0.080
+
+For an irregular layout, provide one base angle per rotor. The phase can then
+rotate the complete pattern without changing its relative geometry::
+
+   Actuator.Drone.arm_angles_degrees = 0.0 90.0 190.0 295.0
+   Actuator.Drone.arm_phase_degrees = 20.0
+
+.. input_param:: Actuator.Drone.num_rotors
+
+   **type:** int, mandatory
+
+   Number of rotors. The value must be positive.
+
+.. input_param:: Actuator.Drone.arm_length
+
+   **type:** Real number or list of real numbers, mandatory
+
+   Body-center to rotor-hub distance in meters. One value is applied to every
+   rotor. Alternatively, provide exactly ``num_rotors`` positive values in
+   rotor order.
+
+.. input_param:: Actuator.Drone.arm_phase_degrees
+
+   **type:** Real number, optional, default = 0.0
+
+   Body-frame angle of the first rotor, measured in degrees from body +x toward
+   body +y when the default uniform pattern is used. More generally, this
+   value is added to every base angle in ``arm_angles_degrees``, rotating the
+   complete arm pattern without changing the relative angles.
+
+.. input_param:: Actuator.Drone.arm_angles_degrees
+
+   **type:** List of real numbers, optional
+
+   Base body-frame arm angles in degrees. The list must contain ``num_rotors``
+   distinct angles. If omitted, the base pattern is uniformly spaced at
+   ``0, 360 / num_rotors, ...``. ``arm_phase_degrees`` is added to every base
+   angle.
+
+.. input_param:: Actuator.Drone.center
+
+   **type:** List of 3 real numbers, conditionally mandatory
+
+   Initial body-center position in meters in the global frame. It is required
+   unless ``position_timetable`` supplies the complete position history, and
+   cannot be combined with that history.
+
+.. input_param:: Actuator.Drone.body_orientation_degrees
+
+   **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0
+
+   Initial body-to-global roll, pitch, and yaw angles in degrees. The rotations
+   are composed as ``Rz(yaw) Ry(pitch) Rx(roll)``. This input cannot be combined
+   with ``orientation_timetable``.
+
+.. input_param:: Actuator.Drone.translation_velocity
+
+   **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0
+
+   Constant body translational velocity in m/s in the global frame. This input
+   is mutually exclusive with ``position_timetable`` and
+   ``velocity_timetable``.
+
+.. input_param:: Actuator.Drone.angular_velocity
+
+   **type:** List of 3 real numbers, optional, default = 0.0 0.0 0.0
+
+   Constant body angular velocity in rad/s in the global frame. This input
+   is mutually exclusive with ``orientation_timetable`` and
+   ``angular_velocity_timetable``.
+
+.. input_param:: Actuator.Drone.rotor_omegas
+
+   **type:** List of real numbers, conditionally mandatory
+
+   Signed rotor angular speeds in rad/s, with exactly one value per rotor.
+   Specify exactly one of ``rotor_omegas`` and ``rotor_speed_timetable``.
+   Positive values follow the clockwise-positive blade convention used by
+   ``ActuatorSector``.
+
+.. input_param:: Actuator.Drone.initial_azimuth_degrees
+
+   **type:** Real number or list of real numbers, optional, default = 0.0
+
+   Initial blade azimuth in degrees. One value is applied to every rotor, or
+   exactly ``num_rotors`` values may be supplied.
+
+.. input_param:: Actuator.Drone.mirror_blades
+
+   **type:** List of Boolean values, optional, default = ``false`` for every
+   rotor
+
+   One value per rotor. A mirrored rotor reverses the sign of its blade twist.
+   This permits geometrically paired counter-rotating propellers without a
+   separate rotation-direction input.
+
+The Drone accepts the prescribed-motion parameters described above using
+either the ``Actuator.Drone`` defaults or the individual
+``Actuator.<label>`` namespace. In addition, the following
+``ActuatorSector`` inputs are shared by all rotors:
+``rotor_diameter``, ``root_radius_fraction``, ``num_blades``, ``epsilon``,
+``epsilon_chord``, ``epsilon_min``, ``epsilon_dr``, ``epsilon_dl``,
+``min_chord_dr``, ``span_locs``, ``chord``, ``twist``, ``airfoil_table``,
+``airfoil_type``, quadrature options, Gaussian options, and
+``output_frequency``. See the ``ActuatorSector`` parameter descriptions for
+their units and defaults.
+
+At ``output_frequency``, the Drone writes the standard actuator NetCDF file
+``<label>.nc``. The ``force`` and ``moment`` variables contain the total
+aerodynamic load acting on the vehicle in global-frame coordinates. These loads
+are equal and opposite to the force applied to the fluid, and the moment is
+taken about the instantaneous drone center. Within the Drone group, the
+``R1``, ``R2``, and subsequent rotor subgroups contain the same blade-load
+variables as a standalone ``ActuatorSector`` output. All Drone and rotor data
+are therefore contained in the single ``<label>.nc`` file.
+
+
 TurbineFastLine
 """""""""""""""
 
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index 4a7929e076..1d8717a415 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -10,4 +10,5 @@ target_sources(${kynema_sgf_lib_name}
   ViewField.cpp
   MLMGOptions.cpp
   MeshMap.cpp
+  vs/quaternion.cpp
   )
diff --git a/src/core/vs/quaternion.H b/src/core/vs/quaternion.H
new file mode 100644
index 0000000000..713dd9f55e
--- /dev/null
+++ b/src/core/vs/quaternion.H
@@ -0,0 +1,104 @@
+#ifndef VS_QUATERNION_H
+#define VS_QUATERNION_H
+
+#include "src/core/vs/tensor.H"
+#include "src/utilities/trig_ops.H"
+
+#include <cmath>
+
+#include "AMReX.H"
+
+namespace kynema_sgf::vs {
+
+using namespace amrex::literals;
+
+/** Scalar-first unit quaternion used to interpolate and compose rotations. */
+struct Quaternion
+{
+    amrex::Real w{1.0_rt};
+    amrex::Real x{0.0_rt};
+    amrex::Real y{0.0_rt};
+    amrex::Real z{0.0_rt};
+
+    //! Normalize this quaternion in place
+    Quaternion& normalize();
+
+    //! Return a normalized copy without modifying this quaternion
+    [[nodiscard]] Quaternion normalized() const;
+
+    //! Return the quaternion conjugate
+    [[nodiscard]] Quaternion conjugate() const;
+};
+
+amrex::Real dot(const Quaternion& a, const Quaternion& b);
+
+//! Rotation composition: apply rhs first, followed by lhs
+Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs);
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor tensor_unit(const Quaternion& q)
+{
+    return {1.0_rt - 2.0_rt * (q.y * q.y + q.z * q.z),
+            2.0_rt * (q.x * q.y - q.z * q.w),
+            2.0_rt * (q.x * q.z + q.y * q.w),
+            2.0_rt * (q.x * q.y + q.z * q.w),
+            1.0_rt - 2.0_rt * (q.x * q.x + q.z * q.z),
+            2.0_rt * (q.y * q.z - q.x * q.w),
+            2.0_rt * (q.x * q.z - q.y * q.w),
+            2.0_rt * (q.y * q.z + q.x * q.w),
+            1.0_rt - 2.0_rt * (q.x * q.x + q.y * q.y)};
+}
+
+Tensor tensor(const Quaternion& input);
+
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Quaternion
+from_axis_angle(const Vector& axis, const amrex::Real angle_degrees)
+{
+    const amrex::Real half_angle = -0.5_rt * utils::radians(angle_degrees);
+    const amrex::Real axis_magnitude = mag(axis);
+    const amrex::Real scale = std::sin(half_angle) / axis_magnitude;
+    return {
+        std::cos(half_angle), scale * axis.x(), scale * axis.y(),
+        scale * axis.z()};
+}
+
+//! Compatibility interface for existing axis-angle tensor rotations
+AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor
+quaternion(const Vector& axis, const amrex::Real angle_degrees)
+{
+    // Preserve the original operation sequence used by existing actuator and
+    // sampling models. Although tensor_unit(from_axis_angle(...)) is
+    // algebraically equivalent, it changes floating-point rounding and
+    // therefore existing regression results.
+    const amrex::Real angle = -1.0_rt * utils::radians(angle_degrees);
+    const amrex::Real cval = std::cos(0.5_rt * angle);
+    const amrex::Real sval = std::sin(0.5_rt * angle);
+    const auto axis_magnitude = mag(axis);
+    const amrex::Real q0 = cval;
+    const amrex::Real q1 = sval * axis.x() / axis_magnitude;
+    const amrex::Real q2 = sval * axis.y() / axis_magnitude;
+    const amrex::Real q3 = sval * axis.z() / axis_magnitude;
+
+    Tensor result;
+    result.xx() = (q0 * q0) + (q1 * q1) - (q2 * q2) - (q3 * q3);
+    result.xy() = 2.0_rt * (q1 * q2 - q0 * q3);
+    result.xz() = 2.0_rt * (q0 * q2 + q1 * q3);
+
+    result.yx() = 2.0_rt * (q1 * q2 + q0 * q3);
+    result.yy() = (q0 * q0) - (q1 * q1) + (q2 * q2) - (q3 * q3);
+    result.yz() = 2.0_rt * (q2 * q3 - q0 * q1);
+
+    result.zx() = 2.0_rt * (q1 * q3 - q0 * q2);
+    result.zy() = 2.0_rt * (q0 * q1 + q2 * q3);
+    result.zz() = (q0 * q0) - (q1 * q1) - (q2 * q2) + (q3 * q3);
+
+    return result;
+}
+
+/** Construct Rz(yaw) Ry(pitch) Rx(roll), with angles in degrees. */
+Quaternion from_roll_pitch_yaw(const Vector& angles);
+
+Quaternion slerp(Quaternion a, Quaternion b, amrex::Real fraction);
+
+} // namespace kynema_sgf::vs
+
+#endif
diff --git a/src/core/vs/quaternion.cpp b/src/core/vs/quaternion.cpp
new file mode 100644
index 0000000000..abed755e2b
--- /dev/null
+++ b/src/core/vs/quaternion.cpp
@@ -0,0 +1,98 @@
+#include "src/core/vs/quaternion.H"
+#include "src/utilities/constants.H"
+
+#include <algorithm>
+#include <cmath>
+
+namespace kynema_sgf::vs {
+
+using namespace amrex::literals;
+
+Quaternion& Quaternion::normalize()
+{
+    const amrex::Real norm = std::sqrt(w * w + x * x + y * y + z * z);
+    if (norm <= constants::EPS) {
+        amrex::Abort("Cannot normalize a zero quaternion");
+    }
+    w /= norm;
+    x /= norm;
+    y /= norm;
+    z /= norm;
+    return *this;
+}
+
+Quaternion Quaternion::normalized() const
+{
+    Quaternion result = *this;
+    result.normalize();
+    return result;
+}
+
+Quaternion Quaternion::conjugate() const { return {w, -x, -y, -z}; }
+
+amrex::Real dot(const Quaternion& a, const Quaternion& b)
+{
+    return a.w * b.w + a.x * b.x + a.y * b.y + a.z * b.z;
+}
+
+Quaternion operator*(const Quaternion& lhs, const Quaternion& rhs)
+{
+    return {
+        lhs.w * rhs.w - lhs.x * rhs.x - lhs.y * rhs.y - lhs.z * rhs.z,
+        lhs.w * rhs.x + lhs.x * rhs.w + lhs.y * rhs.z - lhs.z * rhs.y,
+        lhs.w * rhs.y - lhs.x * rhs.z + lhs.y * rhs.w + lhs.z * rhs.x,
+        lhs.w * rhs.z + lhs.x * rhs.y - lhs.y * rhs.x + lhs.z * rhs.w};
+}
+
+Tensor tensor(const Quaternion& input)
+{
+    return tensor_unit(input.normalized());
+}
+
+Quaternion from_roll_pitch_yaw(const Vector& angles)
+{
+    // The existing vs rotation tensors use the opposite sign from the
+    // conventional active quaternion formula.
+    const amrex::Real roll = -0.5_rt * utils::radians(angles.x());
+    const amrex::Real pitch = -0.5_rt * utils::radians(angles.y());
+    const amrex::Real yaw = -0.5_rt * utils::radians(angles.z());
+    const amrex::Real cr = std::cos(roll);
+    const amrex::Real sr = std::sin(roll);
+    const amrex::Real cp = std::cos(pitch);
+    const amrex::Real sp = std::sin(pitch);
+    const amrex::Real cy = std::cos(yaw);
+    const amrex::Real sy = std::sin(yaw);
+    return Quaternion{
+        cr * cp * cy + sr * sp * sy, sr * cp * cy - cr * sp * sy,
+        cr * sp * cy + sr * cp * sy, cr * cp * sy - sr * sp * cy}
+        .normalized();
+}
+
+Quaternion slerp(Quaternion a, Quaternion b, const amrex::Real fraction)
+{
+    a.normalize();
+    b.normalize();
+    amrex::Real cosine = dot(a, b);
+    // Flip one representation when needed so interpolation follows the
+    // shortest path between orientations.
+    if (cosine < 0.0_rt) {
+        b = {-b.w, -b.x, -b.y, -b.z};
+        cosine = -cosine;
+    }
+    // Scale the near-parallel cutoff with the configured Real precision.
+    if (cosine > 1.0_rt - std::sqrt(constants::EPS)) {
+        return Quaternion{
+            a.w + fraction * (b.w - a.w), a.x + fraction * (b.x - a.x),
+            a.y + fraction * (b.y - a.y), a.z + fraction * (b.z - a.z)}
+            .normalized();
+    }
+    const amrex::Real angle = std::acos(std::clamp(cosine, -1.0_rt, 1.0_rt));
+    const amrex::Real denom = std::sin(angle);
+    const amrex::Real wa = std::sin((1.0_rt - fraction) * angle) / denom;
+    const amrex::Real wb = std::sin(fraction * angle) / denom;
+    return {
+        wa * a.w + wb * b.w, wa * a.x + wb * b.x, wa * a.y + wb * b.y,
+        wa * a.z + wb * b.z};
+}
+
+} // namespace kynema_sgf::vs
diff --git a/src/core/vs/tensor.H b/src/core/vs/tensor.H
index 540140da75..b4e1c1945d 100644
--- a/src/core/vs/tensor.H
+++ b/src/core/vs/tensor.H
@@ -160,5 +160,6 @@ using Tensor = TensorT<amrex::Real>;
 } // namespace kynema_sgf::vs
 
 #include "src/core/vs/tensorI.H"
+#include "src/core/vs/quaternion.H"
 
 #endif /* VS_TENSOR_H */
diff --git a/src/core/vs/tensorI.H b/src/core/vs/tensorI.H
index 7abc72ac9b..d9d71d5852 100644
--- a/src/core/vs/tensorI.H
+++ b/src/core/vs/tensorI.H
@@ -214,34 +214,6 @@ AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor zrot(const amrex::Real angle)
                   0.0_rt, 0.0_rt, 0.0_rt, 1.0_rt};
 }
 
-AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Tensor
-quaternion(const Vector& axis, const amrex::Real angle)
-{
-    const amrex::Real ang = -1.0_rt * utils::radians(angle);
-    const amrex::Real cval = std::cos(0.5_rt * ang);
-    const amrex::Real sval = std::sin(0.5_rt * ang);
-    const auto vmag = mag(axis);
-    const amrex::Real q0 = cval;
-    const amrex::Real q1 = sval * axis.x() / vmag;
-    const amrex::Real q2 = sval * axis.y() / vmag;
-    const amrex::Real q3 = sval * axis.z() / vmag;
-
-    Tensor t;
-    t.xx() = (q0 * q0) + (q1 * q1) - (q2 * q2) - (q3 * q3);
-    t.xy() = 2.0_rt * (q1 * q2 - q0 * q3);
-    t.xz() = 2.0_rt * (q0 * q2 + q1 * q3);
-
-    t.yx() = 2.0_rt * (q1 * q2 + q0 * q3);
-    t.yy() = (q0 * q0) - (q1 * q1) + (q2 * q2) - (q3 * q3);
-    t.yz() = 2.0_rt * (q2 * q3 - q0 * q1);
-
-    t.zx() = 2.0_rt * (q1 * q3 - q0 * q2);
-    t.zy() = 2.0_rt * (q0 * q1 + q2 * q3);
-    t.zz() = (q0 * q0) - (q1 * q1) - (q2 * q2) + (q3 * q3);
-
-    return t;
-}
-
 } // namespace kynema_sgf::vs
 
 #endif /* VS_TENSORI_H */
diff --git a/src/wind_energy/actuator/ActuatorModel.H b/src/wind_energy/actuator/ActuatorModel.H
index 33c0ce82d3..b67436876f 100644
--- a/src/wind_energy/actuator/ActuatorModel.H
+++ b/src/wind_energy/actuator/ActuatorModel.H
@@ -125,6 +125,9 @@ public:
     //! Return the meta info object for this actuator instance
     [[nodiscard]] const auto& meta() const { return m_data.meta(); }
 
+    //! Return the actuator grid for inspection and diagnostics
+    [[nodiscard]] const auto& grid() const { return m_data.grid(); }
+
     void read_inputs(const utils::ActParser& pp) override
     {
         ops::ReadInputsOp<ActTrait, SrcTrait>()(m_data, pp);
diff --git a/src/wind_energy/actuator/CMakeLists.txt b/src/wind_energy/actuator/CMakeLists.txt
index e4307f0337..9aadfa4b48 100644
--- a/src/wind_energy/actuator/CMakeLists.txt
+++ b/src/wind_energy/actuator/CMakeLists.txt
@@ -8,7 +8,9 @@ target_sources(${kynema_sgf_lib_name}
   )
 
 add_subdirectory(aero)
+add_subdirectory(motion)
 add_subdirectory(wing)
 add_subdirectory(sector)
+add_subdirectory(drone)
 add_subdirectory(turbine)
 add_subdirectory(disk)
diff --git a/src/wind_energy/actuator/actuator_types.H b/src/wind_energy/actuator/actuator_types.H
index fdc13c5d9c..0b360aebde 100644
--- a/src/wind_energy/actuator/actuator_types.H
+++ b/src/wind_energy/actuator/actuator_types.H
@@ -78,6 +78,16 @@ struct ActSrcSector : ActSrcType
     static constexpr bool is_sector = true;
 };
 
+/** Composite source representation used by a drone containing sector rotors. */
+struct ActSrcDrone : ActSrcType
+{
+    static std::string identifier() { return ""; }
+
+    static constexpr bool is_line = false;
+    static constexpr bool is_disk = false;
+    static constexpr bool is_sector = false;
+};
+
 using RealList = amrex::Vector<amrex::Real>;
 using RealSlice = ::kynema_sgf::utils::Slice<amrex::Real>;
 using VecList = amrex::Vector<kynema_sgf::vs::Vector>;
diff --git a/src/wind_energy/actuator/drone/CMakeLists.txt b/src/wind_energy/actuator/drone/CMakeLists.txt
new file mode 100644
index 0000000000..612fab6c00
--- /dev/null
+++ b/src/wind_energy/actuator/drone/CMakeLists.txt
@@ -0,0 +1,6 @@
+target_sources(${kynema_sgf_lib_name}
+  PRIVATE
+
+  Drone.cpp
+  drone_ops.cpp
+  )
diff --git a/src/wind_energy/actuator/drone/Drone.H b/src/wind_energy/actuator/drone/Drone.H
new file mode 100644
index 0000000000..602cef4074
--- /dev/null
+++ b/src/wind_energy/actuator/drone/Drone.H
@@ -0,0 +1,69 @@
+#ifndef DRONE_H
+#define DRONE_H
+
+#include "src/wind_energy/actuator/actuator_types.H"
+#include "src/wind_energy/actuator/sector/ActuatorSector.H"
+#include "src/wind_energy/actuator/sector/actuator_sector_ops.H"
+#include "src/wind_energy/actuator/motion/RigidBodyMotion.H"
+#include "src/wind_energy/actuator/motion/RotorMotion.H"
+
+#include <memory>
+#include <string>
+
+namespace kynema_sgf::actuator {
+
+//! One child actuator sector and its fixed offset in the drone body frame
+struct DroneRotor
+{
+    using SectorData = ActuatorSector::DataType;
+
+    DroneRotor(CFDSim& sim, const std::string& label, int id);
+
+    SectorData data;
+    ops::ActSrcOp<ActuatorSector, ActSrcSector> source;
+    ops::ProcessOutputsOp<ActuatorSector, ActSrcSector> output;
+    vs::Vector body_offset{vs::Vector::zero()};
+};
+
+//! Composite actuator that gives all child rotors one rigid-body motion
+struct DroneData
+{
+    int num_rotors{0};
+    RealList arm_lengths;
+    RealList arm_angles_degrees;
+    amrex::Real arm_phase_degrees{0.0_rt};
+    vs::Vector center{vs::Vector::zero()};
+    vs::Vector translation_velocity{vs::Vector::zero()};
+    vs::Vector body_orientation_degrees{vs::Vector::zero()};
+    vs::Tensor body_orientation{vs::Tensor::identity()};
+    std::shared_ptr<motion::RigidBodyMotion> body_motion;
+    std::shared_ptr<motion::RotorMotion> rotor_motion;
+    amrex::Vector<int> mirror_blades;
+    vs::Vector total_force{vs::Vector::zero()};
+    vs::Vector total_moment{vs::Vector::zero()};
+    amrex::Vector<std::unique_ptr<DroneRotor>> rotors;
+};
+
+struct Drone : public ActuatorType
+{
+    using InfoType = ActInfo;
+    using GridType = ActGrid;
+    using MetaType = DroneData;
+    using DataType = ActDataHolder<Drone>;
+
+    static std::string identifier() { return "Drone"; }
+};
+
+namespace drone {
+
+vs::Tensor body_rotation(const vs::Vector& roll_pitch_yaw_degrees);
+
+RealList uniform_arm_angles(int num_rotors, amrex::Real phase_degrees);
+
+VecList
+rotor_body_offsets(const RealList& arm_lengths, const RealList& angles_degrees);
+
+} // namespace drone
+} // namespace kynema_sgf::actuator
+
+#endif
diff --git a/src/wind_energy/actuator/drone/Drone.cpp b/src/wind_energy/actuator/drone/Drone.cpp
new file mode 100644
index 0000000000..f535bfbe5d
--- /dev/null
+++ b/src/wind_energy/actuator/drone/Drone.cpp
@@ -0,0 +1,54 @@
+#include "src/wind_energy/actuator/drone/Drone.H"
+#include "src/wind_energy/actuator/drone/drone_ops.H"
+#include "src/wind_energy/actuator/ActuatorModel.H"
+#include "src/utilities/trig_ops.H"
+
+#include <cmath>
+
+namespace kynema_sgf::actuator {
+
+DroneRotor::DroneRotor(CFDSim& sim, const std::string& label, const int id)
+    : data(sim, label, id), source(data), output(data)
+{}
+
+namespace drone {
+
+vs::Tensor body_rotation(const vs::Vector& roll_pitch_yaw_degrees)
+{
+    // Apply intrinsic body x, then y, then z rotations to body-frame vectors.
+    return vs::zrot(roll_pitch_yaw_degrees.z()) &
+           vs::yrot(roll_pitch_yaw_degrees.y()) &
+           vs::xrot(roll_pitch_yaw_degrees.x());
+}
+
+RealList uniform_arm_angles(const int num_rotors, const amrex::Real phase)
+{
+    RealList angles(num_rotors);
+    for (int i = 0; i < num_rotors; ++i) {
+        angles[i] = phase + 360.0_rt * static_cast<amrex::Real>(i) /
+                                static_cast<amrex::Real>(num_rotors);
+    }
+    return angles;
+}
+
+VecList
+rotor_body_offsets(const RealList& lengths, const RealList& angles_degrees)
+{
+    AMREX_ALWAYS_ASSERT(lengths.size() == angles_degrees.size());
+    // Arms lie in the body x-y plane; rigid-body motion maps these fixed
+    // offsets into the CFD frame at runtime.
+    VecList offsets(lengths.size());
+    for (int i = 0; i < static_cast<int>(lengths.size()); ++i) {
+        const amrex::Real angle =
+            ::kynema_sgf::utils::radians(angles_degrees[i]);
+        offsets[i] = {
+            lengths[i] * std::cos(angle), lengths[i] * std::sin(angle), 0.0_rt};
+    }
+    return offsets;
+}
+
+} // namespace drone
+
+template class ActModel<Drone, ActSrcDrone>;
+
+} // namespace kynema_sgf::actuator
diff --git a/src/wind_energy/actuator/drone/drone_ops.H b/src/wind_energy/actuator/drone/drone_ops.H
new file mode 100644
index 0000000000..089623a1fb
--- /dev/null
+++ b/src/wind_energy/actuator/drone/drone_ops.H
@@ -0,0 +1,72 @@
+#ifndef DRONE_OPS_H
+#define DRONE_OPS_H
+
+#include "src/wind_energy/actuator/drone/Drone.H"
+#include "src/wind_energy/actuator/actuator_ops.H"
+
+namespace kynema_sgf::actuator::ops {
+
+template <>
+struct ReadInputsOp<Drone, ActSrcDrone>
+{
+    void operator()(Drone::DataType& data, const utils::ActParser& pp);
+};
+
+template <>
+struct InitDataOp<Drone, ActSrcDrone>
+{
+    void operator()(Drone::DataType& data);
+};
+
+template <>
+struct UpdatePosOp<Drone, ActSrcDrone>
+{
+    void operator()(Drone::DataType& data);
+};
+
+template <>
+struct UpdateVelOp<Drone, ActSrcDrone>
+{
+    void operator()(Drone::DataType& data);
+};
+
+template <>
+struct ComputeForceOp<Drone, ActSrcDrone>
+{
+    void operator()(Drone::DataType& data);
+};
+
+template <>
+struct ProcessOutputsOp<Drone, ActSrcDrone>
+{
+    explicit ProcessOutputsOp(Drone::DataType& data) : m_data(data) {}
+    void read_io_options(const utils::ActParser& pp)
+    {
+        pp.query("output_frequency", m_output_frequency);
+    }
+    void prepare_outputs(const std::string& out_dir);
+    void write_outputs();
+
+private:
+    Drone::DataType& m_data;
+    std::string m_nc_filename;
+    int m_output_frequency{1};
+};
+
+template <>
+class ActSrcOp<Drone, ActSrcDrone>
+{
+public:
+    explicit ActSrcOp(Drone::DataType& data) : m_data(data) {}
+    void initialize();
+    void setup_op();
+    void
+    operator()(int lev, const amrex::MFIter& mfi, const amrex::Geometry& geom);
+
+private:
+    Drone::DataType& m_data;
+};
+
+} // namespace kynema_sgf::actuator::ops
+
+#endif
diff --git a/src/wind_energy/actuator/drone/drone_ops.cpp b/src/wind_energy/actuator/drone/drone_ops.cpp
new file mode 100644
index 0000000000..3e3331e565
--- /dev/null
+++ b/src/wind_energy/actuator/drone/drone_ops.cpp
@@ -0,0 +1,353 @@
+#include "src/wind_energy/actuator/drone/drone_ops.H"
+
+#include "src/CFDSim.H"
+#include "src/utilities/constants.H"
+#include "src/utilities/io_utils.H"
+#include "src/utilities/ncutils/nc_interface.H"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "AMReX_Print.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf::actuator::ops {
+namespace {
+
+void input_error(const std::string& label, const std::string& message)
+{
+    amrex::Abort("Drone '" + label + "': " + message);
+}
+
+void validate_distinct_angles(const std::string& label, const RealList& angles)
+{
+    for (int i = 0; i < static_cast<int>(angles.size()); ++i) {
+        for (int j = i + 1; j < static_cast<int>(angles.size()); ++j) {
+            amrex::Real delta = std::remainder(angles[i] - angles[j], 360.0_rt);
+            if (std::abs(delta) < constants::TIGHT_TOL) {
+                input_error(label, "arm angles must be distinct");
+            }
+        }
+    }
+}
+
+void prepare_netcdf_file(
+    const std::string& filename, const Drone::DataType& data)
+{
+#ifdef KYNEMA_SGF_USE_NETCDF
+    auto ncf = ncutils::NCFile::create(filename, NC_CLOBBER | NC_NETCDF4);
+    const std::string nt_name{"num_time_steps"};
+
+    ncf.enter_def_mode();
+    ncf.put_attr("title", "Kynema-SGF Drone actuator-load output");
+    ncf.put_attr("version", ioutils::kynema_sgf_version());
+    ncf.put_attr("created_on", ioutils::timestamp());
+    ncf.def_dim(nt_name, NC_UNLIMITED);
+    ncf.def_dim("ndim", AMREX_SPACEDIM);
+
+    auto grp = ncf.def_group(data.info().label);
+    grp.def_var("time", NC_DOUBLE, {nt_name});
+    grp.def_var("time_index", NC_INT, {nt_name});
+    grp.def_var("force", NC_DOUBLE, {nt_name, "ndim"});
+    grp.def_var("moment", NC_DOUBLE, {nt_name, "ndim"});
+    grp.var("force").put_attr(
+        "description",
+        "total density-normalized aerodynamic force acting on the vehicle");
+    grp.var("force").put_attr("units", "m^4 s^-2");
+    grp.var("moment").put_attr(
+        "description",
+        "total density-normalized aerodynamic moment about the drone center");
+    grp.var("moment").put_attr("units", "m^5 s^-2");
+    for (int i = 0; i < static_cast<int>(data.meta().rotors.size()); ++i) {
+        sector::prepare_netcdf_group(
+            grp, "R" + std::to_string(i + 1),
+            data.meta().rotors[i]->data.meta());
+    }
+    ncf.exit_def_mode();
+    for (int i = 0; i < static_cast<int>(data.meta().rotors.size()); ++i) {
+        sector::write_netcdf_group_metadata(
+            grp, "R" + std::to_string(i + 1),
+            data.meta().rotors[i]->data.meta());
+    }
+    ncf.close();
+#else
+    amrex::ignore_unused(filename, data);
+#endif
+}
+
+void write_netcdf(
+    const std::string& filename,
+    const Drone::DataType& data,
+    const amrex::Real time,
+    const int time_index)
+{
+#ifdef KYNEMA_SGF_USE_NETCDF
+    auto ncf = ncutils::NCFile::open(filename, NC_WRITE);
+    const size_t nt = ncf.dim("num_time_steps").len();
+    auto grp = ncf.group(data.info().label);
+    const auto& force = data.meta().total_force;
+    const auto& moment = data.meta().total_moment;
+
+    grp.var("time").put(&time, {nt}, {1});
+    grp.var("time_index").put(&time_index, {nt}, {1});
+    grp.var("force").put(force.data(), {nt, 0}, {1, AMREX_SPACEDIM});
+    grp.var("moment").put(moment.data(), {nt, 0}, {1, AMREX_SPACEDIM});
+    for (int i = 0; i < static_cast<int>(data.meta().rotors.size()); ++i) {
+        const auto& rotor = data.meta().rotors[i]->data;
+        sector::write_netcdf_group(
+            grp, "R" + std::to_string(i + 1), rotor.meta(), rotor.grid(), time,
+            time_index, nt);
+    }
+    ncf.close();
+#else
+    amrex::ignore_unused(filename, data, time, time_index);
+#endif
+}
+
+} // namespace
+
+void ReadInputsOp<Drone, ActSrcDrone>::operator()(
+    Drone::DataType& data, const utils::ActParser& pp)
+{
+    auto& meta = data.meta();
+    const auto& label = data.info().label;
+
+    pp.get("num_rotors", meta.num_rotors);
+    if (meta.num_rotors <= 0) {
+        input_error(label, "num_rotors must be positive");
+    }
+    if (!pp.contains("position_timetable")) {
+        pp.get("center", meta.center);
+    }
+    pp.query("translation_velocity", meta.translation_velocity);
+    pp.query("body_orientation_degrees", meta.body_orientation_degrees);
+    if (pp.contains("orientation_timetable") &&
+        pp.contains("body_orientation_degrees")) {
+        input_error(
+            label,
+            "body_orientation_degrees cannot be combined with "
+            "orientation_timetable");
+    }
+    meta.body_orientation = drone::body_rotation(meta.body_orientation_degrees);
+    // All child sectors share these histories, keeping the rigid-body pose and
+    // rotor clocks synchronized across the composite actuator.
+    meta.body_motion = std::make_shared<motion::RigidBodyMotion>();
+    meta.body_motion->read_inputs(
+        pp, meta.center, meta.body_orientation, "angular_velocity");
+    meta.rotor_motion = std::make_shared<motion::RotorMotion>();
+    meta.rotor_motion->read_drone_inputs(pp, meta.num_rotors);
+
+    // One arm length applies to every rotor; num_rotors values permit an
+    // asymmetric layout. MultiParser preserves instance-over-default priority.
+    if (!pp.contains("arm_length")) {
+        input_error(label, "arm_length is required");
+    }
+    pp.getarr("arm_length", meta.arm_lengths);
+    if (meta.arm_lengths.size() == 1) {
+        const auto arm_length = meta.arm_lengths.front();
+        meta.arm_lengths.assign(meta.num_rotors, arm_length);
+    }
+    if (static_cast<int>(meta.arm_lengths.size()) != meta.num_rotors) {
+        input_error(
+            label,
+            "arm_length must contain either one value or num_rotors values");
+    }
+    if (std::ranges::any_of(meta.arm_lengths, [](const amrex::Real length) {
+            return length <= 0.0_rt;
+        })) {
+        input_error(label, "all arm lengths must be positive");
+    }
+
+    // Explicit angles define the base arm pattern. If they are omitted, use a
+    // uniform pattern. Phase then rotates the complete pattern in the body
+    // x-y plane without changing its relative geometry.
+    if (pp.contains("arm_angles_degrees")) {
+        pp.getarr("arm_angles_degrees", meta.arm_angles_degrees);
+    } else {
+        meta.arm_angles_degrees =
+            drone::uniform_arm_angles(meta.num_rotors, 0.0_rt);
+    }
+    if (static_cast<int>(meta.arm_angles_degrees.size()) != meta.num_rotors) {
+        input_error(label, "arm_angles_degrees must contain num_rotors values");
+    }
+    pp.query("arm_phase_degrees", meta.arm_phase_degrees);
+    for (auto& angle : meta.arm_angles_degrees) {
+        angle += meta.arm_phase_degrees;
+    }
+    validate_distinct_angles(label, meta.arm_angles_degrees);
+
+    meta.mirror_blades.assign(meta.num_rotors, 0);
+    if (pp.contains("mirror_blades")) {
+        amrex::Vector<std::string> mirror_inputs;
+        pp.getarr("mirror_blades", mirror_inputs);
+        if (static_cast<int>(mirror_inputs.size()) != meta.num_rotors) {
+            input_error(label, "mirror_blades must contain num_rotors values");
+        }
+        for (int i = 0; i < meta.num_rotors; ++i) {
+            const auto value = amrex::toLower(mirror_inputs[i]);
+            if ((value != "true") && (value != "false")) {
+                input_error(
+                    label, "mirror_blades values must be true or false");
+            }
+            meta.mirror_blades[i] = (value == "true") ? 1 : 0;
+        }
+    }
+
+    const auto offsets =
+        drone::rotor_body_offsets(meta.arm_lengths, meta.arm_angles_degrees);
+    meta.rotors.reserve(meta.num_rotors);
+    for (int i = 0; i < meta.num_rotors; ++i) {
+        const std::string rotor_label = label + ".R" + std::to_string(i + 1);
+        auto rotor = std::make_unique<DroneRotor>(
+            data.sim(), rotor_label, data.info().id * 100000 + i);
+        rotor->body_offset = offsets[i];
+
+        // A drone owns its rotor configuration. Reuse the same default and
+        // instance namespaces for both drone geometry and sector aerodynamics;
+        // each reader simply ignores inputs outside its responsibility.
+        utils::ActParser sector_pp("Actuator.Drone", "Actuator." + label);
+        auto& rotor_meta = rotor->data.meta();
+        rotor_meta.body_motion = meta.body_motion;
+        rotor_meta.rotor_motion = meta.rotor_motion;
+        rotor_meta.rotor_index = i;
+        rotor_meta.body_offset = offsets[i];
+        rotor_meta.omega = meta.rotor_motion->omega(i, 0.0_rt);
+        rotor_meta.user_omega = true;
+        ReadInputsOp<ActuatorSector, ActSrcSector>()(rotor->data, sector_pp);
+        if (meta.mirror_blades[i] != 0) {
+            // Mirroring the blade geometry reverses twist without introducing
+            // a separate rotation-direction convention.
+            for (auto& twist : rotor->data.meta().twist_inp) {
+                twist = -twist;
+            }
+        }
+        const auto center =
+            meta.body_motion->position(0.0_rt) +
+            (meta.body_motion->orientation(0.0_rt) & rotor->body_offset);
+        const auto rotor_normal =
+            meta.body_motion->orientation(0.0_rt) & vs::Vector::khat();
+        sector::set_placement(
+            rotor->data, center, rotor_normal,
+            meta.body_motion->translation_velocity(0.0_rt));
+        rotor->output.read_io_options(sector_pp);
+        meta.rotors.emplace_back(std::move(rotor));
+    }
+
+    // The composite search region is the union of its child-sector regions.
+    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> lo;
+    amrex::GpuArray<amrex::Real, AMREX_SPACEDIM> hi;
+    for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+        lo[n] = std::numeric_limits<amrex::Real>::max();
+        hi[n] = std::numeric_limits<amrex::Real>::lowest();
+    }
+    for (const auto& rotor : meta.rotors) {
+        const auto& box = rotor->data.info().bound_box;
+        for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+            lo[n] = std::min(lo[n], box.lo(n));
+            hi[n] = std::max(hi[n], box.hi(n));
+        }
+    }
+    data.info().bound_box = amrex::RealBox(lo.data(), hi.data());
+}
+
+void InitDataOp<Drone, ActSrcDrone>::operator()(Drone::DataType& data)
+{
+    auto& meta = data.meta();
+    int total_points = 0;
+    for (auto& rotor : meta.rotors) {
+        InitDataOp<ActuatorSector, ActSrcSector>()(rotor->data);
+        total_points += static_cast<int>(rotor->data.grid().vel_pos.size());
+    }
+    data.grid().resize(0, total_points);
+    UpdatePosOp<Drone, ActSrcDrone>()(data);
+}
+
+void UpdatePosOp<Drone, ActSrcDrone>::operator()(Drone::DataType& data)
+{
+    auto& grid = data.grid();
+    int offset = 0;
+    for (auto& rotor : data.meta().rotors) {
+        UpdatePosOp<ActuatorSector, ActSrcSector>()(rotor->data);
+        const auto& positions = rotor->data.grid().vel_pos;
+        std::copy(
+            positions.begin(), positions.end(), grid.vel_pos.begin() + offset);
+        offset += static_cast<int>(positions.size());
+    }
+}
+
+void UpdateVelOp<Drone, ActSrcDrone>::operator()(Drone::DataType& data)
+{
+    auto& grid = data.grid();
+    int offset = 0;
+    for (auto& rotor : data.meta().rotors) {
+        auto& rotor_grid = rotor->data.grid();
+        const int npts = static_cast<int>(rotor_grid.vel.size());
+        std::copy_n(grid.vel.begin() + offset, npts, rotor_grid.vel.begin());
+        std::copy_n(
+            grid.density.begin() + offset, npts, rotor_grid.density.begin());
+        offset += npts;
+    }
+}
+
+void ComputeForceOp<Drone, ActSrcDrone>::operator()(Drone::DataType& data)
+{
+    auto& meta = data.meta();
+    const auto& time = data.sim().time();
+    const amrex::Real midpoint_time =
+        time.current_time() + 0.5_rt * sector::timestep_width(data.sim());
+    const auto drone_center = meta.body_motion->position(midpoint_time);
+    meta.total_force = vs::Vector::zero();
+    meta.total_moment = vs::Vector::zero();
+    for (auto& rotor : meta.rotors) {
+        ComputeForceOp<ActuatorSector, ActSrcSector>()(rotor->data);
+        const auto& rotor_meta = rotor->data.meta();
+        // Sector forces act on the fluid; report equal-and-opposite vehicle
+        // load.
+        meta.total_force = meta.total_force - rotor_meta.integrated_force;
+        meta.total_moment =
+            meta.total_moment - rotor_meta.integrated_moment -
+            ((rotor_meta.center - drone_center) ^ rotor_meta.integrated_force);
+    }
+}
+
+void ProcessOutputsOp<Drone, ActSrcDrone>::prepare_outputs(
+    const std::string& out_dir)
+{
+    m_nc_filename = out_dir + "/" + m_data.info().label + ".nc";
+    prepare_netcdf_file(m_nc_filename, m_data);
+}
+
+void ProcessOutputsOp<Drone, ActSrcDrone>::write_outputs()
+{
+    const auto& time = m_data.sim().time();
+    if ((m_output_frequency > 0) &&
+        (time.time_index() % m_output_frequency == 0)) {
+        write_netcdf(m_nc_filename, m_data, time.new_time(), time.time_index());
+    }
+}
+
+void ActSrcOp<Drone, ActSrcDrone>::initialize()
+{
+    for (auto& rotor : m_data.meta().rotors) {
+        rotor->source.initialize();
+    }
+}
+
+void ActSrcOp<Drone, ActSrcDrone>::setup_op()
+{
+    for (auto& rotor : m_data.meta().rotors) {
+        rotor->source.setup_op();
+    }
+}
+
+void ActSrcOp<Drone, ActSrcDrone>::operator()(
+    const int lev, const amrex::MFIter& mfi, const amrex::Geometry& geom)
+{
+    for (auto& rotor : m_data.meta().rotors) {
+        rotor->source(lev, mfi, geom);
+    }
+}
+
+} // namespace kynema_sgf::actuator::ops
diff --git a/src/wind_energy/actuator/motion/CMakeLists.txt b/src/wind_energy/actuator/motion/CMakeLists.txt
new file mode 100644
index 0000000000..5bb1fb22d8
--- /dev/null
+++ b/src/wind_energy/actuator/motion/CMakeLists.txt
@@ -0,0 +1,7 @@
+target_sources(${kynema_sgf_lib_name}
+  PRIVATE
+
+  TimeTable.cpp
+  RigidBodyMotion.cpp
+  RotorMotion.cpp
+  )
diff --git a/src/wind_energy/actuator/motion/RigidBodyMotion.H b/src/wind_energy/actuator/motion/RigidBodyMotion.H
new file mode 100644
index 0000000000..6a861047a8
--- /dev/null
+++ b/src/wind_energy/actuator/motion/RigidBodyMotion.H
@@ -0,0 +1,54 @@
+#ifndef ACTUATOR_RIGID_BODY_MOTION_H
+#define ACTUATOR_RIGID_BODY_MOTION_H
+
+#include "src/core/vs/quaternion.H"
+#include "src/wind_energy/actuator/ActParser.H"
+#include "src/wind_energy/actuator/motion/TimeTable.H"
+
+#include <cstdint>
+
+namespace kynema_sgf::actuator::motion {
+
+/** Prescribed translation and orientation of a rigid body.
+ *
+ * Positions and velocities are in the global CFD frame. The orientation is
+ * the active rotation that maps body-frame vectors into that frame.
+ */
+class RigidBodyMotion
+{
+public:
+    void read_inputs(
+        const utils::ActParser& pp,
+        const vs::Vector& initial_position,
+        const vs::Tensor& initial_orientation,
+        const std::string& angular_velocity_key = "angular_velocity",
+        const vs::Vector& default_angular_velocity = vs::Vector::zero());
+
+    [[nodiscard]] vs::Vector position(amrex::Real time) const;
+    [[nodiscard]] vs::Vector translation_velocity(amrex::Real time) const;
+    [[nodiscard]] vs::Tensor orientation(amrex::Real time) const;
+    [[nodiscard]] vs::Vector angular_velocity(amrex::Real time) const;
+
+    [[nodiscard]] bool moves() const;
+
+private:
+    enum class OrientationFormat : std::uint8_t { RollPitchYaw, Quaternion };
+
+    [[nodiscard]] vs::Quaternion orientation_quaternion(int row) const;
+    void
+    warn_ambiguous_orientation_intervals(const std::string& filename) const;
+
+    vs::Vector m_initial_position{vs::Vector::zero()};
+    vs::Vector m_constant_velocity{vs::Vector::zero()};
+    vs::Tensor m_initial_orientation{vs::Tensor::identity()};
+    vs::Vector m_constant_angular_velocity{vs::Vector::zero()};
+    TimeTable m_position_table;
+    TimeTable m_velocity_table;
+    TimeTable m_orientation_table;
+    TimeTable m_angular_velocity_table;
+    OrientationFormat m_orientation_format{OrientationFormat::RollPitchYaw};
+};
+
+} // namespace kynema_sgf::actuator::motion
+
+#endif
diff --git a/src/wind_energy/actuator/motion/RigidBodyMotion.cpp b/src/wind_energy/actuator/motion/RigidBodyMotion.cpp
new file mode 100644
index 0000000000..7a04e7004c
--- /dev/null
+++ b/src/wind_energy/actuator/motion/RigidBodyMotion.cpp
@@ -0,0 +1,274 @@
+#include "src/wind_energy/actuator/motion/RigidBodyMotion.H"
+
+#include "src/core/vs/quaternion.H"
+#include "src/utilities/constants.H"
+#include "src/wind_energy/actuator/sector/actuator_sector_ops.H"
+
+#include <algorithm>
+#include <cmath>
+
+#include "AMReX.H"
+#include "AMReX_Print.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf::actuator::motion {
+namespace {
+
+vs::Vector vector3(const RealList& values)
+{
+    return {values[0], values[1], values[2]};
+}
+
+vs::Tensor integrate_global_angular_velocity(
+    const TimeTable& table,
+    const vs::Tensor& initial_orientation,
+    const amrex::Real time)
+{
+    const auto& times = table.times();
+    vs::Tensor result = initial_orientation;
+    amrex::Real start = times.front();
+    const amrex::Real direction = (time >= start) ? 1.0_rt : -1.0_rt;
+    // Compose global-frame rotations one table interval at a time. Midpoint
+    // sampling gives a second-order update for a varying angular velocity.
+    while (direction * (time - start) > 0.0_rt) {
+        amrex::Real end = time;
+        if (direction > 0.0_rt) {
+            const auto next =
+                std::upper_bound(times.begin(), times.end(), start);
+            if (next != times.end()) {
+                end = std::min(time, *next);
+            }
+        } else {
+            end = time;
+        }
+        const amrex::Real midpoint = 0.5_rt * (start + end);
+        const auto omega = vector3(table.value(midpoint));
+        result =
+            sector::rotation_matrix_from_vector(omega * (end - start)) & result;
+        start = end;
+    }
+    return result;
+}
+
+} // namespace
+
+void RigidBodyMotion::read_inputs(
+    const utils::ActParser& pp,
+    const vs::Vector& initial_position,
+    const vs::Tensor& initial_orientation,
+    const std::string& angular_velocity_key,
+    const vs::Vector& default_angular_velocity)
+{
+    m_initial_position = initial_position;
+    m_initial_orientation = initial_orientation;
+    m_constant_angular_velocity = default_angular_velocity;
+    std::string extrapolation{"hold"};
+    pp.query("timetable_extrapolation", extrapolation);
+
+    // Translation may be prescribed by position or velocity, but never both.
+    const bool has_position = pp.contains("position_timetable");
+    const bool has_velocity_table = pp.contains("velocity_timetable");
+    const bool has_velocity = pp.contains("translation_velocity");
+    if (static_cast<int>(has_position) + static_cast<int>(has_velocity_table) +
+            static_cast<int>(has_velocity) >
+        1) {
+        amrex::Abort(
+            "Specify only one of position_timetable, velocity_timetable, and "
+            "translation_velocity");
+    }
+    if (has_position && pp.contains("center")) {
+        amrex::Abort("center cannot be combined with position_timetable");
+    }
+    if (has_velocity_table && !pp.contains("center")) {
+        amrex::Abort("center is required with velocity_timetable");
+    }
+    if (has_position) {
+        std::string filename;
+        pp.get("position_timetable", filename);
+        m_position_table.read(filename, 3, extrapolation);
+    } else if (has_velocity_table) {
+        std::string filename;
+        pp.get("velocity_timetable", filename);
+        m_velocity_table.read(filename, 3, extrapolation);
+    } else {
+        pp.query("translation_velocity", m_constant_velocity);
+    }
+
+    // Rotation follows the same exclusive-source rule as translation.
+    const bool has_orientation = pp.contains("orientation_timetable");
+    const bool has_angular_table = pp.contains("angular_velocity_timetable");
+    const bool has_angular = pp.contains(angular_velocity_key);
+    if (static_cast<int>(has_orientation) +
+            static_cast<int>(has_angular_table) +
+            static_cast<int>(has_angular) >
+        1) {
+        amrex::Abort(
+            "Specify only one of orientation_timetable, "
+            "angular_velocity_timetable, and constant angular velocity");
+    }
+    std::string angular_frame{"global"};
+    pp.query("angular_velocity_frame", angular_frame);
+    if (amrex::toLower(angular_frame) != "global") {
+        amrex::Abort(
+            "The first motion implementation supports only global "
+            "angular_velocity_frame");
+    }
+    if (has_orientation) {
+        std::string filename;
+        pp.get("orientation_timetable", filename);
+        std::string format{"roll_pitch_yaw"};
+        pp.query("orientation_format", format);
+        format = amrex::toLower(format);
+        if (format == "roll_pitch_yaw") {
+            m_orientation_format = OrientationFormat::RollPitchYaw;
+            m_orientation_table.read(filename, 3, extrapolation);
+        } else if (format == "quaternion") {
+            m_orientation_format = OrientationFormat::Quaternion;
+            m_orientation_table.read(filename, 4, extrapolation);
+        } else {
+            amrex::Abort(
+                "orientation_format must be roll_pitch_yaw or quaternion");
+        }
+        warn_ambiguous_orientation_intervals(filename);
+    } else if (has_angular_table) {
+        std::string filename;
+        pp.get("angular_velocity_timetable", filename);
+        m_angular_velocity_table.read(filename, 3, extrapolation);
+    } else {
+        pp.query(angular_velocity_key, m_constant_angular_velocity);
+    }
+}
+
+vs::Vector RigidBodyMotion::position(const amrex::Real time) const
+{
+    if (!m_position_table.empty()) {
+        return vector3(m_position_table.value(time));
+    }
+    if (!m_velocity_table.empty()) {
+        return m_initial_position + vector3(m_velocity_table.integral(time));
+    }
+    return m_initial_position + m_constant_velocity * time;
+}
+
+vs::Vector RigidBodyMotion::translation_velocity(const amrex::Real time) const
+{
+    if (!m_position_table.empty()) {
+        return vector3(m_position_table.derivative(time));
+    }
+    if (!m_velocity_table.empty()) {
+        return vector3(m_velocity_table.value(time));
+    }
+    return m_constant_velocity;
+}
+
+vs::Quaternion RigidBodyMotion::orientation_quaternion(const int row) const
+{
+    const auto values = m_orientation_table.row(row);
+    if (m_orientation_format == OrientationFormat::Quaternion) {
+        return vs::Quaternion{values[0], values[1], values[2], values[3]}
+            .normalized();
+    }
+    return vs::from_roll_pitch_yaw({values[0], values[1], values[2]});
+}
+
+void RigidBodyMotion::warn_ambiguous_orientation_intervals(
+    const std::string& filename) const
+{
+    const auto& times = m_orientation_table.times();
+    for (int upper = 1; upper < static_cast<int>(times.size()); ++upper) {
+        const auto lower_orientation = orientation_quaternion(upper - 1);
+        const auto upper_orientation = orientation_quaternion(upper);
+        if (std::abs(vs::dot(lower_orientation, upper_orientation)) <=
+            constants::TIGHT_TOL) {
+            amrex::Print()
+                << "WARNING: Consecutive orientations in actuator motion "
+                   "timetable '"
+                << filename << "' at times " << times[upper - 1] << " and "
+                << times[upper]
+                << " are approximately 180 degrees apart. The interpolation "
+                   "direction is ambiguous; add an intermediate orientation "
+                   "to prescribe the rotation direction.\n";
+        }
+    }
+}
+
+vs::Tensor RigidBodyMotion::orientation(const amrex::Real time) const
+{
+    if (!m_orientation_table.empty()) {
+        const auto& times = m_orientation_table.times();
+        if (time <= times.front()) {
+            static_cast<void>(m_orientation_table.value(time));
+            return vs::tensor(orientation_quaternion(0));
+        }
+        if (time >= times.back()) {
+            static_cast<void>(m_orientation_table.value(time));
+            return vs::tensor(
+                orientation_quaternion(static_cast<int>(times.size()) - 1));
+        }
+        const int upper = static_cast<int>(
+            std::upper_bound(times.begin(), times.end(), time) - times.begin());
+        const int lower = upper - 1;
+        const amrex::Real fraction =
+            (time - times[lower]) / (times[upper] - times[lower]);
+        return vs::tensor(
+            vs::slerp(
+                orientation_quaternion(lower), orientation_quaternion(upper),
+                fraction));
+    }
+    if (!m_angular_velocity_table.empty()) {
+        return integrate_global_angular_velocity(
+            m_angular_velocity_table, m_initial_orientation, time);
+    }
+    return sector::rotation_matrix_from_vector(
+               m_constant_angular_velocity * time) &
+           m_initial_orientation;
+}
+
+vs::Vector RigidBodyMotion::angular_velocity(const amrex::Real time) const
+{
+    if (!m_orientation_table.empty()) {
+        const auto& times = m_orientation_table.times();
+        if (times.size() == 1 || time < times.front() || time > times.back()) {
+            static_cast<void>(m_orientation_table.value(time));
+            return vs::Vector::zero();
+        }
+        int upper = static_cast<int>(
+            std::upper_bound(times.begin(), times.end(), time) - times.begin());
+        upper = std::min(upper, static_cast<int>(times.size()) - 1);
+        const int lower = upper - 1;
+        auto a = orientation_quaternion(lower);
+        auto b = orientation_quaternion(upper);
+        if (vs::dot(a, b) < 0.0_rt) {
+            b = {-b.w, -b.x, -b.y, -b.z};
+        }
+        const auto relative = (b * a.conjugate()).normalized();
+        // Convert the relative quaternion over this interval into a constant
+        // body-frame angular velocity, then express it in the CFD frame.
+        const amrex::Real half_angle =
+            std::acos(std::clamp(relative.w, -1.0_rt, 1.0_rt));
+        const amrex::Real sine = std::sin(half_angle);
+        if (std::abs(sine) <= constants::EPS) {
+            return vs::Vector::zero();
+        }
+        const amrex::Real scale =
+            2.0_rt * half_angle / (sine * (times[upper] - times[lower]));
+        const vs::Vector body_omega{
+            relative.x * scale, relative.y * scale, relative.z * scale};
+        return orientation(time) & body_omega;
+    }
+    if (!m_angular_velocity_table.empty()) {
+        return vector3(m_angular_velocity_table.value(time));
+    }
+    return m_constant_angular_velocity;
+}
+
+bool RigidBodyMotion::moves() const
+{
+    return !m_position_table.empty() || !m_velocity_table.empty() ||
+           !m_orientation_table.empty() || !m_angular_velocity_table.empty() ||
+           vs::mag(m_constant_velocity) > 0.0_rt ||
+           vs::mag(m_constant_angular_velocity) > 0.0_rt;
+}
+
+} // namespace kynema_sgf::actuator::motion
diff --git a/src/wind_energy/actuator/motion/RotorMotion.H b/src/wind_energy/actuator/motion/RotorMotion.H
new file mode 100644
index 0000000000..e73bdbdcd2
--- /dev/null
+++ b/src/wind_energy/actuator/motion/RotorMotion.H
@@ -0,0 +1,37 @@
+#ifndef ACTUATOR_ROTOR_MOTION_H
+#define ACTUATOR_ROTOR_MOTION_H
+
+#include "src/wind_energy/actuator/ActParser.H"
+#include "src/wind_energy/actuator/motion/TimeTable.H"
+
+namespace kynema_sgf::actuator::motion {
+
+/** Signed rotor-speed histories and their corresponding azimuths. */
+class RotorMotion
+{
+public:
+    void read_drone_inputs(const utils::ActParser& pp, int num_rotors);
+    void read_sector_inputs(
+        const utils::ActParser& pp,
+        amrex::Real preset_omega,
+        bool omega_is_preset);
+
+    [[nodiscard]] int num_rotors() const
+    {
+        return static_cast<int>(m_constant_omegas.size());
+    }
+    [[nodiscard]] amrex::Real omega(int rotor, amrex::Real time) const;
+    [[nodiscard]] amrex::Real azimuth(int rotor, amrex::Real time) const;
+
+private:
+    void read_initial_azimuths(const utils::ActParser& pp, int num_rotors);
+
+    RealList m_constant_omegas;
+    //! Initial rotor angles stored internally in radians
+    RealList m_initial_azimuths;
+    TimeTable m_speed_table;
+};
+
+} // namespace kynema_sgf::actuator::motion
+
+#endif
diff --git a/src/wind_energy/actuator/motion/RotorMotion.cpp b/src/wind_energy/actuator/motion/RotorMotion.cpp
new file mode 100644
index 0000000000..f94d7dcb22
--- /dev/null
+++ b/src/wind_energy/actuator/motion/RotorMotion.cpp
@@ -0,0 +1,102 @@
+#include "src/wind_energy/actuator/motion/RotorMotion.H"
+
+#include "src/utilities/trig_ops.H"
+
+#include "AMReX.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf::actuator::motion {
+
+void RotorMotion::read_initial_azimuths(
+    const utils::ActParser& pp, const int num_rotors)
+{
+    m_initial_azimuths.assign(num_rotors, 0.0_rt);
+    if (pp.contains("initial_azimuth_degrees")) {
+        RealList input_azimuths;
+        pp.getarr("initial_azimuth_degrees", input_azimuths);
+        if (input_azimuths.size() == 1) {
+            const auto initial_azimuth = input_azimuths.front();
+            input_azimuths.assign(num_rotors, initial_azimuth);
+        }
+        if (static_cast<int>(input_azimuths.size()) != num_rotors) {
+            amrex::Abort(
+                "initial_azimuth_degrees must contain either one value or "
+                "num_rotors values");
+        }
+        for (auto& value : input_azimuths) {
+            value = ::kynema_sgf::utils::radians(value);
+        }
+        m_initial_azimuths = input_azimuths;
+    }
+}
+
+void RotorMotion::read_drone_inputs(
+    const utils::ActParser& pp, const int num_rotors)
+{
+    const bool constant = pp.contains("rotor_omegas");
+    const bool timetable = pp.contains("rotor_speed_timetable");
+    if (constant == timetable) {
+        amrex::Abort(
+            "Specify exactly one of rotor_omegas and rotor_speed_timetable");
+    }
+    m_constant_omegas.assign(num_rotors, 0.0_rt);
+    if (constant) {
+        pp.getarr("rotor_omegas", m_constant_omegas);
+        if (static_cast<int>(m_constant_omegas.size()) != num_rotors) {
+            amrex::Abort("rotor_omegas must contain num_rotors values");
+        }
+    } else {
+        std::string filename;
+        std::string extrapolation{"hold"};
+        pp.get("rotor_speed_timetable", filename);
+        pp.query("timetable_extrapolation", extrapolation);
+        m_speed_table.read(filename, num_rotors, extrapolation);
+    }
+    read_initial_azimuths(pp, num_rotors);
+}
+
+void RotorMotion::read_sector_inputs(
+    const utils::ActParser& pp,
+    const amrex::Real preset_omega,
+    const bool omega_is_preset)
+{
+    const bool constant = pp.contains("omega") || omega_is_preset;
+    const bool timetable = pp.contains("rotor_speed_timetable");
+    if (constant == timetable) {
+        amrex::Abort(
+            "Specify exactly one of omega and rotor_speed_timetable for an "
+            "ActuatorSector");
+    }
+    m_constant_omegas.assign(1, preset_omega);
+    if (!omega_is_preset && pp.contains("omega")) {
+        pp.get("omega", m_constant_omegas[0]);
+    } else if (timetable) {
+        std::string filename;
+        std::string extrapolation{"hold"};
+        pp.get("rotor_speed_timetable", filename);
+        pp.query("timetable_extrapolation", extrapolation);
+        m_speed_table.read(filename, 1, extrapolation);
+    }
+    read_initial_azimuths(pp, 1);
+}
+
+amrex::Real RotorMotion::omega(const int rotor, const amrex::Real time) const
+{
+    AMREX_ALWAYS_ASSERT(rotor >= 0 && rotor < num_rotors());
+    return m_speed_table.empty() ? m_constant_omegas[rotor]
+                                 : m_speed_table.value(time)[rotor];
+}
+
+amrex::Real RotorMotion::azimuth(const int rotor, const amrex::Real time) const
+{
+    AMREX_ALWAYS_ASSERT(rotor >= 0 && rotor < num_rotors());
+    // Integrating omega preserves phase when speed varies; omega(t) * time
+    // would only be correct for a constant rotor speed.
+    const amrex::Real angle = m_speed_table.empty()
+                                  ? m_constant_omegas[rotor] * time
+                                  : m_speed_table.integral(time)[rotor];
+    return m_initial_azimuths[rotor] + angle;
+}
+
+} // namespace kynema_sgf::actuator::motion
diff --git a/src/wind_energy/actuator/motion/TimeTable.H b/src/wind_energy/actuator/motion/TimeTable.H
new file mode 100644
index 0000000000..8da3c759c0
--- /dev/null
+++ b/src/wind_energy/actuator/motion/TimeTable.H
@@ -0,0 +1,50 @@
+#ifndef ACTUATOR_MOTION_TIME_TABLE_H
+#define ACTUATOR_MOTION_TIME_TABLE_H
+
+#include "src/wind_energy/actuator/actuator_types.H"
+
+#include <string>
+
+namespace kynema_sgf::actuator::motion {
+
+/** Piecewise-linear, multi-component time history.
+ *
+ * Values, derivatives, and integrals use the same interpolation intervals so
+ * that motion reconstructed from a velocity history remains consistent.
+ */
+class TimeTable
+{
+public:
+    void read(
+        const std::string& filename,
+        int num_values,
+        const std::string& extrapolation = "hold");
+
+    [[nodiscard]] bool empty() const { return m_time.empty(); }
+    [[nodiscard]] int num_values() const { return m_num_values; }
+    [[nodiscard]] const RealList& times() const { return m_time; }
+    [[nodiscard]] RealList row(int index) const;
+
+    [[nodiscard]] RealList value(amrex::Real time) const;
+    [[nodiscard]] RealList derivative(amrex::Real time) const;
+    [[nodiscard]] RealList integral(amrex::Real time) const;
+
+private:
+    [[nodiscard]] int interval(amrex::Real time) const;
+    void warn_hold(amrex::Real time) const;
+
+    std::string m_filename;
+    int m_num_values{0};
+    RealList m_time;
+    //! Row-major values used directly by interp::linear's ncomp interface
+    RealList m_values;
+    //! Exact trapezoidal integral of each completed interpolation interval
+    amrex::Vector<RealList> m_prefix_integral;
+    mutable bool m_warned_below{false};
+    mutable bool m_warned_above{false};
+    bool m_error_on_extrapolation{false};
+};
+
+} // namespace kynema_sgf::actuator::motion
+
+#endif
diff --git a/src/wind_energy/actuator/motion/TimeTable.cpp b/src/wind_energy/actuator/motion/TimeTable.cpp
new file mode 100644
index 0000000000..7d6cb2d2f6
--- /dev/null
+++ b/src/wind_energy/actuator/motion/TimeTable.cpp
@@ -0,0 +1,182 @@
+#include "src/wind_energy/actuator/motion/TimeTable.H"
+
+#include "src/utilities/linear_interpolation.H"
+
+#include <fstream>
+#include <limits>
+
+#include "AMReX.H"
+#include "AMReX_Print.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf::actuator::motion {
+
+void TimeTable::read(
+    const std::string& filename,
+    const int num_values,
+    const std::string& extrapolation)
+{
+    m_filename = filename;
+    m_num_values = num_values;
+    m_time.clear();
+    m_values.clear();
+    m_prefix_integral.clear();
+    m_warned_below = false;
+    m_warned_above = false;
+    const auto extrapolation_lower = amrex::toLower(extrapolation);
+    if ((extrapolation_lower != "hold") && (extrapolation_lower != "error")) {
+        amrex::Abort("Actuator timetable extrapolation must be hold or error");
+    }
+    m_error_on_extrapolation = (extrapolation_lower == "error");
+    std::ifstream stream(filename);
+    if (!stream.good()) {
+        amrex::Abort("Cannot open actuator motion timetable: " + filename);
+    }
+
+    stream.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
+    amrex::Real time;
+    while (stream >> time) {
+        RealList row(num_values);
+        for (auto& item : row) {
+            if (!(stream >> item)) {
+                amrex::Abort(
+                    "Invalid actuator motion timetable row in: " + filename);
+            }
+        }
+        if (!m_time.empty() && time <= m_time.back()) {
+            amrex::Abort(
+                "Actuator motion timetable times must be strictly "
+                "increasing: " +
+                filename);
+        }
+        m_time.push_back(time);
+        m_values.insert(m_values.end(), row.begin(), row.end());
+    }
+    if (m_time.empty()) {
+        amrex::Abort("Actuator motion timetable contains no data: " + filename);
+    }
+
+    m_prefix_integral.assign(m_time.size(), RealList(num_values, 0.0_rt));
+    for (int i = 1; i < static_cast<int>(m_time.size()); ++i) {
+        const amrex::Real dt = m_time[i] - m_time[i - 1];
+        for (int n = 0; n < num_values; ++n) {
+            m_prefix_integral[i][n] =
+                m_prefix_integral[i - 1][n] +
+                0.5_rt * dt *
+                    (m_values[(i - 1) * m_num_values + n] +
+                     m_values[i * m_num_values + n]);
+        }
+    }
+}
+
+RealList TimeTable::row(const int index) const
+{
+    AMREX_ALWAYS_ASSERT(index >= 0);
+    const auto row_index = static_cast<decltype(m_time.size())>(index);
+    AMREX_ALWAYS_ASSERT(row_index < m_time.size());
+    const auto offset = static_cast<RealList::difference_type>(row_index) *
+                        static_cast<RealList::difference_type>(m_num_values);
+    const auto begin = m_values.begin() + offset;
+    return {begin, begin + m_num_values};
+}
+
+void TimeTable::warn_hold(const amrex::Real time) const
+{
+    if (m_error_on_extrapolation &&
+        ((time < m_time.front()) || (time > m_time.back()))) {
+        amrex::Abort(
+            "Actuator motion timetable evaluated outside its time range: " +
+            m_filename);
+    }
+    if (time < m_time.front() && !m_warned_below) {
+        amrex::Print() << "WARNING: Holding first value of actuator motion "
+                          "timetable '"
+                       << m_filename << "' before time " << m_time.front()
+                       << ".\n";
+        m_warned_below = true;
+    }
+    if (time > m_time.back() && !m_warned_above) {
+        amrex::Print() << "WARNING: Holding last value of actuator motion "
+                          "timetable '"
+                       << m_filename << "' after time " << m_time.back()
+                       << ".\n";
+        m_warned_above = true;
+    }
+}
+
+int TimeTable::interval(const amrex::Real time) const
+{
+    if (m_time.size() == 1 || time <= m_time.front()) {
+        return 0;
+    }
+    if (time >= m_time.back()) {
+        return static_cast<int>(m_time.size()) - 2;
+    }
+    return ::kynema_sgf::interp::bisection_search(
+               m_time.begin(), m_time.end(), time)
+        .idx;
+}
+
+RealList TimeTable::value(const amrex::Real time) const
+{
+    warn_hold(time);
+    RealList result(m_num_values);
+    for (int n = 0; n < m_num_values; ++n) {
+        result[n] = ::kynema_sgf::interp::linear(
+            m_time, m_values, time, m_num_values, n);
+    }
+    return result;
+}
+
+RealList TimeTable::derivative(const amrex::Real time) const
+{
+    warn_hold(time);
+    RealList result(m_num_values, 0.0_rt);
+    if (m_time.size() == 1 || time < m_time.front() || time > m_time.back()) {
+        return result;
+    }
+    const int i = interval(time);
+    const amrex::Real dt = m_time[i + 1] - m_time[i];
+    for (int n = 0; n < m_num_values; ++n) {
+        result[n] = (m_values[(i + 1) * m_num_values + n] -
+                     m_values[i * m_num_values + n]) /
+                    dt;
+    }
+    return result;
+}
+
+RealList TimeTable::integral(const amrex::Real time) const
+{
+    warn_hold(time);
+    RealList result(m_num_values, 0.0_rt);
+    if (time <= m_time.front()) {
+        const amrex::Real dt = time - m_time.front();
+        for (int n = 0; n < m_num_values; ++n) {
+            result[n] = dt * m_values[n];
+        }
+        return result;
+    }
+    if (time >= m_time.back()) {
+        result = m_prefix_integral.back();
+        const amrex::Real dt = time - m_time.back();
+        for (int n = 0; n < m_num_values; ++n) {
+            result[n] += dt * m_values[(m_time.size() - 1) * m_num_values + n];
+        }
+        return result;
+    }
+    const int i = interval(time);
+    result = m_prefix_integral[i];
+    const amrex::Real dt = time - m_time[i];
+    const amrex::Real interval_dt = m_time[i + 1] - m_time[i];
+    for (int n = 0; n < m_num_values; ++n) {
+        const amrex::Real slope = (m_values[(i + 1) * m_num_values + n] -
+                                   m_values[i * m_num_values + n]) /
+                                  interval_dt;
+        result[n] +=
+            m_values[i * m_num_values + n] * dt + 0.5_rt * slope * dt * dt;
+    }
+    return result;
+}
+
+} // namespace kynema_sgf::actuator::motion
diff --git a/src/wind_energy/actuator/sector/ActuatorSector.H b/src/wind_energy/actuator/sector/ActuatorSector.H
index 0f3228c8bc..34984a470c 100644
--- a/src/wind_energy/actuator/sector/ActuatorSector.H
+++ b/src/wind_energy/actuator/sector/ActuatorSector.H
@@ -3,6 +3,8 @@
 
 #include "src/wind_energy/actuator/actuator_types.H"
 #include "src/wind_energy/actuator/aero/AirfoilTable.H"
+#include "src/wind_energy/actuator/motion/RigidBodyMotion.H"
+#include "src/wind_energy/actuator/motion/RotorMotion.H"
 
 #include <memory>
 #include <string>
@@ -24,8 +26,22 @@ namespace kynema_sgf::actuator {
  *      e_r = (cos(theta), -sin(theta), 0)
  *      e_theta = (-sin(theta), -cos(theta), 0)
  */
+// Fields are grouped by purpose instead of minimizing structure padding.
+// NOLINTNEXTLINE(clang-analyzer-optin.performance.Padding)
 struct ActuatorSectorData
 {
+    //! Prescribed rigid-body motion shared with an owning Drone when present
+    std::shared_ptr<motion::RigidBodyMotion> body_motion;
+
+    //! Prescribed rotor-speed motion shared with an owning Drone when present
+    std::shared_ptr<motion::RotorMotion> rotor_motion;
+
+    //! Rotor index within rotor_motion
+    int rotor_index{0};
+
+    //! Hub offset from the rigid-body center in body coordinates [m]
+    vs::Vector body_offset{vs::Vector::zero()};
+
     //! Number of rotor blades [-]
     int num_blades{2};
 
@@ -43,6 +59,7 @@ struct ActuatorSectorData
 
     //! Rotor spin angular speed [rad/s]
     amrex::Real omega{0.0_rt};
+    bool user_omega{false};
 
     //! Initial rotor hub center in the CFD frame [m]
     vs::Vector center0{0.0_rt, 0.0_rt, 0.0_rt};
@@ -178,6 +195,14 @@ struct ActuatorSectorData
 
     //! Integrated torque about the rotor normal [m^5/s^2]
     amrex::Real torque{0.0_rt};
+
+    //! Integrated density-normalized force applied to the fluid in the CFD
+    //! frame [m^4/s^2]
+    vs::Vector integrated_force{0.0_rt, 0.0_rt, 0.0_rt};
+
+    //! Integrated density-normalized moment applied to the fluid about the
+    //! rotor hub [m^5/s^2]
+    vs::Vector integrated_moment{0.0_rt, 0.0_rt, 0.0_rt};
 };
 
 struct ActuatorSector : public ActuatorType
diff --git a/src/wind_energy/actuator/sector/actuator_sector_ops.H b/src/wind_energy/actuator/sector/actuator_sector_ops.H
index 544c5515bd..6748101807 100644
--- a/src/wind_energy/actuator/sector/actuator_sector_ops.H
+++ b/src/wind_energy/actuator/sector/actuator_sector_ops.H
@@ -6,10 +6,15 @@
 #include "src/wind_energy/actuator/actuator_utils.H"
 #include "src/core/FieldRepo.H"
 
+#include <cstddef>
 #include <cmath>
 #include <numbers>
 #include <string>
 
+namespace ncutils {
+class NCGroup;
+}
+
 namespace kynema_sgf::actuator::sector {
 
 void build_radial_grid(ActuatorSectorData& meta);
@@ -43,6 +48,13 @@ void blade_basis(
 
 void update_midpoint_sample_points(ActuatorSector::DataType& data);
 
+//! Set the initial placement and prescribed translation of a mounted rotor.
+void set_placement(
+    ActuatorSector::DataType& data,
+    const vs::Vector& center,
+    const vs::Vector& rotor_normal,
+    const vs::Vector& translation_velocity);
+
 void build_gaussian_table(ActuatorSectorData& meta);
 
 amrex::Real timestep_width(const CFDSim& sim);
@@ -52,6 +64,16 @@ void prepare_netcdf_file(
     const ActuatorSectorData& meta,
     const ActInfo& info);
 
+void prepare_netcdf_group(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta);
+
+void write_netcdf_group_metadata(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta);
+
 void write_netcdf(
     const std::string& ncfile,
     const ActuatorSectorData& meta,
@@ -60,6 +82,15 @@ void write_netcdf(
     amrex::Real time,
     int time_index);
 
+void write_netcdf_group(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta,
+    const ActGrid& grid,
+    amrex::Real time,
+    int time_index,
+    std::size_t time_record);
+
 AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real gaussian_norm_3d()
 {
     const amrex::Real pi = std::numbers::pi_v<amrex::Real>;
diff --git a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp
index 3eca97ffc9..b7a7289844 100644
--- a/src/wind_energy/actuator/sector/actuator_sector_ops.cpp
+++ b/src/wind_energy/actuator/sector/actuator_sector_ops.cpp
@@ -3,6 +3,7 @@
 #include "src/CFDSim.H"
 #include "src/core/Field.H"
 #include "src/core/gpu_utils.H"
+#include "src/utilities/constants.H"
 #include "src/utilities/io_utils.H"
 #include "src/utilities/linear_interpolation.H"
 #include "src/utilities/ncutils/nc_interface.H"
@@ -10,7 +11,6 @@
 
 #include <algorithm>
 #include <cmath>
-#include <limits>
 #include <numbers>
 
 #include "AMReX_Print.H"
@@ -148,7 +148,7 @@ vs::Tensor rotation_matrix_from_euler_degrees(const vs::Vector& angles)
 vs::Tensor orientation_matrix_from_normal(const vs::Vector& normal)
 {
     const amrex::Real normal_mag = vs::mag(normal);
-    if (normal_mag <= std::numeric_limits<amrex::Real>::epsilon()) {
+    if (normal_mag <= constants::EPS) {
         amrex::Abort("ActuatorSector: rotor_normal must be nonzero");
     }
 
@@ -163,7 +163,7 @@ vs::Tensor orientation_matrix_from_normal(const vs::Vector& normal)
 vs::Tensor rotation_matrix_from_vector(const vs::Vector& rotation_vector)
 {
     const amrex::Real angle = vs::mag(rotation_vector);
-    if (angle <= std::numeric_limits<amrex::Real>::epsilon()) {
+    if (angle <= constants::EPS) {
         return vs::Tensor::identity();
     }
     return vs::quaternion(
@@ -222,14 +222,16 @@ void update_midpoint_sample_points(ActuatorSector::DataType& data)
     // convention used elsewhere in the code. The resulting sampled velocity is
     // then used by ComputeForceOp for this step.
     const amrex::Real tmid = 0.5_rt * (time.current_time() + time.new_time());
-    const amrex::Real dt = timestep_width(data.sim());
-    meta.center = meta.center0 + meta.translation_velocity * tmid;
-    meta.azimuth = meta.omega * tmid;
-    meta.delta_azimuth = meta.omega * dt;
-
-    const auto orientation = orientation_at_time(
-        tmid, orientation_matrix_from_normal(meta.rotor_normal),
-        meta.rotor_angular_velocity);
+    const auto orientation = meta.body_motion->orientation(tmid);
+    meta.center =
+        meta.body_motion->position(tmid) + (orientation & meta.body_offset);
+    meta.translation_velocity = meta.body_motion->translation_velocity(tmid);
+    meta.rotor_angular_velocity = meta.body_motion->angular_velocity(tmid);
+    meta.omega = meta.rotor_motion->omega(meta.rotor_index, tmid);
+    meta.azimuth = meta.rotor_motion->azimuth(meta.rotor_index, tmid);
+    meta.delta_azimuth =
+        meta.rotor_motion->azimuth(meta.rotor_index, time.new_time()) -
+        meta.rotor_motion->azimuth(meta.rotor_index, time.current_time());
 
     const int nr = static_cast<int>(meta.radius.size());
     const int nvel = meta.num_blades * nr;
@@ -255,6 +257,38 @@ void update_midpoint_sample_points(ActuatorSector::DataType& data)
     }
 }
 
+void set_placement(
+    ActuatorSector::DataType& data,
+    const vs::Vector& center,
+    const vs::Vector& rotor_normal,
+    const vs::Vector& translation_velocity)
+{
+    auto& meta = data.meta();
+    meta.center0 = center;
+    meta.center = center;
+    meta.rotor_normal = rotor_normal;
+    meta.translation_velocity = translation_velocity;
+    meta.rotor_angular_velocity = vs::Vector::zero();
+    meta.user_rotor_angular_velocity = true;
+
+    const amrex::Real max_eps = local_epsilon(meta, max_interp_chord(meta));
+    const amrex::Real search_radius =
+        meta.rotor_radius + meta.support_radius_over_epsilon * max_eps;
+    if ((meta.body_motion && meta.body_motion->moves()) ||
+        vs::mag(translation_velocity) > constants::EPS) {
+        const auto& geom = data.sim().mesh().Geom(0);
+        const auto plo = geom.ProbLoArray();
+        const auto phi = geom.ProbHiArray();
+        data.info().bound_box =
+            amrex::RealBox(plo[0], plo[1], plo[2], phi[0], phi[1], phi[2]);
+    } else {
+        data.info().bound_box = amrex::RealBox(
+            center.x() - search_radius, center.y() - search_radius,
+            center.z() - search_radius, center.x() + search_radius,
+            center.y() + search_radius, center.z() + search_radius);
+    }
+}
+
 void build_gaussian_table(ActuatorSectorData& meta)
 {
     if (meta.gaussian_table_error <= 0.0_rt) {
@@ -279,31 +313,19 @@ void build_gaussian_table(ActuatorSectorData& meta)
     }
 }
 
-void prepare_netcdf_file(
-    const std::string& ncfile,
-    const ActuatorSectorData& meta,
-    const ActInfo& info)
+void prepare_netcdf_group(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta)
 {
 #ifdef KYNEMA_SGF_USE_NETCDF
-    if (info.root_proc != amrex::ParallelDescriptor::MyProc()) {
-        return;
-    }
-
-    auto ncf = ncutils::NCFile::create(ncfile, NC_CLOBBER | NC_NETCDF4);
     const std::string nt_name = "num_time_steps";
     const std::string nb_name = "num_blades";
     const std::string nr_name = "num_radial_stations";
     const std::vector<std::string> tbr{nt_name, nb_name, nr_name};
     const std::vector<std::string> tbr_vec{nt_name, nb_name, nr_name, "ndim"};
 
-    ncf.enter_def_mode();
-    ncf.put_attr("title", "Kynema-SGF actuator-sector blade-load output");
-    ncf.put_attr("version", ioutils::kynema_sgf_version());
-    ncf.put_attr("created_on", ioutils::timestamp());
-    ncf.def_dim(nt_name, NC_UNLIMITED);
-    ncf.def_dim("ndim", AMREX_SPACEDIM);
-
-    auto grp = ncf.def_group(info.label);
+    auto grp = parent.def_group(group_name);
     grp.put_attr("force_units", "kinematic force per span: m^3 s^-2");
     grp.put_attr("section_force_units", "kinematic force: m^4 s^-2");
     grp.put_attr(
@@ -340,38 +362,69 @@ void prepare_netcdf_file(
         .put_attr("description", "relative velocity along rotor normal");
     grp.var("vel_rel_theta").put_attr("units", "m s^-1");
     grp.var("vel_rel_normal").put_attr("units", "m s^-1");
-    ncf.exit_def_mode();
+#else
+    amrex::ignore_unused(parent, group_name, meta);
+#endif
+}
 
+void write_netcdf_group_metadata(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta)
+{
+#ifdef KYNEMA_SGF_USE_NETCDF
+    auto grp = parent.group(group_name);
     const size_t nr = static_cast<size_t>(meta.radius.size());
     grp.var("radius").put(meta.radius.data(), {0}, {nr});
     grp.var("dr").put(meta.dr.data(), {0}, {nr});
     grp.var("chord").put(meta.chord.data(), {0}, {nr});
     grp.var("twist").put(meta.twist.data(), {0}, {nr});
     grp.var("epsilon").put(meta.epsilon_profile.data(), {0}, {nr});
-    ncf.close();
 #else
-    amrex::ignore_unused(ncfile, meta, info);
+    amrex::ignore_unused(parent, group_name, meta);
 #endif
 }
 
-void write_netcdf(
+void prepare_netcdf_file(
     const std::string& ncfile,
     const ActuatorSectorData& meta,
-    const ActInfo& info,
-    const ActGrid& grid,
-    const amrex::Real time,
-    const int time_index)
+    const ActInfo& info)
 {
 #ifdef KYNEMA_SGF_USE_NETCDF
     if (info.root_proc != amrex::ParallelDescriptor::MyProc()) {
         return;
     }
 
-    auto ncf = ncutils::NCFile::open(ncfile, NC_WRITE);
-    const size_t nt = ncf.dim("num_time_steps").len();
+    auto ncf = ncutils::NCFile::create(ncfile, NC_CLOBBER | NC_NETCDF4);
+    ncf.enter_def_mode();
+    ncf.put_attr("title", "Kynema-SGF actuator-sector blade-load output");
+    ncf.put_attr("version", ioutils::kynema_sgf_version());
+    ncf.put_attr("created_on", ioutils::timestamp());
+    ncf.def_dim("num_time_steps", NC_UNLIMITED);
+    ncf.def_dim("ndim", AMREX_SPACEDIM);
+    prepare_netcdf_group(ncf, info.label, meta);
+    ncf.exit_def_mode();
+    write_netcdf_group_metadata(ncf, info.label, meta);
+    ncf.close();
+#else
+    amrex::ignore_unused(ncfile, meta, info);
+#endif
+}
+
+void write_netcdf_group(
+    const ncutils::NCGroup& parent,
+    const std::string& group_name,
+    const ActuatorSectorData& meta,
+    const ActGrid& grid,
+    const amrex::Real time,
+    const int time_index,
+    const std::size_t time_record)
+{
+#ifdef KYNEMA_SGF_USE_NETCDF
+    const size_t nt = time_record;
     const size_t nb = static_cast<size_t>(meta.num_blades);
     const size_t nr = meta.radius.size();
-    auto grp = ncf.group(info.label);
+    auto grp = parent.group(group_name);
 
     grp.var("time").put(&time, {nt}, {1});
     grp.var("time_index").put(&time_index, {nt}, {1});
@@ -410,6 +463,28 @@ void write_netcdf(
     grp.var("aoa").put(meta.aoa.data(), {nt, 0, 0}, {1, nb, nr});
     grp.var("cl").put(meta.cl.data(), {nt, 0, 0}, {1, nb, nr});
     grp.var("cd").put(meta.cd.data(), {nt, 0, 0}, {1, nb, nr});
+#else
+    amrex::ignore_unused(
+        parent, group_name, meta, grid, time, time_index, time_record);
+#endif
+}
+
+void write_netcdf(
+    const std::string& ncfile,
+    const ActuatorSectorData& meta,
+    const ActInfo& info,
+    const ActGrid& grid,
+    const amrex::Real time,
+    const int time_index)
+{
+#ifdef KYNEMA_SGF_USE_NETCDF
+    if (info.root_proc != amrex::ParallelDescriptor::MyProc()) {
+        return;
+    }
+
+    auto ncf = ncutils::NCFile::open(ncfile, NC_WRITE);
+    const size_t nt = ncf.dim("num_time_steps").len();
+    write_netcdf_group(ncf, info.label, meta, grid, time, time_index, nt);
     ncf.close();
 #else
     amrex::ignore_unused(ncfile, meta, info, grid, time, time_index);
@@ -425,7 +500,15 @@ void ReadInputsOp<ActuatorSector, ActSrcSector>::operator()(
 {
     auto& meta = data.meta();
     pp.get("rotor_diameter", meta.rotor_diameter);
-    pp.get("omega", meta.omega);
+    if (meta.user_omega) {
+        if (pp.contains("omega")) {
+            amrex::Abort(
+                "ActuatorSector omega was set by its owning model and must "
+                "not also be specified in the shared sector inputs");
+        }
+    } else if (!pp.contains("rotor_speed_timetable")) {
+        pp.get("omega", meta.omega);
+    }
     pp.get("airfoil_table", meta.airfoil_file);
     pp.query("airfoil_type", meta.airfoil_type);
     pp.query("num_blades", meta.num_blades);
@@ -433,6 +516,10 @@ void ReadInputsOp<ActuatorSector, ActSrcSector>::operator()(
     pp.query("center", meta.center0);
     pp.query("translation_velocity", meta.translation_velocity);
     pp.query("rotor_normal", meta.rotor_normal);
+    if (pp.contains("orientation_timetable") && pp.contains("rotor_normal")) {
+        amrex::Abort(
+            "rotor_normal cannot be combined with orientation_timetable");
+    }
     if (pp.contains("rotor_orientation") && !pp.contains("rotor_normal")) {
         amrex::Print()
             << "WARNING: ActuatorSector input 'rotor_orientation' is "
@@ -458,6 +545,15 @@ void ReadInputsOp<ActuatorSector, ActSrcSector>::operator()(
         pp.get("rotor_angular_velocity", meta.rotor_angular_velocity);
         meta.user_rotor_angular_velocity = true;
     }
+    if (pp.contains("angular_velocity")) {
+        if (pp.contains("rotor_angular_velocity")) {
+            amrex::Abort(
+                "angular_velocity and rotor_angular_velocity are mutually "
+                "exclusive");
+        }
+        pp.get("angular_velocity", meta.rotor_angular_velocity);
+        meta.user_rotor_angular_velocity = true;
+    }
     pp.queryarr("span_locs", meta.span_locs);
     pp.queryarr("chord", meta.chord_inp);
     pp.queryarr("twist", meta.twist_inp);
@@ -515,6 +611,23 @@ void ReadInputsOp<ActuatorSector, ActSrcSector>::operator()(
             meta.rotor_rotation_degrees_per_revolution;
     }
 
+    // Standalone sectors own their histories. Drone sectors arrive with shared
+    // histories already injected so every rotor uses the same body motion.
+    if (!meta.body_motion) {
+        meta.body_motion = std::make_shared<motion::RigidBodyMotion>();
+        const std::string angular_velocity_key = pp.contains("angular_velocity")
+                                                     ? "angular_velocity"
+                                                     : "rotor_angular_velocity";
+        meta.body_motion->read_inputs(
+            pp, meta.center0,
+            sector::orientation_matrix_from_normal(meta.rotor_normal),
+            angular_velocity_key, meta.rotor_angular_velocity);
+    }
+    if (!meta.rotor_motion) {
+        meta.rotor_motion = std::make_shared<motion::RotorMotion>();
+        meta.rotor_motion->read_sector_inputs(pp, meta.omega, meta.user_omega);
+    }
+
     const amrex::Real max_eps =
         sector::local_epsilon(meta, sector::max_interp_chord(meta));
     const amrex::Real search_radius =
@@ -524,11 +637,12 @@ void ReadInputsOp<ActuatorSector, ActSrcSector>::operator()(
     const auto phi = geom.ProbHiArray();
     const auto& c = meta.center0;
 
-    const bool moves = vs::mag(meta.translation_velocity) >
-                       std::numeric_limits<amrex::Real>::epsilon();
-    const bool rotates = vs::mag(meta.rotor_angular_velocity) >
-                         std::numeric_limits<amrex::Real>::epsilon();
+    const bool moves = (meta.body_motion && meta.body_motion->moves()) ||
+                       vs::mag(meta.translation_velocity) > constants::EPS;
+    const bool rotates = vs::mag(meta.rotor_angular_velocity) > constants::EPS;
     if (moves || rotates) {
+        // A prescribed trajectory is not generally bounded by its initial
+        // placement, so retain the actuator throughout the computational box.
         data.info().bound_box = amrex::RealBox(
             moves ? plo[0] : c.x() - search_radius,
             moves ? plo[1] : c.y() - search_radius,
@@ -590,18 +704,24 @@ void ComputeForceOp<ActuatorSector, ActSrcSector>::operator()(
     const amrex::Real t0 = time.current_time();
     const amrex::Real dt = sector::timestep_width(data.sim());
     const amrex::Real tmid = t0 + 0.5_rt * dt;
-    const amrex::Real start_theta = meta.omega * t0;
-    const amrex::Real dtheta = meta.omega * dt;
-    const auto initial_orientation =
-        sector::orientation_matrix_from_normal(meta.rotor_normal);
-    const auto mid_orientation = sector::orientation_at_time(
-        tmid, initial_orientation, meta.rotor_angular_velocity);
+    // Aerodynamic loads use the same midpoint pose as the CFD velocity samples.
+    const amrex::Real mid_theta =
+        meta.rotor_motion->azimuth(meta.rotor_index, tmid);
+    const auto mid_orientation = meta.body_motion->orientation(tmid);
+    const auto body_position = meta.body_motion->position(tmid);
+    const auto body_velocity = meta.body_motion->translation_velocity(tmid);
+    const auto body_angular_velocity = meta.body_motion->angular_velocity(tmid);
+    const auto hub_offset = mid_orientation & meta.body_offset;
+    meta.center = body_position + hub_offset;
+    meta.omega = meta.rotor_motion->omega(meta.rotor_index, tmid);
 
     const int nr = static_cast<int>(meta.radius.size());
     const int nvel = meta.num_blades * nr;
     VecList force_theta_normal(nvel, vs::Vector::zero());
     meta.thrust = 0.0_rt;
     meta.torque = 0.0_rt;
+    meta.integrated_force = vs::Vector::zero();
+    meta.integrated_moment = vs::Vector::zero();
 
     // First compute blade-section loads at the midpoint sample locations. The
     // force here is per unit span [m^3/s^2] in kinematic units and is converted
@@ -617,14 +737,15 @@ void ComputeForceOp<ActuatorSector, ActSrcSector>::operator()(
             vs::Vector e_theta;
             vs::Vector e_normal;
             sector::blade_basis(
-                start_theta + 0.5_rt * dtheta + phase, mid_orientation, e_r,
-                e_theta, e_normal);
+                mid_theta + phase, mid_orientation, e_r, e_theta, e_normal);
 
             const auto rel_pos = e_r * meta.radius[ir];
-            const auto frame_vel = meta.rotor_angular_velocity ^ rel_pos;
+            // Rigid-body rotation acts about the drone center, so the lever arm
+            // includes both the rotor hub offset and blade-section radius.
+            const auto frame_vel =
+                body_angular_velocity ^ (hub_offset + rel_pos);
             const auto spin_vel = e_theta * (meta.omega * meta.radius[ir]);
-            const auto blade_vel =
-                meta.translation_velocity + frame_vel + spin_vel;
+            const auto blade_vel = body_velocity + frame_vel + spin_vel;
             const auto rel_wind = grid.vel[ip] - blade_vel;
             const amrex::Real vtheta = rel_wind & e_theta;
             const amrex::Real vnormal = rel_wind & e_normal;
@@ -649,9 +770,7 @@ void ComputeForceOp<ActuatorSector, ActSrcSector>::operator()(
             const amrex::Real lift = qval * cl;
             const amrex::Real drag = qval * cd;
             const auto drag_dir =
-                (vmag > std::numeric_limits<amrex::Real>::epsilon())
-                    ? vplane.unit()
-                    : e_theta;
+                (vmag > constants::EPS) ? vplane.unit() : e_theta;
             const auto lift_dir = spin_sign * (drag_dir ^ e_r).unit();
             const auto force_on_fluid =
                 -((lift_dir * lift) + (drag_dir * drag));
@@ -676,11 +795,13 @@ void ComputeForceOp<ActuatorSector, ActSrcSector>::operator()(
     // theta_counts is one and the sector naturally reduces to actuator-line
     // behavior.
     int nforce = 0;
+    const amrex::Real body_translation_speed = vs::mag(body_velocity);
+    const amrex::Real body_angular_speed = vs::mag(body_angular_velocity);
+    const amrex::Real hub_offset_distance = vs::mag(meta.body_offset);
     for (int ir = 0; ir < nr; ++ir) {
         const amrex::Real swept_speed =
-            vs::mag(meta.translation_velocity) +
-            std::abs(meta.omega) * meta.radius[ir] +
-            vs::mag(meta.rotor_angular_velocity) * meta.radius[ir];
+            body_translation_speed + std::abs(meta.omega) * meta.radius[ir] +
+            body_angular_speed * (hub_offset_distance + meta.radius[ir]);
         const int ntheta = amrex::max(
             1, static_cast<int>(std::ceil(
                    swept_speed * dt * meta.epsilon_dl /
@@ -705,21 +826,28 @@ void ComputeForceOp<ActuatorSector, ActSrcSector>::operator()(
                 const amrex::Real xi = (static_cast<amrex::Real>(it) + 0.5_rt) /
                                        static_cast<amrex::Real>(ntheta);
                 const amrex::Real t = t0 + xi * dt;
-                const amrex::Real theta = start_theta + xi * dtheta + phase;
-                const auto orientation = sector::orientation_at_time(
-                    t, initial_orientation, meta.rotor_angular_velocity);
+                // Re-evaluate the prescribed pose and azimuth at every swept
+                // quadrature time rather than approximating the trajectory.
+                const amrex::Real theta =
+                    meta.rotor_motion->azimuth(meta.rotor_index, t) + phase;
+                const auto orientation = meta.body_motion->orientation(t);
+                const auto rotor_center = meta.body_motion->position(t) +
+                                          (orientation & meta.body_offset);
                 vs::Vector e_r;
                 vs::Vector e_theta;
                 vs::Vector e_normal;
                 sector::blade_basis(theta, orientation, e_r, e_theta, e_normal);
-                grid.pos[iq] = meta.center0 + meta.translation_velocity * t +
-                               e_r * meta.radius[ir];
+                grid.pos[iq] = rotor_center + e_r * meta.radius[ir];
                 // Split the radial section force uniformly across the swept
                 // quadrature points so the integrated force remains unchanged.
                 const amrex::Real wt =
                     meta.dr[ir] / static_cast<amrex::Real>(ntheta);
                 grid.force[iq] = wt * ((e_theta * force_theta_normal[ip].x()) +
                                        (e_normal * force_theta_normal[ip].y()));
+                meta.integrated_force = meta.integrated_force + grid.force[iq];
+                meta.integrated_moment =
+                    meta.integrated_moment +
+                    ((grid.pos[iq] - rotor_center) ^ grid.force[iq]);
                 grid.epsilon[iq] = vs::Vector::one() * meta.epsilon_profile[ir];
                 ++iq;
             }
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index aa387d35b4..7c7b005e2b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -247,6 +247,7 @@ add_test_r(act_moving_wing)
 add_test_r(act_pitching_wing_2D)
 add_test_r(act_sector)
 add_test_r(act_sector_quad)
+add_test_r(act_drone_quad)
 add_test_r(box_refinement)
 add_test_r(burggraf_flow)
 add_test_r(channel_builder_multiphase_contraction)
diff --git a/test/test_files/act_drone_quad/act_drone_quad.inp b/test/test_files/act_drone_quad/act_drone_quad.inp
new file mode 100644
index 0000000000..3a1d735932
--- /dev/null
+++ b/test/test_files/act_drone_quad/act_drone_quad.inp
@@ -0,0 +1,83 @@
+# Time controls
+time.max_step = 4
+time.fixed_dt = 1.7453292519943296e-4  # 25 deg rotor sweep per step [s]
+time.use_force_cfl = false
+time.plot_interval = -1
+time.checkpoint_interval = -1
+
+# Domain extents are Cartesian coordinates [m].
+geometry.prob_lo = -0.1 -0.2 -0.2
+geometry.prob_hi = 0.2 0.2 0.2
+geometry.is_periodic = 0 0 0
+amr.n_cell = 48 64 64
+amr.max_level = 0
+amr.max_grid_size = 16
+
+# Constant initial/background fields: density [kg/m^3], velocity [m/s].
+ConstValue.density.value = 1.0
+ConstValue.velocity.value = 1.0 0.0 0.0
+
+incflo.use_godunov = 1
+incflo.diffusion_type = 2
+incflo.do_initial_proj = 1
+incflo.initial_iterations = 3
+transport.viscosity = 1.0e-5
+transport.laminar_prandtl = 0.7
+transport.turbulent_prandtl = 0.3333
+turbulence.model = Laminar
+
+incflo.physics = FreeStream Actuator
+ICNS.source_terms = ActuatorForcing
+
+Actuator.labels = D1
+Actuator.D1.type = Drone
+
+# Shared drone layout and actuator-sector rotor definition. The arm length is
+# the body-center to rotor-hub distance. A 45-degree phase creates an X layout.
+Actuator.Drone.num_rotors = 4
+Actuator.Drone.arm_length = 0.10606601717798213
+Actuator.Drone.arm_phase_degrees = 45.0
+Actuator.Drone.rotor_omegas = 2500.0 -2500.0 2500.0 -2500.0
+Actuator.Drone.mirror_blades = false true false true
+
+Actuator.Drone.rotor_diameter = 0.10
+Actuator.Drone.root_radius_fraction = 0.18
+Actuator.Drone.num_blades = 2
+
+# Rotate the default body xy rotor plane into the CFD yz plane. This maps the
+# body +z rotor normal to CFD +x and reproduces act_sector_quad placement.
+Actuator.D1.center = 0.0 0.0 0.0
+Actuator.D1.body_orientation_degrees = -90.0 0.0 -90.0
+
+# Actuator-sector Gaussian width and generated sector grid
+Actuator.Drone.epsilon_chord = 0.5
+Actuator.Drone.epsilon_dr = 1.5
+Actuator.Drone.epsilon_dl = 1.5
+Actuator.Drone.min_chord_dr = 2.0
+
+# Actuator-sector blade section inputs
+Actuator.Drone.span_locs = 0.0 1.0
+Actuator.Drone.chord = 0.01 0.006
+Actuator.Drone.twist = 12.0 4.0
+Actuator.Drone.airfoil_table = ../act_fixed_wing/DU21_A17.txt
+Actuator.Drone.airfoil_type = openfast
+
+# Actuator-sector force projection and output
+Actuator.Drone.gaussian_type = table
+Actuator.Drone.gaussian_table_error = 1.0e-4
+Actuator.Drone.output_frequency = 1
+
+# Boundary conditions: +x inflow, +x pressure outflow, transverse slip walls.
+xlo.type = "mass_inflow"
+xlo.density = 1.0
+xlo.velocity = 1.0 0.0 0.0
+xhi.type = "pressure_outflow"
+
+ylo.type = "slip_wall"
+yhi.type = "slip_wall"
+
+zlo.type = "slip_wall"
+zhi.type = "slip_wall"
+
+incflo.verbose = 0
+nodal_proj.verbose = 0
diff --git a/unit_tests/wind_energy/actuator/CMakeLists.txt b/unit_tests/wind_energy/actuator/CMakeLists.txt
index baa158acc0..0d74dc2a73 100644
--- a/unit_tests/wind_energy/actuator/CMakeLists.txt
+++ b/unit_tests/wind_energy/actuator/CMakeLists.txt
@@ -5,6 +5,8 @@ target_sources(${kynema_sgf_unit_test_exe_name} PRIVATE
   test_airfoil.cpp
   test_actuator_free_functions.cpp
   test_actuator_sector.cpp
+  test_actuator_motion.cpp
+  test_drone.cpp
   test_disk_uniform_ct.cpp
   test_FLLC.cpp
   test_actuator_joukowsky_disk.cpp
diff --git a/unit_tests/wind_energy/actuator/test_actuator_motion.cpp b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp
new file mode 100644
index 0000000000..0aee280452
--- /dev/null
+++ b/unit_tests/wind_energy/actuator/test_actuator_motion.cpp
@@ -0,0 +1,253 @@
+#include "src/wind_energy/actuator/motion/RigidBodyMotion.H"
+#include "src/wind_energy/actuator/motion/RotorMotion.H"
+#include "src/wind_energy/actuator/motion/TimeTable.H"
+#include "src/utilities/constants.H"
+
+#include "gtest/gtest.h"
+
+#include <cmath>
+#include <cstdio>
+#include <fstream>
+#include <numbers>
+
+#include "AMReX_ParmParse.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf_tests {
+namespace {
+
+using kynema_sgf::actuator::motion::RigidBodyMotion;
+using kynema_sgf::actuator::motion::RotorMotion;
+using kynema_sgf::actuator::motion::TimeTable;
+using kynema_sgf::actuator::utils::ActParser;
+constexpr amrex::Real test_tol = kynema_sgf::constants::TIGHT_TOL;
+
+void write_table(const std::string& filename, const std::string& contents)
+{
+    std::ofstream stream(filename);
+    stream << contents;
+}
+
+TEST(ActuatorMotion, quaternion_normalization)
+{
+    kynema_sgf::vs::Quaternion quaternion{2.0_rt, 0.0_rt, 0.0_rt, 0.0_rt};
+    const auto normalized = quaternion.normalized();
+
+    EXPECT_NEAR(quaternion.w, 2.0_rt, test_tol);
+    EXPECT_NEAR(normalized.w, 1.0_rt, test_tol);
+
+    quaternion.normalize();
+    EXPECT_NEAR(quaternion.w, 1.0_rt, test_tol);
+}
+
+TEST(ActuatorMotion, slerp_flips_any_negative_dot_product)
+{
+    const kynema_sgf::vs::Quaternion start{};
+    const amrex::Real negative_w = -0.5_rt * kynema_sgf::constants::EPS;
+    const kynema_sgf::vs::Quaternion end{
+        negative_w, 0.0_rt, 0.0_rt,
+        std::sqrt(1.0_rt - negative_w * negative_w)};
+
+    ASSERT_LT(kynema_sgf::vs::dot(start, end), 0.0_rt);
+    const auto midpoint = kynema_sgf::vs::slerp(start, end, 0.5_rt);
+    EXPECT_LT(midpoint.z, 0.0_rt);
+}
+
+TEST(ActuatorMotion, timetable_interpolation_derivative_and_integral)
+{
+    const std::string filename{"motion_scalar_table.txt"};
+    write_table(filename, "Time Value\n0 0\n2 2\n");
+
+    TimeTable table;
+    table.read(filename, 1);
+    EXPECT_NEAR(table.value(1.0_rt)[0], 1.0_rt, test_tol);
+    EXPECT_NEAR(table.derivative(1.0_rt)[0], 1.0_rt, test_tol);
+    EXPECT_NEAR(table.integral(1.0_rt)[0], 0.5_rt, test_tol);
+    EXPECT_NEAR(table.value(3.0_rt)[0], 2.0_rt, test_tol);
+    EXPECT_NEAR(table.integral(3.0_rt)[0], 4.0_rt, test_tol);
+
+    std::remove(filename.c_str());
+}
+
+TEST(ActuatorMotion, timetable_error_extrapolation)
+{
+    const std::string filename{"motion_error_table.txt"};
+    write_table(filename, "Time Value\n0 0\n1 1\n");
+    TimeTable table;
+    table.read(filename, 1, "error");
+    EXPECT_THROW(static_cast<void>(table.value(2.0_rt)), amrex::RuntimeError);
+    std::remove(filename.c_str());
+}
+
+TEST(ActuatorMotion, position_and_quaternion_histories)
+{
+    constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL;
+    const std::string position_file{"motion_position_table.txt"};
+    const std::string orientation_file{"motion_orientation_table.txt"};
+    write_table(position_file, "Time X Y Z\n0 0 0 0\n1 2 0 0\n");
+    write_table(
+        orientation_file,
+        "Time Qw Qx Qy Qz\n0 1 0 0 0\n1 0.5 0 0 0.8660254037844386\n");
+
+    amrex::ParmParse pp("MotionHistory");
+    pp.add("position_timetable", position_file);
+    pp.add("orientation_timetable", orientation_file);
+    pp.add("orientation_format", "quaternion");
+    RigidBodyMotion motion;
+    motion.read_inputs(
+        ActParser("UnusedMotionDefaults", "MotionHistory"),
+        kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity());
+
+    const auto position = motion.position(0.5_rt);
+    const auto velocity = motion.translation_velocity(0.5_rt);
+    const auto rotated =
+        motion.orientation(0.5_rt) & kynema_sgf::vs::Vector::ihat();
+    EXPECT_NEAR(position.x(), 1.0_rt, tol);
+    EXPECT_NEAR(velocity.x(), 2.0_rt, tol);
+    EXPECT_NEAR(rotated.x(), 0.5_rt, tol);
+    EXPECT_NEAR(rotated.y(), std::sqrt(3.0_rt) / 2.0_rt, tol);
+
+    std::remove(position_file.c_str());
+    std::remove(orientation_file.c_str());
+}
+
+TEST(ActuatorMotion, orientation_history_warns_for_ambiguous_interval)
+{
+    const std::string filename{"motion_ambiguous_orientation_table.txt"};
+    write_table(filename, "Time Qw Qx Qy Qz\n0 1 0 0 0\n1 0 0 0 1\n");
+
+    amrex::ParmParse pp("MotionAmbiguousOrientation");
+    pp.add("orientation_timetable", filename);
+    pp.add("orientation_format", "quaternion");
+    RigidBodyMotion motion;
+    testing::internal::CaptureStdout();
+    motion.read_inputs(
+        ActParser("UnusedMotionDefaults", "MotionAmbiguousOrientation"),
+        kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity());
+    const auto warning = testing::internal::GetCapturedStdout();
+
+    EXPECT_NE(
+        warning.find("approximately 180 degrees apart"), std::string::npos);
+    std::remove(filename.c_str());
+}
+
+TEST(ActuatorMotion, roll_pitch_yaw_and_quaternion_histories_match)
+{
+    constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL;
+    const std::string rpy_file{"motion_rpy_table.txt"};
+    const std::string quaternion_file{"motion_quaternion_table.txt"};
+    write_table(rpy_file, "Time Roll Pitch Yaw\n0 0 0 0\n1 0 0 -120\n");
+    write_table(
+        quaternion_file,
+        "Time Qw Qx Qy Qz\n0 1 0 0 0\n1 0.5 0 0 0.8660254037844386\n");
+
+    amrex::ParmParse rpy_pp("MotionRpy");
+    rpy_pp.add("orientation_timetable", rpy_file);
+    RigidBodyMotion rpy_motion;
+    rpy_motion.read_inputs(
+        ActParser("UnusedRpyDefaults", "MotionRpy"),
+        kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity());
+
+    amrex::ParmParse quaternion_pp("MotionQuaternion");
+    quaternion_pp.add("orientation_timetable", quaternion_file);
+    quaternion_pp.add("orientation_format", "quaternion");
+    RigidBodyMotion quaternion_motion;
+    quaternion_motion.read_inputs(
+        ActParser("UnusedQuaternionDefaults", "MotionQuaternion"),
+        kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity());
+
+    for (const amrex::Real time : {0.0_rt, 0.25_rt, 0.5_rt, 1.0_rt}) {
+        const auto rpy = rpy_motion.orientation(time);
+        const auto quaternion = quaternion_motion.orientation(time);
+        for (int n = 0; n < kynema_sgf::vs::Tensor::ncomp; ++n) {
+            EXPECT_NEAR(rpy[n], quaternion[n], tol);
+        }
+    }
+
+    std::remove(rpy_file.c_str());
+    std::remove(quaternion_file.c_str());
+}
+
+TEST(ActuatorMotion, velocity_history_integrates_position)
+{
+    const std::string filename{"motion_velocity_table.txt"};
+    write_table(filename, "Time Ux Uy Uz\n0 0 0 0\n2 2 0 0\n");
+
+    amrex::ParmParse pp("MotionVelocity");
+    pp.addarr("center", amrex::Vector<amrex::Real>{1.0_rt, 0.0_rt, 0.0_rt});
+    pp.add("velocity_timetable", filename);
+    RigidBodyMotion motion;
+    motion.read_inputs(
+        ActParser("UnusedVelocityDefaults", "MotionVelocity"),
+        {1.0_rt, 0.0_rt, 0.0_rt}, kynema_sgf::vs::Tensor::identity());
+    EXPECT_NEAR(motion.position(1.0_rt).x(), 1.5_rt, test_tol);
+    EXPECT_NEAR(motion.translation_velocity(1.0_rt).x(), 1.0_rt, test_tol);
+
+    std::remove(filename.c_str());
+}
+
+TEST(ActuatorMotion, position_and_velocity_histories_are_exclusive)
+{
+    const std::string position_file{"motion_exclusive_position.txt"};
+    const std::string velocity_file{"motion_exclusive_velocity.txt"};
+    write_table(position_file, "Time X Y Z\n0 0 0 0\n1 1 0 0\n");
+    write_table(velocity_file, "Time Ux Uy Uz\n0 0 0 0\n1 1 0 0\n");
+
+    amrex::ParmParse pp("MotionExclusive");
+    pp.add("position_timetable", position_file);
+    pp.add("velocity_timetable", velocity_file);
+    RigidBodyMotion motion;
+    EXPECT_THROW(
+        motion.read_inputs(
+            ActParser("UnusedExclusiveDefaults", "MotionExclusive"),
+            kynema_sgf::vs::Vector::zero(), kynema_sgf::vs::Tensor::identity()),
+        amrex::RuntimeError);
+
+    std::remove(position_file.c_str());
+    std::remove(velocity_file.c_str());
+}
+
+TEST(ActuatorMotion, rotor_speed_integrates_azimuth_and_holds_rate)
+{
+    const std::string filename{"motion_rotor_table.txt"};
+    write_table(filename, "Time R1 R2\n0 0 0\n2 2 -2\n");
+
+    amrex::ParmParse pp("RotorHistory");
+    pp.add("rotor_speed_timetable", filename);
+    pp.add("initial_azimuth_degrees", 90.0_rt);
+    RotorMotion motion;
+    motion.read_drone_inputs(
+        ActParser("UnusedRotorDefaults", "RotorHistory"), 2);
+
+    EXPECT_NEAR(motion.omega(0, 1.0_rt), 1.0_rt, test_tol);
+    EXPECT_NEAR(
+        motion.azimuth(0, 1.0_rt),
+        0.5_rt + 0.5_rt * std::numbers::pi_v<amrex::Real>, test_tol);
+    EXPECT_NEAR(
+        motion.azimuth(0, 3.0_rt),
+        4.0_rt + 0.5_rt * std::numbers::pi_v<amrex::Real>, test_tol);
+    EXPECT_NEAR(
+        motion.azimuth(1, 3.0_rt),
+        -4.0_rt + 0.5_rt * std::numbers::pi_v<amrex::Real>, test_tol);
+
+    std::remove(filename.c_str());
+}
+
+TEST(ActuatorMotion, per_rotor_initial_azimuths)
+{
+    amrex::ParmParse pp("RotorInitialAzimuths");
+    pp.addarr("rotor_omegas", amrex::Vector<amrex::Real>{1.0_rt, -1.0_rt});
+    pp.addarr(
+        "initial_azimuth_degrees",
+        amrex::Vector<amrex::Real>{0.0_rt, 180.0_rt});
+    RotorMotion motion;
+    motion.read_drone_inputs(
+        ActParser("UnusedAzimuthDefaults", "RotorInitialAzimuths"), 2);
+    EXPECT_NEAR(motion.azimuth(0, 0.0_rt), 0.0_rt, test_tol);
+    EXPECT_NEAR(
+        motion.azimuth(1, 0.0_rt), std::numbers::pi_v<amrex::Real>, test_tol);
+}
+
+} // namespace
+} // namespace kynema_sgf_tests
diff --git a/unit_tests/wind_energy/actuator/test_drone.cpp b/unit_tests/wind_energy/actuator/test_drone.cpp
new file mode 100644
index 0000000000..a01fc408d3
--- /dev/null
+++ b/unit_tests/wind_energy/actuator/test_drone.cpp
@@ -0,0 +1,368 @@
+#include "src/wind_energy/actuator/drone/Drone.H"
+#include "src/wind_energy/actuator/drone/drone_ops.H"
+#include "src/wind_energy/actuator/Actuator.H"
+#include "src/wind_energy/actuator/ActuatorContainer.H"
+#include "src/wind_energy/actuator/ActuatorModel.H"
+#include "src/utilities/constants.H"
+#include "src/utilities/ncutils/nc_interface.H"
+#include "ks_test_utils/MeshTest.H"
+
+#include "gtest/gtest.h"
+
+#include <cmath>
+#include <filesystem>
+#include <fstream>
+
+#include "AMReX_REAL.H"
+
+using namespace amrex::literals;
+
+namespace kynema_sgf_tests {
+namespace {
+
+using kynema_sgf::actuator::drone::rotor_body_offsets;
+using kynema_sgf::actuator::drone::uniform_arm_angles;
+constexpr amrex::Real test_tol = kynema_sgf::constants::TIGHT_TOL;
+
+class DroneActuatorTest : public MeshTest
+{
+protected:
+    void populate_parameters() override
+    {
+        MeshTest::populate_parameters();
+        amrex::ParmParse pp_amr("amr");
+        pp_amr.add("max_level", 0);
+        pp_amr.add("max_grid_size", 16);
+        pp_amr.addarr("n_cell", amrex::Vector<int>{16, 16, 16});
+
+        amrex::ParmParse pp_time("time");
+        pp_time.add("fixed_dt", 1.0e-4_rt);
+
+        amrex::ParmParse pp_geom("geometry");
+        pp_geom.addarr(
+            "prob_lo", amrex::Vector<amrex::Real>{-0.2_rt, -0.2_rt, -0.2_rt});
+        pp_geom.addarr(
+            "prob_hi", amrex::Vector<amrex::Real>{0.2_rt, 0.2_rt, 0.2_rt});
+    }
+
+    void initialize_domain()
+    {
+        initialize_mesh();
+        sim().repo().declare_field("actuator_src_term", 3, 0);
+        auto& vel = sim().repo().declare_field("velocity", 3, 3);
+        auto& density = sim().repo().declare_field("density", 1, 3);
+        vel.setVal(0.0_rt);
+        density.setVal(1.0_rt);
+        kynema_sgf::actuator::ActuatorContainer::ParticleType::NextID(1U);
+    }
+
+    void populate_inputs()
+    {
+        amrex::ParmParse pp_a("Actuator");
+        pp_a.add("labels", std::string("D1"));
+
+        amrex::ParmParse pp_d("Actuator.Drone");
+        pp_d.add("num_rotors", 4);
+        pp_d.add("arm_length", 0.075_rt);
+        pp_d.add("arm_phase_degrees", 45.0_rt);
+        pp_d.add("rotor_diameter", 0.05_rt);
+        pp_d.add("num_blades", 2);
+        pp_d.add("root_radius_fraction", 0.18_rt);
+        pp_d.add("epsilon_chord", 0.5_rt);
+        pp_d.add("airfoil_table", m_airfoil_file);
+        pp_d.add("airfoil_type", std::string("openfast"));
+        pp_d.addarr("span_locs", amrex::Vector<amrex::Real>{0.0_rt, 1.0_rt});
+        pp_d.addarr("chord", amrex::Vector<amrex::Real>{0.01_rt, 0.006_rt});
+        pp_d.addarr("twist", amrex::Vector<amrex::Real>{12.0_rt, 4.0_rt});
+
+        amrex::ParmParse pp_i("Actuator.D1");
+        pp_i.add("type", std::string("Drone"));
+        pp_i.addarr(
+            "arm_length",
+            amrex::Vector<amrex::Real>{0.075_rt, 0.075_rt, 0.075_rt, 0.075_rt});
+        pp_i.addarr(
+            "arm_angles_degrees",
+            amrex::Vector<amrex::Real>{0.0_rt, 90.0_rt, 180.0_rt, 270.0_rt});
+        pp_i.addarr(
+            "center", amrex::Vector<amrex::Real>{0.0_rt, 0.0_rt, 0.0_rt});
+        pp_i.addarr(
+            "translation_velocity",
+            amrex::Vector<amrex::Real>{0.1_rt, -0.2_rt, 0.3_rt});
+        pp_i.addarr(
+            "rotor_omegas", amrex::Vector<amrex::Real>{
+                                2500.0_rt, -2500.0_rt, 2500.0_rt, -2500.0_rt});
+        pp_i.addarr(
+            "mirror_blades",
+            amrex::Vector<std::string>{"false", "true", "false", "true"});
+
+        amrex::ParmParse pp_s("Actuator.ActuatorSector");
+        pp_s.add("rotor_diameter", 0.05_rt);
+        pp_s.add("num_blades", 2);
+        pp_s.add("root_radius_fraction", 0.18_rt);
+        pp_s.add("epsilon_chord", 0.5_rt);
+        pp_s.add("airfoil_table", m_airfoil_file);
+        pp_s.add("airfoil_type", std::string("openfast"));
+        pp_s.addarr("span_locs", amrex::Vector<amrex::Real>{0.0_rt, 1.0_rt});
+        pp_s.addarr("chord", amrex::Vector<amrex::Real>{0.01_rt, 0.006_rt});
+        pp_s.addarr("twist", amrex::Vector<amrex::Real>{12.0_rt, 4.0_rt});
+
+        const amrex::Real offset = 0.075_rt / std::sqrt(2.0_rt);
+        const amrex::Vector<amrex::Vector<amrex::Real>> centers{
+            {offset, offset, 0.0_rt},
+            {-offset, offset, 0.0_rt},
+            {-offset, -offset, 0.0_rt},
+            {offset, -offset, 0.0_rt}};
+        for (int i = 0; i < 4; ++i) {
+            amrex::ParmParse pp_r("Actuator.R" + std::to_string(i + 1));
+            pp_r.add("omega", (i % 2 == 0) ? 2500.0_rt : -2500.0_rt);
+            if (i % 2 != 0) {
+                pp_r.addarr(
+                    "twist", amrex::Vector<amrex::Real>{-12.0_rt, -4.0_rt});
+            }
+            pp_r.addarr("center", centers[i]);
+            pp_r.addarr(
+                "translation_velocity",
+                amrex::Vector<amrex::Real>{0.1_rt, -0.2_rt, 0.3_rt});
+            pp_r.addarr(
+                "rotor_normal",
+                amrex::Vector<amrex::Real>{0.0_rt, 0.0_rt, 1.0_rt});
+        }
+    }
+
+    void write_airfoil() const
+    {
+        std::ofstream os(m_airfoil_file);
+        os << "! test polar\n5 NumAlf\n! Alpha Cl Cd Cm\n! deg - - -\n"
+              "-180 0 0.04 0\n-10 -0.5 0.02 0\n0 0 0.01 0\n"
+              "10 0.5 0.02 0\n180 0 0.04 0\n";
+    }
+
+    const std::string m_airfoil_file{"drone_airfoil.txt"};
+};
+
+class DronePhysicsTest : public kynema_sgf::actuator::Actuator
+{
+public:
+    explicit DronePhysicsTest(kynema_sgf::CFDSim& sim) : Actuator(sim) {}
+
+protected:
+    void prepare_outputs() override {}
+};
+
+TEST(DroneGeometry, plus_layout)
+{
+    const auto angles = uniform_arm_angles(4, 0.0_rt);
+    const auto offsets =
+        rotor_body_offsets({2.0_rt, 2.0_rt, 2.0_rt, 2.0_rt}, angles);
+
+    ASSERT_EQ(offsets.size(), 4U);
+    EXPECT_NEAR(offsets[0].x(), 2.0_rt, test_tol);
+    EXPECT_NEAR(offsets[0].y(), 0.0_rt, test_tol);
+    EXPECT_NEAR(offsets[1].x(), 0.0_rt, test_tol);
+    EXPECT_NEAR(offsets[1].y(), 2.0_rt, test_tol);
+    EXPECT_NEAR(offsets[2].x(), -2.0_rt, test_tol);
+    EXPECT_NEAR(offsets[3].y(), -2.0_rt, test_tol);
+}
+
+TEST(DroneGeometry, unequal_irregular_arms)
+{
+    const kynema_sgf::actuator::RealList lengths{1.0_rt, 2.0_rt, 3.0_rt};
+    const kynema_sgf::actuator::RealList angles{0.0_rt, 90.0_rt, 225.0_rt};
+    const auto offsets = rotor_body_offsets(lengths, angles);
+
+    EXPECT_NEAR(offsets[0].x(), 1.0_rt, test_tol);
+    EXPECT_NEAR(offsets[1].y(), 2.0_rt, test_tol);
+    EXPECT_NEAR(offsets[2].x(), -3.0_rt / std::sqrt(2.0_rt), test_tol);
+    EXPECT_NEAR(offsets[2].y(), -3.0_rt / std::sqrt(2.0_rt), test_tol);
+}
+
+TEST(DroneGeometry, body_orientation_defaults_to_identity)
+{
+    const auto rotation = kynema_sgf::actuator::drone::body_rotation(
+        kynema_sgf::vs::Vector::zero());
+    const auto value =
+        rotation & kynema_sgf::vs::Vector{1.0_rt, 2.0_rt, 3.0_rt};
+    EXPECT_NEAR(value.x(), 1.0_rt, test_tol);
+    EXPECT_NEAR(value.y(), 2.0_rt, test_tol);
+    EXPECT_NEAR(value.z(), 3.0_rt, test_tol);
+}
+
+TEST_F(DroneActuatorTest, composite_lifecycle)
+{
+    write_airfoil();
+    initialize_domain();
+    populate_inputs();
+
+    DronePhysicsTest actuator(sim());
+    actuator.pre_init_actions();
+    auto* drone = dynamic_cast<kynema_sgf::actuator::ActModel<
+        kynema_sgf::actuator::Drone, kynema_sgf::actuator::ActSrcDrone>*>(
+        &actuator.get_act(0));
+    ASSERT_NE(drone, nullptr);
+    ASSERT_EQ(drone->meta().rotors.size(), 4U);
+    EXPECT_NEAR(
+        drone->meta().rotors[0]->data.meta().center.x(),
+        0.075_rt / std::sqrt(2.0_rt), test_tol);
+    EXPECT_NEAR(
+        drone->meta().rotors[0]->data.meta().center.y(),
+        0.075_rt / std::sqrt(2.0_rt), test_tol);
+
+    actuator.post_init_actions();
+    actuator.pre_advance_work();
+
+    remove(m_airfoil_file.c_str());
+}
+
+TEST_F(DroneActuatorTest, matches_equivalent_standalone_sectors)
+{
+    namespace act = kynema_sgf::actuator;
+    constexpr amrex::Real tol = kynema_sgf::constants::TIGHT_TOL;
+
+    write_airfoil();
+    initialize_domain();
+    populate_inputs();
+
+    act::ActModel<act::Drone, act::ActSrcDrone> drone(sim(), "D1", 0);
+    drone.read_inputs(act::utils::ActParser("Actuator.Drone", "Actuator.D1"));
+    drone.init_actuator_source();
+    ASSERT_EQ(drone.meta().rotors.size(), 4U);
+
+    for (int i = 0; i < 4; ++i) {
+        const std::string label = "R" + std::to_string(i + 1);
+        act::ActModel<act::ActuatorSector, act::ActSrcSector> sector(
+            sim(), label, i);
+        sector.read_inputs(
+            act::utils::ActParser(
+                "Actuator.ActuatorSector", "Actuator." + label));
+        sector.init_actuator_source();
+
+        const auto& drone_data = drone.meta().rotors[i]->data;
+        const auto& drone_meta = drone_data.meta();
+        const auto& sector_meta = sector.meta();
+        const auto& drone_grid = drone_data.grid();
+        const auto& sector_grid = sector.grid();
+
+        EXPECT_EQ(drone_meta.num_blades, sector_meta.num_blades);
+        EXPECT_NEAR(drone_meta.rotor_diameter, sector_meta.rotor_diameter, tol);
+        EXPECT_NEAR(drone_meta.rotor_radius, sector_meta.rotor_radius, tol);
+        EXPECT_NEAR(drone_meta.root_radius, sector_meta.root_radius, tol);
+        EXPECT_NEAR(drone_meta.omega, sector_meta.omega, tol);
+
+        for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+            EXPECT_NEAR(drone_meta.center[n], sector_meta.center[n], tol);
+            EXPECT_NEAR(
+                drone_meta.translation_velocity[n],
+                sector_meta.translation_velocity[n], tol);
+            EXPECT_NEAR(
+                drone_meta.rotor_normal[n], sector_meta.rotor_normal[n], tol);
+        }
+
+        ASSERT_EQ(drone_meta.radius.size(), sector_meta.radius.size());
+        EXPECT_EQ(drone_meta.dr.size(), sector_meta.dr.size());
+        EXPECT_EQ(drone_meta.chord.size(), sector_meta.chord.size());
+        EXPECT_EQ(drone_meta.twist.size(), sector_meta.twist.size());
+        EXPECT_EQ(
+            drone_meta.epsilon_profile.size(),
+            sector_meta.epsilon_profile.size());
+        EXPECT_EQ(drone_meta.theta_counts, sector_meta.theta_counts);
+        for (int j = 0; j < static_cast<int>(drone_meta.radius.size()); ++j) {
+            EXPECT_NEAR(drone_meta.radius[j], sector_meta.radius[j], tol);
+            EXPECT_NEAR(drone_meta.dr[j], sector_meta.dr[j], tol);
+            EXPECT_NEAR(drone_meta.chord[j], sector_meta.chord[j], tol);
+            EXPECT_NEAR(drone_meta.twist[j], sector_meta.twist[j], tol);
+            EXPECT_NEAR(
+                drone_meta.epsilon_profile[j], sector_meta.epsilon_profile[j],
+                tol);
+        }
+
+        EXPECT_EQ(drone_grid.pos.size(), sector_grid.pos.size());
+        EXPECT_EQ(drone_grid.force.size(), sector_grid.force.size());
+        EXPECT_EQ(drone_grid.epsilon.size(), sector_grid.epsilon.size());
+        ASSERT_EQ(drone_grid.vel_pos.size(), sector_grid.vel_pos.size());
+        EXPECT_EQ(drone_grid.vel.size(), sector_grid.vel.size());
+        EXPECT_EQ(drone_grid.density.size(), sector_grid.density.size());
+        for (int j = 0; j < static_cast<int>(drone_grid.vel_pos.size()); ++j) {
+            for (int n = 0; n < AMREX_SPACEDIM; ++n) {
+                EXPECT_NEAR(
+                    drone_grid.vel_pos[j][n], sector_grid.vel_pos[j][n], tol);
+            }
+        }
+    }
+
+    remove(m_airfoil_file.c_str());
+}
+
+TEST_F(
+    DroneActuatorTest, counter_rotating_rotors_have_same_axial_force_direction)
+{
+    namespace act = kynema_sgf::actuator;
+
+    write_airfoil();
+    initialize_domain();
+    populate_inputs();
+
+    act::ActModel<act::Drone, act::ActSrcDrone> drone(sim(), "D1", 0);
+    drone.read_inputs(act::utils::ActParser("Actuator.Drone", "Actuator.D1"));
+    drone.init_actuator_source();
+    ASSERT_EQ(drone.meta().rotors.size(), 4U);
+
+    amrex::Real reference_force_normal = 0.0_rt;
+    for (int i = 0; i < 4; ++i) {
+        auto& rotor_data = drone.meta().rotors[i]->data;
+        act::ops::ComputeForceOp<act::ActuatorSector, act::ActSrcSector>()(
+            rotor_data);
+        const auto& rotor_meta = rotor_data.meta();
+        const amrex::Real force_normal =
+            rotor_meta.integrated_force & rotor_meta.rotor_normal;
+        ASSERT_NE(force_normal, 0.0_rt);
+        if (i == 0) {
+            reference_force_normal = force_normal;
+        } else {
+            EXPECT_GT(force_normal * reference_force_normal, 0.0_rt);
+        }
+    }
+
+    remove(m_airfoil_file.c_str());
+}
+
+#ifdef KYNEMA_SGF_USE_NETCDF
+TEST_F(DroneActuatorTest, writes_rotors_as_nested_netcdf_groups)
+{
+    namespace act = kynema_sgf::actuator;
+
+    write_airfoil();
+    initialize_domain();
+    populate_inputs();
+
+    act::ActModel<act::Drone, act::ActSrcDrone> drone(sim(), "D1", 0);
+    drone.read_inputs(act::utils::ActParser("Actuator.Drone", "Actuator.D1"));
+    drone.init_actuator_source();
+    drone.prepare_outputs(".");
+    drone.write_outputs();
+
+    auto ncf = ncutils::NCFile::open("D1.nc");
+    EXPECT_EQ(ncf.dim("num_time_steps").len(), 1U);
+    ASSERT_TRUE(ncf.has_group("D1"));
+    auto drone_group = ncf.group("D1");
+    EXPECT_TRUE(drone_group.has_var("force"));
+    EXPECT_TRUE(drone_group.has_var("moment"));
+    for (int i = 0; i < 4; ++i) {
+        const std::string rotor_name = "R" + std::to_string(i + 1);
+        ASSERT_TRUE(drone_group.has_group(rotor_name));
+        auto rotor_group = drone_group.group(rotor_name);
+        EXPECT_TRUE(rotor_group.has_var("thrust"));
+        EXPECT_TRUE(rotor_group.has_var("torque"));
+        EXPECT_TRUE(rotor_group.has_var("blade_force"));
+        EXPECT_TRUE(rotor_group.has_var("aoa"));
+        EXPECT_EQ(rotor_group.var("time").shape().front(), 1U);
+    }
+    ncf.close();
+
+    EXPECT_FALSE(std::filesystem::exists("D1.R1.nc"));
+    remove("D1.nc");
+    remove(m_airfoil_file.c_str());
+}
+#endif
+
+} // namespace
+} // namespace kynema_sgf_tests