diff --git a/pycircstat2/base.py b/pycircstat2/base.py
index 34b4d69..e72e335 100644
--- a/pycircstat2/base.py
+++ b/pycircstat2/base.py
@@ -1,4 +1,4 @@
-from typing import Union
+from typing import Optional, Union
 
 import numpy as np
 
@@ -15,7 +15,7 @@
     circ_std,
 )
 from .hypothesis import rayleigh_test
-from .utils import data2rad, rad2data, significance_code
+from .utils import data2rad, rad2data, rotate_data, significance_code
 from .visualization import circ_plot
 
 __names__ = ["Circular", "Axial"]
@@ -146,13 +146,14 @@ class Circular:
     def __init__(
         self,
         data: Union[np.ndarray, list],  # angle
-        w: Union[np.ndarray, list, None] = None,  # frequency
-        bins: Union[int, np.ndarray, None] = None,
+        w: Optional[Union[np.ndarray, list]] = None,  # frequency
+        bins: Optional[Union[int, np.ndarray]] = None,
         unit: str = "degree",
-        n_intervals: Union[
-            int, float, None
-        ] = None,  # number of intervals in the full cycle
+        n_intervals: Optional[Union[
+            int, float
+        ]] = None,  # number of intervals in the full cycle
         n_clusters_max: int = 1,  # number of clusters to be tested for mixture of von Mises
+        rotate: Optional[float] = None, # in rad
         **kwargs,
     ):
         # meta
