diff --git a/docs/changes/824.bugfix.rst b/docs/changes/824.bugfix.rst
new file mode 100644
index 000000000..67c52da32
--- /dev/null
+++ b/docs/changes/824.bugfix.rst
@@ -0,0 +1 @@
+Substitute np.asarray with np.asanyarray everywhere, to avoid copying memory maps into memory if possible
diff --git a/stingray/base.py b/stingray/base.py
index a70fed444..7e8a35335 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -104,7 +104,7 @@ def __init__(cls, *args, **kwargs) -> None:
     def main_array_length(self):
         if getattr(self, self.main_array_attr, None) is None:
             return 0
-        return np.shape(np.asarray(getattr(self, self.main_array_attr)))[0]
+        return np.shape(np.asanyarray(getattr(self, self.main_array_attr)))[0]
 
     def data_attributes(self) -> list[str]:
         """Clean up the list of attributes, only giving out those pointing to data.
@@ -130,7 +130,7 @@ def data_attributes(self) -> list[str]:
                 and not isinstance(getattr(self.__class__, attr, None), property)
                 and not callable(value := getattr(self, attr))
                 and not isinstance(value, StingrayObject)
-                and not np.asarray(value).dtype == "O"
+                and not np.asanyarray(value).dtype == "O"
             )
         ]
 
@@ -368,7 +368,7 @@ def to_astropy_table(self, no_longdouble=False) -> Table:
         array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            vals = np.asarray(getattr(self, attr))
+            vals = np.asanyarray(getattr(self, attr))
             if no_longdouble:
                 vals = reduce_precision_if_extended(vals)
             data[attr] = vals
@@ -455,7 +455,7 @@ def to_xarray(self) -> Dataset:
         array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            new_data = np.asarray(getattr(self, attr))
+            new_data = np.asanyarray(getattr(self, attr))
             ndim = len(np.shape(new_data))
             if ndim > 1:
                 new_data = ([attr + f"_dim{i}" for i in range(ndim)], new_data)
@@ -520,7 +520,7 @@ def to_pandas(self) -> DataFrame:
         array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            values = np.asarray(getattr(self, attr))
+            values = np.asanyarray(getattr(self, attr))
             ndim = len(np.shape(values))
             if ndim > 1:
                 local_data = make_nd_into_arrays(values, attr)
@@ -758,13 +758,13 @@ def apply_mask(self, mask: npt.ArrayLike, inplace: bool = False, filtered_attrs:
             setattr(
                 new_ts,
                 "_" + self.main_array_attr,
-                copy.deepcopy(np.asarray(getattr(self, self.main_array_attr))[mask]),
+                copy.deepcopy(np.asanyarray(getattr(self, self.main_array_attr))[mask]),
             )
         else:
             setattr(
                 new_ts,
                 self.main_array_attr,
-                copy.deepcopy(np.asarray(getattr(self, self.main_array_attr))[mask]),
+                copy.deepcopy(np.asanyarray(getattr(self, self.main_array_attr))[mask]),
             )
 
         for attr in all_attrs:
@@ -772,7 +772,7 @@ def apply_mask(self, mask: npt.ArrayLike, inplace: bool = False, filtered_attrs:
                 # Eliminate all unfiltered attributes
                 setattr(new_ts, attr, None)
             else:
-                setattr(new_ts, attr, copy.deepcopy(np.asarray(getattr(self, attr))[mask]))
+                setattr(new_ts, attr, copy.deepcopy(np.asanyarray(getattr(self, attr))[mask]))
         return new_ts
 
     def _operation_with_other_obj(
@@ -1030,7 +1030,7 @@ def __neg__(self):
 
         ts_new = copy.deepcopy(self)
         for attr in self._default_operated_attrs():
-            setattr(ts_new, attr, -np.asarray(getattr(self, attr)))
+            setattr(ts_new, attr, -np.asanyarray(getattr(self, attr)))
 
         return ts_new
 
@@ -1215,7 +1215,7 @@ def __init__(
         for kw in other_kw:
             setattr(self, kw, other_kw[kw])
         for kw in array_attrs:
-            new_arr = np.asarray(array_attrs[kw])
+            new_arr = np.asanyarray(array_attrs[kw])
             if self.time.shape[0] != new_arr.shape[0]:
                 raise ValueError(f"Lengths of time and {kw} must be equal.")
             setattr(self, kw, new_arr)
@@ -1246,7 +1246,7 @@ def gti(self):
                 dt1 = self.dt[-1]
             else:
                 dt0 = dt1 = self.dt
-            self._gti = np.asarray([[self._time[0] - dt0 / 2, self._time[-1] + dt1 / 2]])
+            self._gti = np.asanyarray([[self._time[0] - dt0 / 2, self._time[-1] + dt1 / 2]])
         return self._gti
 
     @gti.setter
@@ -1254,7 +1254,7 @@ def gti(self, value):
         if value is None:
             self._gti = None
             return
-        value = np.asarray(value)
+        value = np.asanyarray(value)
         self._gti = value
         self._mask = None
 
@@ -1278,15 +1278,15 @@ def _set_times(self, time, high_precision=False):
             return
         time, _ = interpret_times(time, self.mjdref)
         if not high_precision:
-            self._time = np.asarray(time)
+            self._time = np.asanyarray(time)
         else:
-            self._time = np.asarray(time, dtype=np.longdouble)
+            self._time = np.asanyarray(time, dtype=np.longdouble)
 
     def __str__(self) -> str:
         """Return a string representation of the object."""
         return self.pretty_print(
             attrs_to_apply=["gti", "time", "tstart", "tseg", "tstop"],
-            func_to_apply=lambda x: (np.asarray(x) / 86400 + self.mjdref, "MJD"),
+            func_to_apply=lambda x: (np.asanyarray(x) / 86400 + self.mjdref, "MJD"),
             attrs_to_discard=["_mask", "header"],
         )
 
@@ -1318,7 +1318,7 @@ def _validate_and_format(self, value, attr_name, compare_to_attr):
         """
         if value is None:
             return None
-        value = np.asarray(value)
+        value = np.asanyarray(value)
         if len(value.shape) < 1:
             raise ValueError(f"{attr_name} array must be at least 1D")
         # If the attribute we compare it with is the same and it is currently None, we assign it
@@ -1446,7 +1446,7 @@ def to_astropy_timeseries(self) -> TimeSeries:
         for attr in array_attrs:
             if attr == "time":
                 continue
-            data[attr] = np.asarray(getattr(self, attr))
+            data[attr] = np.asanyarray(getattr(self, attr))
 
         if data == {}:
             data = None
@@ -1489,7 +1489,7 @@ def from_astropy_timeseries(cls, ts: TimeSeries) -> StingrayTimeseries:
 
         new_cls = cls()
         time, mjdref = interpret_times(time, mjdref)
-        new_cls.time = np.asarray(time)  # type: ignore
+        new_cls.time = np.asanyarray(time)  # type: ignore
 
         array_attrs = ts.colnames
         for key, val in ts.meta.items():
@@ -1498,7 +1498,7 @@ def from_astropy_timeseries(cls, ts: TimeSeries) -> StingrayTimeseries:
         for attr in array_attrs:
             if attr == "time":
                 continue
-            setattr(new_cls, attr, np.asarray(ts[attr]))
+            setattr(new_cls, attr, np.asanyarray(ts[attr]))
 
         return new_cls
 
@@ -1553,11 +1553,11 @@ def shift(self, time_shift: float, inplace=False) -> StingrayTimeseries:
             ts = self
         else:
             ts = copy.deepcopy(self)
-        ts.time = np.asarray(ts.time) + time_shift  # type: ignore
+        ts.time = np.asanyarray(ts.time) + time_shift  # type: ignore
         # Pay attention here: if the GTIs are created dynamically while we
         # access the property,
         if ts._gti is not None:
-            ts._gti = np.asarray(ts._gti) + time_shift  # type: ignore
+            ts._gti = np.asanyarray(ts._gti) + time_shift  # type: ignore
 
         return ts
 
@@ -1718,7 +1718,9 @@ def __getitem__(self, index):
             delta_gti_start = new_ts.dt[0] * 0.5
             delta_gti_stop = new_ts.dt[-1] * 0.5
 
-        new_gti = np.asarray([[new_ts.time[0] - delta_gti_start, new_ts.time[-1] + delta_gti_stop]])
+        new_gti = np.asanyarray(
+            [[new_ts.time[0] - delta_gti_start, new_ts.time[-1] + delta_gti_stop]]
+        )
         if step > 1 and delta_gti_start > 0:
             new_gt1 = np.array(list(zip(new_ts.time - new_ts.dt / 2, new_ts.time + new_ts.dt / 2)))
             new_gti = cross_two_gtis(new_gti, new_gt1)
