diff --git a/src/swiftest/swiftest_io.f90 b/src/swiftest/swiftest_io.f90
index e7cf51c21..70f57aebd 100644
--- a/src/swiftest/swiftest_io.f90
+++ b/src/swiftest/swiftest_io.f90
@@ -1531,7 +1531,7 @@ module function swiftest_io_netcdf_read_frame_system(self, nc, param) result(ier
             if (status == NF90_NOERR) then
                call netcdf_io_check( nf90_get_var(nc%id, nc%rotphase_varid, rtemp, start=[tslot]), &
                                   "netcdf_io_read_frame_system nf90_getvar rotphase_varid"  )
-               cb%rotphase = rtemp * DEG2RAD
+               cb%rotphase = rtemp(1) * DEG2RAD
             else
                cb%rotphase = 0.0_DP
             end if
diff --git a/swiftest/simulation_class.py b/swiftest/simulation_class.py
index 723e1c690..0511f0a50 100644
--- a/swiftest/simulation_class.py
+++ b/swiftest/simulation_class.py
@@ -2422,7 +2422,8 @@ def add_body(self,
                  Ip: List[float] | npt.NDArray[np.float_] | None=None,
                  J2: float | List[float] | npt.NDArray[np.float_] | None=None,
                  J4: float | List[float] | npt.NDArray[np.float_] | None=None,
-                 c_lm: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None
+                 c_lm: List[float] | List[npt.NDArray[np.float_]] | npt.NDArray[np.float_] | None = None,
+                 rotphase: float | List[float] | npt.NDArray[np.float_] | None=None
                  ) -> None:
         """
         Adds a body (test particle or massive body) to the internal DataSet given a set up 6 vectors (orbital elements
@@ -2465,6 +2466,8 @@ def add_body(self,
             Rotation rate vectors if these are massive bodies with rotation enabled.
         Ip : (3) or (n,3) array-like of float, optional
             Principal axes moments of inertia vectors if these are massive bodies with rotation enabled.
+        rotphase : float, optional
+            rotation phase angle of the central body in degrees
 
         Returns
         -------
@@ -2569,6 +2572,7 @@ def input_to_clm_array(val, n):
         vh,nbodies = input_to_array_3d(vh,nbodies)
         rot,nbodies = input_to_array_3d(rot,nbodies)
         Ip,nbodies = input_to_array_3d(Ip,nbodies)
+        rotphase, nbodies = input_to_array(rotphase, "f", nbodies)
 
         c_lm, nbodies = input_to_clm_array(c_lm, nbodies)