@@ -184,7 +185,7 @@ def __init__(
 
         # data
         self.data = data = np.array(data) if isinstance(data, list) else data
-        self.alpha = alpha = data2rad(data, n_intervals)
+        self.alpha = alpha = data2rad(data, n_intervals) if rotate is None else rotate_data(data2rad(data, n_intervals), rotate, unit="radian")
 
         # data preprocessing
         if bins is None:
diff --git a/pycircstat2/descriptive.py b/pycircstat2/descriptive.py
index 579e80b..43d2e4a 100644
--- a/pycircstat2/descriptive.py
+++ b/pycircstat2/descriptive.py
@@ -1,4 +1,4 @@
-from typing import Tuple, Union
+from typing import Optional, Tuple, Union
 
 import numpy as np
 from scipy.stats import chi2, norm
@@ -7,10 +7,10 @@
 
 
 def circ_r(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    Cbar: Union[float, None] = None,
-    Sbar: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    Cbar: Optional[float] = None,
+    Sbar: Optional[float] = None,
 ) -> float:
     r"""
     Circular mean resultant vector length (r).
@@ -54,7 +54,7 @@ def circ_r(
 
 def circ_mean(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
+    w: Optional[np.ndarray] = None,
 ) -> Union[np.ndarray, float]:
     r"""
     Circular mean (m).
@@ -106,7 +106,7 @@ def circ_mean(
 
 def circ_mean_and_r(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
+    w: Optional[np.ndarray] = None,
 ) -> Tuple[Union[float, np.ndarray], float]:
     """
     Circular mean (m) and resultant vector length (r).
@@ -148,8 +148,8 @@ def circ_mean_and_r(
 
 def circ_mean_and_r_of_means(
     circs: Union[list, None] = None,
-    ms: Union[np.ndarray, None] = None,
-    rs: Union[np.ndarray, None] = None,
+    ms: Optional[np.ndarray] = None,
+    rs: Optional[np.ndarray] = None,
 ) -> Tuple[float, float]:
     """The Mean of a set of Mean Angles
 
@@ -194,7 +194,7 @@ def circ_mean_and_r_of_means(
 
 def circ_moment(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
+    w: Optional[np.ndarray] = None,
     p: int = 1,
     mean: Union[float, np.ndarray, None] = None,
     centered: bool = False,
@@ -241,7 +241,7 @@ def circ_moment(
 
 def circ_dispersion(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
+    w: Optional[np.ndarray] = None,
     mean=None,
 ) -> float:
     r"""
@@ -281,7 +281,7 @@ def circ_dispersion(
     return dispersion
 
 
-def circ_skewness(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float:
+def circ_skewness(alpha: np.ndarray, w: Optional[np.ndarray] = None) -> float:
     r"""
     Circular skewness, as defined by Equation 2.29 (Fisher, 1993):
 
@@ -317,7 +317,7 @@ def circ_skewness(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float
     return skewness
 
 
-def circ_kurtosis(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float:
+def circ_kurtosis(alpha: np.ndarray, w: Optional[np.ndarray] = None) -> float:
     r"""
     Circular kurtosis, as defined by Equation 2.30 (Fisher, 1993):
 
@@ -354,10 +354,10 @@ def circ_kurtosis(alpha: np.ndarray, w: Union[np.ndarray, None] = None) -> float
 
 
 def angular_var(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    r: Union[float, None] = None,
-    bin_size: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    r: Optional[float] = None,
+    bin_size: Optional[float] = None,
 ) -> float:
     r"""
     Angular variance
@@ -391,10 +391,10 @@ def angular_var(
 
 
 def angular_std(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    r: Union[float, None] = None,
-    bin_size: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    r: Optional[float] = None,
+    bin_size: Optional[float] = None,
 ) -> float:
     r"""
     Angular (standard) deviation
@@ -430,10 +430,10 @@ def angular_std(
 
 
 def circ_var(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    r: Union[float, None] = None,
-    bin_size: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    r: Optional[float] = None,
+    bin_size: Optional[float] = None,
 ) -> float:
     r"""
     Circular variance
@@ -494,10 +494,10 @@ def circ_var(
 
 
 def circ_std(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    r: Union[float, None] = None,
-    bin_size: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    r: Optional[float] = None,
+    bin_size: Optional[float] = None,
 ) -> tuple:
     r"""
     Circular standard deviation (s).
@@ -535,7 +535,7 @@ def circ_std(
 
 def circ_median(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
+    w: Optional[np.ndarray] = None,
     method: str = "deviation",
     return_average: bool = True,
     average_method: str = "all",
@@ -817,10 +817,10 @@ def circ_mean_deviation(
 
 
 def circ_mean_ci(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    mean: Union[float, None] = None,
-    r: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    mean: Optional[float] = None,
+    r: Optional[float] = None,
     n: Union[int, None] = None,
     ci: float = 0.95,
     method: str = "approximate",
@@ -921,8 +921,8 @@ def circ_mean_ci(
 
 def _circ_mean_ci_dispersion(
     alpha: np.ndarray,
-    w: Union[np.ndarray, None] = None,
-    mean: Union[float, None] = None,
+    w: Optional[np.ndarray] = None,
+    mean: Optional[float] = None,
     ci: float = 0.95,
 ) -> tuple[float, float]:
     r"""Confidence intervals based on circular dispersion.
@@ -978,10 +978,10 @@ def _circ_mean_ci_dispersion(
 
 
 def _circ_mean_ci_approximate(
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
-    mean: Union[float, None] = None,
-    r: Union[float, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
+    mean: Optional[float] = None,
+    r: Optional[float] = None,
     n: Union[int, None] = None,
     ci: float = 0.95,
 ) -> tuple:
@@ -1146,8 +1146,8 @@ def _circ_mean_resample(alpha, z0, v0):
 
 def circ_median_ci(
     median: float = None,
-    alpha: Union[np.ndarray, None] = None,
-    w: Union[np.ndarray, None] = None,
+    alpha: Optional[np.ndarray] = None,
+    w: Optional[np.ndarray] = None,
     method: str = "deviation",
     ci: float = 0.95,
 ) -> tuple:
@@ -1201,10 +1201,18 @@ def circ_median_ci(
 
         offset = int(1 + np.floor(0.5 * np.sqrt(n) * z))  # fisher:eq(4.19)
 
-        idx_median = np.where(alpha.round(5) < median.round(5))[0][-1]
+        # idx_median = np.where(alpha.round(5) < np.round(median, 5))[0][-1]
+        arr = np.where(alpha.round(5) < np.round(median, 5))[0]
+        if len(arr) == 0:
+            # That means median is smaller than alpha[0] (to 5 decimals).
+            # In a circular sense, the “closest index below” is alpha[-1].
+            idx_median = len(alpha) - 1
+        else:
+            idx_median = arr[-1]
+
         idx_lb = idx_median - offset + 1
         idx_ub = idx_median + offset
-        if median.round(5) in alpha.round(5):  # don't count the median per se
+        if np.round(median, 5) in alpha.round(5):  # don't count the median per se
             idx_ub += 1
 
         if idx_ub > n:
diff --git a/pycircstat2/utils.py b/pycircstat2/utils.py
index a5a6d30..2e15887 100644
--- a/pycircstat2/utils.py
+++ b/pycircstat2/utils.py
@@ -160,7 +160,7 @@ def load_data(
         return csv_data
 
 
-def is_within_circular_range(value, lb, ub):
+def is_within_circular_range(value: float, lb: float, ub: float) -> bool:
     """
     Check if a value lies within the circular range [lb, ub].
 
@@ -188,3 +188,42 @@ def is_within_circular_range(value, lb, ub):
     else:
         # Wrapping range
         return value >= lb or value <= ub
+
+def rotate_data(alpha: np.ndarray, angle: float, unit: str = "radian") -> np.ndarray:
+    """
+    Rotate a circular dataset by a given angle, supporting degrees, radians, and hours.
+
+    Parameters
+    ----------
+    alpha : np.ndarray
+        Angles in the specified unit.
+    angle : float
+        Rotation angle in the specified unit.
+    unit : str, optional
+        Unit of measurement ("degree", "radian", or "hour"). Default is "radian".
+
+    Returns
+    -------
+    rotated_alpha : np.ndarray
+        Rotated angles, normalized within the unit's full cycle.
+    """
+    if unit == "degree":
+        n_intervals = 360
+    elif unit == "radian":
+        n_intervals = 2 * np.pi
+    elif unit == "hour":
+        n_intervals = 24
+    else:
+        raise ValueError("Unit must be 'degree', 'radian', or 'hour'.")
+
+    # Convert to radians for consistent computation
+    alpha_rad = data2rad(alpha, k=n_intervals)
+    angle_rad = data2rad(angle, k=n_intervals)
+
+    # Perform rotation and normalize in radians
+    rotated_alpha_rad = angmod(alpha_rad + angle_rad, bounds=[0, 2 * np.pi])
+
+    # Convert back to the original unit
+    rotated_alpha = rad2data(rotated_alpha_rad, k=n_intervals)
+
+    return rotated_alpha
diff --git a/pyproject.toml b/pyproject.toml
index eb0042f..6fc6c9f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -20,9 +20,7 @@ readme = {file = "README.md", content-type = "text/markdown"}
 
 [project.optional-dependencies]
 dev = [
-    "black",
-    "isort",
-    "jupyter",
+    "ruff",
     "pytest",
 ]
 
diff --git a/tests/test_utils.py b/tests/test_utils.py
index 58dfb92..96148f8 100644
--- a/tests/test_utils.py
+++ b/tests/test_utils.py
@@ -1,6 +1,13 @@
 import numpy as np
 
-from pycircstat2.utils import angmod, angular_distance, data2rad, rad2data, time2float
+from pycircstat2.utils import (
+    angmod,
+    angular_distance,
+    data2rad,
+    rad2data,
+    rotate_data,
+    time2float,
+)
 
 
 def test_data2rad():
@@ -84,3 +91,40 @@ def test_angular_distance():
         np.rad2deg(angular_distance(a=np.deg2rad(190), b=np.deg2rad(5))).round(3),
         175,
     )
+
+def test_rotate_data():
+    # Test case 1: Rotating angles in radians
+    alpha = np.array([0, np.pi / 2, np.pi, 3 * np.pi / 2])
+    rotated = rotate_data(alpha, np.pi / 4, unit="radian")
+    expected = np.array([np.pi / 4, 3 * np.pi / 4, 5 * np.pi / 4, 7 * np.pi / 4])
+    np.testing.assert_allclose(rotated, expected, atol=1e-6)
+
+    # Test case 2: Rotating angles in degrees
+    alpha_deg = np.array([0., 90., 180., 270.])
+    rotated_deg = rotate_data(alpha_deg, 45, unit="degree")
+    expected_deg = np.array([45., 135., 225., 315.])
+    np.testing.assert_allclose(rotated_deg, expected_deg, atol=1e-6)
+
+    # Test case 3: Rotating angles in hours
+    alpha_hour = np.array([0., 6., 12., 18.])
+    rotated_hour = rotate_data(alpha_hour, 3, unit="hour")
+    expected_hour = np.array([3., 9., 15., 21.])
+    np.testing.assert_allclose(rotated_hour, expected_hour, atol=1e-6)
+
+    # Test case 4: Rotation that wraps around the cycle
+    rotated_wrap = rotate_data(np.array([350.]), 20, unit="degree")
+    expected_wrap = np.array([10])
+    np.testing.assert_allclose(rotated_wrap, expected_wrap, atol=1e-6)
+
+    rotated_wrap_hour = rotate_data(np.array([22.]), 5, unit="hour")
+    expected_wrap_hour = np.array([3.])
+    np.testing.assert_allclose(rotated_wrap_hour, expected_wrap_hour, atol=1e-6)
+
+    # Test case 5: Rotating an empty array should return an empty array
+    np.testing.assert_array_equal(rotate_data(np.array([]), 45., unit="degree"), np.array([]))
+
+    # Test case 6: Rotating by 0 should return the same values
+    alpha_no_rotation = np.array([10., 50., 120.])
+    np.testing.assert_allclose(
+        rotate_data(alpha_no_rotation, 0, unit="degree"), alpha_no_rotation
+    )