@@ -1797,7 +1799,8 @@ def _truncate_by_index(self, start, stop):
             dtstop = self.dt[-1]
 
         gti = cross_two_gtis(
-            self.gti, np.asarray([[new_ts.time[0] - 0.5 * dtstart, new_ts.time[-1] + 0.5 * dtstop]])
+            self.gti,
+            np.asanyarray([[new_ts.time[0] - 0.5 * dtstart, new_ts.time[-1] + 0.5 * dtstop]]),
         )
 
         new_ts.gti = gti
@@ -2108,7 +2111,7 @@ def rebin(self, dt_new=None, f=None, method="sum"):
         elif f is not None:
             dt_new = f * self.dt
 
-        if np.any(dt_new < np.asarray(self.dt)):
+        if np.any(dt_new < np.asanyarray(self.dt)):
             raise ValueError("The new time resolution must be larger than the old one!")
 
         gti_new = []
@@ -2150,7 +2153,7 @@ def rebin(self, dt_new=None, f=None, method="sum"):
 
         if len(gti_new) == 0:
             raise ValueError("No valid GTIs after rebin.")
-        new_ts.gti = np.asarray(gti_new)
+        new_ts.gti = np.asanyarray(gti_new)
 
         for attr in self.meta_attrs():
             if attr == "dt":
@@ -2654,7 +2657,7 @@ def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs):
                 res = np.nan
             else:
                 lc_filt = self[st:sp]
-                lc_filt.gti = np.asarray([[tst, tsp]])
+                lc_filt.gti = np.asanyarray([[tst, tsp]])
 
                 res = func(lc_filt, **kwargs)
             results.append(res)
