We are currently solving the following issue within our repository. Here is the issue text:
--- BEGIN ISSUE ---
ExcessStatistics class throws an error for `time_step` with an astropy unit 
### Bug description

The current implementation of the `ExcessStatistics` class doesn't allow astropy units for the parameter `time_step`. See example below.

### Expected outcome

Ideally, `total_time_above_threshold`, `average_times` and `rms_times` should have the same units as `time_step`.

### Minimal complete verifiable example

```Python
from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
import astropy.units as u

signal = [0, 0, 2, 2, 0, 4]
thresholds = 1
time_step = 1 * u.s

excess_statistics = ExcessStatistics(signal, thresholds, time_step)
```


### Package versions

v2023.5.1

### Additional context

This is also relevant for PR #2275. One could also add a check whether time units are used for `time_step` or would this be to cumbersome? 

Since I implemented the `ExcessStatistics` class I would be happy to be assigned to this issue.


--- END ISSUE ---

Below are some code segments, each from a relevant file. One or more of these files may contain bugs.
--- BEGIN FILES ---
Path: `plasmapy/analysis/time_series/excess_statistics.py`
Content:
```
1 """
2 Functionality to calculate excess statistics of time series.
3 
4 .. attention::
5 
6    |expect-api-changes|
7 """
8 
9 __all__ = ["ExcessStatistics"]
10 
11 
12 import numbers
13 import numpy as np
14 
15 from collections.abc import Iterable
16 
17 
18 class ExcessStatistics:
19     """
20     Calculate total time, number of upwards crossings, average time and
21     root-mean-square time above given thresholds of a sequence.
22 
23     Parameters
24     ----------
25     signal : 1D |array_like|
26         Signal to be analyzed.
27 
28     thresholds : 1D |array_like|
29         Threshold values.
30 
31     time_step : int
32         Time step of ``signal``.
33 
34     Raises
35     ------
36     `ValueError`
37         If ``time_step`` ≤ 0.
38 
39     Example
40     -------
41     >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
42     >>> signal = [0, 0, 2, 2, 0, 4]
43     >>> thresholds = [1, 3, 5]
44     >>> time_step = 1
45     >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)
46     >>> excess_statistics.total_time_above_threshold
47     [3, 1, 0]
48     >>> excess_statistics.number_of_crossings
49     [2, 1, 0]
50     >>> excess_statistics.average_times
51     [1.5, 1.0, 0]
52     >>> excess_statistics.rms_times
53     [0.5, 0.0, 0]
54     """
55 
56     def __init__(self, signal, thresholds, time_step):
57         if time_step <= 0:
58             raise ValueError("time_step must be positive")
59 
60         # make sure thresholds is an iterable
61         if not isinstance(thresholds, Iterable):
62             thresholds = [thresholds]
63 
64         self._total_time_above_threshold = []
65         self._number_of_crossings = []
66         self._average_times = []
67         self._rms_times = []
68         self.events_per_threshold = {}
69 
70         self._calculate_excess_statistics(signal, thresholds, time_step)
71 
72     def _calculate_excess_statistics(self, signal, thresholds, time_step):
73         for threshold in thresholds:
74             indices_above_threshold = np.where(np.array(signal) > threshold)[0]
75 
76             if len(indices_above_threshold) == 0:
77                 self._times_above_threshold = []
78                 self._total_time_above_threshold.append(0)
79                 self._number_of_crossings.append(0)
80                 self._average_times.append(0)
81                 self._rms_times.append(0)
82 
83             else:
84                 self._total_time_above_threshold.append(
85                     time_step * len(indices_above_threshold)
86                 )
87 
88                 distances_to_next_index = (
89                     indices_above_threshold[1:] - indices_above_threshold[:-1]
90                 )
91                 split_indices = np.where(distances_to_next_index != 1)[0]
92                 event_lengths = np.split(distances_to_next_index, split_indices)
93 
94                 # set correct length for first event
95                 event_lengths[0] = np.append(event_lengths[0], 1)
96 
97                 self._times_above_threshold = [
98                     time_step * len(event_lengths[i]) for i in range(len(event_lengths))
99                 ]
100 
101                 self._number_of_crossings.append(len(event_lengths))
102                 if indices_above_threshold[0] == 0:
103                     # Don't count the first event if there is no crossing.
104                     self._number_of_crossings[-1] -= 1
105 
106                 self._average_times.append(np.mean(self._times_above_threshold))
107                 self._rms_times.append(np.std(self._times_above_threshold))
108 
109             self.events_per_threshold.update({threshold: self._times_above_threshold})
110 
111     def hist(self, bins=32):
112         """
113         Computes the probability density function of the time above each value
114         in ``thresholds``.
115 
116         Parameters
117         ----------
118         bins : int, default: 32
119             The number of bins in the estimation of the PDF above ``thresholds``.
120 
121         Returns
122         -------
123         hist: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )
124             For each value in ``thresholds``, returns the estimated PDF of time
125             above threshold.
126 
127         bin_centers: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )
128             Bin centers for ``hist``.
129 
130         Raises
131         ------
132         `TypeError`
133             If ``bins`` is not a positive integer.
134 
135         Examples
136         --------
137         >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics
138         >>> signal = [0, 0, 2, 0, 4]
139         >>> thresholds = [1, 3, 5]
140         >>> time_step = 1
141         >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)
142         >>> excess_statistics.hist(2)
143         (array([[0., 2.],
144            [0., 2.],
145            [0., 0.]]), array([[0.75, 1.25],
146            [0.75, 1.25],
147            [0.  , 0.  ]]))
148         """
149 
150         if not isinstance(bins, numbers.Integral):
151             raise TypeError("bins must be an integer")
152 
153         hist = np.zeros((len(self.events_per_threshold), bins))
154         bin_centers = np.zeros((len(self.events_per_threshold), bins))
155 
156         for i, threshold in enumerate(self.events_per_threshold.keys()):
157             if len(self.events_per_threshold[threshold]) >= 1:
158                 hist[i, :], bin_edges = np.histogram(
159                     self.events_per_threshold[threshold], bins=bins, density=True
160                 )
161                 bin_centers[i, :] = (bin_edges[1:] + bin_edges[:-1]) / 2
162         return hist, bin_centers
163 
164     @property
165     def total_time_above_threshold(self):
166         """
167         Total time above threshold(s).
168 
169         Returns
170         -------
171         total_time_above_threshold: 1D |array_like|
172             Total time above threshold for each value in ``thresholds``.
173         """
174 
175         return self._total_time_above_threshold
176 
177     @property
178     def number_of_crossings(self):
179         """
180         Total number of upwards crossings for threshold(s).
181 
182         Returns
183         -------
184         number_of_crossings: 1D |array_like|
185             Total number of upwards crossings for each value in ``thresholds``.
186         """
187 
188         return self._number_of_crossings
189 
190     @property
191     def average_times(self):
192         """
193         Average time above threshold(s).
194 
195         Returns
196         -------
197         average_times: 1D |array_like|
198             Average time above each value in ``thresholds``.
199         """
200 
201         return self._average_times
202 
203     @property
204     def rms_times(self):
205         """
206         Root-mean-square values of time above threshold(s).
207 
208         Returns
209         -------
210         rms_times: 1D |array_like|
211             Root-mean-square values of time above each value in ``thresholds``.
212         """
213 
214         return self._rms_times
215
```
--- END FILES ---

