Repository: StingraySoftware/stingray
Pull request #754 (2023-09-15 10:10:23)

Generalized lightcurves
An attempt at more general time series objects. Not taking for granted counts or count rates, having multiple array attributes, etc.
It might be, for example, spectroscopic or polarimetric light curves with multiple measurements per time point

Documentation for these changes can be found at https://github.com/StingraySoftware/notebooks/pull/73

~(Depends on #774)~ Merged

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) Report
> Merging [#754](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (137d28c) into [main](https://app.codecov.io/gh/StingraySoftware/stingray/commit/5b338dcb6313be5f8420abb2a652052ea08d02c0?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) (5b338dc) will **decrease** coverage by `4.49%`.
> The diff coverage is `85.45%`.

```diff
@@            Coverage Diff             @@
##             main     #754      +/-   ##
==========================================
- Coverage   97.13%   92.65%   -4.49%     
==========================================
  Files          42       41       -1     
  Lines        7932     7600     -332     
==========================================
- Hits         7705     7042     -663     
- Misses        227      558     +331     
```


| [Files Changed](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) | Coverage Δ | |
|---|---|---|
| [stingray/\_\_init\_\_.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvX19pbml0X18ucHk=) | `100.00% <ø> (ø)` | |
| [stingray/base.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvYmFzZS5weQ==) | `79.29% <78.70%> (-17.61%)` | :arrow_down: |
| [stingray/lightcurve.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvbGlnaHRjdXJ2ZS5weQ==) | `96.13% <94.44%> (-1.77%)` | :arrow_down: |
| [stingray/crossspectrum.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvY3Jvc3NzcGVjdHJ1bS5weQ==) | `98.62% <95.23%> (+0.40%)` | :arrow_up: |
| [stingray/utils.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvdXRpbHMucHk=) | `94.65% <97.22%> (-3.85%)` | :arrow_down: |
| [stingray/events.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZXZlbnRzLnB5) | `99.03% <100.00%> (-0.97%)` | :arrow_down: |
| [stingray/fourier.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvZm91cmllci5weQ==) | `99.76% <100.00%> (+0.01%)` | :arrow_up: |
| [stingray/powerspectrum.py](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware#diff-c3RpbmdyYXkvcG93ZXJzcGVjdHJ1bS5weQ==) | `99.67% <100.00%> (+0.22%)` | :arrow_up: |

... and [9 files with indirect coverage changes](https://app.codecov.io/gh/StingraySoftware/stingray/pull/754/indirect-changes?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)

:mega: We’re building smart automated test selection to slash your CI/CD build times. [Learn more](https://about.codecov.io/iterative-testing/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware)

pep8speaks: Hello @matteobachetti! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`stingray/base.py`](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py):

> [Line 77:21](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L77): [E713](https://duckduckgo.com/?q=pep8%20E713) test for membership should be 'not in'
> [Line 552:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L552): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (103 > 100 characters)
> [Line 573:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L573): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (105 > 100 characters)
> [Line 578:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L578): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (103 > 100 characters)
> [Line 604:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L604): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (102 > 100 characters)
> [Line 655:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L655): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (105 > 100 characters)
> [Line 734:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L734): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (105 > 100 characters)
> [Line 1109:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L1109): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (103 > 100 characters)
> [Line 1130:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L1130): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (105 > 100 characters)
> [Line 1135:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/base.py#L1135): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (103 > 100 characters)











* In the file [`stingray/tests/test_base.py`](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/tests/test_base.py):

> [Line 311:101](https://github.com/StingraySoftware/stingray/blob/636a610342bfb0e1770d08939bc5e1e9a2fe1322/stingray/tests/test_base.py#L311): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (102 > 100 characters)










matteobachetti: Fixing things for making them work with the Generalized Light curves produced a monster PR. And there are more modifications to do. I will spin off the fixes to smaller PRs, and try to keep this PR on topic.
matteobachetti: @mgullik thanks for the review! I fixed all the docstrings. Indeed, many  of those methods were moved from `EventList` and `Lightcurve`, and the docstrings were left with the original class names.

In general, all child classes (`EventList`, `Lightcurve`) inherit the methods from `StingrayTimeseries` or re-write them if incompatible. `join` is partially re-written, but e.g. `EventList.join` does call `_join_timeseries`. However, note that there is a check that you _cannot_ join different types. If you start from an event list, you cannot merge it with a light curve, even if both of them are children of `StingrayTimeseries`.

Changes of MJDREF are accounted for. Note that you cannot be a StingrayTimeseries without defining MJDREF, at least as 0 (It is a field automatically defined in the `StingrayTimeseries` class, see line 1085).

Patch:
diff --git a/docs/changes/754.feature.rst b/docs/changes/754.feature.rst
new file mode 100644
index 00000000..0f9f3f00
--- /dev/null
+++ b/docs/changes/754.feature.rst
@@ -0,0 +1,1 @@
+Make StingrayTimeseries into a generalized light curve, with a less strict naming but implementing much of the underlying computing useful for Lightcurve as well.
diff --git a/docs/index.rst b/docs/index.rst
index d2db9152..c6ae9caf 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -123,9 +123,9 @@ For the Gaussian Process modeling in `stingray.modeling.gpmodeling`, you'll need
 + etils
 + typing_extensions
 
-Most of these are installed via ``pip``, but if you have an Nvidia GPU available, you'll want to take special care 
-following the installation instructions for jax and tensorflow(-probability) in order to enable GPU support and 
-take advantage of those speed-ups. 
+Most of these are installed via ``pip``, but if you have an Nvidia GPU available, you'll want to take special care
+following the installation instructions for jax and tensorflow(-probability) in order to enable GPU support and
+take advantage of those speed-ups.
 
 For development work, you will need the following extra libraries:
 
@@ -271,6 +271,7 @@ Using Stingray
 
 Getting started
 ---------------
+
 .. toctree::
    :maxdepth: 2
 
@@ -284,6 +285,7 @@ Advanced
 .. toctree::
    :maxdepth: 2
 
+   timeseries
    modeling
    simulator
    deadtime
diff --git a/docs/notebooks b/docs/notebooks
index 60107493..28e9b776 160000
--- a/docs/notebooks
+++ b/docs/notebooks
@@ -1,1 +1,1 @@
-Subproject commit 601074930f8f32899f0816f5bab4624c4d12592d
+Subproject commit 28e9b776bf76f53fa7ed836c6485e6e92719a0a5
diff --git a/docs/timeseries.rst b/docs/timeseries.rst
new file mode 100644
index 00000000..1a62b586
--- /dev/null
+++ b/docs/timeseries.rst
@@ -0,0 +1,8 @@
+Working with more generic time series
+*************************************
+
+.. toctree::
+   :maxdepth: 2
+
+   notebooks/StingrayTimeseries/StingrayTimeseries Tutorial.ipynb
+   notebooks/StingrayTimeseries/Working with weights and polarization.ipynb
diff --git a/setup.cfg b/setup.cfg
index 342c3916..28e8a3c0 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -108,6 +108,7 @@ filterwarnings =
     ignore:.*HIERARCH card will be created.*:
     ignore:.*FigureCanvasAgg is non-interactive.*:UserWarning
     ignore:.*jax.* deprecated:DeprecationWarning:
+    ignore:.*Converting to lower precision.*:UserWarning
 
 ;addopts = --disable-warnings
 
diff --git a/stingray/__init__.py b/stingray/__init__.py
index 4131ebd4..84c6eeca 100644
--- a/stingray/__init__.py
+++ b/stingray/__init__.py
@@ -9,6 +9,7 @@ from ._astropy_init import *
 
 # For egg_info test builds to pass, put package imports here.
 if not _ASTROPY_SETUP_:
+    from stingray.base import *
     from stingray.events import *
     from stingray.lightcurve import *
     from stingray.utils import *
diff --git a/stingray/base.py b/stingray/base.py
index 28b0ad8f..390fbd67 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -2,15 +2,22 @@
 from __future__ import annotations
 
 from collections.abc import Iterable
+from collections import OrderedDict
+
 import pickle
 import warnings
 import copy
+import os
+import logging
 
 import numpy as np
 from astropy.table import Table
 from astropy.time import Time, TimeDelta
 from astropy.units import Quantity
 
+from stingray.utils import sqsum
+from .io import _can_save_longdouble, _can_serialize_meta
+
 from typing import TYPE_CHECKING, Type, TypeVar, Union
 
 if TYPE_CHECKING:
@@ -24,6 +31,26 @@ if TYPE_CHECKING:
     Tso = TypeVar("Tso", bound="StingrayObject")
 
 
+__all__ = [
+    "convert_table_attrs_to_lowercase",
+    "interpret_times",
+    "reduce_precision_if_extended",
+    "StingrayObject",
+    "StingrayTimeseries",
+]
+
+
+def convert_table_attrs_to_lowercase(table: Table) -> Table:
+    """Convert the column names of an Astropy Table to lowercase."""
+    new_table = Table()
+    for col in table.colnames:
+        new_table[col.lower()] = table[col]
+    for key in table.meta.keys():
+        new_table.meta[key.lower()] = table.meta[key]
+
+    return new_table
+
+
 class StingrayObject(object):
     """This base class defines some general-purpose utilities.
 
@@ -37,24 +64,65 @@ class StingrayObject(object):
     just by defining an attribute called ``main_array_attr``, be able to perform
     the operations above, with no additional effort.
 
-    ``main_array_attr`` is, e.g. ``time`` for :class:`EventList` and
+    ``main_array_attr`` is, e.g. ``time`` for :class:`StingrayTimeseries` and
     :class:`Lightcurve`, ``freq`` for :class:`Crossspectrum`, ``energy`` for
     :class:`VarEnergySpectrum`, and so on. It is the array with wich all other
     attributes are compared: if they are of the same shape, they get saved as
     columns of the table/dataframe, otherwise as metadata.
     """
 
+    not_array_attr: list = []
+
     def __init__(cls, *args, **kwargs) -> None:
         if not hasattr(cls, "main_array_attr"):
             raise RuntimeError(
                 "A StingrayObject needs to have the main_array_attr attribute specified"
             )
 
+    @property
+    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]
+
+    def data_attributes(self) -> list[str]:
+        """Clean up the list of attributes, only giving out those pointing to data.
+
+        List all the attributes that point directly to valid data. This method goes through all the
+        attributes of the class, eliminating methods, properties, and attributes that are complicated
+        to serialize such as other ``StingrayObject``, or arrays of objects.
+
+        This function does not make difference between array-like data and scalar data.
+
+        Returns
+        -------
+        data_attributes : list of str
+            List of attributes pointing to data that are not methods, properties,
+            or other ``StingrayObject`` instances.
+        """
+        return [
+            attr
+            for attr in dir(self)
+            if (
+                not attr.startswith("__")
+                and attr not in ["main_array_attr", "not_array_attr"]
+                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"
+            )
+        ]
+
     def array_attrs(self) -> list[str]:
         """List the names of the array attributes of the Stingray Object.
 
         By array attributes, we mean the ones with the same size and shape as
         ``main_array_attr`` (e.g. ``time`` in ``EventList``)
+
+        Returns
+        -------
+        attributes : list of str
+            List of array attributes.
         """
 
         main_attr = getattr(self, getattr(self, "main_array_attr"))
@@ -63,34 +131,192 @@ class StingrayObject(object):
 
         return [
             attr
-            for attr in dir(self)
+            for attr in self.data_attributes()
             if (
-                isinstance(getattr(self, attr), Iterable)
-                and np.shape(getattr(self, attr)) == np.shape(main_attr)
+                not attr.startswith("_")
+                and not attr == self.main_array_attr
+                and isinstance(getattr(self, attr), Iterable)
+                and attr not in self.not_array_attr
+                and not isinstance(getattr(self, attr), str)
+                and np.shape(getattr(self, attr))[0] == np.shape(main_attr)[0]
             )
         ]
 
+    def internal_array_attrs(self) -> list[str]:
+        """List the names of the internal array attributes of the Stingray Object.
+
+        These are array attributes that can be set by properties, and are generally indicated
+        by an underscore followed by the name of the property that links to it (E.g.
+        ``_counts`` in ``Lightcurve``).
+        By array attributes, we mean the ones with the same size and shape as
+        ``main_array_attr`` (e.g. ``time`` in ``EventList``)
+
+        Returns
+        -------
+        attributes : list of str
+            List of internal array attributes.
+        """
+
+        main_attr = getattr(self, "main_array_attr")
+        main_attr_value = getattr(self, main_attr)
+        if main_attr_value is None:
+            return []
+
+        all_attrs = []
+        for attr in self.data_attributes():
+            if (
+                not attr == "_" + self.main_array_attr  # e.g. _time in lightcurve
+                and attr not in ["_" + a for a in self.not_array_attr]
+                and not np.isscalar(value := getattr(self, attr))
+                and value is not None
+                and not np.size(value) == 0
+                and attr.startswith("_")
+                and np.shape(value)[0] == np.shape(main_attr_value)[0]
+            ):
+                all_attrs.append(attr)
+
+        return all_attrs
+
     def meta_attrs(self) -> list[str]:
         """List the names of the meta attributes of the Stingray Object.
 
         By array attributes, we mean the ones with a different size and shape
         than ``main_array_attr`` (e.g. ``time`` in ``EventList``)
+
+        Returns
+        -------
+        attributes : list of str
+            List of meta attributes.
         """
-        array_attrs = self.array_attrs()
-        return [
+        array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
+
+        all_meta_attrs = [
             attr
-            for attr in dir(self)
-            if (
-                attr not in array_attrs
-                and not attr.startswith("_")
-                # Use new assignment expression (PEP 572). I'm testing that
-                # self.attribute is not callable, and assigning its value to
-                # the variable attr_value for further checks
-                and not callable(attr_value := getattr(self, attr))
-                # a way to avoid EventLists, Lightcurves, etc.
-                and not hasattr(attr_value, "meta_attrs")
-            )
+            for attr in self.data_attributes()
+            if (attr not in array_attrs and not attr.startswith("_"))
         ]
+        if self.not_array_attr is not None and len(self.not_array_attr) >= 1:
+            all_meta_attrs += self.not_array_attr
+        return all_meta_attrs
+
+    def dict(self) -> dict:
+        """Return a dictionary representation of the object."""
+
+        main_attr = self.main_array_attr
+        meta_attrs = self.meta_attrs()
+        array_attrs = self.array_attrs()
+        internal_array_attrs = self.internal_array_attrs()
+
+        results = OrderedDict()
+        results[main_attr] = getattr(self, main_attr)
+
+        for attr in internal_array_attrs:
+            if isinstance(getattr(self.__class__, attr.lstrip("_"), None), property):
+                attr = attr.lstrip("_")
+            results[attr] = getattr(self, attr)
+
+        for attr in array_attrs:
+            results[attr] = getattr(self, attr)
+
+        for attr in meta_attrs:
+            results[attr] = getattr(self, attr)
+
+        return results
+
+    def pretty_print(self, func_to_apply=None, attrs_to_apply=[], attrs_to_discard=[]) -> str:
+        """Return a pretty-printed string representation of the object.
+
+        This is useful for debugging, and for interactive use.
+
+        Optional parameters
+        -------------------
+        func_to_apply : function
+            A function that modifies the attributes listed in ``attrs_to_apply``.
+            It must return the modified attributes and a label to be printed.
+            If ``None``, no function is applied.
+        attrs_to_apply : list of str
+            Attributes to be modified by ``func_to_apply``.
+        attrs_to_discard : list of str
+            Attributes to be discarded from the output.
+        """
+        print(self.__class__.__name__)
+        print("_" * len(self.__class__.__name__))
+        items = self.dict()
+        results = ""
+        np.set_printoptions(threshold=3, edgeitems=1)
+        for attr in items.keys():
+            if attr in attrs_to_discard:
+                continue
+            value = items[attr]
+            label = f"{attr:<15}: {items[attr]}"
+
+            if isinstance(value, Iterable) and not isinstance(value, str):
+                size = np.shape(value)
+                if len(size) == 1:
+                    label += f" (size {size[0]})"
+                else:
+                    label += f" (shape {size})"
+
+            if func_to_apply is not None and attr in attrs_to_apply:
+                new_value, new_label = func_to_apply(items[attr])
+                label += f"\n{attr + ' (' +new_label + ')':<15}: {new_value}"
+
+            results += label + "\n"
+        return results
+
+    def __str__(self) -> str:
+        """Return a string representation of the object."""
+        return self.pretty_print()
+
+    def __eq__(self, other_ts):
+        """Compare two :class:`StingrayObject` instances with ``==``.
+
+        All attributes (internal, array, meta) are compared.
+        """
+
+        if not isinstance(other_ts, type(self)):
+            raise ValueError(f"{type(self)} can only be compared with a {type(self)} Object")
+
+        self_arr_attrs = self.array_attrs()
+        other_arr_attrs = other_ts.array_attrs()
+
+        if not set(self_arr_attrs) == set(other_arr_attrs):
+            return False
+
+        self_meta_attrs = self.meta_attrs()
+        other_meta_attrs = other_ts.meta_attrs()
+
+        if not set(self_meta_attrs) == set(other_meta_attrs):
+            return False
+
+        for attr in self.meta_attrs():
+            # They are either both scalar or arrays
+            if np.isscalar(getattr(self, attr)) != np.isscalar(getattr(other_ts, attr)):
+                return False
+
+            if np.isscalar(getattr(self, attr)):
+                if not getattr(self, attr, None) == getattr(other_ts, attr, None):
+                    return False
+            else:
+                if not np.array_equal(getattr(self, attr, None), getattr(other_ts, attr, None)):
+                    return False
+
+        for attr in self.array_attrs():
+            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr)):
+                return False
+
+        for attr in self.internal_array_attrs():
+            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr)):
+                return False
+
+        return True
+
+    def _default_operated_attrs(self):
+        operated_attrs = [attr for attr in self.array_attrs() if not attr.endswith("_err")]
+        return operated_attrs
+
+    def _default_error_attrs(self):
+        return [attr for attr in self.array_attrs() if attr.endswith("_err")]
 
     def get_meta_dict(self) -> dict:
         """Give a dictionary with all non-None meta attrs of the object."""
@@ -102,22 +328,37 @@ class StingrayObject(object):
                 meta_dict[key] = val
         return meta_dict
 
-    def to_astropy_table(self) -> Table:
+    def to_astropy_table(self, no_longdouble=False) -> Table:
         """Create an Astropy Table from a ``StingrayObject``
 
         Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
         ``EventList``) are converted into columns, while meta attributes
         (``mjdref``, ``gti``, etc.) are saved into the ``meta`` dictionary.
+
+        Other Parameters
+        ----------------
+        no_longdouble : bool
+            If True, reduce the precision of longdouble arrays to double precision.
+            This needs to be done in some cases, e.g. when the table is to be saved
+            in an architecture not supporting extended precision (e.g. ARM), but can
+            also be useful when an extended precision is not needed.
         """
         data = {}
-        array_attrs = self.array_attrs()
+        array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
+            vals = np.asarray(getattr(self, attr))
+            if no_longdouble:
+                vals = reduce_precision_if_extended(vals)
+            data[attr] = vals
 
         ts = Table(data)
+        meta_dict = self.get_meta_dict()
+        for attr in meta_dict.keys():
+            if no_longdouble:
+                meta_dict[attr] = reduce_precision_if_extended(meta_dict[attr])
 
-        ts.meta.update(self.get_meta_dict())
+        ts.meta.update(meta_dict)
 
         return ts
 
@@ -131,9 +372,6 @@ class StingrayObject(object):
         new object, while the attributes in ds.attrs will
         form the new meta attributes of the object.
 
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of the wanted StingrayObject (e.g.
-        ``time``, ``pi``, etc. for ``EventList``)
         """
         cls = cls()
 
@@ -144,11 +382,14 @@ class StingrayObject(object):
         array_attrs = ts.colnames
 
         # Set the main attribute first
-        mainarray = np.array(ts[cls.main_array_attr])  # type: ignore
-        setattr(cls, cls.main_array_attr, mainarray)  # type: ignore
+        for attr in array_attrs:
+            if attr.lower() == cls.main_array_attr:  # type: ignore
+                mainarray = np.array(ts[attr])  # type: ignore
+                setattr(cls, cls.main_array_attr, mainarray)  # type: ignore
+                break
 
         for attr in array_attrs:
-            if attr == cls.main_array_attr:  # type: ignore
+            if attr.lower() == cls.main_array_attr:  # type: ignore
                 continue
             setattr(cls, attr.lower(), np.array(ts[attr]))
 
@@ -167,10 +408,14 @@ class StingrayObject(object):
         from xarray import Dataset
 
         data = {}
-        array_attrs = self.array_attrs()
+        array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
+            new_data = np.asarray(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)
+            data[attr] = new_data
 
         ts = Dataset(data)
 
@@ -187,10 +432,6 @@ class StingrayObject(object):
         The rest of columns will form the array attributes of the
         new object, while the attributes in ds.attrs will
         form the new meta attributes of the object.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of the wanted StingrayObject (e.g.
-        ``time``, ``pi``, etc. for ``EventList``)
         """
         cls = cls()
 
@@ -198,19 +439,20 @@ class StingrayObject(object):
             # return an empty object
             return cls
 
-        array_attrs = ts.coords
-
         # Set the main attribute first
         mainarray = np.array(ts[cls.main_array_attr])  # type: ignore
         setattr(cls, cls.main_array_attr, mainarray)  # type: ignore
 
-        for attr in array_attrs:
-            if attr == cls.main_array_attr:  # type: ignore
-                continue
-            setattr(cls, attr, np.array(ts[attr]))
+        all_array_attrs = []
+        for array_attrs in [ts.coords, ts.data_vars]:
+            for attr in array_attrs:
+                all_array_attrs.append(attr)
+                if attr == cls.main_array_attr:  # type: ignore
+                    continue
+                setattr(cls, attr, np.array(ts[attr]))
 
         for key, val in ts.attrs.items():
-            if key not in array_attrs:
+            if key not in all_array_attrs:
                 setattr(cls, key, val)
 
         return cls
@@ -221,14 +463,27 @@ class StingrayObject(object):
         Array attributes (e.g. ``time``, ``pi``, ``energy``, etc. for
         ``EventList``) are converted into columns, while meta attributes
         (``mjdref``, ``gti``, etc.) are saved into the ``ds.attrs`` dictionary.
+
+        Since pandas does not support n-D data, multi-dimensional arrays are
+        converted into columns before the conversion, with names ``<colname>_dimN_M_K`` etc.
+
+        See documentation of `make_nd_into_arrays` for details.
+
         """
         from pandas import DataFrame
+        from .utils import make_nd_into_arrays
 
         data = {}
-        array_attrs = self.array_attrs()
+        array_attrs = self.array_attrs() + [self.main_array_attr] + self.internal_array_attrs()
 
         for attr in array_attrs:
-            data[attr] = np.asarray(getattr(self, attr))
+            values = np.asarray(getattr(self, attr))
+            ndim = len(np.shape(values))
+            if ndim > 1:
+                local_data = make_nd_into_arrays(values, attr)
+            else:
+                local_data = {attr: values}
+            data.update(local_data)
 
         ts = DataFrame(data)
 
@@ -246,11 +501,15 @@ class StingrayObject(object):
         new object, while the attributes in ds.attrs will
         form the new meta attributes of the object.
 
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of the wanted StingrayObject (e.g.
-        ``time``, ``pi``, etc. for ``EventList``)
+        Since pandas does not support n-D data, multi-dimensional arrays can be
+        specified as ``<colname>_dimN_M_K`` etc.
+
+        See documentation of `make_1d_arrays_into_nd` for details.
 
         """
+        import re
+        from .utils import make_1d_arrays_into_nd
+
         cls = cls()
 
         if len(ts) == 0:
@@ -263,10 +522,17 @@ class StingrayObject(object):
         mainarray = np.array(ts[cls.main_array_attr])  # type: ignore
         setattr(cls, cls.main_array_attr, mainarray)  # type: ignore
 
+        nd_attrs = []
         for attr in array_attrs:
             if attr == cls.main_array_attr:  # type: ignore
                 continue
-            setattr(cls, attr, np.array(ts[attr]))
+            if "_dim" in attr:
+                nd_attrs.append(re.sub("_dim[0-9].*", "", attr))
+            else:
+                setattr(cls, attr, np.array(ts[attr]))
+
+        for attr in list(set(nd_attrs)):
+            setattr(cls, attr, make_1d_arrays_into_nd(ts, attr))
 
         for key, val in ts.attrs.items():
             if key not in array_attrs:
@@ -323,7 +589,7 @@ class StingrayObject(object):
         elif fmt.lower() == "ascii":
             fmt = "ascii.ecsv"
 
-        ts = Table.read(filename, format=fmt)
+        ts = convert_table_attrs_to_lowercase(Table.read(filename, format=fmt))
 
         # For specific formats, and in any case when the format is not
         # specified, make sure that complex values are treated correctly.
@@ -380,7 +646,6 @@ class StingrayObject(object):
             The file format to store the data in.
             Available options are ``pickle``, ``hdf5``, ``ascii``, ``fits``
         """
-
         if fmt is None:
             pass
         elif fmt.lower() == "pickle":
@@ -390,21 +655,694 @@ class StingrayObject(object):
         elif fmt.lower() == "ascii":
             fmt = "ascii.ecsv"
 
-        ts = self.to_astropy_table()
+        probe_file = "probe.bu.bu." + filename[-7:]
+
+        CAN_SAVE_LONGD = _can_save_longdouble(probe_file, fmt)
+        CAN_SERIALIZE_META = _can_serialize_meta(probe_file, fmt)
+
+        to_be_saved = self
+
+        ts = to_be_saved.to_astropy_table(no_longdouble=not CAN_SAVE_LONGD)
+
         if fmt is None or "ascii" in fmt:
             for col in ts.colnames:
-                if np.iscomplex(ts[col][0]):
+                if np.iscomplex(ts[col].flatten()[0]):
                     ts[f"{col}.real"] = ts[col].real
                     ts[f"{col}.imag"] = ts[col].imag
                     ts.remove_column(col)
 
-        try:
+        if CAN_SERIALIZE_META:
             ts.write(filename, format=fmt, overwrite=True, serialize_meta=True)
-        except TypeError:
+        else:
             ts.write(filename, format=fmt, overwrite=True)
 
+    def apply_mask(self, mask: npt.ArrayLike, inplace: bool = False, filtered_attrs: list = None):
+        """Apply a mask to all array attributes of the object
+
+        Parameters
+        ----------
+        mask : array of ``bool``
+            The mask. Has to be of the same length as ``self.time``
+
+        Other parameters
+        ----------------
+        inplace : bool
+            If True, overwrite the current object. Otherwise, return a new one.
+        filtered_attrs : list of str or None
+            Array attributes to be filtered. Defaults to all array attributes if ``None``.
+            The other array attributes will be set to ``None``. The main array attr is always
+            included.
+
+        Returns
+        -------
+        ts_new : StingrayObject object
+            The new object with the mask applied if ``inplace`` is ``False``, otherwise the
+            same object.
+        """
+        all_attrs = self.internal_array_attrs() + self.array_attrs()
+        if filtered_attrs is None:
+            filtered_attrs = all_attrs
+
+        if inplace:
+            new_ts = self
+        else:
+            new_ts = type(self)()
+            for attr in self.meta_attrs():
+                setattr(new_ts, attr, copy.deepcopy(getattr(self, attr)))
+
+        # If the main array attr is managed through an internal attr
+        # (e.g. lightcurve), set the internal attr instead.
+        if hasattr(self, "_" + self.main_array_attr):
+            setattr(
+                new_ts,
+                "_" + self.main_array_attr,
+                copy.deepcopy(np.asarray(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]),
+            )
+
+        for attr in all_attrs:
+            if attr not in filtered_attrs:
+                # Eliminate all unfiltered attributes
+                setattr(new_ts, attr, None)
+            else:
+                setattr(new_ts, attr, copy.deepcopy(np.asarray(getattr(self, attr))[mask]))
+        return new_ts
+
+    def _operation_with_other_obj(
+        self,
+        other,
+        operation,
+        operated_attrs=None,
+        error_attrs=None,
+        error_operation=None,
+        inplace=False,
+    ):
+        """
+        Helper method to codify an operation of one time series with another (e.g. add, subtract).
+        Takes into account the GTIs, and returns a new :class:`StingrayTimeseries` object.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` object
+            A second time series object
+
+        operation : function
+            An operation between the :class:`StingrayTimeseries` object calling this method, and
+            ``other``, operating on all the specified array attributes.
+
+        Other parameters
+        ----------------
+        operated_attrs : list of str or None
+            Array attributes to be operated on. Defaults to all array attributes not ending in
+            ``_err``.
+            The other array attributes will be discarded from the time series to avoid
+            inconsistencies.
+
+        error_attrs : list of str or None
+            Array attributes to be operated on with ``error_operation``. Defaults to all array
+            attributes ending with ``_err``.
+
+        error_operation : function
+            The function used for error propagation. Defaults to the sum of squares.
+
+        Returns
+        -------
+        ts_new : StingrayTimeseries object
+            The new time series calculated in ``operation``
+        """
+
+        if operated_attrs is None:
+            operated_attrs = self._default_operated_attrs()
+
+        if error_attrs is None:
+            error_attrs = self._default_error_attrs()
+
+        if not isinstance(other, type(self)):
+            raise TypeError(
+                f"{type(self)} objects can only be operated with other {type(self)} objects."
+            )
+
+        this_time = getattr(self, self.main_array_attr)
+        # ValueError is raised by Numpy while asserting np.equal over arrays
+        # with different dimensions.
+        try:
+            assert np.array_equal(this_time, getattr(other, self.main_array_attr))
+        except (ValueError, AssertionError):
+            raise ValueError(
+                f"The values of {self.main_array_attr} are different in the two {type(self)} "
+                "objects."
+            )
+
+        if inplace:
+            ts_new = self
+        else:
+            ts_new = type(self)()
+        setattr(ts_new, self.main_array_attr, this_time)
+        for attr in self.meta_attrs():
+            setattr(ts_new, attr, copy.deepcopy(getattr(self, attr)))
+
+        for attr in operated_attrs:
+            setattr(
+                ts_new,
+                attr,
+                operation(getattr(self, attr), getattr(other, attr)),
+            )
+
+        for attr in error_attrs:
+            setattr(
+                ts_new,
+                attr,
+                error_operation(getattr(self, attr), getattr(other, attr)),
+            )
+
+        return ts_new
+
+    def add(
+        self, other, operated_attrs=None, error_attrs=None, error_operation=sqsum, inplace=False
+    ):
+        """Add two :class:`StingrayObject` instances.
+
+        Add the array values of two :class:`StingrayObject` instances element by element, assuming
+        the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` object
+            A second time series object
+
+        Other parameters
+        ----------------
+        operated_attrs : list of str or None
+            Array attributes to be operated on. Defaults to all array attributes not ending in
+            ``_err``.
+            The other array attributes will be discarded from the time series to avoid
+            inconsistencies.
+        error_attrs : list of str or None
+            Array attributes to be operated on with ``error_operation``. Defaults to all array
+            attributes ending with ``_err``.
+        error_operation : function
+            Function to be called to propagate the errors
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
+        """
+        return self._operation_with_other_obj(
+            other,
+            np.add,
+            operated_attrs=operated_attrs,
+            error_attrs=error_attrs,
+            error_operation=error_operation,
+            inplace=inplace,
+        )
+
+    def __add__(self, other):
+        """Operation that gets called with the ``+`` operator.
+
+        Add the array values of two :class:`StingrayObject` instances element by element, assuming
+        the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+        """
+
+        return self._operation_with_other_obj(
+            other,
+            np.add,
+            error_operation=sqsum,
+        )
+
+    def __iadd__(self, other):
+        """Operation that gets called with the ``+=`` operator.
+
+        Add the array values of two :class:`StingrayObject` instances element by element, assuming
+        the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+        """
+
+        return self._operation_with_other_obj(
+            other,
+            np.add,
+            error_operation=sqsum,
+            inplace=True,
+        )
+
+    def sub(
+        self, other, operated_attrs=None, error_attrs=None, error_operation=sqsum, inplace=False
+    ):
+        """
+        Subtract *all the array attrs* of two :class:`StingrayObject` instances element by element, assuming the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` object
+            A second time series object
+
+        Other parameters
+        ----------------
+        operated_attrs : list of str or None
+            Array attributes to be operated on. Defaults to all array attributes not ending in
+            ``_err``.
+            The other array attributes will be discarded from the time series to avoid
+            inconsistencies.
+        error_attrs : list of str or None
+            Array attributes to be operated on with ``error_operation``. Defaults to all array
+            attributes ending with ``_err``.
+        error_operation : function
+            Function to be called to propagate the errors
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
+        """
+        return self._operation_with_other_obj(
+            other,
+            np.subtract,
+            operated_attrs=operated_attrs,
+            error_attrs=error_attrs,
+            error_operation=error_operation,
+            inplace=inplace,
+        )
+
+    def __sub__(self, other):
+        """Operation that gets called with the ``-`` operator.
+
+        Subtract *all the array attrs* of two :class:`StingrayObject` instances element by element, assuming the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+        """
+
+        return self._operation_with_other_obj(
+            other,
+            np.subtract,
+            error_operation=sqsum,
+        )
+
+    def __isub__(self, other):
+        """Operation that gets called with the ``-=`` operator.
+
+        Subtract *all the array attrs* of two :class:`StingrayObject` instances element by element, assuming the main array attributes of the instances match exactly.
+
+        All array attrs ending with ``_err`` are treated as error bars and propagated with the
+        sum of squares.
+
+        GTIs are crossed, so that only common intervals are saved.
+        """
+
+        return self._operation_with_other_obj(
+            other,
+            np.subtract,
+            error_operation=sqsum,
+            inplace=True,
+        )
+
+    def __neg__(self):
+        """
+        Implement the behavior of negation of the array attributes of a :class:`StingrayObject`
+        Error attrs are left alone.
+
+        The negation operator ``-`` is supposed to invert the sign of all array attributes of a
+        time series object, leaving out the ones ending with ``_err``.
+
+        """
+
+        ts_new = copy.deepcopy(self)
+        for attr in self._default_operated_attrs():
+            setattr(ts_new, attr, -np.asarray(getattr(self, attr)))
+
+        return ts_new
+
+    def __len__(self):
+        """
+        Return the number of bins of a the main array attributes
+
+        This method overrides the ``len`` function for a :class:`StingrayObject`
+        object and returns the length of the array attributes (using the main array attribute
+        as probe).
+        """
+        return np.size(getattr(self, self.main_array_attr))
+
+    def __getitem__(self, index):
+        """
+        Return an element or a slice of the :class:`StingrayObject`.
+
+        Parameters
+        ----------
+        index : int or slice instance
+            Index value of the time array or a slice object.
+
+        Returns
+        -------
+        ts_new : :class:`StingrayObject` object
+            The new :class:`StingrayObject` object with the set of selected data.
+        """
+        from .utils import assign_value_if_none
+
+        if isinstance(index, (int, np.integer)):
+            start = index
+            stop = index + 1
+            step = 1
+        elif isinstance(index, slice):
+            start = assign_value_if_none(index.start, 0)
+            stop = assign_value_if_none(index.stop, len(self))
+            step = assign_value_if_none(index.step, 1)
+        else:
+            raise IndexError("The index must be either an integer or a slice " "object !")
+
+        new_ts = type(self)()
+        for attr in self.meta_attrs():
+            setattr(new_ts, attr, copy.deepcopy(getattr(self, attr)))
+
+        for attr in self.array_attrs() + [self.main_array_attr]:
+            setattr(new_ts, attr, getattr(self, attr)[start:stop:step])
+
+        return new_ts
+
 
 class StingrayTimeseries(StingrayObject):
+    """Basic class for time series data.
+
+    This can be events, binned light curves, unevenly sampled light curves, etc. The only
+    requirement is that the data (which can be any quantity, related or not to an electromagnetic
+    measurement) are associated with a time measurement.
+    We make a distinction between the *array* attributes, which have the same length of the
+    ``time`` array, and the *meta* attributes, which can be scalars or arrays of different
+    size. The array attributes can be multidimensional (e.g. a spectrum for each time bin),
+    but their first dimension (``array.shape[0]``) must have same length of the ``time`` array.
+
+    Array attributes are singled out automatically depending on their shape. All filtering
+    operations (e.g. ``apply_gtis``, ``rebin``, etc.) are applied to array attributes only.
+    For this reason, it is advisable to specify whether a given attribute should *not* be
+    considered as an array attribute by adding it to the ``not_array_attr`` list.
+
+    Parameters
+    ----------
+    time: iterable
+        A list or array of time stamps
+
+    Other Parameters
+    ----------------
+    array_attrs : dict
+        Array attributes to be set (e.g. ``{"flux": flux_array, "flux_err": flux_err_array}``).
+        In principle, they could be specified as simple keyword arguments. But this way, we
+        will run a check on the length of the arrays, and raise an error if they are not of a
+        shape compatible with the ``time`` array.
+
+    dt: float
+        The time resolution of the time series. Can be a scalar or an array attribute (useful
+        for non-uniformly sampled data or events from different instruments)
+
+    mjdref : float
+        The MJD used as a reference for the time array.
+
+    gtis: ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
+        Good Time Intervals
+
+    high_precision : bool
+        Change the precision of self.time to float128. Useful while dealing with fast pulsars.
+
+    timeref : str
+        The time reference, as recorded in the FITS file (e.g. SOLARSYSTEM)
+
+    timesys : str
+        The time system, as recorded in the FITS file (e.g. TDB)
+
+    ephem : str
+        The JPL ephemeris used to barycenter the data, if any (e.g. DE430)
+
+    **other_kw :
+        Used internally. Any other keyword arguments will be set as attributes of the object.
+
+    Attributes
+    ----------
+    time: numpy.ndarray
+        The array of time stamps, in seconds from the reference
+        MJD defined in ``mjdref``
+
+    not_array_attr: list
+        List of attributes that are never to be considered as array attributes. For example, GTIs
+        are not array attributes.
+
+    dt: float
+        The time resolution of the measurements. Can be a scalar or an array attribute (useful
+        for non-uniformly sampled data or events from different instruments)
+
+    mjdref : float
+        The MJD used as a reference for the time array.
+
+    gtis: ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
+        Good Time Intervals
+
+    high_precision : bool
+        Change the precision of self.time to float128. Useful while dealing with fast pulsars.
+
+    """
+
+    main_array_attr: str = "time"
+    not_array_attr: list = ["gti"]
+    _time: TTime = None
+    high_precision: bool = False
+    mjdref: TTime = 0.0
+    dt: float = 0.0
+
+    def __init__(
+        self,
+        time: TTime = None,
+        array_attrs: dict = {},
+        mjdref: TTime = 0,
+        notes: str = "",
+        gti: npt.ArrayLike = None,
+        high_precision: bool = False,
+        ephem: str = None,
+        timeref: str = None,
+        timesys: str = None,
+        **other_kw,
+    ):
+        StingrayObject.__init__(self)
+
+        self.notes = notes
+        self.mjdref = mjdref
+        self.gti = gti
+        self.ephem = ephem
+        self.timeref = timeref
+        self.timesys = timesys
+        self._mask = None
+        self.high_precision = high_precision
+        self.dt = other_kw.pop("dt", 0)
+
+        self._set_times(time, high_precision=high_precision)
+        for kw in other_kw:
+            setattr(self, kw, other_kw[kw])
+        for kw in array_attrs:
+            new_arr = np.asarray(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)
+
+    @property
+    def time(self):
+        return self._time
+
+    @time.setter
+    def time(self, value):
+        value = self._validate_and_format(value, "time", "time")
+        if value is None:
+            for attr in self.internal_array_attrs() + self.array_attrs():
+                setattr(self, attr, None)
+        self._set_times(value, high_precision=self.high_precision)
+
+    @property
+    def gti(self):
+        if self._gti is None and self._time is not None:
+            if isinstance(self.dt, Iterable):
+                dt0 = self.dt[0]
+                dt1 = self.dt[-1]
+            else:
+                dt0 = dt1 = self.dt
+            self._gti = np.asarray([[self._time[0] - dt0 / 2, self._time[-1] + dt1 / 2]])
+        return self._gti
+
+    @gti.setter
+    def gti(self, value):
+        if value is None:
+            self._gti = None
+            return
+        value = np.asarray(value)
+        self._gti = value
+        self._mask = None
+
+    @property
+    def mask(self):
+        from .gti import create_gti_mask
+
+        if self._mask is None:
+            self._mask = create_gti_mask(self.time, self.gti, dt=self.dt)
+        return self._mask
+
+    @property
+    def n(self):
+        return self.main_array_length
+
+    def _set_times(self, time, high_precision=False):
+        if time is None or np.size(time) == 0:
+            self._time = None
+            return
+        time, _ = interpret_times(time, self.mjdref)
+        if not high_precision:
+            self._time = np.asarray(time)
+        else:
+            self._time = np.asarray(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"),
+            attrs_to_discard=["_mask", "header"],
+        )
+
+    def _validate_and_format(self, value, attr_name, compare_to_attr):
+        """Check if the size of a value is compatible with the size of another attribute.
+
+        Different cases are possible:
+
+        - If the value is None, we return None
+        - If the value is a scalar, we fail
+        - If the value is an array, we check if it has the correct shape by comparing it with
+          another attribute. In the special case where the attribute is the same, if it is None
+          we assign the new value. Otherwise, the first dimension of the value and the current
+          value of the attribute being compared with has to be the same.
+
+        Parameters
+        ----------
+        value : array-like or None
+            The value to check.
+        attr_name : str
+            The name of the attribute being checked.
+        compare_to_attr : str
+            The name of the attribute to compare with.
+
+        Returns
+        -------
+        value : array-like or None
+            The value to check wrapped in a class:`np.array`, if it is not None. Otherwise None
+        """
+        if value is None:
+            return None
+        value = np.asarray(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
+        # This can happen, e.g. with the time array.
+        compare_with = getattr(self, compare_to_attr, None)
+        if attr_name == compare_to_attr and compare_with is None:
+            return value
+
+        # In the special case where the current value of the attribute being compared
+        # is None, this also has to fail.
+        if compare_with is None:
+            raise ValueError(
+                f"Can only assign new {attr_name} if the {compare_to_attr} array is not None"
+            )
+        if value.shape[0] != compare_with.shape[0]:
+            raise ValueError(
+                f"Can only assign new {attr_name} of the same shape as the {compare_to_attr} array"
+            )
+        return value
+
+    def __eq__(self, other_ts):
+        return super().__eq__(other_ts)
+
+    def apply_gtis(self, new_gti=None, inplace: bool = True):
+        """
+        Apply Good Time Intervals (GTIs) to a time series. Filters all the array attributes, only
+        keeping the bins that fall into GTIs.
+
+        Parameters
+        ----------
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
+
+        """
+        # I import here to avoid the risk of circular imports
+        from .gti import check_gtis, create_gti_mask
+
+        if new_gti is None:
+            new_gti = self.gti
+
+        check_gtis(new_gti)
+
+        # This will automatically be recreated from GTIs once I set it to None
+        good = create_gti_mask(self.time, new_gti, dt=self.dt)
+        newts = self.apply_mask(good, inplace=inplace)
+        # Important, otherwise addition/subtraction ops will go into an infinite loop
+        if inplace:
+            newts.gti = new_gti
+        return newts
+
+    def split_by_gti(self, gti=None, min_points=2):
+        """
+        Split the current :class:`StingrayTimeseries` object into a list of
+        :class:`StingrayTimeseries` objects, one for each continuous GTI segment
+        as defined in the ``gti`` attribute.
+
+        Parameters
+        ----------
+        min_points : int, default 1
+            The minimum number of data points in each time series. Light
+            curves with fewer data points will be ignored.
+
+        Returns
+        -------
+        list_of_tss : list
+            A list of :class:`StingrayTimeseries` objects, one for each GTI segment
+        """
+        from .gti import gti_border_bins, create_gti_mask
+
+        if gti is None:
+            gti = self.gti
+
+        list_of_tss = []
+
+        start_bins, stop_bins = gti_border_bins(gti, self.time, self.dt)
+        for i in range(len(start_bins)):
+            start = start_bins[i]
+            stop = stop_bins[i]
+
+            if (stop - start) < min_points:
+                continue
+
+            new_gti = np.array([gti[i]])
+            mask = create_gti_mask(self.time, new_gti)
+
+            # Note: GTIs are consistent with default in this case!
+            new_ts = self.apply_mask(mask)
+            new_ts.gti = new_gti
+
+            list_of_tss.append(new_ts)
+
+        return list_of_tss
+
     def to_astropy_timeseries(self) -> TimeSeries:
         """Save the ``StingrayTimeseries`` to an ``Astropy`` timeseries.
 
@@ -449,11 +1387,8 @@ class StingrayTimeseries(StingrayObject):
 
         The timeseries has to define at least a column called time,
         the rest of columns will form the array attributes of the
-        new event list, while the attributes in table.meta will
-        form the new meta attributes of the event list.
-
-        It is strongly advisable to define such attributes and columns
-        using the standard attributes of EventList: time, pi, energy, gti etc.
+        new time series, while the attributes in table.meta will
+        form the new meta attributes of the time series.
 
         Parameters
         ----------
@@ -472,21 +1407,22 @@ class StingrayTimeseries(StingrayObject):
         if "mjdref" in ts.meta:
             mjdref = ts.meta["mjdref"]
 
+        new_cls = cls()
         time, mjdref = interpret_times(time, mjdref)
-        cls.time = np.asarray(time)  # type: ignore
+        new_cls.time = np.asarray(time)  # type: ignore
 
         array_attrs = ts.colnames
         for key, val in ts.meta.items():
-            setattr(cls, key, val)
+            setattr(new_cls, key, val)
 
         for attr in array_attrs:
             if attr == "time":
                 continue
-            setattr(cls, attr, ts[attr])
+            setattr(new_cls, attr, np.asarray(ts[attr]))
 
-        return cls
+        return new_cls
 
-    def change_mjdref(self, new_mjdref: float) -> StingrayTimeseries:
+    def change_mjdref(self, new_mjdref: float, inplace=False) -> StingrayTimeseries:
         """Change the MJD reference time (MJDREF) of the time series
 
         The times of the time series will be shifted in order to be referred to
@@ -497,25 +1433,35 @@ class StingrayTimeseries(StingrayObject):
         new_mjdref : float
             New MJDREF
 
+        Other parameters
+        ----------------
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
+
         Returns
         -------
-        new_lc : :class:`StingrayTimeseries` object
+        new_ts : :class:`StingrayTimeseries` object
             The new time series, shifted by MJDREF
         """
         time_shift = (self.mjdref - new_mjdref) * 86400  # type: ignore
 
-        ts = self.shift(time_shift)
+        ts = self.shift(time_shift, inplace=inplace)
         ts.mjdref = new_mjdref  # type: ignore
         return ts
 
-    def shift(self, time_shift: float) -> StingrayTimeseries:
+    def shift(self, time_shift: float, inplace=False) -> StingrayTimeseries:
         """Shift the time and the GTIs by the same amount
 
         Parameters
         ----------
         time_shift: float
-            The time interval by which the light curve will be shifted (in
-            the same units as the time array in :class:`Lightcurve`
+            The time interval by which the time series will be shifted (in
+            the same units as the time array in :class:`StingrayTimeseries`
+
+        Other parameters
+        ----------------
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
 
         Returns
         -------
@@ -523,13 +1469,761 @@ class StingrayTimeseries(StingrayObject):
             The new time series shifted by ``time_shift``
 
         """
-        ts = copy.deepcopy(self)
+        if inplace:
+            ts = self
+        else:
+            ts = copy.deepcopy(self)
         ts.time = np.asarray(ts.time) + time_shift  # type: ignore
-        if hasattr(ts, "gti"):
-            ts.gti = np.asarray(ts.gti) + 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
 
         return ts
 
+    def _operation_with_other_obj(
+        self, other, operation, operated_attrs=None, error_attrs=None, error_operation=None
+    ):
+        """
+        Helper method to codify an operation of one time series with another (e.g. add, subtract).
+        Takes into account the GTIs correctly, and returns a new :class:`StingrayTimeseries` object.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` object
+            A second time series object
+
+        operation : function
+            An operation between the :class:`StingrayTimeseries` object calling this method, and
+            ``other``, operating on all the specified array attributes.
+
+        Other parameters
+        ----------------
+        operated_attrs : list of str or None
+            Array attributes to be operated on. Defaults to all array attributes not ending in
+            ``_err``.
+            The other array attributes will be discarded from the time series to avoid
+            inconsistencies.
+
+        error_attrs : list of str or None
+            Array attributes to be operated on with ``error_operation``. Defaults to all array
+            attributes ending with ``_err``.
+
+        error_operation : function
+            The function used for error propagation. Defaults to the sum of squares.
+
+        Returns
+        -------
+        ts_new : StingrayTimeseries object
+            The new time series calculated in ``operation``
+        """
+
+        if self.mjdref != other.mjdref:
+            warnings.warn("MJDref is different in the two time series")
+            other = other.change_mjdref(self.mjdref)
+
+        if not np.array_equal(self.gti, other.gti):
+            from .gti import cross_two_gtis
+
+            warnings.warn(
+                "The good time intervals in the two time series are different. Data outside the "
+                "common GTIs will be discarded."
+            )
+            common_gti = cross_two_gtis(self.gti, other.gti)
+            masked_self = self.apply_gtis(common_gti)
+            masked_other = other.apply_gtis(common_gti)
+            return masked_self._operation_with_other_obj(
+                masked_other,
+                operation,
+                operated_attrs=operated_attrs,
+                error_attrs=error_attrs,
+                error_operation=error_operation,
+            )
+
+        return super()._operation_with_other_obj(
+            other,
+            operation,
+            operated_attrs=operated_attrs,
+            error_attrs=error_attrs,
+            error_operation=error_operation,
+        )
+
+    def __add__(self, other):
+        """
+        Add the array values of two time series element by element, assuming they
+        have the same time array.
+
+        This magic method adds two :class:`TimeSeries` objects having the same time
+        array such that the corresponding array arrays get summed up.
+
+        GTIs are crossed, so that only common intervals are saved.
+
+        Examples
+        --------
+        >>> time = [5, 10, 15]
+        >>> count1 = [300, 100, 400]
+        >>> count2 = [600, 1200, 800]
+        >>> gti1 = [[0, 25]]
+        >>> gti2 = [[0, 25]]
+        >>> ts1 = StingrayTimeseries(time, array_attrs=dict(counts=count1), gti=gti1, dt=5)
+        >>> ts2 = StingrayTimeseries(time, array_attrs=dict(counts=count2), gti=gti2, dt=5)
+        >>> ts = ts1 + ts2
+        >>> np.allclose(ts.counts, [ 900, 1300, 1200])
+        True
+        """
+
+        return super().__add__(other)
+
+    def __sub__(self, other):
+        """
+        Subtract the counts/flux of one time series from the counts/flux of another
+        time series element by element, assuming the ``time`` arrays of the time series
+        match exactly.
+
+        This magic method adds two :class:`StingrayTimeSeries` objects having the same
+        ``time`` array and subtracts the ``counts`` of one :class:`StingrayTimeseries` with
+        that of another, while also updating ``countrate``, ``counts_err`` and ``countrate_err``
+        correctly.
+
+        GTIs are crossed, so that only common intervals are saved.
+
+        Examples
+        --------
+        >>> time = [10, 20, 30]
+        >>> count1 = [600, 1200, 800]
+        >>> count2 = [300, 100, 400]
+        >>> gti1 = [[0, 35]]
+        >>> gti2 = [[0, 35]]
+        >>> ts1 = StingrayTimeseries(time, array_attrs=dict(counts=count1), gti=gti1, dt=10)
+        >>> ts2 = StingrayTimeseries(time, array_attrs=dict(counts=count2), gti=gti2, dt=10)
+        >>> ts = ts1 - ts2
+        >>> np.allclose(ts.counts, [ 300, 1100,  400])
+        True
+        """
+
+        return super().__sub__(other)
+
+    def __getitem__(self, index):
+        """
+        Return the corresponding count value at the index or a new :class:`StingrayTimeseries`
+        object upon slicing.
+
+        This method adds functionality to retrieve the count value at
+        a particular index. This also can be used for slicing and generating
+        a new :class:`StingrayTimeseries` object. GTIs are recalculated based on the new light
+        curve segment
+
+        If the slice object is of kind ``start:stop:step`` and ``dt`` is not 0, GTIs are also
+        sliced, by crossing with ``zip(time - self.dt /2, time + self.dt / 2)``
+
+        Parameters
+        ----------
+        index : int or slice instance
+            Index value of the time array or a slice object.
+
+        Examples
+        --------
+        >>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        >>> count = [11, 22, 33, 44, 55, 66, 77, 88, 99]
+        >>> ts = StingrayTimeseries(time, array_attrs=dict(counts=count), dt=1)
+        >>> np.allclose(ts[2].counts, [33])
+        True
+        >>> np.allclose(ts[:2].counts, [11, 22])
+        True
+        """
+        from .utils import assign_value_if_none
+        from .gti import cross_two_gtis
+
+        new_ts = super().__getitem__(index)
+        step = 1
+        if isinstance(index, slice):
+            step = assign_value_if_none(index.step, 1)
+
+        dt = self.dt
+        if np.isscalar(dt):
+            delta_gti_start = delta_gti_stop = dt * 0.5
+        else:
+            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]])
+        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)
+        new_gti = cross_two_gtis(self.gti, new_gti)
+
+        new_ts.gti = new_gti
+        return new_ts
+
+    def truncate(self, start=0, stop=None, method="index"):
+        """
+        Truncate a :class:`StingrayTimeseries` object.
+
+        This method takes a ``start`` and a ``stop`` point (either as indices,
+        or as times in the same unit as those in the ``time`` attribute, and truncates
+        all bins before ``start`` and after ``stop``, then returns a new
+        :class:`StingrayTimeseries` object with the truncated time series.
+
+        Parameters
+        ----------
+        start : int, default 0
+            Index (or time stamp) of the starting point of the truncation. If no value is set
+            for the start point, then all points from the first element in the ``time`` array
+            are taken into account.
+
+        stop : int, default ``None``
+            Index (or time stamp) of the ending point (exclusive) of the truncation. If no
+            value of stop is set, then points including the last point in
+            the counts array are taken in count.
+
+        method : {``index`` | ``time``}, optional, default ``index``
+            Type of the start and stop values. If set to ``index`` then
+            the values are treated as indices of the counts array, or
+            if set to ``time``, the values are treated as actual time values.
+
+        Returns
+        -------
+        ts_new: :class:`StingrayTimeseries` object
+            The :class:`StingrayTimeseries` object with truncated time and arrays.
+
+        Examples
+        --------
+        >>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        >>> count = [10, 20, 30, 40, 50, 60, 70, 80, 90]
+        >>> ts = StingrayTimeseries(time, array_attrs={"counts": count}, dt=1)
+        >>> ts_new = ts.truncate(start=2, stop=8)
+        >>> np.allclose(ts_new.counts, [30, 40, 50, 60, 70, 80])
+        True
+        >>> ts_new.time
+        array([3, 4, 5, 6, 7, 8])
+        >>> # Truncation can also be done by time values
+        >>> ts_new = ts.truncate(start=6, method='time')
+        >>> ts_new.time
+        array([6, 7, 8, 9])
+        >>> np.allclose(ts_new.counts, [60, 70, 80, 90])
+        True
+        """
+
+        if not isinstance(method, str):
+            raise TypeError("The method keyword argument is not a string !")
+
+        if method.lower() not in ["index", "time"]:
+            raise ValueError("Unknown method type " + method + ".")
+
+        if method.lower() == "index":
+            new_ts = self._truncate_by_index(start, stop)
+        else:
+            new_ts = self._truncate_by_time(start, stop)
+        new_ts.tstart = new_ts.gti[0, 0]
+        new_ts.tseg = new_ts.gti[-1, 1] - new_ts.gti[0, 0]
+        return new_ts
+
+    def _truncate_by_index(self, start, stop):
+        """Private method for truncation using index values."""
+        from .gti import cross_two_gtis
+
+        new_ts = self.apply_mask(slice(start, stop))
+
+        dtstart = dtstop = new_ts.dt
+        if isinstance(self.dt, Iterable):
+            dtstart = self.dt[0]
+            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]])
+        )
+
+        new_ts.gti = gti
+
+        return new_ts
+
+    def _truncate_by_time(self, start, stop):
+        """Helper method for truncation using time values.
+
+        Parameters
+        ----------
+        start : float
+            start time for new time series; all time bins before this time will be discarded
+
+        stop : float
+            stop time for new time series; all time bins after this point will be discarded
+
+        Returns
+        -------
+            new_ts : StingrayTimeseries
+                A new :class:`StingrayTimeseries` object with the truncated time bins
+
+        """
+
+        if stop is not None:
+            if start > stop:
+                raise ValueError("start time must be less than stop time!")
+
+        if not start == 0:
+            start = self.time.searchsorted(start)
+
+        if stop is not None:
+            stop = self.time.searchsorted(stop)
+
+        return self._truncate_by_index(start, stop)
+
+    def concatenate(self, other, check_gti=True):
+        """
+        Concatenate two :class:`StingrayTimeseries` objects.
+
+        This method concatenates two or more :class:`StingrayTimeseries` objects along the time
+        axis. GTIs are recalculated by merging all the GTIs together. GTIs should not overlap at
+        any point.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` object or list of :class:`StingrayTimeseries` objects
+            A second time series object, or a list of objects to be concatenated
+
+        Other parameters
+        ----------------
+        check_gti : bool
+            Check if the GTIs are overlapping or not. Default: True
+            If this is True and GTIs overlap, an error is raised.
+        """
+
+        if check_gti:
+            treatment = "append"
+        else:
+            treatment = "none"
+        new_ts = self._join_timeseries(other, strategy=treatment)
+        return new_ts
+
+    def _join_timeseries(self, others, strategy="intersection", ignore_meta=[]):
+        """Helper method to join two or more :class:`StingrayTimeseries` objects.
+
+        This is a helper method that can be called by other user-facing methods, such as
+        :class:`StingrayTimeseries().join()`.
+
+        Standard attributes such as ``pi`` and ``energy`` remain ``None`` if they are ``None``
+        in both. Otherwise, ``np.nan`` is used as a default value for the missing values.
+        Arbitrary array attributes are created and joined using the same convention.
+
+        Multiple checks are done on the joined time series. If the time array of the series
+        being joined is empty, it is ignored (and a copy of the original time series is returned
+        instead). If the time resolution is different, the final time series will associate
+        different time resolutions to different time bins.
+        If the MJDREF is different (including being 0), the time reference will be changed to
+        the one of the first time series. An empty time series will be ignored.
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` or class:`list` of :class:`StingrayTimeseries`
+            The other :class:`StingrayTimeseries` object which is supposed to be joined with.
+            If ``other`` is a list, it is assumed to be a list of :class:`StingrayTimeseries`
+            and they are all joined, one by one.
+
+        Other parameters
+        ----------------
+        strategy : {"intersection", "union", "append", "infer", "none"}
+            Method to use to merge the GTIs. If "intersection", the GTIs are merged
+            using the intersection of the GTIs. If "union", the GTIs are merged
+            using the union of the GTIs. If "none", a single GTI with the minimum and
+            the maximum time stamps of all GTIs is returned. If "infer", the strategy
+            is decided based on the GTIs. If there are no overlaps, "union" is used,
+            otherwise "intersection" is used. If "append", the GTIs are simply appended
+            but they must be mutually exclusive.
+
+        Returns
+        -------
+        `ts_new` : :class:`StingrayTimeseries` object
+            The resulting :class:`StingrayTimeseries` object.
+        """
+        from .gti import check_separate, cross_gtis, append_gtis
+
+        new_ts = type(self)()
+
+        if not (
+            isinstance(others, Iterable)
+            and not isinstance(others, str)
+            and not isinstance(others, StingrayObject)
+        ):
+            others = [others]
+        else:
+            others = others
+
+        # First of all, check if there are empty objects
+        for obj in others:
+            if not isinstance(obj, type(self)):
+                raise TypeError(
+                    f"{type(self)} objects can only be merged with other {type(self)} objects."
+                )
+            if getattr(obj, "time", None) is None or np.size(obj.time) == 0:
+                warnings.warn("One of the time series you are joining is empty.")
+                others.remove(obj)
+
+        if len(others) == 0:
+            return copy.deepcopy(self)
+
+        for i, other in enumerate(others):
+            # Tolerance for MJDREF:1 microsecond
+            if not np.isclose(self.mjdref, other.mjdref, atol=1e-6 / 86400):
+                warnings.warn("Attribute mjdref is different in the time series being merged.")
+                others[i] = other.change_mjdref(self.mjdref)
+
+        all_objs = [self] + others
+
+        from .gti import merge_gtis
+
+        # Check if none of the GTIs was already initialized.
+        all_gti = [obj._gti for obj in all_objs if obj._gti is not None]
+
+        if len(all_gti) == 0 or strategy == "none":
+            new_gti = None
+        else:
+            # For this, initialize the GTIs
+            new_gti = merge_gtis([obj.gti for obj in all_objs], strategy=strategy)
+
+        all_time_arrays = [obj.time for obj in all_objs if obj.time is not None]
+
+        new_ts.time = np.concatenate(all_time_arrays)
+        order = np.argsort(new_ts.time)
+        new_ts.time = new_ts.time[order]
+
+        new_ts.gti = new_gti
+
+        dts = list(set([getattr(obj, "dt", None) for obj in all_objs]))
+        if len(dts) != 1:
+            warnings.warn("The time resolution is different. Transforming in array")
+
+            new_dt = np.concatenate([np.zeros_like(obj.time) + obj.dt for obj in all_objs])
+            new_ts.dt = new_dt[order]
+        else:
+            new_ts.dt = dts[0]
+
+        def _get_set_from_many_lists(lists):
+            """Make a single set out of many lists."""
+            all_vals = []
+            for ls in lists:
+                all_vals += ls
+            return set(all_vals)
+
+        def _get_all_array_attrs(objs):
+            """Get all array attributes from the time series being merged. Do not include time."""
+            return _get_set_from_many_lists(
+                [obj.array_attrs() + obj.internal_array_attrs() for obj in objs]
+            )
+
+        for attr in _get_all_array_attrs(all_objs):
+            # if it's here, it means that it's an array attr in at least one object.
+            # So, everywhere it's None, it needs to be set to 0s of the same length as time
+            new_attr_values = []
+            for obj in all_objs:
+                if getattr(obj, attr, None) is None:
+                    warnings.warn(
+                        f"The {attr} array is empty in one of the time series being merged. "
+                        "Setting it to NaN for the affected events"
+                    )
+                    new_attr_values.append(np.zeros_like(obj.time) + np.nan)
+                else:
+                    new_attr_values.append(getattr(obj, attr))
+
+            new_attr = np.concatenate(new_attr_values)[order]
+            setattr(new_ts, attr, new_attr)
+
+        all_meta_attrs = _get_set_from_many_lists([obj.meta_attrs() for obj in all_objs])
+        # The attributes being treated separately are removed from the standard treatment
+        # When energy, pi etc. are None, they might appear in the meta_attrs, so we
+        # also add them to the list of attributes to be removed if present.
+        to_remove = ["gti", "dt"] + new_ts.array_attrs()
+        for attr in to_remove:
+            if attr in all_meta_attrs:
+                all_meta_attrs.remove(attr)
+
+        for attr in ignore_meta:
+            logging.info(f"The {attr} attribute will be removed from the output ")
+            if attr in all_meta_attrs:
+                all_meta_attrs.remove(attr)
+
+        def _safe_concatenate(a, b):
+            if isinstance(a, str) and isinstance(b, str):
+                return a + "," + b
+            else:
+                if isinstance(a, tuple):
+                    return a + (b,)
+                return (a, b)
+
+        for attr in all_meta_attrs:
+            self_attr = getattr(self, attr, None)
+            new_val = self_attr
+            for other in others:
+                other_attr = getattr(other, attr, None)
+                if self_attr != other_attr:
+                    warnings.warn(
+                        "Attribute " + attr + " is different in the time series being merged."
+                    )
+                    new_val = _safe_concatenate(new_val, other_attr)
+            setattr(new_ts, attr, new_val)
+
+        new_ts.mjdref = self.mjdref
+
+        return new_ts
+
+    def join(self, *args, **kwargs):
+        """
+        Join other :class:`StingrayTimeseries` objects with the current one.
+
+        If both are empty, an empty :class:`StingrayTimeseries` is returned.
+
+        Standard attributes such as ``pi`` and ``energy`` remain ``None`` if they are ``None``
+        in both. Otherwise, ``np.nan`` is used as a default value for the missing values.
+        Arbitrary array attributes are created and joined using the same convention.
+
+        Multiple checks are done on the joined time series. If the time array of the series
+        being joined is empty, it is ignored. If the time resolution is different, the final
+        time series will have the rougher time resolution. If the MJDREF is different, the time
+        reference will be changed to the one of the first time series. An empty time series will
+        be ignored.
+
+        Note: ``join`` is not equivalent to ``concatenate``. ``concatenate`` is used to join
+        multiple **non-overlapping** time series along the time axis, while ``join`` is more
+        general, and can be used to join multiple time series with different strategies (see
+        parameter ``strategy`` below).
+
+        Parameters
+        ----------
+        other : :class:`StingrayTimeseries` or class:`list` of :class:`StingrayTimeseries`
+            The other :class:`StingrayTimeseries` object which is supposed to be joined with.
+            If ``other`` is a list, it is assumed to be a list of :class:`StingrayTimeseries`
+            and they are all joined, one by one.
+
+        Other parameters
+        ----------------
+        strategy : {"intersection", "union", "append", "infer", "none"}
+            Method to use to merge the GTIs. If "intersection", the GTIs are merged
+            using the intersection of the GTIs. If "union", the GTIs are merged
+            using the union of the GTIs. If "none", a single GTI with the minimum and
+            the maximum time stamps of all GTIs is returned. If "infer", the strategy
+            is decided based on the GTIs. If there are no overlaps, "union" is used,
+            otherwise "intersection" is used. If "append", the GTIs are simply appended
+            but they must be mutually exclusive.
+
+        Returns
+        -------
+        `ts_new` : :class:`StingrayTimeseries` object
+            The resulting :class:`StingrayTimeseries` object.
+        """
+        return self._join_timeseries(*args, **kwargs)
+
+    def rebin(self, dt_new=None, f=None, method="sum"):
+        """
+        Rebin the time series to a new time resolution. While the new
+        resolution need not be an integer multiple of the previous time
+        resolution, be aware that if it is not, the last bin will be cut
+        off by the fraction left over by the integer division.
+
+        Parameters
+        ----------
+        dt_new: float
+            The new time resolution of the time series. Must be larger than
+            the time resolution of the old time series!
+
+        method: {``sum`` | ``mean`` | ``average``}, optional, default ``sum``
+            This keyword argument sets whether the counts in the new bins
+            should be summed or averaged.
+
+        Other Parameters
+        ----------------
+        f: float
+            the rebin factor. If specified, it substitutes ``dt_new`` with
+            ``f*self.dt``
+
+        Returns
+        -------
+        ts_new: :class:`StingrayTimeseries` object
+            The :class:`StingrayTimeseries` object with the new, binned time series.
+        """
+        from .utils import rebin_data
+
+        if f is None and dt_new is None:
+            raise ValueError("You need to specify at least one between f and " "dt_new")
+        elif f is not None:
+            dt_new = f * self.dt
+
+        if np.any(dt_new < np.asarray(self.dt)):
+            raise ValueError("The new time resolution must be larger than the old one!")
+
+        gti_new = []
+
+        new_ts = type(self)()
+
+        for attr in self.array_attrs() + self.internal_array_attrs():
+            bin_time, bin_counts, bin_err = [], [], []
+            if attr.endswith("_err"):
+                continue
+            e_temp = None
+            for g in self.gti:
+                if g[1] - g[0] < dt_new:
+                    continue
+                else:
+                    # find start and end of GTI segment in data
+                    start_ind = self.time.searchsorted(g[0])
+                    end_ind = self.time.searchsorted(g[1])
+
+                    t_temp = self.time[start_ind:end_ind]
+                    c_temp = getattr(self, attr)[start_ind:end_ind]
+
+                    if hasattr(self, attr + "_err"):
+                        e_temp = getattr(self, attr + "_err")[start_ind:end_ind]
+
+                    bin_t, bin_c, bin_e, _ = rebin_data(
+                        t_temp, c_temp, dt_new, yerr=e_temp, method=method, dx=self.dt
+                    )
+
+                    bin_time.extend(bin_t)
+                    bin_counts.extend(bin_c)
+                    bin_err.extend(bin_e)
+                    gti_new.append(g)
+            if new_ts.time is None:
+                new_ts.time = np.array(bin_time)
+            setattr(new_ts, attr, bin_counts)
+            if e_temp is not None:
+                setattr(new_ts, attr + "_err", bin_err)
+
+        if len(gti_new) == 0:
+            raise ValueError("No valid GTIs after rebin.")
+        new_ts.gti = np.asarray(gti_new)
+
+        for attr in self.meta_attrs():
+            if attr == "dt":
+                continue
+            setattr(new_ts, attr, copy.deepcopy(getattr(self, attr)))
+        new_ts.dt = dt_new
+        return new_ts
+
+    def sort(self, reverse=False, inplace=False):
+        """
+        Sort a ``StingrayTimeseries`` object by time.
+
+        A ``StingrayTimeseries`` can be sorted in either increasing or decreasing order
+        using this method. The time array gets sorted and the counts array is
+        changed accordingly.
+
+        Parameters
+        ----------
+        reverse : boolean, default False
+            If True then the object is sorted in reverse order.
+        inplace : bool
+            If True, overwrite the current time series. Otherwise, return a new one.
+
+        Examples
+        --------
+        >>> time = [2, 1, 3]
+        >>> count = [200, 100, 300]
+        >>> ts = StingrayTimeseries(time, array_attrs={"counts": count}, dt=1)
+        >>> ts_new = ts.sort()
+        >>> ts_new.time
+        array([1, 2, 3])
+        >>> np.allclose(ts_new.counts, [100, 200, 300])
+        True
+
+        Returns
+        -------
+        ts_new: :class:`StingrayTimeseries` object
+            The :class:`StingrayTimeseries` object with sorted time and counts
+            arrays.
+        """
+
+        mask = np.argsort(self.time)
+        if reverse:
+            mask = mask[::-1]
+        return self.apply_mask(mask, inplace=inplace)
+
+    def plot(
+        self,
+        attr,
+        witherrors=False,
+        labels=None,
+        ax=None,
+        title=None,
+        marker="-",
+        save=False,
+        filename=None,
+        plot_btis=True,
+    ):
+        """
+        Plot the time series using ``matplotlib``.
+
+        Plot the time series object on a graph ``self.time`` on x-axis and
+        ``self.counts`` on y-axis with ``self.counts_err`` optionally
+        as error bars.
+
+        Parameters
+        ----------
+        attr: str
+            Attribute to plot.
+
+        Other parameters
+        ----------------
+        witherrors: boolean, default False
+            Whether to plot the StingrayTimeseries with errorbars or not
+        labels : iterable, default ``None``
+            A list or tuple with ``xlabel`` and ``ylabel`` as strings. E.g.
+            if the attribute is ``'counts'``, the list of labels
+            could be ``['Time (s)', 'Counts (s^-1)']``
+        ax : ``matplotlib.pyplot.axis`` object
+            Axis to be used for plotting. Defaults to creating a new one.
+        title : str, default ``None``
+            The title of the plot.
+        marker : str, default '-'
+            Line style and color of the plot. Line styles and colors are
+            combined in a single format string, as in ``'bo'`` for blue
+            circles. See ``matplotlib.pyplot.plot`` for more options.
+        save : boolean, optional, default ``False``
+            If ``True``, save the figure with specified filename.
+        filename : str
+            File name of the image to save. Depends on the boolean ``save``.
+        plot_btis : bool
+            Plot the bad time intervals as red areas on the plot
+        """
+        import matplotlib.pyplot as plt
+        from .gti import get_btis
+
+        if ax is None:
+            plt.figure(attr)
+            ax = plt.gca()
+
+        if labels is None:
+            labels = ["Time (s)"] + [attr]
+
+        ylabel = labels[1]
+        xlabel = labels[0]
+
+        ax.plot(self.time, getattr(self, attr), marker, ds="steps-mid", label=attr, zorder=10)
+
+        if witherrors and attr + "_err" in self.array_attrs():
+            ax.errorbar(
+                self.time,
+                getattr(self, attr),
+                yerr=getattr(self, attr + "_err"),
+                fmt="o",
+            )
+
+        ax.set_ylabel(ylabel)
+        ax.set_xlabel(xlabel)
+
+        if title is not None:
+            ax.set_title(title)
+
+        if save:
+            if filename is None:
+                ax.figure.savefig("out.png")
+            else:
+                ax.figure.savefig(filename)
+
+        if plot_btis and self.gti is not None and len(self.gti) > 1:
+            tstart = min(self.time[0] - self.dt / 2, self.gti[0, 0])
+            tend = max(self.time[-1] + self.dt / 2, self.gti[-1, 1])
+            btis = get_btis(self.gti, tstart, tend)
+            for bti in btis:
+                plt.axvspan(bti[0], bti[1], alpha=0.5, color="r", zorder=10)
+        return ax
+
 
 def interpret_times(time: TTime, mjdref: float = 0) -> tuple[npt.ArrayLike, float]:
     """Understand the format of input times, and return seconds from MJDREF
@@ -555,6 +2249,9 @@ def interpret_times(time: TTime, mjdref: float = 0) -> tuple[npt.ArrayLike, floa
     Examples
     --------
     >>> import astropy.units as u
+    >>> newt, mjdref = interpret_times(None)
+    >>> newt is None
+    True
     >>> time = Time(57483, format='mjd')
     >>> newt, mjdref = interpret_times(time)
     >>> newt == 0
@@ -595,6 +2292,9 @@ def interpret_times(time: TTime, mjdref: float = 0) -> tuple[npt.ArrayLike, floa
     ...
     ValueError: Unknown time format: ...
     """
+    if time is None:
+        return None, mjdref
+
     if isinstance(time, TimeDelta):
         out_times = time.to("s").value
         return out_times, mjdref
@@ -626,3 +2326,47 @@ def interpret_times(time: TTime, mjdref: float = 0) -> tuple[npt.ArrayLike, floa
             pass
 
     raise ValueError(f"Unknown time format: {type(time)}")
+
+
+def reduce_precision_if_extended(
+    x, probe_types=["float128", "float96", "float80", "longdouble"], destination=float
+):
+    """Reduce a number to a standard float if extended precision.
+
+    Ignore all non-float types.
+
+    Parameters
+    ----------
+    x : float
+        The number to be reduced
+
+    Returns
+    -------
+    x_red : same type of input
+        The input, only reduce to ``float`` precision if ``np.float128``
+
+    Examples
+    --------
+    >>> x = 1.0
+    >>> val = reduce_precision_if_extended(x, probe_types=["float64"])
+    >>> val is x
+    True
+    >>> x = np.asanyarray(1.0).astype(int)
+    >>> val = reduce_precision_if_extended(x, probe_types=["float64"])
+    >>> val is x
+    True
+    >>> x = np.asanyarray([1.0]).astype(int)
+    >>> val = reduce_precision_if_extended(x, probe_types=["float64"])
+    >>> val is x
+    True
+    >>> x = np.asanyarray(1.0).astype(np.float64)
+    >>> reduce_precision_if_extended(x, probe_types=["float64"], destination=np.float32) is x
+    False
+    >>> x = np.asanyarray([1.0]).astype(np.float64)
+    >>> reduce_precision_if_extended(x, probe_types=["float64"], destination=np.float32) is x
+    False
+    """
+    if any([t in str(np.obj2sctype(x)) for t in probe_types]):
+        x_ret = x.astype(destination)
+        return x_ret
+    return x
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index a0cc3d77..e9e606e1 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -1406,6 +1406,76 @@ class Crossspectrum(StingrayObject):
             gti=gti,
         )
 
+    @staticmethod
+    def from_stingray_timeseries(
+        ts1,
+        ts2,
+        flux_attr,
+        error_flux_attr=None,
+        segment_size=None,
+        norm="none",
+        power_type="all",
+        silent=False,
+        fullspec=False,
+        use_common_mean=True,
+        gti=None,
+    ):
+        """Calculate AveragedCrossspectrum from two light curves
+
+        Parameters
+        ----------
+        ts1 : `stingray.Timeseries`
+            Time series from channel 1
+        ts2 : `stingray.Timeseries`
+            Time series from channel 2
+        flux_attr : `str`
+            What attribute of the time series will be used.
+
+        Other parameters
+        ----------------
+        error_flux_attr : `str`
+            What attribute of the time series will be used as error bar.
+        segment_size : float
+            The length, in seconds, of the light curve segments that will be averaged.
+            Only relevant (and required) for AveragedCrossspectrum
+        norm : str, default "frac"
+            The normalization of the periodogram. "abs" is absolute rms, "frac" is
+            fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+            unnormalized periodogram
+        use_common_mean : bool, default True
+            The mean of the light curve can be estimated in each interval, or on
+            the full light curve. This gives different results (Alston+2013).
+            Here we assume the mean is calculated on the full light curve, but
+            the user can set ``use_common_mean`` to False to calculate it on a
+            per-segment basis.
+        fullspec : bool, default False
+            Return the full periodogram, including negative frequencies
+        silent : bool, default False
+            Silence the progress bars
+        power_type : str, default 'all'
+            If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+            the real part
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals. Defaults to the common GTIs from the two input
+            objects. Could throw errors if these GTIs have overlaps with the
+            input object GTIs! If you're getting errors regarding your GTIs,
+            don't  use this and only give GTIs to the input objects before
+            making the cross spectrum.
+        """
+        return crossspectrum_from_timeseries(
+            ts1,
+            ts2,
+            flux_attr=flux_attr,
+            error_flux_attr=error_flux_attr,
+            segment_size=segment_size,
+            norm=norm,
+            power_type=power_type,
+            silent=silent,
+            fullspec=fullspec,
+            use_common_mean=use_common_mean,
+            gti=gti,
+        )
+
     @staticmethod
     def from_lc_iterable(
         iter_lc1,
@@ -2329,6 +2399,88 @@ def crossspectrum_from_lightcurve(
         This is likely to fill up your RAM on medium-sized datasets, and to
         slow down the computation when rebinning.
 
+    Returns
+    -------
+    spec : `AveragedCrossspectrum` or `Crossspectrum`
+        The output cross spectrum.
+    """
+    error_flux_attr = None
+
+    if lc1.err_dist == "gauss":
+        error_flux_attr = "_counts_err"
+
+    return crossspectrum_from_timeseries(
+        lc1,
+        lc2,
+        "_counts",
+        error_flux_attr=error_flux_attr,
+        segment_size=segment_size,
+        norm=norm,
+        power_type=power_type,
+        silent=silent,
+        fullspec=fullspec,
+        use_common_mean=use_common_mean,
+        gti=gti,
+        save_all=save_all,
+    )
+
+
+def crossspectrum_from_timeseries(
+    ts1,
+    ts2,
+    flux_attr,
+    error_flux_attr=None,
+    segment_size=None,
+    norm="none",
+    power_type="all",
+    silent=False,
+    fullspec=False,
+    use_common_mean=True,
+    gti=None,
+    save_all=False,
+):
+    """Calculate AveragedCrossspectrum from two time series
+
+    Parameters
+    ----------
+    ts1 : `stingray.StingrayTimeseries`
+        Time series from channel 1
+    ts2 : `stingray.StingrayTimeseries`
+        Time series from channel 2
+    flux_attr : `str`
+        What attribute of the time series will be used.
+
+    Other parameters
+    ----------------
+    error_flux_attr : `str`
+        What attribute of the time series will be used as error bar.
+    segment_size : float, default None
+        The length, in seconds, of the light curve segments that will be averaged
+    norm : str, default "frac"
+        The normalization of the periodogram. "abs" is absolute rms, "frac" is
+        fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+        unnormalized periodogram
+    use_common_mean : bool, default True
+        The mean of the light curve can be estimated in each interval, or on
+        the full light curve. This gives different results (Alston+2013).
+        Here we assume the mean is calculated on the full light curve, but
+        the user can set ``use_common_mean`` to False to calculate it on a
+        per-segment basis.
+    fullspec : bool, default False
+        Return the full periodogram, including negative frequencies
+    silent : bool, default False
+        Silence the progress bars
+    power_type : str, default 'all'
+        If 'all', give complex powers. If 'abs', the absolute value; if 'real',
+        the real part
+    gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+        Good Time intervals. Defaults to the common GTIs from the two input
+        objects
+    save_all : bool, default False
+        Save all intermediate spectra used for the final average. Use with care.
+        This is likely to fill up your RAM on medium-sized datasets, and to
+        slow down the computation when rebinning.
+
     Returns
     -------
     spec : `AveragedCrossspectrum` or `Crossspectrum`
@@ -2338,26 +2490,26 @@ def crossspectrum_from_lightcurve(
     # Suppress progress bar for single periodogram
     silent = silent or (segment_size is None)
     if gti is None:
-        gti = cross_two_gtis(lc1.gti, lc2.gti)
+        gti = cross_two_gtis(ts1.gti, ts2.gti)
 
     err1 = err2 = None
-    if lc1.err_dist == "gauss":
-        err1 = lc1._counts_err
-        err2 = lc2._counts_err
+    if error_flux_attr is not None:
+        err1 = getattr(ts1, error_flux_attr)
+        err2 = getattr(ts2, error_flux_attr)
 
     results = avg_cs_from_events(
-        lc1.time,
-        lc2.time,
+        ts1.time,
+        ts2.time,
         gti,
         segment_size,
-        lc1.dt,
+        ts1.dt,
         norm=norm,
         use_common_mean=use_common_mean,
         fullspec=fullspec,
         silent=silent,
         power_type=power_type,
-        fluxes1=lc1.counts,
-        fluxes2=lc2.counts,
+        fluxes1=getattr(ts1, flux_attr),
+        fluxes2=getattr(ts2, flux_attr),
         errors1=err1,
         errors2=err2,
         return_auxil=True,
diff --git a/stingray/events.py b/stingray/events.py
index 0f2e1811..90ce22ab 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -14,24 +14,44 @@ import numpy as np
 import numpy.random as ra
 from astropy.table import Table
 
+from stingray.utils import _int_sum_non_zero
+
 from .base import StingrayObject, StingrayTimeseries
 from .filters import get_deadtime_mask
 from .gti import append_gtis, check_separate, cross_gtis, generate_indices_of_boundaries
 from .io import load_events_and_gtis
 from .lightcurve import Lightcurve
-from .utils import assign_value_if_none, simon, interpret_times, njit
+from .utils import assign_value_if_none, simon, njit
+from .utils import histogram
 
 __all__ = ["EventList"]
 
 
 @njit
 def _from_lc_numba(times, counts, empty_times):
+    """Create a rough event list from a light curve.
+
+    This function creates as many events as the counts in each time bin of the light curve,
+    with event times equal to the light curve time stamps.
+
+    Parameters
+    ----------
+    times : array-like
+        Array of time stamps
+    counts : array-like
+        Array of counts
+    empty_times : array-like
+        Empty array to be filled with time stamps
+    """
     last = 0
     for t, c in zip(times, counts):
+        if c <= 0:
+            continue
         val = c + last
         empty_times[last:val] = t
         last = val
-    return empty_times
+    # If c < 0 in some cases, some times will be empty
+    return empty_times[:val]
 
 
 def simple_events_from_lc(lc):
@@ -39,6 +59,8 @@ def simple_events_from_lc(lc):
     Create an :class:`EventList` from a :class:`stingray.Lightcurve` object. Note that all
     events in a given time bin will have the same time stamp.
 
+    Bins with negative counts will be ignored.
+
     Parameters
     ----------
     lc: :class:`stingray.Lightcurve` object
@@ -52,14 +74,14 @@ def simple_events_from_lc(lc):
     Examples
     --------
     >>> from stingray import Lightcurve
-    >>> lc = Lightcurve([0, 1], [2, 3], dt=1)
+    >>> lc = Lightcurve([0, 1, 2], [2, 3, -1], dt=1)
     >>> ev = simple_events_from_lc(lc)
     >>> np.allclose(ev.time, [0, 0, 1, 1, 1])
     True
     """
-    times = _from_lc_numba(
-        lc.time, lc.counts.astype(int), np.zeros(np.sum(lc.counts).astype(int), dtype=float)
-    )
+    counts = lc.counts.astype(int)
+    allcounts = _int_sum_non_zero(counts)
+    times = _from_lc_numba(lc.time, counts, np.zeros(allcounts, dtype=float))
     return EventList(time=times, gti=lc.gti)
 
 
@@ -88,7 +110,7 @@ class EventList(StingrayTimeseries):
         The MJD used as a reference for the time array.
 
     ncounts: int
-        Number of desired data points in event list.
+        Number of desired data points in event list. Deprecated
 
     gtis: ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
         Good Time Intervals
@@ -191,40 +213,44 @@ class EventList(StingrayTimeseries):
         timesys=None,
         **other_kw,
     ):
-        StingrayObject.__init__(self)
-
-        self.energy = None if energy is None else np.asarray(energy)
-        self.notes = notes
-        self.dt = dt
-        self.mjdref = mjdref
-        self.gti = np.asarray(gti) if gti is not None else None
-        self.pi = None if pi is None else np.asarray(pi)
-        self.ncounts = ncounts
-        self.mission = mission
-        self.instr = instr
-        self.detector_id = detector_id
-        self.header = header
-        self.ephem = ephem
-        self.timeref = timeref
-        self.timesys = timesys
+        if ncounts is not None:
+            warnings.warn(
+                "The ncounts keyword does nothing, and is maintained for backwards compatibility.",
+                DeprecationWarning,
+            )
+
+        StingrayTimeseries.__init__(
+            self,
+            time=time,
+            energy=None if energy is None else np.asarray(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),
+            high_precision=high_precision,
+            mission=mission,
+            instr=instr,
+            header=header,
+            detector_id=detector_id,
+            ephem=ephem,
+            timeref=timeref,
+            timesys=timesys,
+            **other_kw,
+        )
 
         if other_kw != {}:
             warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")
 
-        if time is not None:
-            time, mjdref = interpret_times(time, mjdref)
-            if not high_precision:
-                self.time = np.asarray(time)
-            else:
-                self.time = np.asarray(time, dtype=np.longdouble)
-            self.ncounts = self.time.size
-        else:
-            self.time = None
-
         if (self.time is not None) and (self.energy is not None):
-            if self.time.size != self.energy.size:
+            if np.size(self.time) != np.size(self.energy):
                 raise ValueError("Lengths of time and energy must be equal.")
 
+    @property
+    def ncounts(self):
+        """Number of events in the event list."""
+        return self.n
+
     def to_lc(self, dt, tstart=None, tseg=None):
         """
         Convert event list to a :class:`stingray.Lightcurve` object.
@@ -246,14 +272,59 @@ class EventList(StingrayTimeseries):
         -------
         lc: :class:`stingray.Lightcurve` object
         """
-        if tstart is None and self.gti is not None:
-            tstart = self.gti[0][0]
-            tseg = self.gti[-1][1] - tstart
-
         return Lightcurve.make_lightcurve(
-            self.time, dt, tstart=tstart, gti=self.gti, tseg=tseg, mjdref=self.mjdref
+            self.time, dt, tstart=tstart, gti=self._gti, tseg=tseg, mjdref=self.mjdref
         )
 
+    def to_binned_timeseries(self, dt, array_attrs=None):
+        """Convert the event list to a binned :class:`stingray.StingrayTimeseries` object.
+
+        The result will be something similar to a light curve, but with arbitrary
+        attributes corresponding to a weighted sum of each specified attribute of
+        the event list.
+
+        E.g. if the event list has a ``q`` attribute, the final time series will
+        have a ``q`` attribute, which is the sum of all ``q`` values in each time bin.
+
+        Parameters
+        ----------
+        dt: float
+            Binning time of the light curve
+
+        Other Parameters
+        ----------------
+        array_attrs: list of str
+            List of attributes to be converted to light curve arrays. If None,
+            all array attributes will be converted.
+
+        Returns
+        -------
+        lc: :class:`stingray.Lightcurve` object
+        """
+        if array_attrs is None:
+            array_attrs = self.array_attrs()
+
+        ranges = [self.gti[0, 0], self.gti[-1, 1]]
+        nbins = int((ranges[1] - ranges[0]) / dt)
+        ranges = [ranges[0], ranges[0] + nbins * dt]
+        times = np.arange(ranges[0] + dt * 0.5, ranges[1], dt)
+
+        counts = histogram(self.time, range=ranges, bins=nbins)
+
+        attr_dict = dict(counts=counts)
+
+        for attr in array_attrs:
+            if getattr(self, attr, None) is not None:
+                logging.info(f"Creating the {attr} array")
+
+                attr_dict[attr] = histogram(
+                    self.time, bins=nbins, weights=getattr(self, attr), range=ranges
+                )
+        meta_attrs = dict((attr, getattr(self, attr)) for attr in self.meta_attrs())
+        new_ts = StingrayTimeseries(times, array_attrs=attr_dict, **meta_attrs)
+        new_ts.dt = dt
+        return new_ts
+
     def to_lc_iter(self, dt, segment_size=None):
         """Convert event list to a generator of Lightcurves.
 
@@ -317,6 +388,8 @@ class EventList(StingrayTimeseries):
         Create an :class:`EventList` from a :class:`stingray.Lightcurve` object. Note that all
         events in a given time bin will have the same time stamp.
 
+        Bins with negative counts will be ignored.
+
         Parameters
         ----------
         lc: :class:`stingray.Lightcurve` object
@@ -362,9 +435,9 @@ class EventList(StingrayTimeseries):
         if bin_time is not None:
             warnings.warn("Bin time will be ignored in simulate_times", DeprecationWarning)
 
-        self.time = simulate_times(lc, use_spline=use_spline)
-        self.gti = lc.gti
-        self.ncounts = len(self.time)
+        vals = simulate_times(lc, use_spline=use_spline)
+        self.time = vals
+        self._gti = lc.gti
 
     def simulate_energies(self, spectrum, use_spline=False):
         """
@@ -384,8 +457,8 @@ class EventList(StingrayTimeseries):
         """
         from .simulator.base import simulate_with_inverse_cdf
 
-        if self.ncounts is None:
-            simon("Either set time values or explicity provide counts.")
+        if self.ncounts is None or self.ncounts == 0:
+            simon("Simulating on an empty event list")
             return
 
         if isinstance(spectrum, list) or isinstance(spectrum, np.ndarray):
@@ -449,7 +522,7 @@ class EventList(StingrayTimeseries):
         order = np.argsort(self.time)
         return self.apply_mask(order, inplace=inplace)
 
-    def join(self, other):
+    def join(self, other, strategy="infer"):
         """
         Join two :class:`EventList` objects into one.
 
@@ -476,143 +549,24 @@ class EventList(StingrayTimeseries):
             If ``other`` is a list, it is assumed to be a list of :class:`EventList` objects
             and they are all joined, one by one.
 
+        Other parameters
+        ----------------
+        strategy : {"intersection", "union", "append", "infer", "none"}
+            Method to use to merge the GTIs. If "intersection", the GTIs are merged
+            using the intersection of the GTIs. If "union", the GTIs are merged
+            using the union of the GTIs. If "none", a single GTI with the minimum and
+            the maximum time stamps of all GTIs is returned. If "infer", the strategy
+            is decided based on the GTIs. If there are no overlaps, "union" is used,
+            otherwise "intersection" is used. If "append", the GTIs are simply appended
+            but they must be mutually exclusive.
+
         Returns
         -------
         `ev_new` : :class:`EventList` object
             The resulting :class:`EventList` object.
         """
 
-        ev_new = EventList()
-
-        if isinstance(other, EventList):
-            others = [other]
-        else:
-            others = other
-
-        # First of all, check if there are empty event lists
-        for obj in others:
-            if getattr(obj, "time", None) is None or np.size(obj.time) == 0:
-                warnings.warn("One of the event lists you are joining is empty.")
-                others.remove(obj)
-
-        if len(others) == 0:
-            return copy.deepcopy(self)
-
-        for i, other in enumerate(others):
-            # Tolerance for MJDREF:1 microsecond
-            if not np.isclose(self.mjdref, other.mjdref, atol=1e-6 / 86400):
-                warnings.warn("Attribute mjdref is different in the event lists being merged.")
-                others[i] = other.change_mjdref(self.mjdref)
-
-        all_objs = [self] + others
-
-        dts = list(set([getattr(obj, "dt", None) for obj in all_objs]))
-        if len(dts) != 1:
-            warnings.warn("The time resolution is different. Using the rougher by default")
-
-            ev_new.dt = np.max(dts)
-
-        all_time_arrays = [obj.time for obj in all_objs if obj.time is not None]
-
-        ev_new.time = np.concatenate(all_time_arrays)
-        order = np.argsort(ev_new.time)
-        ev_new.time = ev_new.time[order]
-
-        def _get_set_from_many_lists(lists):
-            """Make a single set out of many lists."""
-            all_vals = []
-            for l in lists:
-                all_vals += l
-            return set(all_vals)
-
-        def _get_all_array_attrs(objs):
-            """Get all array attributes from the event lists being merged. Do not include time."""
-            all_attrs = []
-            for obj in objs:
-                if obj.time is not None and len(obj.time) > 0:
-                    all_attrs += obj.array_attrs()
-
-            all_attrs = list(set(all_attrs))
-            if "time" in all_attrs:
-                all_attrs.remove("time")
-            return all_attrs
-
-        for attr in _get_all_array_attrs(all_objs):
-            # if it's here, it means that it's an array attr in at least one object.
-            # So, everywhere it's None, it needs to be set to 0s of the same length as time
-            new_attr_values = []
-            for obj in all_objs:
-                if getattr(obj, attr, None) is None:
-                    warnings.warn(
-                        f"The {attr} array is empty in one of the event lists being merged. Setting it to NaN for the affected events"
-                    )
-                    new_attr_values.append(np.zeros_like(obj.time) + np.nan)
-                else:
-                    new_attr_values.append(getattr(obj, attr))
-            new_attr = np.concatenate(new_attr_values)[order]
-            setattr(ev_new, attr, new_attr)
-
-        if np.all([obj.gti is None for obj in all_objs]):
-            ev_new.gti = None
-        else:
-            all_gti_lists = []
-
-            for obj in all_objs:
-                if obj.gti is None and len(obj.time) > 0:
-                    obj.gti = assign_value_if_none(
-                        obj.gti,
-                        np.asarray([[obj.time[0] - obj.dt / 2, obj.time[-1] + obj.dt / 2]]),
-                    )
-                if obj.gti is not None:
-                    all_gti_lists.append(obj.gti)
-
-            new_gtis = all_gti_lists[0]
-            for gti in all_gti_lists[1:]:
-                if check_separate(new_gtis, gti):
-                    new_gtis = append_gtis(new_gtis, gti)
-                    warnings.warn(
-                        "GTIs in these two event lists do not overlap at all."
-                        "Merging instead of returning an overlap."
-                    )
-                else:
-                    new_gtis = cross_gtis([new_gtis, gti])
-            ev_new.gti = new_gtis
-
-        all_meta_attrs = _get_set_from_many_lists([obj.meta_attrs() for obj in all_objs])
-        # The attributes being treated separately are removed from the standard treatment
-        # When energy, pi etc. are None, they might appear in the meta_attrs, so we
-        # also add them to the list of attributes to be removed if present.
-        to_remove = ["gti", "header", "ncounts", "dt"] + ev_new.array_attrs()
-
-        for attrs in to_remove:
-            if attrs in all_meta_attrs:
-                all_meta_attrs.remove(attrs)
-
-        logging.info("The header attribute will be removed from the output event list.")
-
-        def _safe_concatenate(a, b):
-            if isinstance(a, str) and isinstance(b, str):
-                return a + "," + b
-            else:
-                if isinstance(a, tuple):
-                    return a + (b,)
-                return (a, b)
-
-        for attr in all_meta_attrs:
-            self_attr = getattr(self, attr, None)
-            new_val = self_attr
-            for other in others:
-                other_attr = getattr(other, attr, None)
-                if self_attr != other_attr:
-                    warnings.warn(
-                        "Attribute " + attr + " is different in the event lists being merged."
-                    )
-                    new_val = _safe_concatenate(new_val, other_attr)
-            setattr(ev_new, attr, new_val)
-
-        ev_new.mjdref = self.mjdref
-
-        return ev_new
+        return self._join_timeseries(other, strategy=strategy, ignore_meta=["header", "ncounts"])
 
     @classmethod
     def read(cls, filename, fmt=None, **kwargs):
@@ -719,50 +673,6 @@ class EventList(StingrayTimeseries):
 
         return self.apply_mask(mask, inplace=inplace)
 
-    def apply_mask(self, mask, inplace=False):
-        """Apply a mask to all array attributes of the event list
-
-        Parameters
-        ----------
-        mask : array of ``bool``
-            The mask. Has to be of the same length as ``self.time``
-
-        Other parameters
-        ----------------
-        inplace : bool
-            If True, overwrite the current event list. Otherwise, return a new one.
-
-        Examples
-        --------
-        >>> evt = EventList(time=[0, 1, 2], mission="nustar")
-        >>> evt.bubuattr = [222, 111, 333]
-        >>> newev0 = evt.apply_mask([True, True, False], inplace=False);
-        >>> newev1 = evt.apply_mask([True, True, False], inplace=True);
-        >>> newev0.mission == "nustar"
-        True
-        >>> np.allclose(newev0.time, [0, 1])
-        True
-        >>> np.allclose(newev0.bubuattr, [222, 111])
-        True
-        >>> np.allclose(newev1.time, [0, 1])
-        True
-        >>> evt is newev1
-        True
-        """
-        array_attrs = self.array_attrs()
-
-        if inplace:
-            new_ev = self
-        else:
-            new_ev = EventList()
-            for attr in self.meta_attrs():
-                setattr(new_ev, attr, copy.deepcopy(getattr(self, attr)))
-
-        for attr in array_attrs:
-            if hasattr(self, attr) and getattr(self, attr) is not None:
-                setattr(new_ev, attr, copy.deepcopy(np.asarray(getattr(self, attr))[mask]))
-        return new_ev
-
     def apply_deadtime(self, deadtime, inplace=False, **kwargs):
         """Apply deadtime filter to this event list.
 
@@ -788,7 +698,7 @@ class EventList(StingrayTimeseries):
         Examples
         --------
         >>> events = np.array([1, 1.05, 1.07, 1.08, 1.1, 2, 2.2, 3, 3.1, 3.2])
-        >>> events = EventList(events)
+        >>> events = EventList(events, gti=[[0, 3.3]])
         >>> events.pi=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
         >>> events.energy=np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
         >>> events.mjdref = 10
diff --git a/stingray/gti.py b/stingray/gti.py
index 7c9e97bf..4844de66 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -436,11 +436,11 @@ def create_gti_mask(
     new_gtis : ``Nx2`` array
         An array of new GTIs created by this function.
     """
-    gtis = np.array(gtis, dtype=np.longdouble)
-    if time.size == 0:
+    if time is None or np.size(time) == 0:
         raise ValueError("Passing an empty time array to create_gti_mask")
-    if gtis.size == 0:
+    if gtis is None or np.size(gtis) == 0:
         raise ValueError("Passing an empty GTI array to create_gti_mask")
+    gtis = np.array(gtis, dtype=np.longdouble)
 
     mask = np.zeros(len(time), dtype=bool)
 
@@ -968,7 +968,7 @@ def join_equal_gti_boundaries(gti, threshold=0.0):
     return np.asarray(ng)
 
 
-def merge_gtis(gti_list, gti_treatment):
+def merge_gtis(gti_list, strategy):
     """Merge a list of GTIs using the specified method.
 
     Invalid GTI lists (None or empty) are ignored.
@@ -980,7 +980,7 @@ def merge_gtis(gti_list, gti_treatment):
 
     Other parameters
     ----------------
-    gti_treatment : {"intersection", "union", "append", "infer", "none"}
+    strategy : {"intersection", "union", "append", "infer", "none"}
         Method to use to merge the GTIs. If "intersection", the GTIs are merged
         using the intersection of the GTIs. If "union", the GTIs are merged
         using the union of the GTIs. If "none", a single GTI with the minimum and
@@ -1021,7 +1021,7 @@ def merge_gtis(gti_list, gti_treatment):
     if len(all_gti_lists) == 0:
         return None
 
-    if gti_treatment == "none":
+    if strategy == "none":
         return np.asarray([[global_min, global_max]])
 
     if len(all_gti_lists) == 1:
@@ -1030,20 +1030,20 @@ def merge_gtis(gti_list, gti_treatment):
     cross = cross_gtis(all_gti_lists)
     if len(cross) == 0:
         cross = None
-    if gti_treatment == "infer":
+    if strategy == "infer":
         if cross is None:
-            gti_treatment = "union"
+            strategy = "union"
         else:
-            gti_treatment = "intersection"
+            strategy = "intersection"
 
-    if gti_treatment == "intersection":
+    if strategy == "intersection":
         return cross
 
     gti0 = all_gti_lists[0]
     for gti in all_gti_lists[1:]:
-        if gti_treatment == "union":
+        if strategy == "union":
             gti0 = join_gtis(gti0, gti)
-        elif gti_treatment == "append":
+        elif strategy == "append":
             gti0 = append_gtis(gti0, gti)
     return gti0
 
diff --git a/stingray/io.py b/stingray/io.py
index 06f57c37..2bf8ec59 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -28,6 +28,13 @@ except ImportError:
     _H5PY_INSTALLED = False
 
 
+HAS_128 = True
+try:
+    np.float128
+except AttributeError:  # pragma: no cover
+    HAS_128 = False
+
+
 def rough_calibration(pis, mission):
     """Make a rough conversion betwenn PI channel and energy.
 
@@ -916,3 +923,83 @@ def savefig(filename, **kwargs):
         )
 
     plt.savefig(filename, **kwargs)
+
+
+def _can_save_longdouble(probe_file: str, fmt: str) -> bool:
+    """Check if a given file format can save tables with longdoubles.
+
+    Try to save a table with a longdouble column, and if it doesn't work, catch the exception.
+    If the exception is related to longdouble, return False (otherwise just raise it, this
+    would mean there are larger problems that need to be solved). In this case, also warn that
+    probably part of the data will not be saved.
+
+    If no exception is raised, return True.
+
+    Parameters
+    ----------
+    probe_file : str
+        The name of the file to be used for probing
+    fmt : str
+        The format to be used for probing, in the ``format`` argument of ``Table.write``
+
+    Returns
+    -------
+    yes_it_can : bool
+        Whether the format can serialize the metadata
+    """
+    if not HAS_128:  # pragma: no cover
+        # There are no known issues with saving longdoubles where numpy.float128 is not defined
+        return True
+
+    try:
+        Table({"a": np.arange(0, 3, 1.212314).astype(np.float128)}).write(
+            probe_file, format=fmt, overwrite=True
+        )
+        yes_it_can = True
+        os.unlink(probe_file)
+    except ValueError as e:
+        if "float128" not in str(e):  # pragma: no cover
+            raise
+        warnings.warn(
+            f"{fmt} output does not allow saving metadata at maximum precision. "
+            "Converting to lower precision"
+        )
+        yes_it_can = False
+    return yes_it_can
+
+
+def _can_serialize_meta(probe_file: str, fmt: str) -> bool:
+    """
+    Try to save a table with meta to be serialized, and if it doesn't work, catch the exception.
+    If the exception is related to serialization, return False (otherwise just raise it, this
+    would mean there are larger problems that need to be solved). In this case, also warn that
+    probably part of the data will not be saved.
+
+    If no exception is raised, return True.
+
+    Parameters
+    ----------
+    probe_file : str
+        The name of the file to be used for probing
+    fmt : str
+        The format to be used for probing, in the ``format`` argument of ``Table.write``
+
+    Returns
+    -------
+    yes_it_can : bool
+        Whether the format can serialize the metadata
+    """
+    try:
+        Table({"a": [3]}).write(probe_file, overwrite=True, format=fmt, serialize_meta=True)
+
+        os.unlink(probe_file)
+        yes_it_can = True
+    except TypeError as e:
+        if "serialize_meta" not in str(e):  # pragma: no cover
+            raise
+        warnings.warn(
+            f"{fmt} output does not serialize the metadata at the moment. "
+            "Some attributes will be lost."
+        )
+        yes_it_can = False
+    return yes_it_can
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 9b76965f..82f69c74 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -17,7 +17,7 @@ from astropy.table import Table
 from astropy.time import TimeDelta, Time
 from astropy import units as u
 
-from stingray.base import StingrayTimeseries
+from stingray.base import StingrayTimeseries, reduce_precision_if_extended
 import stingray.utils as utils
 from stingray.exceptions import StingrayError
 from stingray.gti import (
@@ -33,12 +33,12 @@ from stingray.utils import (
     baseline_als,
     poisson_symmetrical_errors,
     simon,
-    interpret_times,
     is_sorted,
     check_isallfinite,
 )
 from stingray.io import lcurve_from_fits
 from stingray import bexvar
+from stingray.base import interpret_times
 
 __all__ = ["Lightcurve"]
 
@@ -229,11 +229,6 @@ class Lightcurve(StingrayTimeseries):
         header=None,
         **other_kw,
     ):
-        StingrayTimeseries.__init__(self)
-
-        if other_kw != {}:
-            warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")
-
         self._time = None
         self._mask = None
         self._counts = None
@@ -245,6 +240,11 @@ class Lightcurve(StingrayTimeseries):
         self._bin_lo = None
         self._bin_hi = None
         self._n = None
+        StingrayTimeseries.__init__(self)
+
+        if other_kw != {}:
+            warnings.warn(f"Unrecognized keywords: {list(other_kw.keys())}")
+
         self.mission = mission
         self.instr = instr
         self.header = header
@@ -256,7 +256,6 @@ class Lightcurve(StingrayTimeseries):
         self.mjdref = mjdref
 
         if time is None or len(time) == 0:
-            warnings.warn("No time values passed to Lightcurve object!")
             return
 
         if counts is None or np.size(time) != np.size(counts):
@@ -327,6 +326,13 @@ class Lightcurve(StingrayTimeseries):
         if gti is not None:
             self._gti = np.asarray(gti)
 
+        if os.name == "nt":
+            warnings.warn(
+                "On Windows, the size of an integer is 32 bits. "
+                "To avoid integer overflow, I'm converting the input array to float"
+            )
+            counts = counts.astype(float)
+
         if input_counts:
             self._counts = np.asarray(counts)
             self._counts_err = err
@@ -349,12 +355,6 @@ class Lightcurve(StingrayTimeseries):
 
         if not skip_checks:
             self.check_lightcurve()
-        if os.name == "nt":
-            warnings.warn(
-                "On Windows, the size of an integer is 32 bits. "
-                "To avoid integer overflow, I'm converting the input array to float"
-            )
-            counts = counts.astype(float)
 
     @property
     def time(self):
@@ -362,37 +362,14 @@ class Lightcurve(StingrayTimeseries):
 
     @time.setter
     def time(self, value):
-        value = np.asarray(value)
-        if not value.shape == self.time.shape:
-            raise ValueError("Can only assign new times of the same shape as " "the original array")
+        value = self._validate_and_format(value, "time", "time")
+        if value is None:
+            for attr in self.internal_array_attrs():
+                setattr(self, attr, None)
         self._time = value
         self._bin_lo = None
         self._bin_hi = None
 
-    @property
-    def gti(self):
-        if self._gti is None:
-            self._gti = np.asarray([[self.tstart, self.tstart + self.tseg]])
-        return self._gti
-
-    @gti.setter
-    def gti(self, value):
-        value = np.asarray(value)
-        self._gti = value
-        self._mask = None
-
-    @property
-    def mask(self):
-        if self._mask is None:
-            self._mask = create_gti_mask(self.time, self.gti, dt=self.dt)
-        return self._mask
-
-    @property
-    def n(self):
-        if self._n is None:
-            self._n = self.counts.shape[0]
-        return self._n
-
     @property
     def meanrate(self):
         if self._meanrate is None:
@@ -408,7 +385,7 @@ class Lightcurve(StingrayTimeseries):
     @property
     def counts(self):
         counts = self._counts
-        if self._counts is None:
+        if self._counts is None and self._countrate is not None:
             counts = self._countrate * self.dt
             # If not in low-memory regime, cache the values
             if not self.low_memory or self.input_counts:
@@ -418,11 +395,7 @@ class Lightcurve(StingrayTimeseries):
 
     @counts.setter
     def counts(self, value):
-        value = np.asarray(value)
-        if not value.shape == self.counts.shape:
-            raise ValueError(
-                "Can only assign new counts array of the same " "shape as the original array"
-            )
+        value = self._validate_and_format(value, "counts", "time")
         self._counts = value
         self._countrate = None
         self._meancounts = None
@@ -450,18 +423,15 @@ class Lightcurve(StingrayTimeseries):
 
     @counts_err.setter
     def counts_err(self, value):
-        value = np.asarray(value)
-        if not value.shape == self.counts.shape:
-            raise ValueError(
-                "Can only assign new error array of the same " "shape as the original array"
-            )
+        value = self._validate_and_format(value, "counts_err", "counts")
+
         self._counts_err = value
         self._countrate_err = None
 
     @property
     def countrate(self):
         countrate = self._countrate
-        if countrate is None:
+        if countrate is None and self._counts is not None:
             countrate = self._counts / self.dt
             # If not in low-memory regime, cache the values
             if not self.low_memory or not self.input_counts:
@@ -471,11 +441,8 @@ class Lightcurve(StingrayTimeseries):
 
     @countrate.setter
     def countrate(self, value):
-        value = np.asarray(value)
-        if not value.shape == self.countrate.shape:
-            raise ValueError(
-                "Can only assign new countrate array of the same " "shape as the original array"
-            )
+        value = self._validate_and_format(value, "countrate", "time")
+
         self._countrate = value
         self._counts = None
         self._meancounts = None
@@ -500,11 +467,8 @@ class Lightcurve(StingrayTimeseries):
 
     @countrate_err.setter
     def countrate_err(self, value):
-        value = np.asarray(value)
-        if not value.shape == self.countrate.shape:
-            raise ValueError(
-                "Can only assign new error array of the same " "shape as the original array"
-            )
+        value = self._validate_and_format(value, "countrate_err", "countrate")
+
         self._countrate_err = value
         self._counts_err = None
 
@@ -594,7 +558,7 @@ class Lightcurve(StingrayTimeseries):
                 "Only use with LombScargleCrossspectrum, LombScarglePowerspectrum and QPO using GPResult"
             )
 
-    def _operation_with_other_lc(self, other, operation):
+    def _operation_with_other_obj(self, other, operation):
         """
         Helper method to codify an operation of one light curve with another (e.g. add, subtract, ...).
         Takes into account the GTIs correctly, and returns a new :class:`Lightcurve` object.
@@ -618,6 +582,11 @@ class Lightcurve(StingrayTimeseries):
             other = other.change_mjdref(self.mjdref)
 
         common_gti = cross_two_gtis(self.gti, other.gti)
+        if not np.array_equal(self.gti, common_gti):
+            warnings.warn(
+                "The good time intervals in the two time series are different. Data outside the "
+                "common GTIs will be discarded."
+            )
         mask_self = create_gti_mask(self.time, common_gti, dt=self.dt)
         mask_other = create_gti_mask(other.time, common_gti, dt=other.dt)
 
@@ -689,7 +658,7 @@ class Lightcurve(StingrayTimeseries):
         True
         """
 
-        return self._operation_with_other_lc(other, np.add)
+        return self._operation_with_other_obj(other, np.add)
 
     def __sub__(self, other):
         """
@@ -710,7 +679,7 @@ class Lightcurve(StingrayTimeseries):
         >>> count1 = [600, 1200, 800]
         >>> count2 = [300, 100, 400]
         >>> gti1 = [[0, 35]]
-        >>> gti2 = [[5, 40]]
+        >>> gti2 = [[0, 35]]
         >>> lc1 = Lightcurve(time, count1, gti=gti1, dt=10)
         >>> lc2 = Lightcurve(time, count2, gti=gti2, dt=10)
         >>> lc = lc1 - lc2
@@ -718,7 +687,7 @@ class Lightcurve(StingrayTimeseries):
         True
         """
 
-        return self._operation_with_other_lc(other, np.subtract)
+        return self._operation_with_other_obj(other, np.subtract)
 
     def __neg__(self):
         """
@@ -750,24 +719,6 @@ class Lightcurve(StingrayTimeseries):
 
         return lc_new
 
-    def __len__(self):
-        """
-        Return the number of time bins of a light curve.
-
-        This method implements overrides the ``len`` function for a :class:`Lightcurve`
-        object and returns the length of the ``time`` array (which should be equal to the
-        length of the ``counts`` and ``countrate`` arrays).
-
-        Examples
-        --------
-        >>> time = [1, 2, 3]
-        >>> count = [100, 200, 300]
-        >>> lc = Lightcurve(time, count, dt=1)
-        >>> len(lc)
-        3
-        """
-        return self.n
-
     def __getitem__(self, index):
         """
         Return the corresponding count value at the index or a new :class:`Lightcurve`
@@ -828,28 +779,6 @@ class Lightcurve(StingrayTimeseries):
         else:
             raise IndexError("The index must be either an integer or a slice " "object !")
 
-    def __eq__(self, other_lc):
-        """
-        Compares two :class:`Lightcurve` objects.
-
-        Light curves are equal only if their counts as well as times at which those counts occur equal.
-
-        Examples
-        --------
-        >>> time = [1, 2, 3]
-        >>> count1 = [100, 200, 300]
-        >>> count2 = [100, 200, 300]
-        >>> lc1 = Lightcurve(time, count1, dt=1)
-        >>> lc2 = Lightcurve(time, count2, dt=1)
-        >>> lc1 == lc2
-        True
-        """
-        if not isinstance(other_lc, Lightcurve):
-            raise ValueError("Lightcurve can only be compared with a Lightcurve Object")
-        if np.allclose(self.time, other_lc.time) and np.allclose(self.counts, other_lc.counts):
-            return True
-        return False
-
     def baseline(self, lam, p, niter=10, offset_correction=False):
         """Calculate the baseline of the light curve, accounting for GTIs.
 
@@ -1237,85 +1166,7 @@ class Lightcurve(StingrayTimeseries):
         True
         """
 
-        if not isinstance(method, str):
-            raise TypeError("method key word argument is not " "a string !")
-
-        if method.lower() not in ["index", "time"]:
-            raise ValueError("Unknown method type " + method + ".")
-
-        if method.lower() == "index":
-            new_lc = self._truncate_by_index(start, stop)
-        else:
-            new_lc = self._truncate_by_time(start, stop)
-        new_lc.tstart = new_lc.gti[0, 0]
-        new_lc.tseg = new_lc.gti[-1, 1] - new_lc.gti[0, 0]
-        return new_lc
-
-    def _truncate_by_index(self, start, stop):
-        """Private method for truncation using index values."""
-
-        new_lc = self.apply_mask(slice(start, stop))
-
-        dtstart = dtstop = new_lc.dt
-        if isinstance(self.dt, Iterable):
-            dtstart = self.dt[0]
-            dtstop = self.dt[-1]
-
-        gti = cross_two_gtis(
-            self.gti, np.asarray([[new_lc.time[0] - 0.5 * dtstart, new_lc.time[-1] + 0.5 * dtstop]])
-        )
-
-        new_lc.gti = gti
-
-        return new_lc
-
-    def _truncate_by_time(self, start, stop):
-        """Helper method for truncation using time values.
-
-        Parameters
-        ----------
-        start : float
-            start time for new light curve; all time bins before this time will be discarded
-
-        stop : float
-            stop time for new light curve; all time bins after this point will be discarded
-
-        Returns
-        -------
-            new_lc : Lightcurve
-                A new :class:`Lightcurve` object with the truncated time bins
-
-        """
-
-        if stop is not None:
-            if start > stop:
-                raise ValueError("start time must be less than stop time!")
-
-        if not start == 0:
-            start = self.time.searchsorted(start)
-
-        if stop is not None:
-            stop = self.time.searchsorted(stop)
-
-        return self._truncate_by_index(start, stop)
-
-    def meta_attrs(self):
-        """Extends StingrayObject.meta_attrs to the specifics of Lightcurve."""
-        attrs = super().meta_attrs()
-        sure_array = ["counts", "counts_err", "countrate", "countrate_err"]
-        for attr in sure_array:
-            if attr in attrs:
-                attrs.remove(attr)
-        return attrs
-
-    def array_attrs(self):
-        """Extends StingrayObject.array_attrs to the specifics of Lightcurve."""
-        attrs = super().array_attrs()
-        sure_array = ["counts", "counts_err", "countrate", "countrate_err"]
-        for attr in sure_array:
-            if attr not in attrs:
-                attrs.append(attr)
-        return attrs
+        return super().truncate(start=start, stop=stop, method=method)
 
     def split(self, min_gap, min_points=1):
         """
@@ -1625,13 +1476,48 @@ class Lightcurve(StingrayTimeseries):
             skip_checks=skip_checks,
         )
 
-    def to_astropy_timeseries(self):
-        return self._to_astropy_object(kind="timeseries")
+    def to_astropy_timeseries(self, **kwargs):
+        """Save the light curve to an :class:`astropy.timeseries.TimeSeries` object.
+
+        The time array and all the array attributes become columns. The meta attributes become
+        metadata of the :class:`astropy.timeseries.TimeSeries` object.
+        The time array is saved as a TimeDelta object.
+
+        Other Parameters
+        ----------------
+        no_longdouble : bool, default False
+            If True, the data are converted to double precision before being saved.
+            This is useful, e.g., for saving to FITS files, which do not support long double precision.
+        """
+        return self._to_astropy_object(kind="timeseries", **kwargs)
+
+    def to_astropy_table(self, **kwargs):
+        """Save the light curve to an :class:`astropy.table.Table` object.
 
-    def to_astropy_table(self):
-        return self._to_astropy_object(kind="table")
+        The time array and all the array attributes become columns. The meta attributes become
+        metadata of the :class:`astropy.table.Table` object.
 
-    def _to_astropy_object(self, kind="table"):
+        Other Parameters
+        ----------------
+        no_longdouble : bool, default False
+            If True, the data are converted to double precision before being saved.
+            This is useful, e.g., for saving to FITS files, which do not support long double precision.
+        """
+        return self._to_astropy_object(kind="table", **kwargs)
+
+    def _to_astropy_object(self, kind="table", no_longdouble=False):
+        """Save the light curve to an :class:`astropy.table.Table` or :class:`astropy.timeseries.TimeSeries` object.
+
+        If ``kind`` is ``timeseries``, the time array and all the array attributes become columns.
+
+        Other Parameters
+        ----------------
+        kind : str, default ``table``
+            The type of object to return. Accepted values are ``table`` or ``timeseries``.
+        no_longdouble : bool, default False
+            If True, the data are converted to double precision before being saved.
+            This is useful, e.g., for saving to FITS files, which do not support long double precision.
+        """
         data = {}
 
         for attr in [
@@ -1643,15 +1529,22 @@ class Lightcurve(StingrayTimeseries):
             "_bin_hi",
         ]:
             if hasattr(self, attr) and getattr(self, attr) is not None:
-                data[attr.lstrip("_")] = np.asarray(getattr(self, attr))
+                vals = np.asarray(getattr(self, attr))
+                if no_longdouble:
+                    vals = reduce_precision_if_extended(vals)
+                data[attr.lstrip("_")] = vals
+
+        time_array = self.time
+        if no_longdouble:
+            time_array = reduce_precision_if_extended(time_array)
 
         if kind.lower() == "table":
-            data["time"] = self.time
+            data["time"] = time_array
             ts = Table(data)
         elif kind.lower() == "timeseries":
             from astropy.timeseries import TimeSeries
 
-            ts = TimeSeries(data=data, time=TimeDelta(self.time * u.s))
+            ts = TimeSeries(data=data, time=TimeDelta(time_array * u.s))
         else:  # pragma: no cover
             raise ValueError("Invalid kind (accepted: table or timeseries)")
 
@@ -1666,6 +1559,9 @@ class Lightcurve(StingrayTimeseries):
             "err_dist",
         ]:
             if hasattr(self, attr) and getattr(self, attr) is not None:
+                vals = getattr(self, attr)
+                if no_longdouble:
+                    vals = reduce_precision_if_extended(vals)
                 ts.meta[attr.lstrip("_")] = getattr(self, attr)
 
         return ts
@@ -1854,130 +1750,6 @@ class Lightcurve(StingrayTimeseries):
 
         return super().read(filename=filename, fmt=fmt)
 
-    def split_by_gti(self, gti=None, min_points=2):
-        """
-        Split the current :class:`Lightcurve` object into a list of :class:`Lightcurve` objects, one
-        for each continuous GTI segment as defined in the ``gti`` attribute.
-
-        Parameters
-        ----------
-        min_points : int, default 1
-            The minimum number of data points in each light curve. Light
-            curves with fewer data points will be ignored.
-
-        Returns
-        -------
-        list_of_lcs : list
-            A list of :class:`Lightcurve` objects, one for each GTI segment
-        """
-
-        if gti is None:
-            gti = self.gti
-
-        list_of_lcs = []
-
-        start_bins, stop_bins = gti_border_bins(gti, self.time, self.dt)
-        for i in range(len(start_bins)):
-            start = start_bins[i]
-            stop = stop_bins[i]
-
-            if (stop - start) < min_points:
-                continue
-
-            new_gti = np.array([gti[i]])
-            mask = create_gti_mask(self.time, new_gti)
-
-            # Note: GTIs are consistent with default in this case!
-            new_lc = self.apply_mask(mask)
-            new_lc.gti = new_gti
-
-            list_of_lcs.append(new_lc)
-
-        return list_of_lcs
-
-    def apply_mask(self, mask, inplace=False):
-        """Apply a mask to all array attributes of the event list
-
-        Parameters
-        ----------
-        mask : array of ``bool``
-            The mask. Has to be of the same length as ``self.time``
-
-        Other parameters
-        ----------------
-        inplace : bool
-            If True, overwrite the current light curve. Otherwise, return a new one.
-
-        Examples
-        --------
-        >>> lc = Lightcurve(time=[0, 1, 2], counts=[2, 3, 4], mission="nustar")
-        >>> lc.bubuattr = [222, 111, 333]
-        >>> newlc0 = lc.apply_mask([True, True, False], inplace=False);
-        >>> newlc1 = lc.apply_mask([True, True, False], inplace=True);
-        >>> newlc0.mission == "nustar"
-        True
-        >>> np.allclose(newlc0.time, [0, 1])
-        True
-        >>> np.allclose(newlc0.bubuattr, [222, 111])
-        True
-        >>> np.allclose(newlc1.time, [0, 1])
-        True
-        >>> lc is newlc1
-        True
-        """
-        array_attrs = self.array_attrs()
-
-        self._mask = self._n = None
-        if isinstance(self.dt, Iterable):
-            new_dt = self.dt[mask]
-        else:
-            new_dt = self.dt
-        if inplace:
-            new_ev = self
-            # If they don't exist, they get set
-            self.counts, self.counts_err
-            # eliminate possible conflicts
-            self._countrate = self._countrate_err = None
-            # Set time, counts and errors
-            self._time = self._time[mask]
-            self._counts = self._counts[mask]
-            if self._counts_err is not None:
-                self._counts_err = self._counts_err[mask]
-            new_ev.dt = new_dt
-        else:
-            with warnings.catch_warnings():
-                warnings.filterwarnings(
-                    "ignore", message="Some functionalities of Stingray Lightcurve.*"
-                )
-                new_ev = Lightcurve(
-                    time=self.time[mask],
-                    counts=self.counts[mask],
-                    skip_checks=True,
-                    gti=self.gti,
-                    dt=new_dt,
-                )
-            if self._counts_err is not None:
-                new_ev.counts_err = self.counts_err[mask]
-            for attr in self.meta_attrs():
-                try:
-                    setattr(new_ev, attr, copy.deepcopy(getattr(self, attr)))
-                except AttributeError:
-                    continue
-        for attr in array_attrs:
-            if hasattr(self, "_" + attr) or attr in [
-                "time",
-                "counts",
-                "counts_err",
-                "dt",
-                "_time",
-                "_counts",
-                "_counts_err",
-            ]:
-                continue
-            if hasattr(self, attr) and getattr(self, attr) is not None:
-                setattr(new_ev, attr, copy.deepcopy(np.asarray(getattr(self, attr))[mask]))
-        return new_ev
-
     def apply_gtis(self, inplace=True):
         """
         Apply GTIs to a light curve. Filters the ``time``, ``counts``,
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 6fb037bf..cdadaa79 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -549,6 +549,63 @@ class Powerspectrum(Crossspectrum):
             use_common_mean=use_common_mean,
         )
 
+    @staticmethod
+    def from_stingray_timeseries(
+        ts,
+        flux_attr,
+        error_flux_attr=None,
+        segment_size=None,
+        norm="none",
+        silent=False,
+        use_common_mean=True,
+        gti=None,
+    ):
+        """Calculate AveragedPowerspectrum from a time series.
+
+        Parameters
+        ----------
+        ts : `stingray.Timeseries`
+            Input Time Series
+        flux_attr : `str`
+            What attribute of the time series will be used.
+
+        Other parameters
+        ----------------
+        error_flux_attr : `str`
+            What attribute of the time series will be used as error bar.
+        segment_size : float
+            The length, in seconds, of the light curve segments that will be averaged.
+            Only relevant (and required) for AveragedCrossspectrum
+        norm : str, default "frac"
+            The normalization of the periodogram. "abs" is absolute rms, "frac" is
+            fractional rms, "leahy" is Leahy+83 normalization, and "none" is the
+            unnormalized periodogram
+        use_common_mean : bool, default True
+            The mean of the light curve can be estimated in each interval, or on
+            the full light curve. This gives different results (Alston+2013).
+            Here we assume the mean is calculated on the full light curve, but
+            the user can set ``use_common_mean`` to False to calculate it on a
+            per-segment basis.
+        silent : bool, default False
+            Silence the progress bars
+        gti: [[gti0_0, gti0_1], [gti1_0, gti1_1], ...]
+            Good Time intervals. Defaults to the common GTIs from the two input
+            objects. Could throw errors if these GTIs have overlaps with the
+            input object GTIs! If you're getting errors regarding your GTIs,
+            don't  use this and only give GTIs to the input objects before
+            making the cross spectrum.
+        """
+        return powerspectrum_from_timeseries(
+            ts,
+            flux_attr=flux_attr,
+            error_flux_attr=error_flux_attr,
+            segment_size=segment_size,
+            norm=norm,
+            silent=silent,
+            use_common_mean=use_common_mean,
+            gti=gti,
+        )
+
     @staticmethod
     def from_lightcurve(
         lc, segment_size=None, gti=None, norm="frac", silent=False, use_common_mean=True
@@ -1152,11 +1209,8 @@ def powerspectrum_from_lightcurve(
 
     Parameters
     ----------
-    events : `stingray.Lightcurve`
+    lc : `stingray.Lightcurve`
         Light curve to be analyzed.
-    dt : float
-        The time resolution of the intermediate light curves
-        (sets the Nyquist frequency)
 
     Other parameters
     ----------------
@@ -1215,6 +1269,86 @@ def powerspectrum_from_lightcurve(
     return _create_powerspectrum_from_result_table(table, force_averaged=force_averaged)
 
 
+def powerspectrum_from_timeseries(
+    ts,
+    flux_attr,
+    error_flux_attr=None,
+    segment_size=None,
+    norm="none",
+    silent=False,
+    use_common_mean=True,
+    gti=None,
+    save_all=False,
+):
+    """Calculate power spectrum from a time series
+
+    Parameters
+    ----------
+    ts : `stingray.StingrayTimeseries`
+        Input time series
+    flux_attr : `str`
+        What attribute of the time series will be used.
+
+    Other parameters
+    ----------------
+    error_flux_attr : `str`
+        What attribute of the time series will be used as error bar.
+    segment_size : float
+        The length, in seconds, of the light curve segments that will be
+        averaged. Only required (and used) for ``AveragedPowerspectrum``.
+    gti : ``[[gti0_0, gti0_1], [gti1_0, gti1_1], ...]``
+        Additional, optional Good Time intervals that get intersected with
+        the GTIs of the input object. Can cause errors if there are
+        overlaps between these GTIs and the input object GTIs. If that
+        happens, assign the desired GTIs to the input object.
+    norm : str, default "frac"
+        The normalization of the periodogram. `abs` is absolute rms, `frac`
+        is fractional rms, `leahy` is Leahy+83 normalization, and `none` is
+        the unnormalized periodogram.
+    use_common_mean : bool, default True
+        The mean of the light curve can be estimated in each interval, or
+        on the full light curve. This gives different results
+        (Alston+2013). By default, we assume the mean is calculated on the
+        full light curve, but the user can set ``use_common_mean`` to False
+        to calculate it on a per-segment basis.
+    silent : bool, default False
+        Silence the progress bars.
+    save_all : bool, default False
+        Save all intermediate PDSs used for the final average. Use with care.
+        This is likely to fill up your RAM on medium-sized datasets, and to
+        slow down the computation when rebinning.
+
+    Returns
+    -------
+    spec : `AveragedCrossspectrum` or `Crossspectrum`
+        The output cross spectrum.
+    """
+    force_averaged = segment_size is not None
+    # Suppress progress bar for single periodogram
+    silent = silent or (segment_size is None)
+    if gti is None:
+        gti = ts.gti
+
+    err = None
+    if error_flux_attr is not None:
+        err = getattr(ts, error_flux_attr)
+
+    results = avg_pds_from_events(
+        ts.time,
+        gti,
+        segment_size,
+        ts.dt,
+        norm=norm,
+        use_common_mean=use_common_mean,
+        silent=silent,
+        fluxes=getattr(ts, flux_attr),
+        errors=err,
+        return_subcs=save_all,
+    )
+
+    return _create_powerspectrum_from_result_table(results, force_averaged=force_averaged)
+
+
 def powerspectrum_from_lc_iterable(
     iter_lc,
     dt,
diff --git a/stingray/utils.py b/stingray/utils.py
index f9ba9a22..c24ad288 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -1,5 +1,6 @@
 import numbers
 import os
+import re
 import copy
 import random
 import string
@@ -13,7 +14,6 @@ import scipy
 from numpy import histogram as histogram_np
 from numpy import histogram2d as histogram2d_np
 from numpy import histogramdd as histogramdd_np
-from .base import interpret_times
 
 try:
     import pyfftw
@@ -140,6 +140,153 @@ __all__ = [
 ]
 
 
+@njit
+def any_complex_in_array(array):
+    """Check if any element of an array is complex.
+
+    Examples
+    --------
+    >>> any_complex_in_array(np.array([1, 2, 3]))
+    False
+    >>> any_complex_in_array(np.array([1, 2 + 1.j, 3]))
+    True
+    """
+    for a in array:
+        if np.iscomplex(a):
+            return True
+    return False
+
+
+def make_nd_into_arrays(array: np.ndarray, label: str) -> dict:
+    """If an array is n-dimensional, make it into many 1-dimensional arrays.
+
+    Call additional dimensions, e.g. ``_dimN_M``. See examples below.
+
+    Parameters
+    ----------
+    array : `np.ndarray`
+        Input data
+    label : `str`
+        Label for the array
+
+    Returns
+    -------
+    data : `dict`
+        Dictionary of arrays. Defauls to ``{label: array}`` if ``array`` is 1-dimensional,
+        otherwise, e.g.: ``{label_dim1_2_3: array[1, 2, 3], ... }``
+
+    Examples
+    --------
+    >>> a1, a2, a3 = np.arange(3), np.arange(3, 6), np.arange(6, 9)
+    >>> A = np.array([a1, a2, a3]).T
+    >>> data = make_nd_into_arrays(A, "test")
+    >>> np.array_equal(data["test_dim0"], a1)
+    True
+    >>> np.array_equal(data["test_dim1"], a2)
+    True
+    >>> np.array_equal(data["test_dim2"], a3)
+    True
+    >>> A3 = [[[1, 2], [3, 4]], [[5, 6], [7, 8]]]
+    >>> data = make_nd_into_arrays(A3, "test")
+    >>> np.array_equal(data["test_dim0_0"], [1, 5])
+    True
+    """
+    data = {}
+    array = np.asarray(array)
+    shape = np.shape(array)
+    ndim = len(shape)
+    if ndim <= 1:
+        data[label] = array
+    else:
+        for i in range(shape[1]):
+            new_label = f"_dim{i}" if "_dim" not in label else f"_{i}"
+            dumdata = make_nd_into_arrays(array[:, i], label=label + new_label)
+            data.update(dumdata)
+    return data
+
+
+def get_dimensions_from_list_of_column_labels(labels: list, label: str) -> list:
+    """Get the dimensions of a multi-dimensional array from a list of column labels.
+
+    Examples
+    --------
+    >>> labels = ['test_dim0_0', 'test_dim0_1', 'test_dim0_2',
+    ...           'test_dim1_0', 'test_dim1_1', 'test_dim1_2', 'test', 'bu']
+    >>> keys, dimensions = get_dimensions_from_list_of_column_labels(labels, "test")
+    >>> for key0, key1 in zip(labels[:6], keys): assert key0 == key1
+    >>> np.array_equal(dimensions, [2, 3])
+    True
+    """
+    all_keys = []
+    count_dimensions = None
+    for key in labels:
+        if label not in key:
+            continue
+        match = re.search("^" + label + r"_dim([0-9]+(_[0-9]+)*)", key)
+        if match is None:
+            continue
+        all_keys.append(key)
+        new_count_dimensions = [int(val) for val in match.groups()[0].split("_")]
+        if count_dimensions is None:
+            count_dimensions = np.array(new_count_dimensions)
+        else:
+            count_dimensions = np.max([count_dimensions, new_count_dimensions], axis=0)
+
+    return sorted(all_keys), count_dimensions + 1
+
+
+def make_1d_arrays_into_nd(data: dict, label: str) -> np.ndarray:
+    """Literally the opposite of make_nd_into_arrays.
+
+    Call additional dimensions, e.g. ``_dimN_M``
+
+    Parameters
+    ----------
+    data : dict
+        Input data
+    label : `str`
+        Label for the array
+
+    Returns
+    -------
+    array : `np.array`
+        N-dimensional array that was stored in the data.
+
+    Examples
+    --------
+    >>> a1, a2, a3 = np.arange(3), np.arange(3, 6), np.arange(6, 9)
+    >>> A = np.array([a1, a2, a3]).T
+    >>> data = make_nd_into_arrays(A, "test")
+    >>> A_ret = make_1d_arrays_into_nd(data, "test")
+    >>> np.array_equal(A, A_ret)
+    True
+    >>> A = np.array([[[1, 2, 12], [3, 4, 34]],
+    ...               [[5, 6, 56], [7, 8, 78]],
+    ...               [[9, 10, 910], [11, 12, 1112]],
+    ...               [[13, 14, 1314], [15, 16, 1516]]])
+    >>> data = make_nd_into_arrays(A, "_test")
+    >>> A_ret = make_1d_arrays_into_nd(data, "_test")
+    >>> np.array_equal(A, A_ret)
+    True
+    >>> data = make_nd_into_arrays(a1, "_test")
+    >>> A_ret = make_1d_arrays_into_nd(data, "_test")
+    >>> np.array_equal(a1, A_ret)
+    True
+    """
+
+    if label in list(data.keys()):
+        return data[label]
+
+    # Get the dimensionality of the data
+    dim = 0
+    all_keys = []
+
+    all_keys, dimensions = get_dimensions_from_list_of_column_labels(list(data.keys()), label)
+    arrays = np.array([np.array(data[key]) for key in all_keys])
+
+    return arrays.T.reshape([len(arrays[0])] + list(dimensions))
+
+
 @njit
 def _check_isallfinite_numba(array):
     """Check if all elements of an array are finite.
@@ -335,12 +482,12 @@ def rebin_data(x, y, dx_new, yerr=None, method="sum", dx=None):
     else:
         yerr = np.asarray(yerr)
 
-    if not dx:
+    if isinstance(dx, Iterable):
+        dx_old = dx
+    elif dx is None or dx == 0:
         dx_old = np.diff(x)
-    elif np.size(dx) == 1:
-        dx_old = np.array([dx])
     else:
-        dx_old = dx
+        dx_old = np.array([dx])
 
     if np.any(dx_new < dx_old):
         raise ValueError(
@@ -2129,3 +2276,24 @@ def assign_if_not_finite(value, default):
     if not np.isfinite(value):
         return default
     return value
+
+
+def sqsum(array1, array2):
+    """Return the square root of the sum of the squares of two arrays."""
+    return np.sqrt(np.add(np.square(array1), np.square(array2)))
+
+
+@njit
+def _int_sum_non_zero(array):
+    """Sum all positive elements of an array of integers.
+
+    Parameters
+    ----------
+    array : array-like
+        Array of integers
+    """
+    sum = 0
+    for a in array:
+        if a > 0:
+            sum += int(a)
+    return sum
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index 49ecaaed..8445396f 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -227,8 +227,8 @@ class VarEnergySpectrum(StingrayObject, metaclass=ABCMeta):
 
         self._create_empty_spectrum()
 
-        if len(events.time) == 0:
-            simon("There are no events in your event list!" + "Can't make a spectrum!")
+        if events.time is None or len(events.time) == 0:
+            simon("There are no events in your event list! Can't make a spectrum!")
         else:
             self._spectrum_function()
 


Test patch:
diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py
index 488f4cde..c3d10438 100644
--- a/stingray/tests/test_base.py
+++ b/stingray/tests/test_base.py
@@ -2,6 +2,7 @@ import os
 import copy
 import pytest
 import numpy as np
+import matplotlib.pyplot as plt
 from stingray.base import StingrayObject, StingrayTimeseries
 
 _HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = True
@@ -29,6 +30,7 @@ class DummyStingrayObj(StingrayObject):
 
     def __init__(self, dummy=None):
         self.guefus = dummy
+        self._mask = None
         # StingrayObject.__init__(self)
 
 
@@ -47,36 +49,22 @@ class BadStingrayObj(StingrayObject):
         StingrayObject.__init__(self)
 
 
-def _check_equal(so, new_so):
-    for attr in ["time", "guefus", "pardulas", "panesapa"]:
-        if not hasattr(so, attr):
-            assert not hasattr(new_so, attr)
-            continue
-        so_attr = at if (at := getattr(so, attr)) is not None else []
-        new_so_attr = at if (at := getattr(new_so, attr)) is not None else []
-
-        assert np.allclose(so_attr, new_so_attr)
-
-    for attr in ["mjdref", "pirichitus", "parafritus"]:
-        if not hasattr(so, attr):
-            assert not hasattr(new_so, attr)
-            continue
-        so_attr = getattr(so, attr)
-        new_so_attr = getattr(new_so, attr)
-        assert so_attr == new_so_attr
-
-
 class TestStingrayObject:
     @classmethod
     def setup_class(cls):
-        cls.arr = [4, 5, 2]
+        cls.arr = np.asarray([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]]
+        sting_obj._sebadas = [[0, 1], [2, 3], [4, 5]]
         sting_obj.pirichitus = 4
         sting_obj.parafritus = "bonus!"
         sting_obj.panesapa = [[41, 25], [98, 3]]
         cls.sting_obj = sting_obj
 
+    def test_print(self):
+        print(self.sting_obj)
+
     def test_preliminary(self):
         assert np.allclose(self.sting_obj.guefus, self.arr)
 
@@ -84,6 +72,199 @@ class TestStingrayObject:
         with pytest.raises(RuntimeError):
             BadStingrayObj(self.arr)
 
+    def test_equality(self):
+        ts1 = copy.deepcopy(self.sting_obj)
+        ts2 = copy.deepcopy(self.sting_obj)
+        assert ts1 == ts2
+
+    def test_different_array_attributes(self):
+        ts1 = copy.deepcopy(self.sting_obj)
+        ts2 = copy.deepcopy(self.sting_obj)
+        # Add a meta attribute only to ts1. This will fail
+        ts1.blah = 2
+        assert ts1 != ts2
+
+        # Add a non-scalar meta attribute, but with the same name, to ts2.
+        ts2.blah = [2]
+        assert ts1 != ts2
+
+        # Get back to normal
+        del ts1.blah, ts2.blah
+        assert ts1 == ts2
+
+        # Add a non-scalar meta attribute to both, just slightly different
+        ts1.blah = [2]
+        ts2.blah = [3]
+        assert ts1 != ts2
+
+        # Get back to normal
+        del ts1.blah, ts2.blah
+        assert ts1 == ts2
+
+        # Add a meta attribute only to ts2. This will also fail
+        ts2.blah = 3
+        assert ts1 != ts2
+
+    def test_different_meta_attributes(self):
+        ts1 = copy.deepcopy(self.sting_obj)
+        ts2 = copy.deepcopy(self.sting_obj)
+        # Add an array attribute to ts1. This will fail
+        ts1.blah = ts1.guefus
+        assert ts1 != ts2
+
+        # Get back to normal
+        del ts1.blah
+        assert ts1 == ts2
+        # Add an array attribute to ts2. This will fail
+        ts2.blah = ts1.guefus
+        assert ts1 != ts2
+
+        # Get back to normal
+        del ts2.blah
+        assert ts1 == ts2
+
+    @pytest.mark.parametrize("inplace", [True, False])
+    def test_apply_mask(self, inplace):
+        ts = copy.deepcopy(self.sting_obj)
+        obj = ts.apply_mask([True, True, False], inplace=inplace)
+
+        assert obj.parafritus == "bonus!"
+        assert np.array_equal(obj.guefus, [4, 5])
+        assert np.array_equal(obj.panesapa, ts.panesapa)
+        assert np.array_equal(obj.pardulas, [3.0 + 1.0j, 2.0j])
+        assert np.array_equal(obj.sebadas, [[0, 1], [2, 3]])
+        if inplace:
+            # Only if masking in place, the final object will be the same as the starting one.
+            assert ts is obj
+        else:
+            # If not, the objects have to be different
+            assert ts is not obj
+
+    @pytest.mark.parametrize("inplace", [True, False])
+    def test_partial_apply_mask(self, inplace):
+        ts = copy.deepcopy(self.sting_obj)
+        obj = ts.apply_mask([True, True, False], inplace=inplace, filtered_attrs=["pardulas"])
+        assert obj.parafritus == "bonus!"
+        assert np.array_equal(obj.guefus, [4, 5])
+        assert np.array_equal(obj.panesapa, ts.panesapa)
+        assert np.array_equal(obj.pardulas, [3.0 + 1.0j, 2.0j])
+        assert obj.sebadas is None
+
+        if inplace:
+            assert ts is obj
+        else:
+            assert ts is not obj
+
+    def test_operations(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(guefus)
+        ts1.counts = count1
+        ts2.counts = count2
+        lc = ts1 + ts2  # Test __add__
+        assert np.allclose(lc.counts, [900, 1300, 1200])
+        assert np.array_equal(lc.guefus, guefus)
+        lc = ts1 - ts2  # Test __sub__
+        assert np.allclose(lc.counts, [-300, -1100, -400])
+        assert np.array_equal(lc.guefus, guefus)
+        lc = -ts2 + ts1  # Test __neg__
+        assert np.allclose(lc.counts, [-300, -1100, -400])
+        assert np.array_equal(lc.guefus, guefus)
+
+    def test_inplace_add(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(guefus)
+        ts1.counts = count1
+        ts2.counts = count2
+        lc = ts1 + ts2  # Test __add__
+        ts1 += ts2  # Test __iadd__
+        assert np.allclose(ts1.counts, [900, 1300, 1200])
+        assert np.array_equal(ts1.guefus, guefus)
+        assert lc == ts1
+
+    def test_inplace_sub(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(guefus)
+        ts1.counts = count1
+        ts2.counts = count2
+        lc = ts1 - ts2  # Test __sub__
+        ts1 -= ts2  # Test __isub__
+        assert np.allclose(ts1.counts, [-300, -1100, -400])
+        assert np.array_equal(ts1.guefus, guefus)
+        assert lc == ts1
+
+    def test_inplace_add_with_method(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(guefus)
+        ts1.counts = count1
+        ts2.counts = count2
+        lc = ts1.add(ts2)
+        assert lc is not ts1
+        lc_ip = ts2.add(ts1, inplace=True)
+        assert lc == lc_ip
+        assert lc_ip is ts2
+        assert np.allclose(lc.counts, [900, 1300, 1200])
+        assert np.array_equal(lc.guefus, guefus)
+
+    def test_inplace_sub_with_method(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(guefus)
+        ts1.counts = count1
+        ts2.counts = count2
+        lc = ts1.sub(ts2)
+        assert lc is not ts1
+        lc_ip = ts2.sub(ts1, inplace=True)
+        assert lc == -lc_ip
+        assert lc_ip is ts2
+        assert np.allclose(lc.counts, [-300, -1100, -400])
+        assert np.array_equal(lc.guefus, guefus)
+
+    def test_failed_operations(self):
+        guefus = [5, 10, 15]
+        count1 = [300, 100, 400]
+        ts1 = DummyStingrayObj(guefus)
+        ts2 = DummyStingrayObj(np.array(guefus) + 1)
+        ts1.counts = np.array(count1)
+        ts2.counts = np.array(count1)
+        with pytest.raises(TypeError, match=".*objects can only be operated with other.*"):
+            ts1._operation_with_other_obj(float(3), np.add)
+        with pytest.raises(ValueError, match="The values of guefus are different"):
+            ts1 + ts2
+
+    def test_len(self):
+        assert len(self.sting_obj) == 3
+
+    def test_slice(self):
+        ts1 = self.sting_obj
+        ts_filt = ts1[1]
+        assert np.array_equal(ts_filt.guefus, [5])
+        assert ts_filt.parafritus == "bonus!"
+        assert np.array_equal(ts_filt.panesapa, ts1.panesapa)
+        assert np.array_equal(ts_filt.pardulas, [2.0j])
+
+        ts_filt = ts1[:2]
+        assert np.array_equal(ts_filt.guefus, [4, 5])
+        assert ts_filt.parafritus == "bonus!"
+        assert np.array_equal(ts_filt.panesapa, ts1.panesapa)
+        assert np.array_equal(ts_filt.pardulas, [3.0 + 1.0j, 2.0j])
+
+        with pytest.raises(IndexError, match="The index must be either an integer or a slice"):
+            ts1[1.0]
+
     def test_side_effects(self):
         so = copy.deepcopy(self.sting_obj)
         assert np.allclose(so.guefus, [4, 5, 2])
@@ -100,7 +281,7 @@ class TestStingrayObject:
         so.stingattr = DummyStingrayObj([3, 4, 5])
         ts = so.to_astropy_table()
         new_so = DummyStingrayObj.from_astropy_table(ts)
-        _check_equal(so, new_so)
+        assert so == new_so
         assert not hasattr(new_so, "stingattr")
 
     @pytest.mark.skipif("not _HAS_XARRAY")
@@ -110,8 +291,7 @@ class TestStingrayObject:
         so.guefus = np.random.randint(0, 4, 3)
         ts = so.to_xarray()
         new_so = DummyStingrayObj.from_xarray(ts)
-
-        _check_equal(so, new_so)
+        assert so == new_so
 
     @pytest.mark.skipif("not _HAS_PANDAS")
     def test_pandas_roundtrip(self):
@@ -120,8 +300,9 @@ class TestStingrayObject:
         so.guefus = np.random.randint(0, 4, 3)
         ts = so.to_pandas()
         new_so = DummyStingrayObj.from_pandas(ts)
-
-        _check_equal(so, new_so)
+        # assert not hasattr(new_so, "sebadas")
+        # new_so.sebadas = so.sebadas
+        assert so == new_so
 
     def test_astropy_roundtrip_empty(self):
         # Set an attribute to a DummyStingrayObj. It will *not* be saved
@@ -153,22 +334,37 @@ class TestStingrayObject:
         new_so = so.read("dummy.hdf5")
         os.unlink("dummy.hdf5")
 
-        _check_equal(so, new_so)
+        assert so == new_so
 
     def test_file_roundtrip_fits(self):
         so = copy.deepcopy(self.sting_obj)
         so.guefus = np.random.randint(0, 4, 3)
         so.panesapa = np.random.randint(5, 9, (6, 2))
-        so.write("dummy.fits")
+        with pytest.warns(
+            UserWarning, match=".* output does not serialize the metadata at the moment"
+        ):
+            so.write("dummy.fits")
         new_so = DummyStingrayObj.read("dummy.fits")
         os.unlink("dummy.fits")
         # panesapa is invalid for FITS header and got lost
         assert not hasattr(new_so, "panesapa")
         new_so.panesapa = so.panesapa
-        _check_equal(so, new_so)
+        assert so == new_so
 
-    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv"])
+    @pytest.mark.parametrize("fmt", ["ascii", "ascii.ecsv"])
     def test_file_roundtrip(self, fmt):
+        so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = DummyStingrayObj.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        assert so == new_so
+
+    def test_file_roundtrip_pickle(self):
+        fmt = "pickle"
         so = copy.deepcopy(self.sting_obj)
         so.guefus = np.random.randint(0, 4, 3)
         so.panesapa = np.random.randint(5, 9, (6, 2))
@@ -176,10 +372,665 @@ class TestStingrayObject:
         new_so = DummyStingrayObj.read(f"dummy.{fmt}", fmt=fmt)
         os.unlink(f"dummy.{fmt}")
 
-        _check_equal(so, new_so)
+        assert so == new_so
 
 
 class TestStingrayTimeseries:
+    @classmethod
+    def setup_class(cls):
+        cls.time = np.arange(0, 10, 1)
+        cls.arr = cls.time + 2
+        sting_obj = 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]]),
+        )
+        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]]),
+            high_precision=True,
+        )
+        cls.sting_obj = sting_obj
+        cls.sting_obj_highp = sting_obj_highp
+
+    def test_print(self):
+        print(self.sting_obj)
+
+    def test_invalid_instantiation(self):
+        with pytest.raises(ValueError, match="Lengths of time and guefus must be equal"):
+            StingrayTimeseries(time=np.arange(10), array_attrs=dict(guefus=np.arange(11)))
+        with pytest.raises(ValueError, match="Lengths of time and guefus must be equal"):
+            StingrayTimeseries(time=np.arange(10), array_attrs=dict(guefus=np.zeros((5, 2))))
+
+    def test_apply_mask(self):
+        ts = copy.deepcopy(self.sting_obj)
+        mask = [True, True] + 8 * [False]
+        newts0 = ts.apply_mask(mask, inplace=False)
+        newts1 = ts.apply_mask(mask, inplace=True)
+        for obj in [newts1, newts0]:
+            for attr in ["parafritus", "mjdref"]:
+                assert getattr(obj, attr) == getattr(ts, attr)
+            for attr in ["panesapa", "gti"]:
+                assert np.array_equal(getattr(obj, attr), getattr(ts, attr))
+
+            assert np.array_equal(obj.guefus, [2, 3])
+            assert np.array_equal(obj.time, [0, 1])
+        assert ts is newts1
+        assert ts is not newts0
+
+    def test_comparison(self):
+        time = [5, 10, 15]
+        count1 = [300, 100, 400]
+
+        ts1 = StingrayTimeseries(
+            time=time,
+            array_attrs=dict(counts=np.array(count1), _counts=np.array(count1)),
+            mjdref=55000,
+        )
+        ts2 = StingrayTimeseries(
+            time=time,
+            array_attrs=dict(counts=np.array(count1), _counts=np.array(count1)),
+            mjdref=55000,
+        )
+
+        assert ts1 == ts2
+        # Change one attribute, check that they are not equal
+        ts2.counts[0] += 1
+        assert ts1 != ts2
+        # Put back, check there are no side effects
+        ts2.counts = np.array(count1)
+        assert ts1 == ts2
+        # Now check a meta attribute
+        ts2.mjdref += 1
+        assert ts1 != ts2
+        # Put back, check there are no side effects
+        ts2.mjdref = 55000
+        assert ts1 == ts2
+        # Now check an internal array attribute
+        ts2._counts[0] += 1
+        assert ts1 != ts2
+        # Put back, check there are no side effects
+        ts2._counts = np.array(count1)
+        assert ts1 == ts2
+
+    def test_zero_out_timeseries(self):
+        time = [5, 10, 15]
+        count1 = [300, 100, 400]
+
+        ts1 = StingrayTimeseries(
+            time=time,
+            array_attrs=dict(counts=np.array(count1), _bla=np.array(count1)),
+            mjdref=55000,
+        )
+        # All has been set correctly
+        assert np.array_equal(ts1.counts, count1)
+        assert np.array_equal(ts1._bla, count1)
+
+        # Now zero out times and verify that everything else has been zeroed out
+        ts1.time = None
+        assert ts1.counts is None
+        assert ts1.time is None
+        assert ts1._bla is None
+
+    def test_n_property(self):
+        ts = StingrayTimeseries()
+        assert ts.n == 0
+
+        time = [5, 10, 15]
+        count1 = [300, 100, 400]
+
+        ts1 = StingrayTimeseries(
+            time=time,
+            array_attrs=dict(counts=np.array(count1)),
+            mjdref=55000,
+        )
+        # All has been set correctly
+        assert ts1.n == 3
+
+        ts1.time = None
+        assert ts1.n == 0
+
+    def test_what_is_array_and_what_is_not(self):
+        """Test that array_attrs are not confused with other attributes.
+
+        In particular, time, gti and panesapa have the same length. Verify that panesapa
+        is considered an array attribute, but not gti."""
+        ts = StingrayTimeseries(
+            [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]])
+            ),
+            dt=1,
+        )
+        array_attrs = ts.array_attrs()
+        internal_array_attrs = ts.internal_array_attrs()
+        assert "panesapa" in array_attrs
+        assert "_panesapa" not in array_attrs
+        assert "_panesapa" in internal_array_attrs
+        assert "gti" not in array_attrs
+        assert "time" not in array_attrs
+
+    def test_operations(self):
+        time = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+
+        ts1 = StingrayTimeseries(time=time)
+        ts2 = StingrayTimeseries(time=time)
+        ts1.counts = count1
+        ts2.counts = count2
+        ts1.counts_err = np.zeros_like(count1) + 1
+        ts2.counts_err = np.zeros_like(count2) + 1
+
+        lc = ts1 + ts2  # Test __add__
+        assert np.allclose(lc.counts, [900, 1300, 1200])
+        assert np.array_equal(lc.time, time)
+        assert np.allclose(lc.counts_err, np.sqrt(2))
+        lc = ts1 - ts2  # Test __sub__
+        assert np.allclose(lc.counts, [-300, -1100, -400])
+        assert np.array_equal(lc.time, time)
+        assert np.allclose(lc.counts_err, np.sqrt(2))
+        lc = -ts2 + ts1  # Test __neg__
+        assert np.allclose(lc.counts, [-300, -1100, -400])
+        assert np.array_equal(lc.time, time)
+        assert np.allclose(lc.counts_err, np.sqrt(2))
+
+    def test_operations_different_mjdref(self):
+        time = [5, 10, 15]
+        count1 = [300, 100, 400]
+        count2 = [600, 1200, 800]
+
+        ts1 = StingrayTimeseries(time=time, array_attrs=dict(counts=count1), mjdref=55000)
+        # Now I create a second time series, I make sure that the times are the same, but
+        # Then I change the mjdref. From now on, the time array will not have the same
+        # values as `ts1.time`. The sum operation will warn the user about this, but then
+        # change the mjdref of the second time series to match the first one, and the times
+        # will be aligned again.
+        ts2 = StingrayTimeseries(time=time, array_attrs=dict(counts=count2), mjdref=55000)
+        ts2.change_mjdref(54000, inplace=True)
+        with pytest.warns(UserWarning, match="MJDref is different in the two time series"):
+            lc = ts1 + ts2  # Test __add__
+        assert np.allclose(lc.counts, [900, 1300, 1200])
+        assert np.array_equal(lc.time, time)
+        assert np.array_equal(lc.mjdref, ts1.mjdref)
+        with pytest.warns(UserWarning, match="MJDref is different in the two time series"):
+            lc = ts2 + ts1  # Test __add__ of the other curve. The mjdref will be the other one now
+        assert np.allclose(lc.counts, [900, 1300, 1200])
+        assert np.array_equal(lc.time, ts2.time)
+        assert np.array_equal(lc.mjdref, ts2.mjdref)
+
+    def test_operation_with_diff_gti(self):
+        time = [10, 20, 30, 40]
+        count1 = [600, 1200, 800, 400]
+        count2 = [300, 100, 400, 100]
+        gti1 = [[0, 35]]
+        gti2 = [[5, 40]]
+        ts1 = StingrayTimeseries(time, array_attrs=dict(counts=count1), gti=gti1, dt=10)
+        ts2 = StingrayTimeseries(time, array_attrs=dict(counts=count2), gti=gti2, dt=10)
+        with pytest.warns(
+            UserWarning, match="The good time intervals in the two time series are different."
+        ):
+            lc = ts1 - ts2
+        assert np.allclose(lc.counts, [300, 1100, 400])
+        assert np.allclose(lc.time, [10, 20, 30])
+
+    def test_len(self):
+        assert len(self.sting_obj) == 10
+
+    def test_slice(self):
+        ts1 = copy.deepcopy(self.sting_obj)
+
+        ts_filt = ts1[1]
+        assert np.array_equal(ts_filt.guefus, [3])
+        assert ts_filt.parafritus == "bonus!"
+        assert np.array_equal(ts_filt.panesapa, ts1.panesapa)
+
+        ts_filt = ts1[:2]
+        assert np.array_equal(ts_filt.guefus, [2, 3])
+        assert ts_filt.parafritus == "bonus!"
+        assert np.array_equal(ts_filt.panesapa, ts1.panesapa)
+
+        ts_filt = ts1[0:3:2]
+        # If dt >0, gtis are also altered. Otherwise they're left alone
+        ts1.dt = 1
+        ts_filt_dt1 = ts1[0:3:2]
+        # Also try with array dt
+        ts1.dt = np.ones_like(ts1.time)
+        ts_filt_dtarr = ts1[0:3:2]
+        for ts_f in [ts_filt, ts_filt_dt1, ts_filt_dtarr]:
+            assert np.array_equal(ts_f.guefus, [2, 4])
+            assert ts_f.parafritus == "bonus!"
+            assert np.array_equal(ts_f.panesapa, ts1.panesapa)
+        assert np.allclose(ts_filt_dt1.gti, [[-0.5, 0.5], [1.5, 2.5]])
+        assert np.allclose(ts_filt_dtarr.gti, [[-0.5, 0.5], [1.5, 2.5]])
+        assert np.allclose(ts_filt.gti, [[0, 2.0]])
+
+        with pytest.raises(IndexError, match="The index must be either an integer or a slice"):
+            ts1[1.0]
+
+    @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]])
+        so2 = so.apply_gtis()
+        if inplace:
+            assert so2 is so
+
+        assert np.allclose(so2.time, [0, 1, 2])
+        assert np.allclose(so2.guefus, [2, 3, 4])
+        assert np.allclose(so2.gti, [[-0.1, 2.1]])
+        assert np.allclose(so2.mjdref, 59777.000)
+
+    def test_split_ts_by_gtis(self):
+        times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+        counts = [1, 1, 1, 1, 2, 3, 3, 2, 3, 3]
+        bg_counts = [0, 0, 0, 1, 0, 1, 2, 0, 0, 1]
+        bg_ratio = [0.1, 0.1, 0.1, 0.2, 0.1, 0.2, 0.2, 0.1, 0.1, 0.2]
+        frac_exp = [1, 0.5, 1, 1, 1, 0.5, 0.5, 1, 1, 1]
+        gti = [[0.5, 4.5], [5.5, 7.5], [8.5, 9.5]]
+
+        ts = StingrayTimeseries(
+            times,
+            array_attrs=dict(
+                counts=counts, bg_counts=bg_counts, bg_ratio=bg_ratio, frac_exp=frac_exp
+            ),
+            gti=gti,
+        )
+        list_of_tss = ts.split_by_gti(min_points=0)
+        assert len(list_of_tss) == 3
+
+        ts0 = list_of_tss[0]
+        ts1 = list_of_tss[1]
+        ts2 = list_of_tss[2]
+        assert np.allclose(ts0.time, [1, 2, 3, 4])
+        assert np.allclose(ts1.time, [6, 7])
+        assert np.allclose(ts2.time, [9])
+        assert np.allclose(ts0.counts, [1, 1, 1, 1])
+        assert np.allclose(ts1.counts, [3, 3])
+        assert np.allclose(ts1.counts, [3])
+        assert np.allclose(ts0.gti, [[0.5, 4.5]])
+        assert np.allclose(ts1.gti, [[5.5, 7.5]])
+        assert np.allclose(ts2.gti, [[8.5, 9.5]])
+        # Check if new attributes are also splited accordingly
+        assert np.allclose(ts0.bg_counts, [0, 0, 0, 1])
+        assert np.allclose(ts1.bg_counts, [1, 2])
+        assert np.allclose(ts0.bg_ratio, [0.1, 0.1, 0.1, 0.2])
+        assert np.allclose(ts1.bg_ratio, [0.2, 0.2])
+        assert np.allclose(ts0.frac_exp, [1, 0.5, 1, 1])
+        assert np.allclose(ts1.frac_exp, [0.5, 0.5])
+
+    def test_truncate(self):
+        time = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        count = [10, 20, 30, 40, 50, 60, 70, 80, 90]
+        lc = StingrayTimeseries(time, array_attrs={"counts": count}, dt=1)
+        lc_new = lc.truncate(start=2, stop=8, method="index")
+        assert np.allclose(lc_new.counts, [30, 40, 50, 60, 70, 80])
+        assert np.array_equal(lc_new.time, [3, 4, 5, 6, 7, 8])
+
+        # Truncation can also be done by time values
+        lc_new = lc.truncate(start=6, method="time")
+        assert np.array_equal(lc_new.time, [6, 7, 8, 9])
+        assert np.allclose(lc_new.counts, [60, 70, 80, 90])
+
+    def test_truncate_not_str(self):
+        with pytest.raises(TypeError, match="The method keyword argument"):
+            self.sting_obj.truncate(method=1)
+
+    def test_truncate_invalid(self):
+        with pytest.raises(ValueError, match="Unknown method type"):
+            self.sting_obj.truncate(method="ababalksdfja")
+
+    def test_concatenate(self):
+        time0 = [1, 2, 3]
+        time1 = [5, 6, 7, 8]
+        time2 = [10]
+        count0 = [10, 20, 30]
+        count1 = [50, 60, 70, 80]
+        count2 = [100]
+        gti0 = [[0.5, 3.5]]
+        gti1 = [[4.5, 8.5]]
+        gti2 = [[9.5, 10.5]]
+        lc0 = StingrayTimeseries(
+            time0, array_attrs={"counts": count0, "_bla": count0}, dt=1, gti=gti0
+        )
+        lc1 = StingrayTimeseries(
+            time1, array_attrs={"counts": count1, "_bla": count1}, dt=1, gti=gti1
+        )
+        lc2 = StingrayTimeseries(
+            time2, array_attrs={"counts": count2, "_bla": count2}, dt=1, gti=gti2
+        )
+        lc = lc0.concatenate([lc1, lc2])
+        assert np.allclose(lc._bla, count0 + count1 + count2)
+        assert np.allclose(lc.counts, count0 + count1 + count2)
+        assert np.allclose(lc.time, time0 + time1 + time2)
+        assert np.allclose(lc.gti, [[0.5, 3.5], [4.5, 8.5], [9.5, 10.5]])
+
+    def test_concatenate_invalid(self):
+        with pytest.raises(TypeError, match="objects can only be merged with other"):
+            self.sting_obj.concatenate(1)
+
+    def test_concatenate_gtis_overlap(self):
+        time0 = [1, 2, 3, 4]
+        time1 = [5, 6, 7, 8, 9]
+        count0 = [10, 20, 30, 40]
+        count1 = [50, 60, 70, 80, 90]
+        gti0 = [[0.5, 4.5]]
+        gti1 = [[3.5, 9.5]]
+        lc0 = StingrayTimeseries(
+            time0, array_attrs={"counts": count0, "_bla": count0}, dt=1, gti=gti0
+        )
+        lc1 = StingrayTimeseries(
+            time1, array_attrs={"counts": count1, "_bla": count1}, dt=1, gti=gti1
+        )
+        with pytest.raises(ValueError, match="In order to append, GTIs must be mutually"):
+            lc0.concatenate(lc1)
+
+        # Instead, this will work
+        lc0.concatenate(lc1, check_gti=False)
+
+    def test_concatenate_diff_mjdref(self):
+        time0 = [1, 2, 3, 4]
+        time1 = [5, 6, 7, 8, 9]
+        count0 = [10, 20, 30, 40]
+        count1 = [50, 60, 70, 80, 90]
+        gti0 = [[0.5, 4.49]]
+        gti1 = [[4.51, 9.5]]
+        lc0 = StingrayTimeseries(
+            time0, array_attrs={"counts": count0, "_bla": count0}, dt=1, gti=gti0, mjdref=55000
+        )
+        lc1 = StingrayTimeseries(
+            time1, array_attrs={"counts": count1, "_bla": count1}, dt=1, gti=gti1, mjdref=55000
+        )
+        lc1.change_mjdref(50001, inplace=True)
+        with pytest.warns(UserWarning, match="mjdref is different"):
+            lc = lc0.concatenate(lc1)
+        assert lc.mjdref == 55000
+
+    def test_rebin(self):
+        time0 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        count0 = [10, 20, 30, 40, 50, 60, 70, 80, 90]
+        count0_err = [1] * 9
+        gti0 = [[0.5, 9.5]]
+        lc0 = StingrayTimeseries(
+            time0,
+            array_attrs={"counts": count0, "counts_err": count0_err, "_bla": count0},
+            dt=1,
+            gti=gti0,
+        )
+        # With new dt=2
+        lc1 = lc0.rebin(dt_new=2)
+        assert np.allclose(lc1.counts, [30, 70, 110, 150])
+        assert np.allclose(lc1.counts_err, [np.sqrt(2)] * 4)
+        assert np.allclose(lc1._bla, [30, 70, 110, 150])
+        assert np.allclose(lc1.time, [1.5, 3.5, 5.5, 7.5])
+        assert lc1.dt == 2
+        # With a factor of two. Should give the same result
+        lc1 = lc0.rebin(f=2)
+        assert np.allclose(lc1.counts, [30, 70, 110, 150])
+        assert np.allclose(lc1.counts_err, [np.sqrt(2)] * 4)
+        assert np.allclose(lc1._bla, [30, 70, 110, 150])
+        assert np.allclose(lc1.time, [1.5, 3.5, 5.5, 7.5])
+        assert lc1.dt == 2
+
+    def test_rebin_irregular(self):
+        x0 = [-10]
+        x1 = np.linspace(0, 10, 11)
+        x2 = np.linspace(10.33, 20.0, 30)
+        x3 = np.linspace(21, 30, 10)
+        x = np.hstack([x0, x1, x2, x3])
+        dt = np.hstack(
+            [
+                [1],
+                [np.diff(x1).mean()] * x1.size,
+                [np.diff(x2).mean()] * x3.size,
+                [np.diff(x3).mean()] * x3.size,
+            ]
+        )
+
+        counts = 2.0
+        y = np.zeros_like(x) + counts
+
+        yerr = np.sqrt(y)
+        # Note: the first point of x is outside the GTIs
+        lc0 = StingrayTimeseries(
+            x,
+            array_attrs={"counts": y, "counts_err": yerr, "_bla": y},
+            dt=dt,
+            gti=[[x1[0] - 0.5 * dt[1], x3[-1] + 0.5 * dt[-1]]],
+        )
+
+        dx_new = 1.5
+        lc1 = lc0.rebin(dt_new=dx_new)
+        from stingray import utils
+
+        # Verifiy that the rebinning of irregular data is sampled correctly,
+        # Including all data but the first point, which is outsid GTIs
+        xbin, ybin, yerr_bin, step_size = utils.rebin_data(x[1:], y[1:], dx_new, yerr[1:])
+
+        assert np.allclose(lc1.time, xbin)
+        assert np.allclose(lc1.counts, ybin)
+        assert np.allclose(lc1._bla, ybin)
+        assert np.allclose(lc1.counts_err, yerr_bin)
+
+    def test_rebin_no_good_gtis(self):
+        time0 = [1, 2, 3, 4]
+        count0 = [10, 20, 30, 40]
+        gti0 = [[0.5, 4.5]]
+        lc0 = StingrayTimeseries(
+            time0,
+            array_attrs={"counts": count0},
+            dt=1,
+            gti=gti0,
+        )
+        with pytest.raises(ValueError, match="No valid GTIs after rebin."):
+            print(lc0.rebin(dt_new=5).counts)
+
+    def test_rebin_no_input(self):
+        with pytest.raises(ValueError, match="You need to specify at least one between f and"):
+            self.sting_obj.rebin()
+
+    def test_rebin_less_than_dt(self):
+        time0 = [1, 2, 3, 4]
+        count0 = [10, 20, 30, 40]
+        lc0 = StingrayTimeseries(time0, array_attrs={"counts": count0}, dt=1)
+        with pytest.raises(ValueError, match="The new time resolution must be larger than"):
+            lc0.rebin(dt_new=0.1)
+
+    def test_sort(self):
+        times = [2, 1, 3, 4]
+        blah = np.asarray([40, 10, 20, 5])
+        bleh = [4, 1, 2, 0.5]
+        mjdref = 57000
+
+        lc = StingrayTimeseries(
+            times, array_attrs={"blah": blah, "_bleh": bleh}, dt=1, mjdref=mjdref
+        )
+
+        lc_new = lc.sort()
+
+        assert np.allclose(lc_new._bleh, np.array([1, 4, 2, 0.5]))
+        assert np.allclose(lc_new.blah, np.array([10, 40, 20, 5]))
+        assert np.allclose(lc_new.time, np.array([1, 2, 3, 4]))
+        assert lc_new.mjdref == mjdref
+
+        lc_new = lc.sort(reverse=True)
+
+        assert np.allclose(lc_new._bleh, np.array([0.5, 2, 4, 1]))
+        assert np.allclose(lc_new.blah, np.array([5, 20, 40, 10]))
+        assert np.allclose(lc_new.time, np.array([4, 3, 2, 1]))
+        assert lc_new.mjdref == mjdref
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_astropy_roundtrip(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        ts = so.to_astropy_table()
+        new_so = StingrayTimeseries.from_astropy_table(ts)
+        assert so == new_so
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_astropy_ts_roundtrip(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        ts = so.to_astropy_timeseries()
+        new_so = StingrayTimeseries.from_astropy_timeseries(ts)
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_XARRAY")
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_xarray_roundtrip(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        ts = so.to_xarray()
+        new_so = StingrayTimeseries.from_xarray(ts)
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_PANDAS")
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_pandas_roundtrip(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        ts = so.to_pandas()
+        new_so = StingrayTimeseries.from_pandas(ts)
+        # assert not hasattr(new_so, "sebadas")
+        # new_so.sebadas = so.sebadas
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_H5PY")
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_hdf_roundtrip(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.write("dummy.hdf5")
+        new_so = so.read("dummy.hdf5")
+        os.unlink("dummy.hdf5")
+
+        assert so == new_so
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_file_roundtrip_fits(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, self.time.shape)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        with pytest.warns(
+            UserWarning, match=".* output does not serialize the metadata at the moment"
+        ):
+            so.write("dummy.fits")
+        new_so = StingrayTimeseries.read("dummy.fits")
+        os.unlink("dummy.fits")
+        # panesapa is invalid for FITS header and got lost
+        assert not hasattr(new_so, "panesapa")
+        new_so.panesapa = so.panesapa
+        new_so.gti = so.gti
+        new_so == so
+
+    @pytest.mark.parametrize("fmt", ["ascii", "ascii.ecsv"])
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_file_roundtrip(self, fmt, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        with pytest.warns(
+            UserWarning, match=f".* output does not serialize the metadata at the moment"
+        ):
+            so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = StingrayTimeseries.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        assert so == new_so
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_file_roundtrip_pickle(self, highprec):
+        fmt = "pickle"
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        so.guefus = np.random.randint(0, 4, 3)
+        so.panesapa = np.random.randint(5, 9, (6, 2))
+        so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = StingrayTimeseries.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        assert so == new_so
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_shift_time(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        new_so = so.shift(1)
+        assert np.allclose(new_so.time - 1, self.sting_obj.time)
+        assert np.allclose(new_so.gti - 1, self.sting_obj.gti)
+
+    @pytest.mark.parametrize("highprec", [True, False])
+    def test_change_mjdref(self, highprec):
+        if highprec:
+            so = copy.deepcopy(self.sting_obj_highp)
+        else:
+            so = copy.deepcopy(self.sting_obj)
+        new_so = so.change_mjdref(59776.5)
+        assert new_so.mjdref == 59776.5
+        assert np.allclose(new_so.time - 43200, self.sting_obj.time)
+        assert np.allclose(new_so.gti - 43200, self.sting_obj.gti)
+
+    def test_plot_simple(self):
+        time0 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+        count0 = [10, 20, 30, 40, 50, 60, 70, 80, 90]
+        count0_err = [1] * 9
+        gti0 = [[0.5, 3.5], [4.5, 9.5]]
+        lc0 = StingrayTimeseries(
+            time0,
+            array_attrs={"counts": count0, "counts_err": count0_err, "_bla": count0},
+            dt=1,
+            gti=gti0,
+        )
+        plt.close("all")
+        lc0.plot("counts", title="Counts", witherrors=True)
+        assert plt.fignum_exists("counts")
+        plt.close("all")
+
+    def test_plot_default_filename(self):
+        self.sting_obj.plot("guefus", save=True)
+        assert os.path.isfile("out.png")
+        os.unlink("out.png")
+        plt.close("all")
+
+    def test_plot_custom_filename(self):
+        self.sting_obj.plot("guefus", save=True, filename="lc.png")
+        assert os.path.isfile("lc.png")
+        os.unlink("lc.png")
+        plt.close("all")
+
+
+class TestStingrayTimeseriesSubclass:
     @classmethod
     def setup_class(cls):
         cls.arr = [4, 5, 2]
@@ -191,18 +1042,21 @@ class TestStingrayTimeseries:
         sting_obj.gti = np.asarray([[-0.5, 2.5]])
         cls.sting_obj = sting_obj
 
+    def test_print(self):
+        print(self.sting_obj)
+
     def test_astropy_roundtrip(self):
         so = copy.deepcopy(self.sting_obj)
         # Set an attribute to a DummyStingrayObj. It will *not* be saved
         ts = so.to_astropy_table()
         new_so = DummyStingrayTs.from_astropy_table(ts)
-        _check_equal(so, new_so)
+        assert so == new_so
 
     def test_astropy_ts_roundtrip(self):
         so = copy.deepcopy(self.sting_obj)
         ts = so.to_astropy_timeseries()
         new_so = DummyStingrayTs.from_astropy_timeseries(ts)
-        _check_equal(so, new_so)
+        assert so == new_so
 
     def test_shift_time(self):
         new_so = self.sting_obj.shift(1)
@@ -214,3 +1068,220 @@ class TestStingrayTimeseries:
         assert new_so.mjdref == 59776.5
         assert np.allclose(new_so.time - 43200, self.sting_obj.time)
         assert np.allclose(new_so.gti - 43200, self.sting_obj.gti)
+
+
+class TestJoinEvents:
+    def test_join_without_times_simulated(self):
+        """Test if exception is raised when join method is
+        called before first simulating times.
+        """
+        ts = StingrayTimeseries()
+        ts_other = StingrayTimeseries()
+
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            assert ts.join(ts_other, strategy="union").time is None
+
+    def test_join_empty_lists(self):
+        """Test if an empty event list can be concatenated
+        with a non-empty event list.
+        """
+        ts = StingrayTimeseries(time=[1, 2, 3])
+        ts_other = StingrayTimeseries()
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ts_new = ts.join(ts_other, strategy="union")
+        assert np.allclose(ts_new.time, [1, 2, 3])
+
+        ts = StingrayTimeseries()
+        ts_other = StingrayTimeseries(time=[1, 2, 3])
+        ts_new = ts.join(ts_other, strategy="union")
+        assert np.allclose(ts_new.time, [1, 2, 3])
+
+        ts = StingrayTimeseries()
+        ts_other = StingrayTimeseries()
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ts_new = ts.join(ts_other, strategy="union")
+        assert ts_new.time is None
+        assert ts_new.gti is None
+
+        ts = StingrayTimeseries(time=[1, 2, 3])
+        ts_other = StingrayTimeseries([])
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ts_new = ts.join(ts_other, strategy="union")
+        assert np.allclose(ts_new.time, [1, 2, 3])
+
+        ts = StingrayTimeseries([])
+        ts_other = StingrayTimeseries(time=[1, 2, 3])
+        ts_new = ts.join(ts_other, strategy="union")
+        assert np.allclose(ts_new.time, [1, 2, 3])
+
+    def test_join_different_dt(self):
+        ts = StingrayTimeseries(time=[10, 20, 30], dt=1)
+        ts_other = StingrayTimeseries(time=[40, 50, 60], dt=3)
+        with pytest.warns(UserWarning, match="The time resolution is different."):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.array_equal(ts_new.dt, [1, 1, 1, 3, 3, 3])
+        assert np.allclose(ts_new.time, [10, 20, 30, 40, 50, 60])
+
+    def test_join_different_instr(self):
+        ts = StingrayTimeseries(time=[10, 20, 30], instr="fpma")
+        ts_other = StingrayTimeseries(time=[40, 50, 60], instr="fpmb")
+        with pytest.warns(
+            UserWarning,
+            match="Attribute instr is different in the time series being merged.",
+        ):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert ts_new.instr == "fpma,fpmb"
+
+    def test_join_different_meta_attribute(self):
+        ts = StingrayTimeseries(time=[10, 20, 30])
+        ts_other = StingrayTimeseries(time=[40, 50, 60])
+        ts_other.bubu = "settete"
+        ts.whatstheanswer = 42
+        ts.unmovimentopara = "arriba"
+        ts_other.unmovimentopara = "abajo"
+
+        with pytest.warns(
+            UserWarning,
+            match=(
+                "Attribute (bubu|whatstheanswer|unmovimentopara) is different "
+                "in the time series being merged."
+            ),
+        ):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert ts_new.bubu == (None, "settete")
+        assert ts_new.whatstheanswer == (42, None)
+        assert ts_new.unmovimentopara == "arriba,abajo"
+
+    def test_join_without_energy(self):
+        ts = StingrayTimeseries(time=[1, 2, 3], energy=[3, 3, 3])
+        ts_other = StingrayTimeseries(time=[4, 5])
+        with pytest.warns(
+            UserWarning, match="The energy array is empty in one of the time series being merged."
+        ):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.allclose(ts_new.energy, [3, 3, 3, np.nan, np.nan], equal_nan=True)
+
+    def test_join_without_pi(self):
+        ts = StingrayTimeseries(time=[1, 2, 3], pi=[3, 3, 3])
+        ts_other = StingrayTimeseries(time=[4, 5])
+        with pytest.warns(
+            UserWarning, match="The pi array is empty in one of the time series being merged."
+        ):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.allclose(ts_new.pi, [3, 3, 3, np.nan, np.nan], equal_nan=True)
+
+    def test_join_with_arbitrary_attribute(self):
+        ts = StingrayTimeseries(time=[1, 2, 4])
+        ts_other = StingrayTimeseries(time=[3, 5])
+        ts.u = [3, 3, 3]
+        ts_other.q = [1, 2]
+        with pytest.warns(
+            UserWarning, match="The (u|q) array is empty in one of the time series being merged."
+        ):
+            ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.allclose(ts_new.q, [np.nan, np.nan, 1, np.nan, 2], equal_nan=True)
+        assert np.allclose(ts_new.u, [3, 3, np.nan, 3, np.nan], equal_nan=True)
+
+    def test_join_with_gti_none(self):
+        ts = StingrayTimeseries(time=[1, 2, 3])
+        ts_other = StingrayTimeseries(time=[4, 5], gti=[[3.5, 5.5]])
+        ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.allclose(ts_new.gti, [[1, 3], [3.5, 5.5]])
+
+        ts = StingrayTimeseries(time=[1, 2, 3], gti=[[0.5, 3.5]])
+        ts_other = StingrayTimeseries(time=[4, 5])
+        ts_new = ts.join(ts_other, strategy="union")
+
+        assert np.allclose(ts_new.gti, [[0.5, 3.5], [4, 5]])
+
+        ts = StingrayTimeseries(time=[1, 2, 3])
+        ts_other = StingrayTimeseries(time=[4, 5])
+        ts_new = ts.join(ts_other, strategy="union")
+
+        assert ts_new._gti is None
+
+    def test_non_overlapping_join_infer(self):
+        """Join two overlapping event lists."""
+        ts = StingrayTimeseries(
+            time=[1, 1.1, 2, 3, 4], energy=[3, 4, 7, 4, 3], gti=[[1, 2], [3, 4]]
+        )
+        ts_other = StingrayTimeseries(time=[5, 6, 6.1, 7, 10], energy=[4, 3, 8, 1, 2], gti=[[6, 7]])
+        ts_new = ts.join(ts_other, strategy="infer")
+
+        assert (ts_new.time == np.array([1, 1.1, 2, 3, 4, 5, 6, 6.1, 7, 10])).all()
+        assert (ts_new.energy == np.array([3, 4, 7, 4, 3, 4, 3, 8, 1, 2])).all()
+        assert (ts_new.gti == np.array([[1, 2], [3, 4], [6, 7]])).all()
+
+    def test_overlapping_join_infer(self):
+        """Join two non-overlapping event lists."""
+        ts = StingrayTimeseries(
+            time=[1, 1.1, 10, 6, 5], energy=[10, 6, 3, 11, 2], gti=[[1, 3], [5, 6]]
+        )
+        ts_other = StingrayTimeseries(
+            time=[5.1, 7, 6.1, 6.11, 10.1], energy=[2, 3, 8, 1, 2], gti=[[5, 7], [8, 10]]
+        )
+        ts_new = ts.join(ts_other, strategy="infer")
+
+        assert (ts_new.time == np.array([1, 1.1, 5, 5.1, 6, 6.1, 6.11, 7, 10, 10.1])).all()
+        assert (ts_new.energy == np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
+        assert (ts_new.gti == np.array([[5, 6]])).all()
+
+    def test_overlapping_join_change_mjdref(self):
+        """Join two non-overlapping event lists."""
+        ts = StingrayTimeseries(
+            time=[1, 1.1, 10, 6, 5], energy=[10, 6, 3, 11, 2], gti=[[1, 3], [5, 6]], mjdref=57001
+        )
+        ts_other = StingrayTimeseries(
+            time=np.asarray([5.1, 7, 6.1, 6.11, 10.1]) + 86400,
+            energy=[2, 3, 8, 1, 2],
+            gti=np.asarray([[5, 7], [8, 10]]) + 86400,
+            mjdref=57000,
+        )
+        with pytest.warns(UserWarning, match="Attribute mjdref is different"):
+            ts_new = ts.join(ts_other, strategy="intersection")
+
+        assert np.allclose(ts_new.time, np.array([1, 1.1, 5, 5.1, 6, 6.1, 6.11, 7, 10, 10.1]))
+        assert (ts_new.energy == np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
+        assert np.allclose(ts_new.gti, np.array([[5, 6]]))
+
+    def test_multiple_join(self):
+        """Test if multiple event lists can be joined."""
+        ts = StingrayTimeseries(time=[1, 2, 4], instr="a", mission=1)
+        ts_other = StingrayTimeseries(time=[3, 5, 7], instr="b", mission=2)
+        ts_other2 = StingrayTimeseries(time=[6, 8, 9], instr="c", mission=3)
+
+        ts.pibiri = [1, 1, 1]
+        ts_other.pibiri = [2, 2, 2]
+        ts_other2.pibiri = [3, 3, 3]
+
+        with pytest.warns(
+            UserWarning,
+            match="Attribute (instr|mission) is different in the time series being merged.",
+        ):
+            ts_new = ts.join([ts_other, ts_other2], strategy="union")
+        assert np.allclose(ts_new.time, [1, 2, 3, 4, 5, 6, 7, 8, 9])
+        assert np.allclose(ts_new.pibiri, [1, 1, 2, 1, 2, 3, 2, 3, 3])
+        assert ts_new.instr == "a,b,c"
+        assert ts_new.mission == (1, 2, 3)
+
+    def test_join_ignore_attr(self):
+        """Test if multiple event lists can be joined."""
+        ts = StingrayTimeseries(time=[1, 2, 4], instr="a", mission=1)
+        ts_other = StingrayTimeseries(time=[3, 5, 7], instr="b", mission=2)
+
+        with pytest.warns(
+            UserWarning,
+            match="Attribute mission is different in the time series being merged.",
+        ):
+            ts_new = ts.join([ts_other], strategy="union", ignore_meta=["instr"])
+
+        assert np.allclose(ts_new.time, [1, 2, 3, 4, 5, 7])
+        assert not hasattr(ts_new, "instr")
+        assert ts_new.mission == (1, 2)
diff --git a/stingray/tests/test_crossspectrum.py b/stingray/tests/test_crossspectrum.py
index d466b206..08ec13b8 100644
--- a/stingray/tests/test_crossspectrum.py
+++ b/stingray/tests/test_crossspectrum.py
@@ -175,6 +175,8 @@ class TestAveragedCrossspectrumEvents(object):
 
         self.events1 = EventList(times1, gti=gti)
         self.events2 = EventList(times2, gti=gti)
+        self.events1.fake_weights = np.ones_like(self.events1.time)
+        self.events2.fake_weights = np.ones_like(self.events2.time)
 
         self.cs = Crossspectrum(self.events1, self.events2, dt=self.dt, norm="none")
 
@@ -358,6 +360,45 @@ class TestAveragedCrossspectrumEvents(object):
         for attr in ["power", "freq", "m", "n", "nphots1", "nphots2", "segment_size"]:
             assert np.allclose(getattr(pds, attr), getattr(pds_ev, attr))
 
+    @pytest.mark.parametrize("norm", ["frac", "abs", "none", "leahy"])
+    def test_from_timeseries_with_err_works(self, norm):
+        lc1 = self.events1.to_binned_timeseries(self.dt)
+        lc2 = self.events2.to_binned_timeseries(self.dt)
+        lc1.counts_err = np.sqrt(lc1.counts.mean()) + np.zeros_like(lc1.counts)
+        lc2.counts_err = np.sqrt(lc2.counts.mean()) + np.zeros_like(lc2.counts)
+        pds = AveragedCrossspectrum.from_stingray_timeseries(
+            lc1, lc2, "counts", "counts_err", segment_size=self.segment_size, norm=norm, silent=True
+        )
+        pds = AveragedCrossspectrum.from_stingray_timeseries(
+            lc1,
+            lc2,
+            "counts",
+            "counts_err",
+            segment_size=self.segment_size,
+            norm=norm,
+            silent=True,
+        )
+        pds_weight = AveragedCrossspectrum.from_stingray_timeseries(
+            lc1,
+            lc2,
+            "fake_weights",
+            "counts_err",
+            segment_size=self.segment_size,
+            norm=norm,
+            silent=True,
+        )
+        pds_ev = AveragedCrossspectrum.from_events(
+            self.events1,
+            self.events2,
+            segment_size=self.segment_size,
+            dt=self.dt,
+            norm=norm,
+            silent=True,
+        )
+        for attr in ["power", "freq", "m", "n", "nphots1", "nphots2", "segment_size"]:
+            assert np.allclose(getattr(pds, attr), getattr(pds_ev, attr))
+            assert np.allclose(getattr(pds_weight, attr), getattr(pds_ev, attr))
+
     def test_it_works_with_events(self):
         lc1 = self.events1.to_lc(self.dt)
         lc2 = self.events2.to_lc(self.dt)
@@ -1302,14 +1343,14 @@ class TestRoundTrip:
 
         self._check_equal(so, new_so)
 
-    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv", "fits", "hdf5"])
+    @pytest.mark.parametrize("fmt", ["pickle", "hdf5"])
     def test_file_roundtrip(self, fmt):
         so = self.cs
         fname = f"dummy.{fmt}"
         if not _HAS_H5PY and fmt == "hdf5":
             with pytest.raises(Exception) as excinfo:
                 so.write(fname, fmt=fmt)
-                assert h5py in str(excinfo.value)
+                assert "h5py" in str(excinfo.value)
             return
         so.write(fname, fmt=fmt)
         new_so = so.read(fname, fmt=fmt)
@@ -1317,6 +1358,17 @@ class TestRoundTrip:
 
         self._check_equal(so, new_so)
 
+    @pytest.mark.parametrize("fmt", ["ascii", "ascii.ecsv", "fits"])
+    def test_file_roundtrip_lossy(self, fmt):
+        so = self.cs
+        fname = f"dummy.{fmt}"
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            so.write(fname, fmt=fmt)
+        new_so = so.read(fname, fmt=fmt)
+        os.unlink(fname)
+
+        self._check_equal(so, new_so)
+
 
 class TestDynamicalCrossspectrum(object):
     def setup_class(cls):
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index a37ac4bb..f2bc04af 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -1,5 +1,6 @@
 import warnings
 import numpy as np
+import copy
 import os
 import pytest
 from astropy.time import Time
@@ -55,12 +56,22 @@ class TestEvents(object):
             _ = EventList(self.time, self.counts, gti=self.gti, bubu="settete")
         assert np.any(["Unrecognized keywords:" in r.message.args[0] for r in record])
 
+    def test_warn_wrong_keywords(self):
+        with pytest.warns(DeprecationWarning, match="The ncounts keyword does nothing"):
+            _ = EventList(self.time, self.counts, gti=self.gti, ncounts=10)
+
     def test_initiate_from_ndarray(self):
         times = np.sort(np.random.uniform(1e8, 1e8 + 1000, 101).astype(np.longdouble))
         ev = EventList(times, mjdref=54600)
         assert np.allclose(ev.time, times, atol=1e-15)
         assert np.allclose(ev.mjdref, 54600)
 
+    def test_print(self):
+        times = [1.01, 2, 3]
+        ev = EventList(times, mjdref=54600)
+
+        print(ev)
+
     def test_initiate_from_astropy_time(self):
         times = np.sort(np.random.uniform(1e8, 1e8 + 1000, 101).astype(np.longdouble))
         mjdref = 54600
@@ -81,6 +92,17 @@ class TestEvents(object):
         with pytest.raises(ValueError):
             EventList(time=[1, 2, 3], energy=[10, 12])
 
+    def test_concatenate(self):
+        """Join two overlapping event lists."""
+        ev = EventList(time=[1, 1.1, 2, 3, 4], energy=[3, 4, 7, 4, 3], gti=[[1, 2], [3, 4]])
+        ev_other1 = EventList(time=[5, 6, 6.1], energy=[4, 3, 8], gti=[[6, 6.2]])
+        ev_other2 = EventList(time=[7, 10], energy=[1, 2], gti=[[6.5, 7]])
+        ev_new = ev.concatenate([ev_other1, ev_other2])
+
+        assert (ev_new.time == np.array([1, 1.1, 2, 3, 4, 5, 6, 6.1, 7, 10])).all()
+        assert (ev_new.energy == np.array([3, 4, 7, 4, 3, 4, 3, 8, 1, 2])).all()
+        assert (ev_new.gti == np.array([[1, 2], [3, 4], [6, 6.2], [6.5, 7]])).all()
+
     def test_to_lc(self):
         """Create a light curve from event list."""
         ev = EventList(self.time, gti=self.gti)
@@ -88,12 +110,23 @@ class TestEvents(object):
         assert np.allclose(lc.time, [0.5, 1.5, 2.5, 3.5])
         assert (lc.gti == self.gti).all()
 
+    def test_to_timeseries(self):
+        """Create a time series from event list."""
+        ev = EventList(self.time, gti=self.gti)
+        ev.bla = np.zeros_like(ev.time) + 2
+        lc = ev.to_lc(1)
+        ts = ev.to_binned_timeseries(1)
+        assert np.allclose(ts.time, [0.5, 1.5, 2.5, 3.5])
+        assert (ts.gti == self.gti).all()
+        assert np.array_equal(ts.counts, lc.counts)
+        assert np.array_equal(ts.bla, ts.counts * 2)
+
     def test_from_lc(self):
         """Load event list from lightcurve"""
-        lc = Lightcurve(time=[0.5, 1.5, 2.5], counts=[2, 1, 2])
+        lc = Lightcurve(time=[0.5, 1.5, 2.5], counts=[2, -1, 2])
         ev = EventList.from_lc(lc)
 
-        assert (ev.time == np.array([0.5, 0.5, 1.5, 2.5, 2.5])).all()
+        assert np.array_equal(ev.time, np.array([0.5, 0.5, 2.5, 2.5]))
 
     def test_simulate_times_warns_bin_time(self):
         """Simulate photon arrival times for an event list
@@ -120,34 +153,34 @@ class TestEvents(object):
 
     def test_simulate_energies(self):
         """Assign photon energies to an event list."""
-        ev = EventList(ncounts=100)
+        ev = EventList(np.arange(10))
         ev.simulate_energies(self.spectrum)
 
     def test_simulate_energies_with_1d_spectrum(self):
         """Test that simulate_energies() method raises index
         error exception is spectrum is 1-d.
         """
-        ev = EventList(ncounts=100)
+        ev = EventList(np.arange(10))
         with pytest.raises(IndexError):
             ev.simulate_energies(self.spectrum[0])
 
     def test_simulate_energies_with_wrong_spectrum_type(self):
         """Test that simulate_energies() method raises type error
-        exception when wrong sepctrum type is supplied.
+        exception when wrong spectrum type is supplied.
         """
-        ev = EventList(ncounts=100)
+        ev = EventList(np.arange(10))
         with pytest.raises(TypeError):
             ev.simulate_energies(1)
 
     def test_simulate_energies_with_counts_not_set(self):
         ev = EventList()
-        with pytest.warns(UserWarning):
+        with pytest.warns(UserWarning, match="empty event list"):
             ev.simulate_energies(self.spectrum)
 
     def test_compare_energy(self):
         """Compare the simulated energy distribution to actual distribution."""
         fluxes = np.array(self.spectrum[1])
-        ev = EventList(ncounts=1000)
+        ev = EventList(np.arange(1000))
         ev.simulate_energies(self.spectrum)
 
         # Note: I'm passing the edges: when the bin center is 1, the
@@ -164,9 +197,9 @@ class TestEvents(object):
     @pytest.mark.skipif("not (_HAS_YAML)")
     def test_io_with_ascii(self):
         ev = EventList(self.time)
-        ev.write("ascii_ev.ecsv", fmt="ascii")
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            ev.write("ascii_ev.ecsv", fmt="ascii")
         ev = ev.read("ascii_ev.ecsv", fmt="ascii")
-        print(ev.time, self.time)
         assert np.allclose(ev.time, self.time)
         os.remove("ascii_ev.ecsv")
 
@@ -197,7 +230,8 @@ class TestEvents(object):
 
     def test_io_with_fits(self):
         ev = EventList(time=self.time, mjdref=54000)
-        ev.write("ev.fits", fmt="fits")
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            ev.write("ev.fits", fmt="fits")
         ev = ev.read("ev.fits", fmt="fits")
         assert np.allclose(ev.time, self.time)
         os.remove("ev.fits")
@@ -310,8 +344,8 @@ class TestJoinEvents:
         ev = EventList()
         ev_other = EventList()
 
-        with pytest.warns(UserWarning, match="One of the event lists you are joining is empty"):
-            assert ev.join(ev_other).time is None
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            assert ev.join(ev_other, strategy="union").time is None
 
     def test_join_empty_lists(self):
         """Test if an empty event list can be concatenated
@@ -319,51 +353,52 @@ class TestJoinEvents:
         """
         ev = EventList(time=[1, 2, 3])
         ev_other = EventList()
-        with pytest.warns(UserWarning, match="One of the event lists you are joining"):
-            ev_new = ev.join(ev_other)
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ev_new = ev.join(ev_other, strategy="union")
         assert np.allclose(ev_new.time, [1, 2, 3])
 
         ev = EventList()
         ev_other = EventList(time=[1, 2, 3])
-        ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="union")
         assert np.allclose(ev_new.time, [1, 2, 3])
 
         ev = EventList()
         ev_other = EventList()
-        with pytest.warns(UserWarning, match="One of the event lists you are joining"):
-            ev_new = ev.join(ev_other)
-        assert ev_new.time == None
-        assert ev_new.gti == None
-        assert ev_new.pi == None
-        assert ev_new.energy == None
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ev_new = ev.join(ev_other, strategy="union")
+        assert ev_new.time is None
+        assert ev_new.gti is None
+        assert ev_new.pi is None
+        assert ev_new.energy is None
 
         ev = EventList(time=[1, 2, 3])
         ev_other = EventList([])
-        with pytest.warns(UserWarning, match="One of the event lists you are joining"):
-            ev_new = ev.join(ev_other)
+        with pytest.warns(UserWarning, match="One of the time series you are joining is empty."):
+            ev_new = ev.join(ev_other, strategy="union")
         assert np.allclose(ev_new.time, [1, 2, 3])
 
         ev = EventList([])
         ev_other = EventList(time=[1, 2, 3])
-        ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="union")
         assert np.allclose(ev_new.time, [1, 2, 3])
 
     def test_join_different_dt(self):
         ev = EventList(time=[10, 20, 30], dt=1)
         ev_other = EventList(time=[40, 50, 60], dt=3)
         with pytest.warns(UserWarning, match="The time resolution is different."):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
-        assert ev_new.dt == 3
+        assert np.array_equal(ev_new.dt, [1, 1, 1, 3, 3, 3])
         assert np.allclose(ev_new.time, [10, 20, 30, 40, 50, 60])
 
     def test_join_different_instr(self):
         ev = EventList(time=[10, 20, 30], instr="fpma")
         ev_other = EventList(time=[40, 50, 60], instr="fpmb")
         with pytest.warns(
-            UserWarning, match="Attribute instr is different in the event lists being merged."
+            UserWarning,
+            match="Attribute instr is different in the time series being merged.",
         ):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
         assert ev_new.instr == "fpma,fpmb"
 
@@ -377,9 +412,12 @@ class TestJoinEvents:
 
         with pytest.warns(
             UserWarning,
-            match="Attribute (bubu|whatstheanswer|unmovimentopara) is different in the event lists being merged.",
+            match=(
+                "Attribute (bubu|whatstheanswer|unmovimentopara) is different "
+                "in the time series being merged."
+            ),
         ):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
         assert ev_new.bubu == (None, "settete")
         assert ev_new.whatstheanswer == (42, None)
@@ -389,9 +427,9 @@ class TestJoinEvents:
         ev = EventList(time=[1, 2, 3], energy=[3, 3, 3])
         ev_other = EventList(time=[4, 5])
         with pytest.warns(
-            UserWarning, match="The energy array is empty in one of the event lists being merged."
+            UserWarning, match="The energy array is empty in one of the time series being merged."
         ):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
         assert np.allclose(ev_new.energy, [3, 3, 3, np.nan, np.nan], equal_nan=True)
 
@@ -399,9 +437,9 @@ class TestJoinEvents:
         ev = EventList(time=[1, 2, 3], pi=[3, 3, 3])
         ev_other = EventList(time=[4, 5])
         with pytest.warns(
-            UserWarning, match="The pi array is empty in one of the event lists being merged."
+            UserWarning, match="The pi array is empty in one of the time series being merged."
         ):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
         assert np.allclose(ev_new.pi, [3, 3, 3, np.nan, np.nan], equal_nan=True)
 
@@ -411,9 +449,9 @@ class TestJoinEvents:
         ev.u = [3, 3, 3]
         ev_other.q = [1, 2]
         with pytest.warns(
-            UserWarning, match="The (u|q) array is empty in one of the event lists being merged."
+            UserWarning, match="The (u|q) array is empty in one of the time series being merged."
         ):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="union")
 
         assert np.allclose(ev_new.q, [np.nan, np.nan, 1, np.nan, 2], equal_nan=True)
         assert np.allclose(ev_new.u, [3, 3, np.nan, 3, np.nan], equal_nan=True)
@@ -421,42 +459,39 @@ class TestJoinEvents:
     def test_join_with_gti_none(self):
         ev = EventList(time=[1, 2, 3])
         ev_other = EventList(time=[4, 5], gti=[[3.5, 5.5]])
-        with pytest.warns(UserWarning, match="GTIs in these two event lists do not overlap"):
-            ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="union")
 
         assert np.allclose(ev_new.gti, [[1, 3], [3.5, 5.5]])
 
         ev = EventList(time=[1, 2, 3], gti=[[0.5, 3.5]])
         ev_other = EventList(time=[4, 5])
-        with pytest.warns(UserWarning, match="GTIs in these two event lists do not overlap"):
-            ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="union")
 
         assert np.allclose(ev_new.gti, [[0.5, 3.5], [4, 5]])
 
         ev = EventList(time=[1, 2, 3])
         ev_other = EventList(time=[4, 5])
-        ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="union")
 
-        assert ev_new.gti == None
+        assert ev_new._gti is None
 
-    def test_non_overlapping_join(self):
+    def test_non_overlapping_join_infer(self):
         """Join two overlapping event lists."""
         ev = EventList(time=[1, 1.1, 2, 3, 4], energy=[3, 4, 7, 4, 3], gti=[[1, 2], [3, 4]])
         ev_other = EventList(time=[5, 6, 6.1, 7, 10], energy=[4, 3, 8, 1, 2], gti=[[6, 7]])
-        with pytest.warns(UserWarning, match="GTIs in these two event lists do not overlap"):
-            ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="infer")
 
         assert (ev_new.time == np.array([1, 1.1, 2, 3, 4, 5, 6, 6.1, 7, 10])).all()
         assert (ev_new.energy == np.array([3, 4, 7, 4, 3, 4, 3, 8, 1, 2])).all()
         assert (ev_new.gti == np.array([[1, 2], [3, 4], [6, 7]])).all()
 
-    def test_overlapping_join(self):
+    def test_overlapping_join_infer(self):
         """Join two non-overlapping event lists."""
         ev = EventList(time=[1, 1.1, 10, 6, 5], energy=[10, 6, 3, 11, 2], gti=[[1, 3], [5, 6]])
         ev_other = EventList(
             time=[5.1, 7, 6.1, 6.11, 10.1], energy=[2, 3, 8, 1, 2], gti=[[5, 7], [8, 10]]
         )
-        ev_new = ev.join(ev_other)
+        ev_new = ev.join(ev_other, strategy="infer")
 
         assert (ev_new.time == np.array([1, 1.1, 5, 5.1, 6, 6.1, 6.11, 7, 10, 10.1])).all()
         assert (ev_new.energy == np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
@@ -474,7 +509,7 @@ class TestJoinEvents:
             mjdref=57000,
         )
         with pytest.warns(UserWarning, match="Attribute mjdref is different"):
-            ev_new = ev.join(ev_other)
+            ev_new = ev.join(ev_other, strategy="intersection")
 
         assert np.allclose(ev_new.time, np.array([1, 1.1, 5, 5.1, 6, 6.1, 6.11, 7, 10, 10.1]))
         assert (ev_new.energy == np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
@@ -492,10 +527,74 @@ class TestJoinEvents:
 
         with pytest.warns(
             UserWarning,
-            match="Attribute (instr|mission) is different in the event lists being merged.",
+            match="Attribute (instr|mission) is different in the time series being merged.",
         ):
-            ev_new = ev.join([ev_other, ev_other2])
+            ev_new = ev.join([ev_other, ev_other2], strategy="union")
         assert np.allclose(ev_new.time, [1, 2, 3, 4, 5, 6, 7, 8, 9])
         assert np.allclose(ev_new.pibiri, [1, 1, 2, 1, 2, 3, 2, 3, 3])
         assert ev_new.instr == "a,b,c"
         assert ev_new.mission == (1, 2, 3)
+
+
+class TestFilters(object):
+    @classmethod
+    def setup_class(cls):
+        events = np.array([1, 1.05, 1.07, 1.08, 1.1, 2, 2.2, 3, 3.1, 3.2])
+        events = EventList(events, gti=[[0, 3.3]])
+        events.pi = np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
+        events.energy = np.array([1, 2, 2, 2, 2, 1, 1, 1, 2, 1])
+        events.mjdref = 10
+        cls.events = events
+
+    @pytest.mark.parametrize("inplace", [True, False])
+    def test_apply_mask(self, inplace):
+        events = copy.deepcopy(self.events)
+        mask = [True, False, False, False, False, True, True, True, False, True]
+        filt_events = events.apply_mask(mask, inplace=inplace)
+        if inplace:
+            assert filt_events is events
+            assert np.allclose(events.pi, 1)
+        else:
+            assert filt_events is not events
+            assert not np.allclose(events.pi, 1)
+
+        expected = np.array([1, 2, 2.2, 3, 3.2])
+        assert np.allclose(filt_events.time, expected)
+        assert np.allclose(filt_events.pi, 1)
+        assert np.allclose(filt_events.energy, 1)
+
+    @pytest.mark.parametrize("inplace", [True, False])
+    @pytest.mark.parametrize("use_pi", [True, False])
+    def test_filter_energy_range(self, inplace, use_pi):
+        events = copy.deepcopy(self.events)
+
+        filt_events = events.filter_energy_range([0.5, 1.5], use_pi=use_pi, inplace=inplace)
+        if inplace:
+            assert filt_events is events
+            assert np.allclose(events.pi, 1)
+        else:
+            assert filt_events is not events
+            assert not np.allclose(events.pi, 1)
+
+        expected = np.array([1, 2, 2.2, 3, 3.2])
+        assert np.allclose(filt_events.time, expected)
+        assert np.allclose(filt_events.pi, 1)
+        assert np.allclose(filt_events.energy, 1)
+
+    @pytest.mark.parametrize("inplace", [True, False])
+    def test_apply_deadtime(self, inplace):
+        events = copy.deepcopy(self.events)
+        filt_events, _ = events.apply_deadtime(
+            0.11, inplace=inplace, verbose=False, return_all=True
+        )
+        if inplace:
+            assert filt_events is events
+            assert np.allclose(events.pi, 1)
+        else:
+            assert filt_events is not events
+            assert not np.allclose(events.pi, 1)
+
+        expected = np.array([1, 2, 2.2, 3, 3.2])
+        assert np.allclose(filt_events.time, expected)
+        assert np.allclose(filt_events.pi, 1)
+        assert np.allclose(filt_events.energy, 1)
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index bc237107..f385f9b7 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -48,6 +48,8 @@ except ImportError:
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
 
+plt.close("all")
+
 
 def fvar_fun(lc):
     from stingray.utils import excess_variance
@@ -82,11 +84,15 @@ class TestProperties(object):
         cls.lc_lowmem = Lightcurve(times, counts, gti=cls.gti, low_memory=True)
 
     def test_empty_lightcurve(self):
-        with pytest.warns(UserWarning, match="No time values passed to Lightcurve object!"):
-            Lightcurve()
+        lc0 = Lightcurve()
+        lc1 = Lightcurve([], [])
+        assert lc0.time is None
+        assert lc1.time is None
 
-        with pytest.warns(UserWarning, match="No time values passed to Lightcurve object!"):
-            Lightcurve([], [])
+    def test_add_data_to_empty_lightcurve(self):
+        lc0 = Lightcurve()
+        lc0.time = [1, 2, 3]
+        lc0.counts = [1, 2, 3]
 
     def test_bad_counts_lightcurve(self):
         with pytest.raises(StingrayError, match="Empty or invalid counts array. "):
@@ -150,6 +156,12 @@ class TestProperties(object):
         _ = lc.bin_lo
         assert lc._bin_lo is not None
 
+    def test_make_all_none(self):
+        lc = copy.deepcopy(self.lc)
+        lc.time = None
+        assert lc.counts is None
+        assert lc._counts is None
+
     def test_lightcurve_from_astropy_time(self):
         time = Time([57483, 57484], format="mjd")
         counts = np.array([2, 2])
@@ -246,18 +258,36 @@ class TestProperties(object):
         # (because low_memory, and input_counts is now False)
         assert lc._counts_err is None
 
-    @pytest.mark.parametrize(
-        "property", "time,counts,counts_err," "countrate,countrate_err".split(",")
-    )
-    def test_assign_bad_shape_fails(self, property):
+    @pytest.mark.parametrize("attr", "time,counts,countrate".split(","))
+    def test_add_data_to_empty_lightcurve_wrong(self, attr):
+        lc0 = Lightcurve()
+        lc0.time = [1, 2, 3]
+        with pytest.raises(ValueError, match=f".*the same shape as the time array"):
+            setattr(lc0, attr, [1, 2, 3, 4])
+
+    @pytest.mark.parametrize("attr", "counts,countrate".split(","))
+    def test_add_err_data_to_empty_lightcurve_wrong_order(self, attr):
+        lc0 = Lightcurve()
+        lc0.time = [1, 2, 3]
+        with pytest.raises(ValueError, match=f"if the {attr} array is not None"):
+            setattr(lc0, attr + "_err", [1, 2, 3])
+
+    @pytest.mark.parametrize("attr", "counts,countrate".split(","))
+    def test_add_err_data_to_empty_lightcurve_wrong_size(self, attr):
+        lc0 = Lightcurve()
+        lc0.time = [1, 2, 3]
+        setattr(lc0, attr, [1, 2, 1])
+        with pytest.raises(ValueError, match=f"the same shape as the {attr} array"):
+            setattr(lc0, attr + "_err", [1, 2, 3, 4])
+
+    @pytest.mark.parametrize("attr", "time,counts,counts_err,countrate,countrate_err".split(","))
+    def test_assign_scalar_data(self, attr):
         lc = copy.deepcopy(self.lc)
         # Same shape passes
-        setattr(lc, property, np.zeros_like(lc.time))
+        setattr(lc, attr, np.zeros_like(lc.time))
         # Different shape doesn't
-        with pytest.raises(ValueError):
-            setattr(lc, property, 3)
-        with pytest.raises(ValueError):
-            setattr(lc, property, np.arange(2))
+        with pytest.raises(ValueError, match="at least 1D"):
+            setattr(lc, attr, 3)
 
 
 class TestChunks(object):
@@ -322,6 +352,15 @@ class TestLightcurve(object):
         """
         lc = Lightcurve(self.times, self.counts)
 
+    def test_print(self, capsys):
+        lc = Lightcurve(self.times, self.counts, header="TEST")
+
+        print(lc)
+        captured = capsys.readouterr()
+        assert "header" not in captured.out
+        assert "time" in captured.out
+        assert "counts" in captured.out
+
     def test_irregular_time_warning(self):
         """
         Check if inputting an irregularly spaced time iterable throws out
@@ -556,12 +595,12 @@ class TestLightcurve(object):
     def test_add_with_different_time_arrays(self):
         _times = [1.1, 2.1, 3.1, 4.1, 5.1]
         _counts = [2, 2, 2, 2, 2]
+        lc1 = Lightcurve(self.times, self.counts)
+        lc2 = Lightcurve(_times, _counts)
 
-        with pytest.raises(ValueError):
-            lc1 = Lightcurve(self.times, self.counts)
-            lc2 = Lightcurve(_times, _counts)
-
-            lc = lc1 + lc2
+        with pytest.warns(UserWarning, match="The good time intervals in the two time series"):
+            with pytest.raises(ValueError):
+                lc = lc1 + lc2
 
     def test_add_with_different_err_dist(self):
         lc1 = Lightcurve(self.times, self.counts)
@@ -581,16 +620,17 @@ class TestLightcurve(object):
         gti = [[0.0, 3.5]]
         lc1 = Lightcurve(self.times, self.counts, gti=self.gti)
         lc2 = Lightcurve(self.times, self.counts, gti=gti)
-        lc = lc1 + lc2
+        with pytest.warns(UserWarning, match="The good time intervals in the two time series"):
+            lc = lc1 + lc2
         np.testing.assert_almost_equal(lc.gti, [[0.5, 3.5]])
 
     def test_add_with_unequal_time_arrays(self):
         _times = [1, 3, 5, 7]
 
-        with pytest.raises(ValueError):
-            lc1 = Lightcurve(self.times, self.counts)
-            lc2 = Lightcurve(_times, self.counts)
+        lc1 = Lightcurve(self.times, self.counts)
+        lc2 = Lightcurve(_times, self.counts)
 
+        with pytest.raises(ValueError):
             lc = lc1 + lc2
 
     def test_add_with_equal_time_arrays(self):
@@ -609,11 +649,12 @@ class TestLightcurve(object):
         _times = [1.1, 2.1, 3.1, 4.1, 5.1]
         _counts = [2, 2, 2, 2, 2]
 
-        with pytest.raises(ValueError):
-            lc1 = Lightcurve(self.times, self.counts)
-            lc2 = Lightcurve(_times, _counts)
+        lc1 = Lightcurve(self.times, self.counts)
+        lc2 = Lightcurve(_times, _counts)
 
-            _ = lc1 - lc2
+        with pytest.warns(UserWarning, match="The good time intervals in the two time series"):
+            with pytest.raises(ValueError):
+                _ = lc1 - lc2
 
     def test_sub_with_different_err_dist(self):
         lc1 = Lightcurve(self.times, self.counts)
@@ -743,6 +784,23 @@ class TestLightcurve(object):
 
         assert np.allclose(newlc.counts, 0)
 
+    def test_concatenate(self):
+        time0 = [1, 2, 3, 4]
+        time1 = [5, 6, 7, 8, 9]
+        count0 = [10, 20, 30, 40]
+        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)
+        lc1 = Lightcurve(time1, counts=count1, dt=1, gti=gti1)
+        with pytest.warns(UserWarning, match="The _counts_err array is empty in one of the"):
+            lc = lc0.concatenate(lc1)
+        assert np.allclose(lc.counts, count0 + count1)
+        # Errors have been defined inside
+        assert len(lc.counts_err) == len(lc.counts)
+        assert np.allclose(lc.time, time0 + time1)
+        assert np.allclose(lc.gti, [[0.5, 9.5]])
+
     def test_join_disjoint_time_arrays(self):
         _times = [5, 6, 7, 8]
         _counts = [2, 2, 2, 2]
@@ -1079,11 +1137,13 @@ class TestLightcurve(object):
         assert_allclose(sr.counts_err, lc.flux_err)
 
     def test_plot_simple(self):
+        plt.close("all")
         lc = Lightcurve(self.times, self.counts)
         with warnings.catch_warnings():
             warnings.simplefilter("ignore", category=UserWarning)
             lc.plot()
         assert plt.fignum_exists(1)
+        plt.close("all")
 
     def test_plot_wrong_label_type(self):
         lc = Lightcurve(self.times, self.counts)
@@ -1091,6 +1151,7 @@ class TestLightcurve(object):
         with pytest.raises(TypeError):
             with pytest.warns(UserWarning, match="must be either a list or tuple") as w:
                 lc.plot(labels=123)
+        plt.close("all")
 
     def test_plot_labels_index_error(self):
         lc = Lightcurve(self.times, self.counts)
@@ -1098,18 +1159,21 @@ class TestLightcurve(object):
             lc.plot(labels=("x"))
 
             assert np.any(["must have two labels" in str(wi.message) for wi in w])
+        plt.close("all")
 
     def test_plot_default_filename(self):
         lc = Lightcurve(self.times, self.counts)
         lc.plot(save=True)
         assert os.path.isfile("out.png")
         os.unlink("out.png")
+        plt.close("all")
 
     def test_plot_custom_filename(self):
         lc = Lightcurve(self.times, self.counts)
         lc.plot(save=True, filename="lc.png")
         assert os.path.isfile("lc.png")
         os.unlink("lc.png")
+        plt.close("all")
 
     def test_plot_axis(self):
         lc = Lightcurve(self.times, self.counts)
@@ -1117,6 +1181,7 @@ class TestLightcurve(object):
             warnings.simplefilter("ignore", category=UserWarning)
             lc.plot(axis=[0, 1, 0, 100])
         assert plt.fignum_exists(1)
+        plt.close("all")
 
     def test_plot_title(self):
         lc = Lightcurve(self.times, self.counts)
@@ -1124,6 +1189,7 @@ class TestLightcurve(object):
             warnings.simplefilter("ignore", category=UserWarning)
             lc.plot(title="Test Lightcurve")
         assert plt.fignum_exists(1)
+        plt.close("all")
 
     def test_read_from_lcurve_1(self):
         fname = "lcurveA.fits"
@@ -1144,12 +1210,22 @@ class TestLightcurve(object):
     @pytest.mark.skipif("not _HAS_YAML")
     def test_io_with_ascii(self):
         lc = Lightcurve(self.times, self.counts)
-        lc.write("ascii_lc.ecsv", fmt="ascii")
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            lc.write("ascii_lc.ecsv", fmt="ascii")
         lc = lc.read("ascii_lc.ecsv", fmt="ascii")
         assert np.allclose(lc.time, self.times)
         assert np.allclose(lc.counts, self.counts)
         os.remove("ascii_lc.ecsv")
 
+    def test_io_with_fits(self):
+        lc = Lightcurve(self.times, self.counts)
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            lc.write("ascii_lc.fits", fmt="fits")
+        lc = lc.read("ascii_lc.fits", fmt="fits")
+        assert np.allclose(lc.time, self.times)
+        assert np.allclose(lc.counts, self.counts)
+        os.remove("ascii_lc.fits")
+
     def test_io_with_pickle(self):
         lc = Lightcurve(self.times, self.counts)
         lc.write("lc.pickle", fmt="pickle")
diff --git a/stingray/tests/test_lombscargle.py b/stingray/tests/test_lombscargle.py
index d352f707..4a8f76f5 100644
--- a/stingray/tests/test_lombscargle.py
+++ b/stingray/tests/test_lombscargle.py
@@ -59,7 +59,7 @@ class TestLombScargleCrossspectrum:
         ev_lscs = LombScargleCrossspectrum(ev1, ev2, dt=1)
         lc_lscs = LombScargleCrossspectrum(lc1, lc2, dt=1)
 
-        assert np.argmax(lc_lscs) == np.argmax(ev_lscs)
+        assert np.argmax(lc_lscs.power) == np.argmax(ev_lscs.power)
         assert np.all(ev_lscs.freq == lc_lscs.freq)
         assert np.all(ev_lscs.power == lc_lscs.power)
         assert ev_lscs.freq[np.argmax(ev_lscs.power)] == lc_lscs.freq[np.argmax(lc_lscs.power)] != 0
diff --git a/stingray/tests/test_powerspectrum.py b/stingray/tests/test_powerspectrum.py
index 1a4cff65..9cf34fa6 100644
--- a/stingray/tests/test_powerspectrum.py
+++ b/stingray/tests/test_powerspectrum.py
@@ -53,6 +53,7 @@ class TestAveragedPowerspectrumEvents(object):
         gti = np.array([[tstart, tend]])
 
         cls.events = EventList(times, gti=gti)
+        cls.events.fake_weights = np.ones_like(cls.events.time)
 
         cls.lc = cls.events
         cls.leahy_pds = AveragedPowerspectrum(
@@ -244,6 +245,39 @@ class TestAveragedPowerspectrumEvents(object):
         for attr in ["power", "freq", "m", "n", "nphots", "segment_size"]:
             assert np.allclose(getattr(pds, attr), getattr(pds_ev, attr))
 
+    @pytest.mark.parametrize("norm", ["frac", "abs", "none", "leahy"])
+    def test_from_timeseries_with_err_works(self, norm):
+        lc = self.events.to_binned_timeseries(self.dt)
+        lc._counts_err = np.sqrt(lc.counts.mean()) + np.zeros_like(lc.counts)
+        pds = AveragedPowerspectrum.from_stingray_timeseries(
+            lc,
+            "counts",
+            "_counts_err",
+            segment_size=self.segment_size,
+            norm=norm,
+            silent=True,
+            gti=lc.gti,
+        )
+        pds_weight = AveragedPowerspectrum.from_stingray_timeseries(
+            lc,
+            "fake_weights",
+            "_counts_err",
+            segment_size=self.segment_size,
+            norm=norm,
+            silent=True,
+        )
+        pds_ev = AveragedPowerspectrum.from_events(
+            self.events,
+            segment_size=self.segment_size,
+            dt=self.dt,
+            norm=norm,
+            silent=True,
+            gti=self.events.gti,
+        )
+        for attr in ["power", "freq", "m", "n", "nphots", "segment_size"]:
+            assert np.allclose(getattr(pds, attr), getattr(pds_ev, attr))
+            assert np.allclose(getattr(pds_weight, attr), getattr(pds_ev, attr))
+
     def test_init_without_segment(self):
         with pytest.raises(ValueError):
             assert AveragedPowerspectrum(self.lc, dt=self.dt)
@@ -1174,7 +1208,7 @@ class TestRoundTrip:
 
         self._check_equal(so, new_so)
 
-    @pytest.mark.parametrize("fmt", ["pickle", "ascii", "ascii.ecsv", "fits"])
+    @pytest.mark.parametrize("fmt", ["pickle"])
     def test_file_roundtrip(self, fmt):
         so = self.cs
         fname = f"dummy.{fmt}"
@@ -1183,3 +1217,14 @@ class TestRoundTrip:
         # os.unlink(fname)
 
         self._check_equal(so, new_so)
+
+    @pytest.mark.parametrize("fmt", ["ascii", "ascii.ecsv", "fits"])
+    def test_file_roundtrip_lossy(self, fmt):
+        so = self.cs
+        fname = f"dummy.{fmt}"
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            so.write(fname, fmt=fmt)
+        new_so = so.read(fname, fmt=fmt)
+        os.unlink(fname)
+
+        self._check_equal(so, new_so)
diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
index ea6d640f..a46f4e74 100644
--- a/stingray/tests/test_varenergyspectrum.py
+++ b/stingray/tests/test_varenergyspectrum.py
@@ -486,10 +486,11 @@ class TestRoundTrip:
         os.unlink("dummy.hdf5")
         self._check_equal(so, new_so)
 
-    @pytest.mark.parametrize("fmt", ["ascii.ecsv", "fits"])
+    @pytest.mark.parametrize("fmt", ["ascii.ecsv", "ascii", "fits"])
     def test_file_export(self, fmt):
         so = self.vespec
-        so.write("dummy", fmt=fmt)
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            so.write("dummy", fmt=fmt)
         new_so = Table.read("dummy", format=fmt)
         os.unlink("dummy")
         self._check_equal(so, new_so)