diff --git a/stingray/events.py b/stingray/events.py
index ff98ba0eb..6cd11e424 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -232,12 +232,12 @@ def __init__(
         StingrayTimeseries.__init__(
             self,
             time=time,
-            energy=None if energy is None else np.asarray(energy),
+            energy=None if energy is None else np.asanyarray(energy),
             mjdref=mjdref,
             dt=dt,
             notes=notes,
-            gti=np.asarray(gti) if gti is not None else None,
-            pi=None if pi is None else np.asarray(pi),
+            gti=np.asanyarray(gti) if gti is not None else None,
+            pi=None if pi is None else np.asanyarray(pi),
             high_precision=high_precision,
             mission=mission,
             instr=instr,
@@ -367,7 +367,7 @@ def to_lc_iter(self, dt, segment_size=None):
                 self.time[idx_st : idx_end + 1],
                 dt,
                 tstart=st,
-                gti=np.asarray([[st, end]]),
+                gti=np.asanyarray([[st, end]]),
                 tseg=tseg,
                 mjdref=self.mjdref,
                 use_hist=True,
@@ -474,8 +474,8 @@ def simulate_energies(self, spectrum, use_spline=False):
             return
 
         if isinstance(spectrum, list) or isinstance(spectrum, np.ndarray):
-            energy = np.asarray(spectrum)[0]
-            fluxes = np.asarray(spectrum)[1]
+            energy = np.asanyarray(spectrum)[0]
+            fluxes = np.asanyarray(spectrum)[1]
 
             if not isinstance(energy, np.ndarray):
                 raise IndexError("Spectrum must be a 2-d array or list")
diff --git a/stingray/fourier.py b/stingray/fourier.py
index 8c52d3844..613e687eb 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -99,7 +99,7 @@ def integrate_power_in_frequency_range(
     if power_err is None:
         power_err_to_integrate = powers_to_integrate / np.sqrt(m)
     else:
-        power_err_to_integrate = np.asarray(power_err)[frequency_mask]
+        power_err_to_integrate = np.asanyarray(power_err)[frequency_mask]
 
     power_integrated = np.sum((powers_to_integrate - poisson_power) * dfs_to_integrate)
     power_err_integrated = np.sqrt(np.sum((power_err_to_integrate * dfs_to_integrate) ** 2))
@@ -1250,7 +1250,7 @@ def get_average_ctrate(times, gti, segment_size, counts=None):
     Examples
     --------
     >>> times = np.sort(np.random.uniform(0, 1000, 1000))
-    >>> gti = np.asarray([[0, 1000]])
+    >>> gti = np.asanyarray([[0, 1000]])
     >>> counts, _ = np.histogram(times, bins=np.linspace(0, 1000, 11))
     >>> bin_times = np.arange(50, 1000, 100)
     >>> assert get_average_ctrate(bin_times, gti, 1000, counts=counts) == 1.0
@@ -1326,7 +1326,7 @@ def get_flux_iterable_from_segments(
         dt = np.median(np.diff(times[:100]))
 
     if binned:
-        fluxes = np.asarray(fluxes)
+        fluxes = np.asanyarray(fluxes)
         if np.iscomplexobj(fluxes):
             cast_kind = complex
 
@@ -2399,7 +2399,7 @@ def lsft_slow(
         An array of Fourier transformed data.
     """
     y_ = y - np.mean(y)
-    freqs = np.asarray(freqs[np.asarray(freqs) >= 0])
+    freqs = np.asanyarray(freqs[np.asanyarray(freqs) >= 0])
 
     ft_real = np.zeros_like(freqs)
     ft_imag = np.zeros_like(freqs)
diff --git a/stingray/gti.py b/stingray/gti.py
index 007638e66..6a0648f13 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -673,8 +673,8 @@ def cross_two_gtis(gti0, gti1):
     >>> len(newgti)
     0
     """
-    gti0 = join_equal_gti_boundaries(np.asarray(gti0))
-    gti1 = join_equal_gti_boundaries(np.asarray(gti1))
+    gti0 = join_equal_gti_boundaries(np.asanyarray(gti0))
+    gti1 = join_equal_gti_boundaries(np.asanyarray(gti1))
     # Check GTIs
     check_gtis(gti0)
     check_gtis(gti1)
@@ -826,7 +826,7 @@ def get_btis(gtis, start_time=None, stop_time=None):
         if start_time is None or stop_time is None:
             raise ValueError("Empty GTI and no valid start_time " "and stop_time. BAD!")
 
-        return np.asarray([[start_time, stop_time]])
+        return np.asanyarray([[start_time, stop_time]])
     check_gtis(gtis)
 
     start_time = assign_value_if_none(start_time, gtis[0][0])
@@ -844,7 +844,7 @@ def get_btis(gtis, start_time=None, stop_time=None):
     if stop_time > gtis[-1][1]:
         btis.extend([[gtis[-1][1], stop_time]])
 
-    return np.asarray(btis)
+    return np.asanyarray(btis)
 
 
 @jit(nopython=True)
@@ -911,8 +911,8 @@ def check_separate(gti0, gti1):
     >>> assert check_separate(gti0, gti1)
     """
 
-    gti0 = np.asarray(gti0)
-    gti1 = np.asarray(gti1)
+    gti0 = np.asanyarray(gti0)
+    gti1 = np.asanyarray(gti1)
     if len(gti0) == 0 or len(gti1) == 0:
         return True
 
@@ -953,7 +953,7 @@ def join_equal_gti_boundaries(gti, threshold=0.0):
             ng.append(new_gtis[count])
         count += 1
     ng.append(new_gtis[-1])
-    return np.asarray(ng)
+    return np.asanyarray(ng)
 
 
 def merge_gtis(gti_list, strategy):
@@ -1005,7 +1005,7 @@ def merge_gtis(gti_list, strategy):
         return None
 
     if strategy == "none":
-        return np.asarray([[global_min, global_max]])
+        return np.asanyarray([[global_min, global_max]])
 
     if len(all_gti_lists) == 1:
         return all_gti_lists[0]
@@ -1060,8 +1060,8 @@ def append_gtis(gti0, gti1):
     >>> assert np.allclose(append_gtis([[0, 1]], [[1, 3]]), [[0, 3]])
     """
 
-    gti0 = np.asarray(gti0)
-    gti1 = np.asarray(gti1)
+    gti0 = np.asanyarray(gti0)
+    gti1 = np.asanyarray(gti1)
     # Check if independently GTIs are well behaved.
     check_gtis(gti0)
     check_gtis(gti1)
@@ -1123,8 +1123,8 @@ def join_gtis(gti0, gti1):
         The newly created GTI
     """
 
-    gti0 = np.asarray(gti0)
-    gti1 = np.asarray(gti1)
+    gti0 = np.asanyarray(gti0)
+    gti1 = np.asanyarray(gti1)
 
     # Check if independently GTIs are well behaved.
     check_gtis(gti0)
@@ -1135,11 +1135,11 @@ def join_gtis(gti0, gti1):
 
     g0 = gti0.flatten()
     # Opening GTI: type = 1; Closing: type = -1
-    g0_type = np.asarray(
+    g0_type = np.asanyarray(
         list(zip(-np.ones(int(len(g0) / 2), dtype=int), np.ones(int(len(g0) / 2), dtype=int)))
     )
     g1 = gti1.flatten()
-    g1_type = np.asarray(
+    g1_type = np.asanyarray(
         list(zip(-np.ones(int(len(g1) / 2), dtype=int), np.ones(int(len(g1) / 2), dtype=int)))
     )
 
@@ -1332,8 +1332,8 @@ def bin_intervals_from_gtis(gtis, segment_size, time, dt=None, fraction_step=1,
     >>> assert np.allclose(time[start_bins[0]:stop_bins[0]], [0.5, 1.5])
     >>> assert np.allclose(time[start_bins[1]:stop_bins[1]], [2.5, 3.5])
     """
-    time = np.asarray(time)
-    gtis = np.asarray(gtis)
+    time = np.asanyarray(time)
+    gtis = np.asanyarray(gtis)
     if dt is None:
         dt = np.median(np.diff(time))
 
@@ -1433,8 +1433,8 @@ def gti_border_bins(gtis, time, dt=None, epsilon=0.001):
     >>> assert np.allclose(times[start_bins[0]:stop_bins[0]], [0.5, 1.5, 2.5, 3.5, 4.5])
     >>> np.allclose(times[start_bins[1]:stop_bins[1]], [6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5])
     True"""
-    time = np.asarray(time)
-    gtis = np.asarray(gtis)
+    time = np.asanyarray(time)
+    gtis = np.asanyarray(gtis)
     if dt is None:
         dt = np.median(np.diff(time))
 
@@ -1557,8 +1557,8 @@ def generate_indices_of_gti_boundaries(times, gti, dt=0):
     >>> assert np.allclose(v0[:2], gtis[0])
     >>> assert np.allclose(v0[2:], [0, 3])
     """
-    gti = np.asarray(gti)
-    times = np.asarray(times)
+    gti = np.asanyarray(gti)
+    times = np.asanyarray(times)
     startidx, stopidx = gti_border_bins(gti, times, dt=dt)
 
     for s, e, idx0, idx1 in zip(gti[:, 0], gti[:, 1], startidx, stopidx):
@@ -1614,8 +1614,8 @@ def generate_indices_of_segment_boundaries_unbinned(times, gti, segment_size, ch
     >>> # Again: 1.1 is not included in the interval
     >>> assert np.allclose(v1[2:], [3, 4])
     """
-    gti = np.asarray(gti)
-    times = np.asarray(times)
+    gti = np.asanyarray(gti)
+    times = np.asanyarray(times)
 
     start, stop = time_intervals_from_gtis(gti, segment_size)
 
@@ -1632,10 +1632,10 @@ def generate_indices_of_segment_boundaries_unbinned(times, gti, segment_size, ch
         )
     )
 
-    idxs = np.searchsorted(times, all_times)
+    idxs = times.searchsorted(all_times)
     idx_dict = dict([(s, a) for s, a in zip(all_times, idxs)])
-    startidx = np.asarray([idx_dict[s] for s in start])
-    stopidx = np.asarray([idx_dict[s] for s in stop])
+    startidx = np.asanyarray([idx_dict[s] for s in start])
+    stopidx = np.asanyarray([idx_dict[s] for s in stop])
 
     for s, e, idx0, idx1 in zip(start, stop, startidx, stopidx):
         yield s, e, idx0, idx1
@@ -1679,8 +1679,8 @@ def generate_indices_of_segment_boundaries_binned(times, gti, segment_size, dt=N
     >>> assert np.allclose(v0[:2], [0.05, 0.55])
     >>> assert np.allclose(v0[2:], [0, 5])
     """
-    gti = np.asarray(gti)
-    times = np.asarray(times)
+    gti = np.asanyarray(gti)
+    times = np.asanyarray(times)
     startidx, stopidx = bin_intervals_from_gtis(gti, segment_size, times, dt=dt)
 
     if dt is None:
diff --git a/stingray/io.py b/stingray/io.py
index 31c5ac39a..8e3f48419 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -941,7 +941,7 @@ def split_numbers(number, shift=0):
     >>> assert np.allclose(split_numbers(n, -1), (10.0, 2.34))
     """
     if isinstance(number, Iterable):
-        number = np.asarray(number)
+        number = np.asanyarray(number)
         number *= 10**shift
         mods = [math.modf(n) for n in number]
         number_F = [f for f, _ in mods]
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index ad7fed6c6..367082cfc 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -266,11 +266,11 @@ def __init__(
         time, mjdref = interpret_times(time, mjdref=mjdref)
         self.mjdref = mjdref
 
-        time = np.asarray(time)
-        counts = np.asarray(counts)
+        time = np.asanyarray(time)
+        counts = np.asanyarray(counts)
 
         if err is not None:
-            err = np.asarray(err)
+            err = np.asanyarray(err)
 
         if not skip_checks:
             time, counts, err = self.initial_optional_checks(time, counts, err, gti=gti)
@@ -323,7 +323,7 @@ def __init__(
 
         self._gti = None
         if gti is not None:
-            self._gti = np.asarray(gti)
+            self._gti = np.asanyarray(gti)
 
         if os.name == "nt":
             warnings.warn(
@@ -333,22 +333,22 @@ def __init__(
             counts = counts.astype(float)
 
         if input_counts:
-            self._counts = np.asarray(counts)
+            self._counts = np.asanyarray(counts)
             self._counts_err = err
         else:
-            self._countrate = np.asarray(counts)
+            self._countrate = np.asanyarray(counts)
             self._countrate_err = err
 
         if bg_counts is not None:
-            self.bg_counts = np.asarray(bg_counts)
+            self.bg_counts = np.asanyarray(bg_counts)
         else:
             self.bg_counts = None
         if bg_ratio is not None:
-            self.bg_ratio = np.asarray(bg_ratio)
+            self.bg_ratio = np.asanyarray(bg_ratio)
         else:
             self.bg_ratio = None
         if frac_exp is not None:
-            self.frac_exp = np.asarray(frac_exp)
+            self.frac_exp = np.asanyarray(frac_exp)
         else:
             self.frac_exp = None
 
@@ -752,7 +752,7 @@ def __getitem__(self, index):
             new_time = self.time[start:stop:step]
 
             new_gti = [[self.time[start] - 0.5 * self.dt, self.time[stop - 1] + 0.5 * self.dt]]
-            new_gti = np.asarray(new_gti)
+            new_gti = np.asanyarray(new_gti)
             if step > 1:
                 new_gt1 = np.array(list(zip(new_time - self.dt / 2, new_time + self.dt / 2)))
                 new_gti = cross_two_gtis(new_gti, new_gt1)
@@ -860,7 +860,7 @@ def make_lightcurve(toa, dt, tseg=None, tstart=None, gti=None, mjdref=0, use_his
         """
         toa, mjdref = interpret_times(toa, mjdref=mjdref)
 
-        toa = np.sort(np.asarray(toa))
+        toa = np.sort(np.asanyarray(toa))
         # tstart is an optional parameter to set a starting time for
         # the light curve in case this does not coincide with the first photon
         if tstart is None:
@@ -1092,9 +1092,9 @@ def join(self, other, skip_checks=False):
             new_counts = np.concatenate([first_lc.counts, second_lc.counts])
             new_counts_err = np.concatenate([first_lc.counts_err, second_lc.counts_err])
 
-        new_time = np.asarray(new_time)
-        new_counts = np.asarray(new_counts)
-        new_counts_err = np.asarray(new_counts_err)
+        new_time = np.asanyarray(new_time)
+        new_counts = np.asanyarray(new_counts)
+        new_counts_err = np.asanyarray(new_counts_err)
         gti = join_gtis(self.gti, other.gti)
 
         lc_new = Lightcurve(
@@ -1477,7 +1477,7 @@ def _to_astropy_object(self, kind="table", no_longdouble=False):
             "_bin_hi",
         ]:
             if hasattr(self, attr) and getattr(self, attr) is not None:
-                vals = np.asarray(getattr(self, attr))
+                vals = np.asanyarray(getattr(self, attr))
                 if no_longdouble:
                     vals = reduce_precision_if_extended(vals)
                 data[attr.lstrip("_")] = vals
diff --git a/stingray/modeling/gpmodeling.py b/stingray/modeling/gpmodeling.py
index 0b387dfa4..0c5dbe5b1 100644
--- a/stingray/modeling/gpmodeling.py
+++ b/stingray/modeling/gpmodeling.py
@@ -782,10 +782,12 @@ def weighted_posterior_plot(
         )
 
         nbins = max(10, int(jnp.sqrt(self.results.ESS)) + 1)
-        binsx = jnp.linspace(*jnp.percentile(samples_resampled, jnp.asarray([0, 100])), 2 * nbins)
+        binsx = jnp.linspace(
+            *jnp.percentile(samples_resampled, jnp.asanyarray([0, 100])), 2 * nbins
+        )
 
         plt.hist(
-            np.asarray(samples_resampled),
+            np.asanyarray(samples_resampled),
             bins=binsx,
             density=True,
             alpha=1.0,
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
index 73441c016..c17d7a1f7 100644
--- a/stingray/modeling/parameterestimation.py
+++ b/stingray/modeling/parameterestimation.py
@@ -174,7 +174,7 @@ def _compute_covariance(self, lpost, res):
 
         if hasattr(res, "hess_inv"):
             if not isinstance(res.hess_inv, np.ndarray):
-                self.cov = np.asarray(res.hess_inv.todense())
+                self.cov = np.asanyarray(res.hess_inv.todense())
             else:
                 self.cov = res.hess_inv
 
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
index 0fbbece53..891d3709d 100644
--- a/stingray/modeling/tests/test_parameterestimation.py
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -346,7 +346,7 @@ def test_compute_criteria_returns_correct_attributes(self):
     def test_compute_covariance_with_hess_inverse(self):
         self.optres._compute_covariance(self.lpost, self.opt)
 
-        assert np.allclose(self.optres.cov, np.asarray(self.opt.hess_inv))
+        assert np.allclose(self.optres.cov, np.asanyarray(self.opt.hess_inv))
         assert np.allclose(self.optres.err, np.sqrt(np.diag(self.opt.hess_inv)))
 
     @pytest.mark.skipif("comp_hessian")
diff --git a/stingray/modeling/tests/test_scripts.py b/stingray/modeling/tests/test_scripts.py
index 8790a6ed6..67b854805 100644
--- a/stingray/modeling/tests/test_scripts.py
+++ b/stingray/modeling/tests/test_scripts.py
@@ -67,7 +67,7 @@ def setup_class(cls):
         cls.cs.df = cls.cs.freq[1] - cls.cs.freq[0]
         cls.cs.m = 1
 
-        cls.t0 = np.asarray([200.0, 0.5, 0.1, 100.0, 2.0, 1.0, 50.0, 7.5, 0.5, 2.0])
+        cls.t0 = np.asanyarray([200.0, 0.5, 0.1, 100.0, 2.0, 1.0, 50.0, 7.5, 0.5, 2.0])
 
         cls.parest, cls.res = fit_lorentzians(cls.ps, cls.nlor, cls.t0)
 
diff --git a/stingray/power_colors.py b/stingray/power_colors.py
index 71b515e4e..48c7e1961 100644
--- a/stingray/power_colors.py
+++ b/stingray/power_colors.py
@@ -2,6 +2,7 @@
 from scipy.interpolate import interp1d
 import numpy as np
 from .fourier import integrate_power_in_frequency_range
+from .utils import force_array
 from collections.abc import Iterable
 
 DEFAULT_COLOR_CONFIGURATION = {
@@ -177,9 +178,9 @@ def _create_pc_plot(
         ax.set_ylim(yrange)
         return ax
 
-    center = np.log10(np.asarray(configuration["center"]))
-    ax.set_xlim(center[0] + np.asarray(xrange))
-    ax.set_ylim(center[1] + np.asarray(yrange))
+    center = np.log10(np.asanyarray(configuration["center"]))
+    ax.set_xlim(center[0] + np.asanyarray(xrange))
+    ax.set_ylim(center[1] + np.asanyarray(yrange))
 
     for angle in range(0, 360, 20):
         x, y = _hue_line_data(center, np.radians(angle), ref_angle=configuration["ref_angle"])
@@ -264,7 +265,7 @@ def plot_hues(
     plot_spans=False,
     configuration=DEFAULT_COLOR_CONFIGURATION,
 ):
-    hues = hue_from_power_color(pc1, pc2)
+    hues = hue_from_power_color(force_array(pc1), force_array(pc2))
 
     ax = _create_rms_hue_plot(polar=polar, plot_spans=plot_spans, configuration=configuration)
     hues = hues % (np.pi * 2)
@@ -348,12 +349,12 @@ def power_color(
     PC1_err : float
         The error on the second power color
     """
-    freq_edges = np.asarray(freq_edges)
+    freq_edges = np.asanyarray(freq_edges)
     if len(freq_edges) != 5:
         raise ValueError("freq_edges must have 5 elements")
 
-    frequency = np.asarray(frequency)
-    power = np.asarray(power)
+    frequency = np.asanyarray(frequency)
+    power = np.asanyarray(power)
 
     if df is None:
         df = np.median(np.diff(frequency))
@@ -368,7 +369,7 @@ def power_color(
     if power_err is None:
         power_err = power / np.sqrt(m)
     else:
-        power_err = np.asarray(power_err)
+        power_err = np.asanyarray(power_err)
 
     if freqs_to_exclude is not None:
         if len(np.shape(freqs_to_exclude)) == 1:
@@ -461,7 +462,7 @@ def hue_from_power_color(pc0, pc1, center=[4.51920, 0.453724]):
     pc0 = np.log10(pc0)
     pc1 = np.log10(pc1)
 
-    center = np.log10(np.asarray(center))
+    center = np.log10(np.asanyarray(center))
 
     return hue_from_logpower_color(pc0, pc1, center=center)
 
diff --git a/stingray/pulse/accelsearch.py b/stingray/pulse/accelsearch.py
index 8183e7c94..754f61e23 100644
--- a/stingray/pulse/accelsearch.py
+++ b/stingray/pulse/accelsearch.py
@@ -166,7 +166,7 @@ def _convolve_with_response(
     """
     response, j = response_and_j
     r_freqs = np.arange(A.size)
-    if np.asarray(response).size == 1:
+    if np.asanyarray(response).size == 1:
         accel = A
     else:
         accel = convolve(A, response, memout=memout)
@@ -347,10 +347,8 @@ def accelsearch(
         the time and the observation length.
 
     """
-    if not isinstance(times, np.ndarray):
-        times = np.asarray(times)
-    if not isinstance(signal, np.ndarray):
-        signal = np.asarray(signal)
+    times = np.asanyarray(times)
+    signal = np.asanyarray(signal)
 
     dt = times[1] - times[0]
     n_photons = np.sum(signal)
@@ -359,7 +357,7 @@ def accelsearch(
             "Data contain multiple GTIs. Bad time intervals will be "
             "filled with the mean of the signal."
         )
-        gti = np.asarray(gti)
+        gti = np.asanyarray(gti)
         # Fill in the data with a constant outside GTIs
         gti_mask = create_gti_mask(times, gti)
         expo_fraction = np.count_nonzero(gti_mask) / len(gti_mask)
@@ -496,8 +494,8 @@ def interbin_fft(freq, fft):
     """
     import numpy as np
 
-    freq = np.asarray(freq)
-    fft = np.asarray(fft)
+    freq = np.asanyarray(freq)
+    fft = np.asanyarray(fft)
 
     neglast = freq[-1] < 0
     if neglast:
diff --git a/stingray/pulse/modeling.py b/stingray/pulse/modeling.py
index 521aabd26..dc864d33f 100644
--- a/stingray/pulse/modeling.py
+++ b/stingray/pulse/modeling.py
@@ -97,7 +97,7 @@ def sinc_square_deriv(x, amplitude=1.0, mean=0.0, width=1.0):
         * (x * np.cos((x - mean) / width) - np.sin((x - mean) / width))
         / ((x - mean) / width) ** 2
     )
-    d_x = np.asarray(d_x)
+    d_x = np.asanyarray(d_x)
     d_amplitude = sinc((x - mean) / width) ** 2
     d_x[x_is_zero] = 0
 
diff --git a/stingray/pulse/pulsar.py b/stingray/pulse/pulsar.py
index dac17457b..4e84ab8bd 100644
--- a/stingray/pulse/pulsar.py
+++ b/stingray/pulse/pulsar.py
@@ -272,7 +272,7 @@ def fold_events(times, *frequency_derivatives, **opts):
     if gti is None:
         gti = [[times[0], times[-1]]]
     # Be safe if gtis are a list
-    gti = np.asarray(gti)
+    gti = np.asanyarray(gti)
     ref_time = _default_value_if_no_key(opts, "ref_time", 0)
     expocorr = _default_value_if_no_key(opts, "expocorr", False)
 
@@ -648,7 +648,7 @@ def z_n(data, n, datatype="events", err=None, norm=None):
     z2_n : float
         The Z^2_n statistics of the events.
     """
-    data = np.asarray(data)
+    data = np.asanyarray(data)
 
     if norm is not None:
         warnings.warn(
diff --git a/stingray/pulse/search.py b/stingray/pulse/search.py
index 29e779eef..1906cb155 100644
--- a/stingray/pulse/search.py
+++ b/stingray/pulse/search.py
@@ -31,7 +31,7 @@ def _folding_search(
     stat_func, times, frequencies, segment_size=np.inf, use_times=False, fdots=0, **kwargs
 ):
     fgrid, fdgrid = np.meshgrid(
-        np.asarray(frequencies).astype(np.float64), np.asarray(fdots).astype(np.float64)
+        np.asanyarray(frequencies).astype(np.float64), np.asanyarray(fdots).astype(np.float64)
     )
     stats = np.zeros_like(fgrid)
     times = (times - times[0]).astype(np.float64)
@@ -391,8 +391,8 @@ def search_best_peaks(x, stat, threshold):
     []
 
     """
-    stat = np.asarray(stat)
-    x = np.asarray(x)
+    stat = np.asanyarray(stat)
+    x = np.asanyarray(x)
     peaks = stat >= threshold
     regions = contiguous_regions(peaks)
     if len(regions) == 0:
diff --git a/stingray/simulator/base.py b/stingray/simulator/base.py
index 2a4bb8254..db3659337 100644
--- a/stingray/simulator/base.py
+++ b/stingray/simulator/base.py
@@ -111,9 +111,9 @@ def simulate_times_from_count_array(time, counts, gti, dt, use_spline=False):
         ...
     ValueError: simulate_times can only work with...
     """
-    time = np.asarray(time)
-    counts = np.asarray(counts).astype(float)
-    gti = np.asarray(gti)
+    time = np.asanyarray(time)
+    counts = np.asanyarray(counts).astype(float)
+    gti = np.asanyarray(gti)
     kind = "linear"
     if use_spline and time.size > 2:
         kind = "cubic"
@@ -220,7 +220,7 @@ def simulate_with_inverse_cdf(
     >>> assert np.all(np.diff(vals)) >= 0
 
     """
-    binned_pdf = np.asarray(binned_pdf).astype(float)
+    binned_pdf = np.asanyarray(binned_pdf).astype(float)
 
     if x_range is None:
         x_range = [0, 1]
diff --git a/stingray/simulator/transfer.py b/stingray/simulator/transfer.py
index f56711350..afbfe67f2 100644
--- a/stingray/simulator/transfer.py
+++ b/stingray/simulator/transfer.py
@@ -49,7 +49,7 @@ class TransferFunction(object):
     """
 
     def __init__(self, data, dt=1, de=1, tstart=0, estart=0, time=None, energy=None):
-        self.data = np.asarray(data)
+        self.data = np.asanyarray(data)
         self.dt = dt
         self.de = de
         self.tstart = tstart
diff --git a/stingray/spectroscopy.py b/stingray/spectroscopy.py
index 8b395dac8..2f407024d 100644
--- a/stingray/spectroscopy.py
+++ b/stingray/spectroscopy.py
@@ -39,8 +39,8 @@ def load_lc_fits(file, counts_type=True):
     meta = lc_fits.meta
     dt = meta["DT"]
 
-    ref = np.asarray(lc_fits["REF"].T, dtype=np.float64)
-    ci = np.asarray(lc_fits["CI"].T, dtype=np.float64)
+    ref = np.asanyarray(lc_fits["REF"].T, dtype=np.float64)
+    ci = np.asanyarray(lc_fits["CI"].T, dtype=np.float64)
 
     if not counts_type:
         print(dt)
diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py
index 37128486d..145d194a3 100644
--- a/stingray/tests/test_base.py
+++ b/stingray/tests/test_base.py
@@ -40,7 +40,7 @@ def __init__(self, dummy=None):
 class TestStingrayObject:
     @classmethod
     def setup_class(cls):
-        cls.arr = np.asarray([4, 5, 2])
+        cls.arr = np.asanyarray([4, 5, 2])
         sting_obj = DummyStingrayObj(cls.arr)
         sting_obj.pardulas = [3.0 + 1.0j, 2.0j, 1.0 + 0.0j]
         sting_obj.sebadas = [[0, 1], [2, 3], [4, 5]]
@@ -376,16 +376,16 @@ def setup_class(cls):
             mjdref=59777.000,
             array_attrs=dict(guefus=cls.arr),
             parafritus="bonus!",
-            panesapa=np.asarray([[41, 25], [98, 3]]),
-            gti=np.asarray([[-0.5, 10.5]]),
+            panesapa=np.asanyarray([[41, 25], [98, 3]]),
+            gti=np.asanyarray([[-0.5, 10.5]]),
         )
         sting_obj_highp = StingrayTimeseries(
             time=cls.time,
             mjdref=59777.000,
             array_attrs=dict(guefus=cls.arr),
             parafritus="bonus!",
-            panesapa=np.asarray([[41, 25], [98, 3]]),
-            gti=np.asarray([[-0.5, 10.5]]),
+            panesapa=np.asanyarray([[41, 25], [98, 3]]),
+            gti=np.asanyarray([[-0.5, 10.5]]),
             high_precision=True,
         )
         cls.sting_obj = sting_obj
@@ -507,7 +507,8 @@ def test_what_is_array_and_what_is_not(self):
             [0, 3],
             gti=[[0.5, 1.5], [2.5, 3.5]],
             array_attrs=dict(
-                panesapa=np.asarray([[41, 25], [98, 3]]), _panesapa=np.asarray([[41, 25], [98, 3]])
+                panesapa=np.asanyarray([[41, 25], [98, 3]]),
+                _panesapa=np.asanyarray([[41, 25], [98, 3]]),
             ),
             dt=1,
         )
@@ -620,7 +621,7 @@ def test_slice(self):
     @pytest.mark.parametrize("inplace", [True, False])
     def test_apply_gti(self, inplace):
         so = copy.deepcopy(self.sting_obj)
-        so.gti = np.asarray([[-0.1, 2.1]])
+        so.gti = np.asanyarray([[-0.1, 2.1]])
         so2 = so.apply_gtis()
         if inplace:
             assert so2 is so
@@ -847,7 +848,7 @@ def test_rebin_less_than_dt(self):
 
     def test_sort(self):
         times = [2, 1, 3, 4]
-        blah = np.asarray([40, 10, 20, 5])
+        blah = np.asanyarray([40, 10, 20, 5])
         bleh = [4, 1, 2, 0.5]
         mjdref = 57000
 
@@ -1045,11 +1046,11 @@ class TestStingrayTimeseriesSubclass:
     def setup_class(cls):
         cls.arr = [4, 5, 2]
         sting_obj = DummyStingrayTs(cls.arr)
-        sting_obj.time = np.asarray([0, 1, 2])
+        sting_obj.time = np.asanyarray([0, 1, 2])
         sting_obj.mjdref = 59777.000
         sting_obj.parafritus = "bonus!"
-        sting_obj.panesapa = np.asarray([[41, 25], [98, 3]])
-        sting_obj.gti = np.asarray([[-0.5, 2.5]])
+        sting_obj.panesapa = np.asanyarray([[41, 25], [98, 3]])
+        sting_obj.gti = np.asanyarray([[-0.5, 2.5]])
         cls.sting_obj = sting_obj
 
     def test_print(self, capsys):
@@ -1258,9 +1259,9 @@ def test_overlapping_join_change_mjdref(self):
             )
         with pytest.warns(UserWarning, match="The time array is not sorted."):
             ts_other = StingrayTimeseries(
-                time=np.asarray([5.1, 7, 6.1, 6.11, 10.1]) + 86400,
+                time=np.asanyarray([5.1, 7, 6.1, 6.11, 10.1]) + 86400,
                 energy=[2, 3, 8, 1, 2],
-                gti=np.asarray([[5, 7], [8, 10]]) + 86400,
+                gti=np.asanyarray([[5, 7], [8, 10]]) + 86400,
                 mjdref=57000,
             )
         with pytest.warns(UserWarning, match="Attribute mjdref is different"):
@@ -1331,7 +1332,9 @@ def test_no_btis_returns_copy(self):
     def test_event_like(self):
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]])
+        ev_like_filt.gti = np.asanyarray(
+            [[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]]
+        )
         ev_new = ev_like_filt.fill_bad_time_intervals()
 
         assert np.allclose(ev_new.gti, self.gti)
@@ -1347,7 +1350,7 @@ def test_event_like(self):
     def test_no_counts_in_buffer(self):
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 490], [491, 498], [500, 505], [510, 520], [522, 700]])
+        ev_like_filt.gti = np.asanyarray([[0, 490], [491, 498], [500, 505], [510, 520], [522, 700]])
 
         # I empty out two GTIs
         bad = (ev_like_filt.time > 490) & (ev_like_filt.time < 510)
@@ -1359,7 +1362,9 @@ def test_no_counts_in_buffer(self):
     def test_lc_like(self):
         lc_like_filt = copy.deepcopy(self.lc_like)
         # I introduce a small gap in the GTIs
-        lc_like_filt.gti = np.asarray([[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]])
+        lc_like_filt.gti = np.asanyarray(
+            [[0, 498], [500, 520], [522, 700], [702, 900], [950, 1000]]
+        )
         lc_new = lc_like_filt.fill_bad_time_intervals()
         assert np.allclose(lc_new.gti, self.gti)
 
@@ -1378,7 +1383,7 @@ def test_lc_like(self):
     def test_ignore_attrs_ev_like(self):
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        ev_like_filt.gti = np.asanyarray([[0, 498], [500, 900], [950, 1000]])
         ev_new0 = ev_like_filt.fill_bad_time_intervals(seed=1234)
         ev_new1 = ev_like_filt.fill_bad_time_intervals(seed=1234, attrs_to_randomize=["energy"])
         assert np.allclose(ev_new0.gti, ev_new1.gti)
@@ -1391,7 +1396,7 @@ def test_ignore_attrs_ev_like(self):
     def test_ignore_attrs_lc_like(self):
         lc_like_filt = copy.deepcopy(self.lc_like)
         # I introduce a small gap in the GTIs
-        lc_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        lc_like_filt.gti = np.asanyarray([[0, 498], [500, 900], [950, 1000]])
         lc_new0 = lc_like_filt.fill_bad_time_intervals(seed=1234)
         lc_new1 = lc_like_filt.fill_bad_time_intervals(seed=1234, attrs_to_randomize=["counts"])
         assert np.allclose(lc_new0.gti, lc_new1.gti)
@@ -1404,7 +1409,7 @@ def test_ignore_attrs_lc_like(self):
     def test_forcing_non_uniform(self):
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        ev_like_filt.gti = np.asanyarray([[0, 498], [500, 900], [950, 1000]])
         # Results should be exactly the same
         ev_new0 = ev_like_filt.fill_bad_time_intervals(even_sampling=False, seed=201903)
         ev_new1 = ev_like_filt.fill_bad_time_intervals(even_sampling=None, seed=201903)
@@ -1414,7 +1419,7 @@ def test_forcing_non_uniform(self):
     def test_forcing_uniform(self):
         lc_like_filt = copy.deepcopy(self.lc_like)
         # I introduce a small gap in the GTIs
-        lc_like_filt.gti = np.asarray([[0, 498], [500, 900], [950, 1000]])
+        lc_like_filt.gti = np.asanyarray([[0, 498], [500, 900], [950, 1000]])
         # Results should be exactly the same
         lc_new0 = lc_like_filt.fill_bad_time_intervals(even_sampling=True, seed=201903)
         lc_new1 = lc_like_filt.fill_bad_time_intervals(even_sampling=None, seed=201903)
@@ -1424,26 +1429,26 @@ def test_forcing_uniform(self):
     def test_bti_close_to_edge_event_like(self):
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 0.5], [1, 900], [950, 1000]])
+        ev_like_filt.gti = np.asanyarray([[0, 0.5], [1, 900], [950, 1000]])
         ev_new = ev_like_filt.fill_bad_time_intervals()
         assert np.allclose(ev_new.gti, self.gti)
 
         ev_like_filt = copy.deepcopy(self.ev_like)
         # I introduce a small gap in the GTIs
-        ev_like_filt.gti = np.asarray([[0, 900], [950, 999], [999.5, 1000]])
+        ev_like_filt.gti = np.asanyarray([[0, 900], [950, 999], [999.5, 1000]])
         ev_new = ev_like_filt.fill_bad_time_intervals()
         assert np.allclose(ev_new.gti, self.gti)
 
     def test_bti_close_to_edge_lc_like(self):
         lc_like_filt = copy.deepcopy(self.lc_like)
         # I introduce a small gap in the GTIs
-        lc_like_filt.gti = np.asarray([[0, 0.5], [1, 900], [950, 1000]])
+        lc_like_filt.gti = np.asanyarray([[0, 0.5], [1, 900], [950, 1000]])
         lc_new = lc_like_filt.fill_bad_time_intervals()
         assert np.allclose(lc_new.gti, self.gti)
 
         lc_like_filt = copy.deepcopy(self.lc_like)
         # I introduce a small gap in the GTIs
-        lc_like_filt.gti = np.asarray([[0, 900], [950, 999], [999.5, 1000]])
+        lc_like_filt.gti = np.asanyarray([[0, 900], [950, 999], [999.5, 1000]])
         lc_new = lc_like_filt.fill_bad_time_intervals()
         assert np.allclose(lc_new.gti, self.gti)
 
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index 3f909edb6..2ae34d406 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -25,7 +25,7 @@ def setup_class(self):
         self.counts = [3000, 2000, 2200, 3600]
         self.counts_flat = [3000, 3000, 3000, 3000]
         self.spectrum = [[1, 2, 3, 4, 5, 6], [1000, 2040, 1000, 3000, 4020, 2070]]
-        self.gti = np.asarray([[0, 4]])
+        self.gti = np.asanyarray([[0, 4]])
 
     def test_warn_wrong_keywords(self):
         with pytest.warns(UserWarning) as record:
@@ -511,9 +511,9 @@ def test_overlapping_join_change_mjdref(self):
             )
         with pytest.warns(UserWarning, match="The time array is not sorted."):
             ev_other = EventList(
-                time=np.asarray([5.1, 7, 6.1, 6.11, 10.1]) + 86400,
+                time=np.asanyarray([5.1, 7, 6.1, 6.11, 10.1]) + 86400,
                 energy=[2, 3, 8, 1, 2],
-                gti=np.asarray([[5, 7], [8, 10]]) + 86400,
+                gti=np.asanyarray([[5, 7], [8, 10]]) + 86400,
                 mjdref=57000,
             )
         with pytest.warns(UserWarning, match="Attribute mjdref is different"):
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index beffcf1a4..acd3fcfd8 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -75,7 +75,7 @@ def test_flux_iterables(dtype):
     times = np.arange(4)
     fluxes = np.ones(4).astype(dtype)
     errors = np.ones(4).astype(dtype) * np.sqrt(2)
-    gti = np.asarray([[-0.5, 3.5]])
+    gti = np.asanyarray([[-0.5, 3.5]])
     iter = get_flux_iterable_from_segments(times, gti, 2, n_bin=None, fluxes=fluxes, errors=errors)
     cast_kind = float
     if np.iscomplexobj(fluxes):
@@ -90,7 +90,7 @@ def test_flux_iterables(dtype):
 def test_avg_pds_imperfect_lc_size():
     times = np.arange(100)
     fluxes = np.ones(100).astype(float)
-    gti = np.asarray([[-0.5, 99.5]])
+    gti = np.asanyarray([[-0.5, 99.5]])
     segment_size = 5.99
     dt = 1
     res = avg_pds_from_timeseries(times, gti, segment_size, dt, fluxes=fluxes)
@@ -101,7 +101,7 @@ def test_avg_pds_imperfect_lc_size():
 def test_avg_pds_from_events_warns():
     times = np.arange(100)
     fluxes = np.ones(100).astype(float)
-    gti = np.asarray([[-0.5, 99.5]])
+    gti = np.asanyarray([[-0.5, 99.5]])
     segment_size = 5.99
     dt = 1
     with pytest.warns(DeprecationWarning, match="avg_pds_from_events is deprecated"):
@@ -114,7 +114,7 @@ def test_avg_cs_imperfect_lc_size():
     times1 = times2 = np.arange(100)
     fluxes1 = np.ones(100).astype(float)
     fluxes2 = np.ones(100).astype(float)
-    gti = np.asarray([[-0.5, 99.5]])
+    gti = np.asanyarray([[-0.5, 99.5]])
     segment_size = 5.99
     dt = 1
     res = avg_cs_from_timeseries(
@@ -128,7 +128,7 @@ def test_avg_cs_from_events_warns():
     times1 = times2 = np.arange(100)
     fluxes1 = np.ones(100).astype(float)
     fluxes2 = np.ones(100).astype(float)
-    gti = np.asarray([[-0.5, 99.5]])
+    gti = np.asanyarray([[-0.5, 99.5]])
     segment_size = 5.99
     dt = 1
     with pytest.warns(DeprecationWarning, match="avg_cs_from_events is deprecated"):
@@ -224,7 +224,7 @@ def setup_class(cls):
         cls.N = np.rint(cls.length / cls.dt).astype(int)
         cls.dt = cls.length / cls.N
         cls.times = np.sort(rng.uniform(0, cls.length, int(cls.length * cls.ctrate)))
-        cls.gti = np.asarray([[0, cls.length]])
+        cls.gti = np.asanyarray([[0, cls.length]])
         cls.counts, bins = np.histogram(cls.times, bins=np.linspace(0, cls.length, cls.N + 1))
         cls.errs = np.ones_like(cls.counts) * np.sqrt(cls.ctrate)
         cls.bin_times = (bins[:-1] + bins[1:]) / 2
diff --git a/stingray/tests/test_gti.py b/stingray/tests/test_gti.py
index b3cf4647a..1cf29cda1 100644
--- a/stingray/tests/test_gti.py
+++ b/stingray/tests/test_gti.py
@@ -14,7 +14,6 @@
 
 
 class TestGTI(object):
-
     """Real unit tests."""
 
     def test_crossgti1(self):
@@ -107,7 +106,7 @@ def test_bti_empty_valid(self):
         gti = np.array([])
 
         bti = get_btis(gti, start_time=0, stop_time=1)
-        assert np.allclose(bti, np.asarray([[0, 1]]))
+        assert np.allclose(bti, np.asanyarray([[0, 1]]))
 
     def test_bti_fail(self):
         gti = np.array([])
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index b247b4bc1..8dc528b2c 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -808,7 +808,7 @@ def test_concatenate(self):
         count1 = [50, 60, 70, 80, 90]
         gti0 = [[0.5, 4.5]]
         gti1 = [[4.5, 9.5]]
-        lc0 = Lightcurve(time0, counts=count0, err=np.asarray(count0) / 2, dt=1, gti=gti0)
+        lc0 = Lightcurve(time0, counts=count0, err=np.asanyarray(count0) / 2, dt=1, gti=gti0)
         lc1 = Lightcurve(time1, counts=count1, dt=1, gti=gti1)
         with pytest.warns(UserWarning) as record:
             lc = lc0.concatenate(lc1)
diff --git a/stingray/tests/test_power_colors.py b/stingray/tests/test_power_colors.py
index 0f62e499a..4bca1a63f 100644
--- a/stingray/tests/test_power_colors.py
+++ b/stingray/tests/test_power_colors.py
@@ -78,7 +78,7 @@ def test_with_power_err(self):
 
     def test_hue(self):
         center = (4.51920, 0.453724)
-        log_center = np.log10(np.asarray(center))
+        log_center = np.log10(np.asanyarray(center))
         for angle in np.radians(np.arange(0, 380, 20)):
             factor = rng.uniform(0.1, 10)
             x = factor * np.cos(3 / 4 * np.pi - angle) + log_center[0]
diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
index b6c0e309d..e37678942 100644
--- a/stingray/tests/test_varenergyspectrum.py
+++ b/stingray/tests/test_varenergyspectrum.py
@@ -142,7 +142,9 @@ def test_construct_lightcurves_no_exclude(self):
 
     def test_construct_lightcurves_pi(self):
         events = EventList(
-            [0.09, 0.21, 0.23, 0.32, 0.4, 0.54], pi=np.asarray([0, 0, 0, 0, 1, 1]), gti=[[0, 0.65]]
+            [0.09, 0.21, 0.23, 0.32, 0.4, 0.54],
+            pi=np.asanyarray([0, 0, 0, 0, 1, 1]),
+            gti=[[0, 0.65]],
         )
         vespec = DummyVarEnergy(
             events, [0.0, 10000], (0, 1, 2, "lin"), [0.5, 1.1], use_pi=True, bin_time=0.1
@@ -178,7 +180,7 @@ def setup_class(cls):
         # No need for huge count rates
         flux = data / 40
         times = np.arange(data.size) * cls.bin_time
-        gti = np.asarray([[0, data.size * cls.bin_time]])
+        gti = np.asanyarray([[0, data.size * cls.bin_time]])
         test_lc = Lightcurve(
             times, flux, err_dist="gauss", gti=gti, dt=cls.bin_time, skip_checks=True
         )
diff --git a/stingray/utils.py b/stingray/utils.py
index 264fc17da..ee08796a0 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -112,7 +112,7 @@ def mad(data, c=0.6745, axis=None):
             Axis along which to calculate ``mad``. Default is ``0``, can also
             be ``None``
         """
-        data = np.asarray(data)
+        data = np.asanyarray(data)
         if axis is not None:
             center = np.apply_over_axes(np.median, data, axis)
         else:
@@ -144,6 +144,32 @@ def mad(data, c=0.6745, axis=None):
 ]
 
 
+def force_array(x):
+    """Convert an input to a numpy array.
+
+    If it is not iterable, convert to a 1-element array.
+
+    Parameters
+    ----------
+    x : iterable or number
+        The input to be converted
+
+    Returns
+    -------
+    x : numpy.ndarray
+        The input converted to a numpy array
+
+    Examples
+    --------
+    >>> assert isinstance(force_array(1), np.ndarray)
+    >>> assert isinstance(force_array([1]), np.ndarray)
+    """
+    if not isinstance(x, Iterable):
+        x = [x]
+
+    return np.asanyarray(x)
+
+
 @njit()
 def heaviside(x):
     """Heaviside function. Returns 1 if x>0, and 0 otherwise.
@@ -208,7 +234,7 @@ def make_nd_into_arrays(array: np.ndarray, label: str) -> dict:
     >>> assert np.array_equal(data["test_dim0_0"], [1, 5])
     """
     data = {}
-    array = np.asarray(array)
+    array = np.asanyarray(array)
     shape = np.shape(array)
     ndim = len(shape)
     if ndim <= 1:
@@ -336,7 +362,7 @@ def check_isallfinite(array):
         # Numba is very picky about the type of the input array. If an exception
         # occurs in the numba-compiled function, use the default Numpy implementation.
         try:
-            return _check_isallfinite_numba(np.asarray(array))
+            return _check_isallfinite_numba(np.asanyarray(array))
         except Exception:
             pass
     return bool(np.all(np.isfinite(array)))
@@ -359,7 +385,7 @@ def is_sorted(array):
         True if the array is sorted, False otherwise
     """
 
-    array = np.asarray(array)
+    array = np.asanyarray(array)
     # If the array is empty or has only one element, it is sorted
     if array.size <= 1:
         return True
@@ -402,7 +428,7 @@ def _is_sorted_numba(array):
 
 
 def _root_squared_mean(array):
-    array = np.asarray(array)
+    array = np.asanyarray(array)
     return np.sqrt(np.sum(array**2)) / array.size
 
 
@@ -482,11 +508,11 @@ def rebin_data(x, y, dx_new, yerr=None, method="sum", dx=None):
     >>> assert np.allclose(ybinerr, 0.05)
     """
 
-    y = np.asarray(y)
+    y = np.asanyarray(y)
     if yerr is None:
         yerr = np.zeros_like(y)
     else:
-        yerr = np.asarray(yerr)
+        yerr = np.asanyarray(yerr)
 
     if isinstance(dx, Iterable):
         dx_old = dx
@@ -619,9 +645,9 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
     """
 
     dx_init = apply_function_if_none(dx, np.diff(x), np.median)
-    x = np.asarray(x)
-    y = np.asarray(y)
-    y_err = np.asarray(apply_function_if_none(y_err, y, np.zeros_like))
+    x = np.asanyarray(x)
+    y = np.asanyarray(y)
+    y_err = np.asanyarray(apply_function_if_none(y_err, y, np.zeros_like))
 
     if x.shape[0] != y.shape[0]:
         raise ValueError("x and y must be of the same length!")
@@ -639,7 +665,7 @@ def rebin_data_log(x, y, f, y_err=None, dx=None):
         binx_for_stats.append(binx_for_stats[-1] + dx * (1.0 + f))
         dx = binx_for_stats[-1] - binx_for_stats[-2]
 
-    binx_for_stats = np.asarray(binx_for_stats)
+    binx_for_stats = np.asanyarray(binx_for_stats)
 
     real = y.real
     real_err = y_err.real
@@ -1292,10 +1318,10 @@ def poisson_symmetrical_errors(counts):
     >>> from astropy.stats import poisson_conf_interval
     >>> counts = np.random.randint(0, 1000, 100)
     >>> # ---- Do it without the lookup table ----
-    >>> err_low, err_high = poisson_conf_interval(np.asarray(counts),
+    >>> err_low, err_high = poisson_conf_interval(np.asanyarray(counts),
     ...                 interval='frequentist-confidence', sigma=1)
-    >>> err_low -= np.asarray(counts)
-    >>> err_high -= np.asarray(counts)
+    >>> err_low -= np.asanyarray(counts)
+    >>> err_high -= np.asanyarray(counts)
     >>> err = (np.absolute(err_low) + np.absolute(err_high))/2.0
     >>> # Do it with this function
     >>> err_thisfun = poisson_symmetrical_errors(counts)
@@ -1304,14 +1330,14 @@ def poisson_symmetrical_errors(counts):
     """
     from astropy.stats import poisson_conf_interval
 
-    counts_int = np.asarray(counts, dtype=np.int64)
+    counts_int = np.asanyarray(counts, dtype=np.int64)
     count_values = np.nonzero(np.bincount(counts_int))[0]
     err_low, err_high = poisson_conf_interval(
         count_values, interval="frequentist-confidence", sigma=1
     )
     # calculate approximately symmetric uncertainties
-    err_low -= np.asarray(count_values)
-    err_high -= np.asarray(count_values)
+    err_low -= np.asanyarray(count_values)
+    err_high -= np.asanyarray(count_values)
     err = (np.absolute(err_low) + np.absolute(err_high)) / 2.0
 
     idxs = np.searchsorted(count_values, counts_int)
@@ -1462,8 +1488,8 @@ def check_allclose_and_print(
     try:
         assert np.allclose(v1, v2, rtol, atol)
     except Exception as e:
-        v1 = np.asarray(v1)
-        v2 = np.asarray(v2)
+        v1 = np.asanyarray(v1)
+        v2 = np.asanyarray(v2)
         bad = np.abs(v1 - v2) >= (atol + rtol * np.abs(v2))
 
         raise AssertionError(
@@ -1610,7 +1636,7 @@ def hist1d_numba_seq(a, bins, range, use_memmap=False, tmp=None):
     """
     hist_arr = _allocate_array_or_memmap((bins,), a.dtype, use_memmap=use_memmap, tmp=tmp)
 
-    return _hist1d_numba_seq(hist_arr, a, bins, np.asarray(range))
+    return _hist1d_numba_seq(hist_arr, a, bins, np.asanyarray(range))
 
 
 @njit(nogil=True, parallel=False)
@@ -1674,7 +1700,7 @@ def hist2d_numba_seq(x, y, bins, range, use_memmap=False, tmp=None):
     """
 
     H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
-    return _hist2d_numba_seq(H, np.array([x, y]), np.asarray(list(bins)), np.asarray(range))
+    return _hist2d_numba_seq(H, np.array([x, y]), np.asanyarray(list(bins)), np.asanyarray(range))
 
 
 @njit(nogil=True, parallel=False)
@@ -1737,7 +1763,9 @@ def hist3d_numba_seq(tracks, bins, range, use_memmap=False, tmp=None):
     """
     H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
 
-    return _hist3d_numba_seq(H, np.asarray(tracks), np.asarray(list(bins)), np.asarray(range))
+    return _hist3d_numba_seq(
+        H, np.asanyarray(tracks), np.asanyarray(list(bins)), np.asanyarray(range)
+    )
 
 
 @njit(nogil=True, parallel=False)
@@ -1810,7 +1838,7 @@ def hist1d_numba_seq_weight(a, weights, bins, range, use_memmap=False, tmp=None)
     else:
         hist_arr = np.zeros((bins,), dtype=a.dtype)
 
-    return _hist1d_numba_seq_weight(hist_arr, a, weights, bins, np.asarray(range))
+    return _hist1d_numba_seq_weight(hist_arr, a, weights, bins, np.asanyarray(range))
 
 
 @njit(nogil=True, parallel=False)
@@ -1879,8 +1907,8 @@ def hist2d_numba_seq_weight(x, y, weights, bins, range, use_memmap=False, tmp=No
         H,
         np.array([x, y]),
         weights,
-        np.asarray(list(bins)),
-        np.asarray(range),
+        np.asanyarray(list(bins)),
+        np.asanyarray(range),
     )
 
 
@@ -1948,10 +1976,10 @@ def hist3d_numba_seq_weight(tracks, weights, bins, range, use_memmap=False, tmp=
     H = _allocate_array_or_memmap(bins, np.double, use_memmap=use_memmap, tmp=tmp)
     return _hist3d_numba_seq_weight(
         H,
-        np.asarray(tracks),
+        np.asanyarray(tracks),
         weights,
-        np.asarray(list(bins)),
-        np.asarray(range),
+        np.asanyarray(list(bins)),
+        np.asanyarray(range),
     )
 
 
@@ -2038,7 +2066,7 @@ def histnd_numba_seq(tracks, bins, range, use_memmap=False, tmp=None):
     ...                       range=np.array([[0., 1.], [2., 3.], [4., 5.]]))
     >>> assert np.all(H == Hn)
     """
-    tracks = np.asarray(tracks)
+    tracks = np.asanyarray(tracks)
     H = _allocate_array_or_memmap(bins, np.uint64, use_memmap=use_memmap, tmp=tmp)
     slice_int = np.zeros(len(bins), dtype=np.uint64)
 
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index 9b1d7a355..4369d80f6 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -62,12 +62,12 @@ def get_non_overlapping_ref_band(channel_band, ref_band):
     >>> new_ref = get_non_overlapping_ref_band([0, 1], [[2, 3]])
     >>> assert np.allclose(new_ref, [[2, 3]])
     """
-    channel_band = np.asarray(channel_band)
-    ref_band = np.asarray(ref_band)
+    channel_band = np.asanyarray(channel_band)
+    ref_band = np.asanyarray(ref_band)
     if len(ref_band.shape) <= 1:
-        ref_band = np.asarray([ref_band])
+        ref_band = np.asanyarray([ref_band])
     if check_separate(ref_band, [channel_band]):
-        return np.asarray(ref_band)
+        return np.asanyarray(ref_band)
     not_channel_band = [
         [0, channel_band[0]],
         [channel_band[1], np.max([np.max(ref_band), channel_band[1] + 1])],
@@ -211,14 +211,14 @@ def __init__(
         if isinstance(energy_spec, tuple):
             energies = _decode_energy_specification(energy_spec)
         else:
-            energies = np.asarray(energy_spec)
+            energies = np.asanyarray(energy_spec)
 
         self.energy_intervals = list(zip(energies[0:-1], energies[1:]))
 
-        self.ref_band = np.asarray(assign_value_if_none(ref_band, [0, np.inf]))
+        self.ref_band = np.asanyarray(assign_value_if_none(ref_band, [0, np.inf]))
 
         if len(self.ref_band.shape) <= 1:
-            self.ref_band = np.asarray([self.ref_band])
+            self.ref_band = np.asanyarray([self.ref_band])
 
         self.segment_size = self.delta_nu = None
         if segment_size is not None:
@@ -1057,7 +1057,7 @@ def _spectrum_function(self):
 
             # Convert the cross spectrum to a covariance.
             cov, cov_e = cross_to_covariance(
-                np.asarray([Cmean, Ce]), mean_ref_power, ref_power_noise, delta_nu
+                np.asanyarray([Cmean, Ce]), mean_ref_power, ref_power_noise, delta_nu
             )
 
             meanrate = mean / self.bin_time