Please first localize the bug based on the issue statement, and then generate a patch according to the `git diff` format fenced by three backticks.

Here is an example:

```diff
diff --git a/examples/server_async.py b/examples/server_async.py
--- a/examples/server_async.py
+++ b/examples/server_async.py
@@ -313,4 +313,4 @@
 
 
 if __name__ == "__main__":
-    asyncio.run(run_async_server("."), debug=True)
+    asyncio.run(run_async_server(), debug=True)
diff --git a/examples/server_sync.py b/examples/server_sync.py
--- a/examples/server_sync.py
+++ b/examples/server_sync.py
@@ -313,5 +313,5 @@
 
 
 if __name__ == "__main__":
-    server = run_sync_server(".")
+    server = run_sync_server()
     server.shutdown()

```

diff --git a/plasmapy/analysis/time_series/excess_statistics.py b/plasmapy/analysis/time_series/excess_statistics.py
--- a/plasmapy/analysis/time_series/excess_statistics.py
+++ b/plasmapy/analysis/time_series/excess_statistics.py
@@ -9,6 +9,7 @@
 __all__ = ["ExcessStatistics"]
 
 
+import astropy.units as u
 import numbers
 import numpy as np
 
@@ -98,6 +99,9 @@
                     time_step * len(event_lengths[i]) for i in range(len(event_lengths))
                 ]
 
+                if isinstance(time_step, u.Quantity):
+                    self._times_above_threshold *= time_step.unit
+
                 self._number_of_crossings.append(len(event_lengths))
                 if indices_above_threshold[0] == 0:
                     # Don't count the first event if there is no crossing.

{"golden_diff": "diff --git a/plasmapy/analysis/time_series/excess_statistics.py b/plasmapy/analysis/time_series/excess_statistics.py\n--- a/plasmapy/analysis/time_series/excess_statistics.py\n+++ b/plasmapy/analysis/time_series/excess_statistics.py\n@@ -9,6 +9,7 @@\n __all__ = [\"ExcessStatistics\"]\n \n \n+import astropy.units as u\n import numbers\n import numpy as np\n \n@@ -98,6 +99,9 @@\n                     time_step * len(event_lengths[i]) for i in range(len(event_lengths))\n                 ]\n \n+                if isinstance(time_step, u.Quantity):\n+                    self._times_above_threshold *= time_step.unit\n+\n                 self._number_of_crossings.append(len(event_lengths))\n                 if indices_above_threshold[0] == 0:\n                     # Don't count the first event if there is no crossing.\n", "issue": "ExcessStatistics class throws an error for `time_step` with an astropy unit \n### Bug description\n\nThe current implementation of the `ExcessStatistics` class doesn't allow astropy units for the parameter `time_step`. See example below.\n\n### Expected outcome\n\nIdeally, `total_time_above_threshold`, `average_times` and `rms_times` should have the same units as `time_step`.\n\n### Minimal complete verifiable example\n\n```Python\nfrom plasmapy.analysis.time_series.excess_statistics import ExcessStatistics\r\nimport astropy.units as u\r\n\r\nsignal = [0, 0, 2, 2, 0, 4]\r\nthresholds = 1\r\ntime_step = 1 * u.s\r\n\r\nexcess_statistics = ExcessStatistics(signal, thresholds, time_step)\n```\n\n\n### Package versions\n\nv2023.5.1\n\n### Additional context\n\nThis is also relevant for PR #2275. One could also add a check whether time units are used for `time_step` or would this be to cumbersome? \r\n\r\nSince I implemented the `ExcessStatistics` class I would be happy to be assigned to this issue.\r\n\n", "before_files": [{"content": "\"\"\"\nFunctionality to calculate excess statistics of time series.\n\n.. attention::\n\n   |expect-api-changes|\n\"\"\"\n\n__all__ = [\"ExcessStatistics\"]\n\n\nimport numbers\nimport numpy as np\n\nfrom collections.abc import Iterable\n\n\nclass ExcessStatistics:\n    \"\"\"\n    Calculate total time, number of upwards crossings, average time and\n    root-mean-square time above given thresholds of a sequence.\n\n    Parameters\n    ----------\n    signal : 1D |array_like|\n        Signal to be analyzed.\n\n    thresholds : 1D |array_like|\n        Threshold values.\n\n    time_step : int\n        Time step of ``signal``.\n\n    Raises\n    ------\n    `ValueError`\n        If ``time_step`` \u2264 0.\n\n    Example\n    -------\n    >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics\n    >>> signal = [0, 0, 2, 2, 0, 4]\n    >>> thresholds = [1, 3, 5]\n    >>> time_step = 1\n    >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)\n    >>> excess_statistics.total_time_above_threshold\n    [3, 1, 0]\n    >>> excess_statistics.number_of_crossings\n    [2, 1, 0]\n    >>> excess_statistics.average_times\n    [1.5, 1.0, 0]\n    >>> excess_statistics.rms_times\n    [0.5, 0.0, 0]\n    \"\"\"\n\n    def __init__(self, signal, thresholds, time_step):\n        if time_step <= 0:\n            raise ValueError(\"time_step must be positive\")\n\n        # make sure thresholds is an iterable\n        if not isinstance(thresholds, Iterable):\n            thresholds = [thresholds]\n\n        self._total_time_above_threshold = []\n        self._number_of_crossings = []\n        self._average_times = []\n        self._rms_times = []\n        self.events_per_threshold = {}\n\n        self._calculate_excess_statistics(signal, thresholds, time_step)\n\n    def _calculate_excess_statistics(self, signal, thresholds, time_step):\n        for threshold in thresholds:\n            indices_above_threshold = np.where(np.array(signal) > threshold)[0]\n\n            if len(indices_above_threshold) == 0:\n                self._times_above_threshold = []\n                self._total_time_above_threshold.append(0)\n                self._number_of_crossings.append(0)\n                self._average_times.append(0)\n                self._rms_times.append(0)\n\n            else:\n                self._total_time_above_threshold.append(\n                    time_step * len(indices_above_threshold)\n                )\n\n                distances_to_next_index = (\n                    indices_above_threshold[1:] - indices_above_threshold[:-1]\n                )\n                split_indices = np.where(distances_to_next_index != 1)[0]\n                event_lengths = np.split(distances_to_next_index, split_indices)\n\n                # set correct length for first event\n                event_lengths[0] = np.append(event_lengths[0], 1)\n\n                self._times_above_threshold = [\n                    time_step * len(event_lengths[i]) for i in range(len(event_lengths))\n                ]\n\n                self._number_of_crossings.append(len(event_lengths))\n                if indices_above_threshold[0] == 0:\n                    # Don't count the first event if there is no crossing.\n                    self._number_of_crossings[-1] -= 1\n\n                self._average_times.append(np.mean(self._times_above_threshold))\n                self._rms_times.append(np.std(self._times_above_threshold))\n\n            self.events_per_threshold.update({threshold: self._times_above_threshold})\n\n    def hist(self, bins=32):\n        \"\"\"\n        Computes the probability density function of the time above each value\n        in ``thresholds``.\n\n        Parameters\n        ----------\n        bins : int, default: 32\n            The number of bins in the estimation of the PDF above ``thresholds``.\n\n        Returns\n        -------\n        hist: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )\n            For each value in ``thresholds``, returns the estimated PDF of time\n            above threshold.\n\n        bin_centers: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )\n            Bin centers for ``hist``.\n\n        Raises\n        ------\n        `TypeError`\n            If ``bins`` is not a positive integer.\n\n        Examples\n        --------\n        >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics\n        >>> signal = [0, 0, 2, 0, 4]\n        >>> thresholds = [1, 3, 5]\n        >>> time_step = 1\n        >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)\n        >>> excess_statistics.hist(2)\n        (array([[0., 2.],\n           [0., 2.],\n           [0., 0.]]), array([[0.75, 1.25],\n           [0.75, 1.25],\n           [0.  , 0.  ]]))\n        \"\"\"\n\n        if not isinstance(bins, numbers.Integral):\n            raise TypeError(\"bins must be an integer\")\n\n        hist = np.zeros((len(self.events_per_threshold), bins))\n        bin_centers = np.zeros((len(self.events_per_threshold), bins))\n\n        for i, threshold in enumerate(self.events_per_threshold.keys()):\n            if len(self.events_per_threshold[threshold]) >= 1:\n                hist[i, :], bin_edges = np.histogram(\n                    self.events_per_threshold[threshold], bins=bins, density=True\n                )\n                bin_centers[i, :] = (bin_edges[1:] + bin_edges[:-1]) / 2\n        return hist, bin_centers\n\n    @property\n    def total_time_above_threshold(self):\n        \"\"\"\n        Total time above threshold(s).\n\n        Returns\n        -------\n        total_time_above_threshold: 1D |array_like|\n            Total time above threshold for each value in ``thresholds``.\n        \"\"\"\n\n        return self._total_time_above_threshold\n\n    @property\n    def number_of_crossings(self):\n        \"\"\"\n        Total number of upwards crossings for threshold(s).\n\n        Returns\n        -------\n        number_of_crossings: 1D |array_like|\n            Total number of upwards crossings for each value in ``thresholds``.\n        \"\"\"\n\n        return self._number_of_crossings\n\n    @property\n    def average_times(self):\n        \"\"\"\n        Average time above threshold(s).\n\n        Returns\n        -------\n        average_times: 1D |array_like|\n            Average time above each value in ``thresholds``.\n        \"\"\"\n\n        return self._average_times\n\n    @property\n    def rms_times(self):\n        \"\"\"\n        Root-mean-square values of time above threshold(s).\n\n        Returns\n        -------\n        rms_times: 1D |array_like|\n            Root-mean-square values of time above each value in ``thresholds``.\n        \"\"\"\n\n        return self._rms_times\n", "path": "plasmapy/analysis/time_series/excess_statistics.py"}], "after_files": [{"content": "\"\"\"\nFunctionality to calculate excess statistics of time series.\n\n.. attention::\n\n   |expect-api-changes|\n\"\"\"\n\n__all__ = [\"ExcessStatistics\"]\n\n\nimport astropy.units as u\nimport numbers\nimport numpy as np\n\nfrom collections.abc import Iterable\n\n\nclass ExcessStatistics:\n    \"\"\"\n    Calculate total time, number of upwards crossings, average time and\n    root-mean-square time above given thresholds of a sequence.\n\n    Parameters\n    ----------\n    signal : 1D |array_like|\n        Signal to be analyzed.\n\n    thresholds : 1D |array_like|\n        Threshold values.\n\n    time_step : int\n        Time step of ``signal``.\n\n    Raises\n    ------\n    `ValueError`\n        If ``time_step`` \u2264 0.\n\n    Example\n    -------\n    >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics\n    >>> signal = [0, 0, 2, 2, 0, 4]\n    >>> thresholds = [1, 3, 5]\n    >>> time_step = 1\n    >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)\n    >>> excess_statistics.total_time_above_threshold\n    [3, 1, 0]\n    >>> excess_statistics.number_of_crossings\n    [2, 1, 0]\n    >>> excess_statistics.average_times\n    [1.5, 1.0, 0]\n    >>> excess_statistics.rms_times\n    [0.5, 0.0, 0]\n    \"\"\"\n\n    def __init__(self, signal, thresholds, time_step):\n        if time_step <= 0:\n            raise ValueError(\"time_step must be positive\")\n\n        # make sure thresholds is an iterable\n        if not isinstance(thresholds, Iterable):\n            thresholds = [thresholds]\n\n        self._total_time_above_threshold = []\n        self._number_of_crossings = []\n        self._average_times = []\n        self._rms_times = []\n        self.events_per_threshold = {}\n\n        self._calculate_excess_statistics(signal, thresholds, time_step)\n\n    def _calculate_excess_statistics(self, signal, thresholds, time_step):\n        for threshold in thresholds:\n            indices_above_threshold = np.where(np.array(signal) > threshold)[0]\n\n            if len(indices_above_threshold) == 0:\n                self._times_above_threshold = []\n                self._total_time_above_threshold.append(0)\n                self._number_of_crossings.append(0)\n                self._average_times.append(0)\n                self._rms_times.append(0)\n\n            else:\n                self._total_time_above_threshold.append(\n                    time_step * len(indices_above_threshold)\n                )\n\n                distances_to_next_index = (\n                    indices_above_threshold[1:] - indices_above_threshold[:-1]\n                )\n                split_indices = np.where(distances_to_next_index != 1)[0]\n                event_lengths = np.split(distances_to_next_index, split_indices)\n\n                # set correct length for first event\n                event_lengths[0] = np.append(event_lengths[0], 1)\n\n                self._times_above_threshold = [\n                    time_step * len(event_lengths[i]) for i in range(len(event_lengths))\n                ]\n\n                if isinstance(time_step, u.Quantity):\n                    self._times_above_threshold *= time_step.unit\n\n                self._number_of_crossings.append(len(event_lengths))\n                if indices_above_threshold[0] == 0:\n                    # Don't count the first event if there is no crossing.\n                    self._number_of_crossings[-1] -= 1\n\n                self._average_times.append(np.mean(self._times_above_threshold))\n                self._rms_times.append(np.std(self._times_above_threshold))\n\n            self.events_per_threshold.update({threshold: self._times_above_threshold})\n\n    def hist(self, bins=32):\n        \"\"\"\n        Computes the probability density function of the time above each value\n        in ``thresholds``.\n\n        Parameters\n        ----------\n        bins : int, default: 32\n            The number of bins in the estimation of the PDF above ``thresholds``.\n\n        Returns\n        -------\n        hist: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )\n            For each value in ``thresholds``, returns the estimated PDF of time\n            above threshold.\n\n        bin_centers: 2D `~numpy.ndarray`, shape (``thresholds.size``, ``bins`` )\n            Bin centers for ``hist``.\n\n        Raises\n        ------\n        `TypeError`\n            If ``bins`` is not a positive integer.\n\n        Examples\n        --------\n        >>> from plasmapy.analysis.time_series.excess_statistics import ExcessStatistics\n        >>> signal = [0, 0, 2, 0, 4]\n        >>> thresholds = [1, 3, 5]\n        >>> time_step = 1\n        >>> excess_statistics = ExcessStatistics(signal, thresholds, time_step)\n        >>> excess_statistics.hist(2)\n        (array([[0., 2.],\n           [0., 2.],\n           [0., 0.]]), array([[0.75, 1.25],\n           [0.75, 1.25],\n           [0.  , 0.  ]]))\n        \"\"\"\n\n        if not isinstance(bins, numbers.Integral):\n            raise TypeError(\"bins must be an integer\")\n\n        hist = np.zeros((len(self.events_per_threshold), bins))\n        bin_centers = np.zeros((len(self.events_per_threshold), bins))\n\n        for i, threshold in enumerate(self.events_per_threshold.keys()):\n            if len(self.events_per_threshold[threshold]) >= 1:\n                hist[i, :], bin_edges = np.histogram(\n                    self.events_per_threshold[threshold], bins=bins, density=True\n                )\n                bin_centers[i, :] = (bin_edges[1:] + bin_edges[:-1]) / 2\n        return hist, bin_centers\n\n    @property\n    def total_time_above_threshold(self):\n        \"\"\"\n        Total time above threshold(s).\n\n        Returns\n        -------\n        total_time_above_threshold: 1D |array_like|\n            Total time above threshold for each value in ``thresholds``.\n        \"\"\"\n\n        return self._total_time_above_threshold\n\n    @property\n    def number_of_crossings(self):\n        \"\"\"\n        Total number of upwards crossings for threshold(s).\n\n        Returns\n        -------\n        number_of_crossings: 1D |array_like|\n            Total number of upwards crossings for each value in ``thresholds``.\n        \"\"\"\n\n        return self._number_of_crossings\n\n    @property\n    def average_times(self):\n        \"\"\"\n        Average time above threshold(s).\n\n        Returns\n        -------\n        average_times: 1D |array_like|\n            Average time above each value in ``thresholds``.\n        \"\"\"\n\n        return self._average_times\n\n    @property\n    def rms_times(self):\n        \"\"\"\n        Root-mean-square values of time above threshold(s).\n\n        Returns\n        -------\n        rms_times: 1D |array_like|\n            Root-mean-square values of time above each value in ``thresholds``.\n        \"\"\"\n\n        return self._rms_times\n", "path": "plasmapy/analysis/time_series/excess_statistics.py"}]}