Repository: StingraySoftware/stingray
Pull request #796 (2024-02-06 09:36:29)

Update and reorder tutorials
This is an update of the spectral timing tutorial, but also a reorder of the list of tutorials in the main page to give more visibility and centrality to the spectral timing tutorial (that shows off a lot of Stingray's tools, including the new HIDs and power colors).
I set up a build on readthedocs to show the new page: https://docs.stingray.science/en/update_docs/
(Compare to https://docs.stingray.science/en/latest/)

Comments welcome.

codecov[bot]: ## [Codecov](https://app.codecov.io/gh/StingraySoftware/stingray/pull/796?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`5285c23`)](https://app.codecov.io/gh/StingraySoftware/stingray/commit/5285c23221ed81c946007e8ae116c84fdc9947d4?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) 96.35% compared to head [(`0ae8abf`)](https://app.codecov.io/gh/StingraySoftware/stingray/pull/796?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware) 94.35%.

> :exclamation: Current head 0ae8abf differs from pull request most recent head ce2d3e6. Consider uploading reports for the commit ce2d3e6 to get more accurate results


<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main     #796      +/-   ##
==========================================
- Coverage   96.35%   94.35%   -2.01%     
==========================================
  Files          43       43              
  Lines        8784     8784              
==========================================
- Hits         8464     8288     -176     
- Misses        320      496     +176     
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/StingraySoftware/stingray/pull/796?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=StingraySoftware).

matteobachetti: Thanks @mgullik, I will work on your comments. A couple of replies:

0) In the tutorial, I read the DET_ID column and renormalize the intensities for the number of detectors that are in during the observation. However, some colleagues pointed out that there are some issues that switch off the detectors during the observation, so the filtering might be better. I don't know what's the cleanest way to handle this, but the procedure is appropriate for the specific case in the tutorial.
1) The fractional rms spectrum is in (rms/mean)^2 Hz^-1, you take out the Hz^-1 factor when you multiply by the frequency.
2) Stingray does the calibration PI/PHA-> energy automatically for a few missions, including NICER (stingray.io.rough_calibration)
matteobachetti: Changes to spectral timing notebook: https://github.com/StingraySoftware/notebooks/pull/90

matteobachetti: About failing tests: The linkcheck is failing for all stackoverflow links, I asked a question on SO but haven't received an answer yet.
The devdeps test seems to fail when importing Astropy, I'm discussing the issue on Slack and will open an issue on Github if needed.
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/power_colors.py`](https://github.com/StingraySoftware/stingray/blob/d48d3ad0e98939969eca5d4e9ac6287384da4538/stingray/power_colors.py):

> [Line 67:9](https://github.com/StingraySoftware/stingray/blob/d48d3ad0e98939969eca5d4e9ac6287384da4538/stingray/power_colors.py#L67): [E731](https://duckduckgo.com/?q=pep8%20E731) do not assign a lambda expression, use a def
> [Line 69:13](https://github.com/StingraySoftware/stingray/blob/d48d3ad0e98939969eca5d4e9ac6287384da4538/stingray/power_colors.py#L69): [E731](https://duckduckgo.com/?q=pep8%20E731) do not assign a lambda expression, use a def











Patch:
diff --git a/docs/dataexplo.rst b/docs/dataexplo.rst
index df9b98a8..f340c0a4 100644
--- a/docs/dataexplo.rst
+++ b/docs/dataexplo.rst
@@ -1,9 +1,10 @@
-Data Exploration
-****************
+More Data Exploration
+*********************
 
 These notebook tutorials show some ways to explore data with
 Stingray.
 
+
 A quick look at a NuSTAR observation
 ====================================
 
@@ -16,18 +17,6 @@ In this Tutorial, we will show an example quicklook of a NuSTAR observation.
    notebooks/DataQuickLook/Quicklook NuSTAR data with Stingray.ipynb
 
 
-Spectral timing exploration with NICER
-======================================
-
-In this Tutorial, we will show an example spectral timing exploration of a
-black hole binary using NICER data.
-
-.. toctree::
-   :maxdepth: 2
-
-   notebooks/Spectral Timing/Spectral Timing Exploration.ipynb
-
-
 Studying very slow variability with the Lomb-Scargle periodogram
 ================================================================
 
diff --git a/docs/index.rst b/docs/index.rst
index 7aa5ea11..2169dd1d 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -271,15 +271,31 @@ the stingray source directory.
 Using Stingray
 ===============
 
-Getting started
----------------
+A Spectral timing exploration
+-----------------------------
 
+In this Tutorial, we will show an example spectral timing exploration of a
+black hole binary using NICER data. The tutorial includes a hardness-intensity
+diagram, the modeling of the power density spectrum, power colors, lag-frequency,
+lag-energy, and rms/covariance spectra.
+
+.. toctree::
+   :maxdepth: 1
+
+   notebooks/Spectral Timing/Spectral Timing Exploration.ipynb
+
+
+Stingray fundamentals
+---------------------
 .. toctree::
    :maxdepth: 2
 
    core
    dataexplo
    pulsar
+   modeling
+   simulator
+   deadtime
 
 Advanced
 --------
@@ -288,9 +304,6 @@ Advanced
    :maxdepth: 2
 
    timeseries
-   modeling
-   simulator
-   deadtime
    api
 
 Additional information
diff --git a/docs/notebooks b/docs/notebooks
index d3443305..1de8b585 160000
--- a/docs/notebooks
+++ b/docs/notebooks
@@ -1,1 +1,1 @@
-Subproject commit d344330557a96105e83304fe4b931287395201cd
+Subproject commit 1de8b5854a2ee3e1f70fd1abcf6d2e7ad07b1d87
diff --git a/stingray/base.py b/stingray/base.py
index a3d4a4f7..6d9451eb 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -8,7 +8,6 @@ from collections import OrderedDict
 import pickle
 import warnings
 import copy
-import logging
 
 import numpy as np
 from astropy.table import Table
@@ -38,7 +37,6 @@ from .gti import (
     bin_intervals_from_gtis,
     time_intervals_from_gtis,
 )
-
 from typing import TYPE_CHECKING, Type, TypeVar, Union
 
 if TYPE_CHECKING:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 5a6d8486..39bc5f03 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -18,8 +18,9 @@ from .lightcurve import Lightcurve
 from .fourier import avg_cs_from_iterables, error_on_averaged_cross_spectrum
 from .fourier import avg_cs_from_timeseries, poisson_level
 from .fourier import normalize_periodograms, raw_coherence
-from .fourier import get_flux_iterable_from_segments, power_color
+from .fourier import get_flux_iterable_from_segments
 from .fourier import get_rms_from_unnorm_periodogram
+from .power_colors import power_color
 
 from scipy.special import factorial
 
diff --git a/stingray/events.py b/stingray/events.py
index 343bbb25..b1335ea0 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -4,7 +4,6 @@ Definition of :class:`EventList`.
 :class:`EventList` is used to handle photon arrival times.
 """
 
-import logging
 import warnings
 
 import numpy as np
diff --git a/stingray/fourier.py b/stingray/fourier.py
index ce0b2943..fec7500c 100644
--- a/stingray/fourier.py
+++ b/stingray/fourier.py
@@ -12,6 +12,7 @@ from .gti import (
     generate_indices_of_segment_boundaries_binned,
     generate_indices_of_segment_boundaries_unbinned,
 )
+
 from .utils import (
     fft,
     fftfreq,
@@ -105,226 +106,6 @@ def integrate_power_in_frequency_range(
     return power_integrated, power_err_integrated
 
 
-def power_color(
-    frequency,
-    power,
-    power_err=None,
-    freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
-    df=None,
-    m=1,
-    freqs_to_exclude=None,
-    poisson_power=0,
-    return_log=False,
-):
-    """
-    Calculate two power colors from a power spectrum.
-
-    Power colors are an alternative to spectral colors to understand the spectral state of an
-    accreting source. They are defined as the ratio of the power in two frequency ranges,
-    analogously to the colors calculated from electromagnetic spectra.
-
-    This function calculates two power colors, using the four frequency ranges contained
-    between the five frequency edges in ``freq_edges``. Given [f0, f1, f2, f3, f4], the
-    two power colors are calculated as the following ratios of the integrated power
-    (which are variances):
-
-    + PC0 = Var([f0, f1]) / Var([f2, f3])
-
-    + PC1 = Var([f1, f2]) / Var([f3, f4])
-
-    Errors are calculated using simple error propagation from the integrated power errors.
-
-    See Heil et al. 2015, MNRAS, 448, 3348
-
-    Parameters
-    ----------
-    frequency : iterable
-        The frequencies of the power spectrum
-    power : iterable
-        The power at each frequency
-
-    Other Parameters
-    ----------------
-    power_err : iterable
-        The power error bar at each frequency
-    freq_edges : iterable, optional, default ``[0.0039, 0.031, 0.25, 2.0, 16.0]``
-        The five edges defining the four frequency intervals to use to calculate the power color.
-        If empty, the power color is calculated using the frequencies from Heil et al. 2015.
-    df : float or float iterable, optional, default None
-        The frequency resolution of the input data. If None, it is calculated
-        from the median difference of input frequencies.
-    m : int, optional, default 1
-        The number of segments and/or contiguous frequency bins averaged to obtain power
-    freqs_to_exclude : 1-d or 2-d iterable, optional, default None
-        The ranges of frequencies to exclude from the calculation of the power color.
-        For example, the frequencies containing strong QPOs.
-        A 1-d iterable should contain two values for the edges of a single range. (E.g.
-        ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and 3-4 Hz.
-    poisson_power : float or iterable, optional, default 0
-        The Poisson noise level of the power spectrum. If iterable, it should have the same
-        length as ``frequency``. (This might apply to the case of a power spectrum with a
-        strong dead time distortion
-    return_log : bool, optional, default False
-        Return the base-10 logarithm of the power color and the errors
-
-    Returns
-    -------
-    PC0 : float
-        The first power color
-    PC0_err : float
-        The error on the first power color
-    PC1 : float
-        The second power color
-    PC1_err : float
-        The error on the second power color
-    """
-    freq_edges = np.asarray(freq_edges)
-    if len(freq_edges) != 5:
-        raise ValueError("freq_edges must have 5 elements")
-
-    frequency = np.asarray(frequency)
-    power = np.asarray(power)
-
-    if df is None:
-        df = np.median(np.diff(frequency))
-    input_frequency_low_edges = frequency - df / 2
-    input_frequency_high_edges = frequency + df / 2
-
-    if freq_edges.min() < input_frequency_low_edges[0]:
-        raise ValueError("The minimum frequency is larger than the first frequency edge")
-    if freq_edges.max() > input_frequency_high_edges[-1]:
-        raise ValueError("The maximum frequency is lower than the last frequency edge")
-
-    if power_err is None:
-        power_err = power / np.sqrt(m)
-    else:
-        power_err = np.asarray(power_err)
-
-    if freqs_to_exclude is not None:
-        if len(np.shape(freqs_to_exclude)) == 1:
-            freqs_to_exclude = [freqs_to_exclude]
-
-        if (
-            not isinstance(freqs_to_exclude, Iterable)
-            or len(np.shape(freqs_to_exclude)) != 2
-            or np.shape(freqs_to_exclude)[1] != 2
-        ):
-            raise ValueError("freqs_to_exclude must be of format [[f0, f1], [f2, f3], ...]")
-        for f0, f1 in freqs_to_exclude:
-            frequency_mask = (input_frequency_low_edges > f0) & (input_frequency_high_edges < f1)
-            idx0, idx1 = np.searchsorted(frequency, [f0, f1])
-            power[frequency_mask] = np.mean([power[idx0], power[idx1]])
-
-    var00, var00_err = integrate_power_in_frequency_range(
-        frequency,
-        power,
-        freq_edges[:2],
-        power_err=power_err,
-        df=df,
-        m=m,
-        poisson_power=poisson_power,
-    )
-    var01, var01_err = integrate_power_in_frequency_range(
-        frequency,
-        power,
-        freq_edges[2:4],
-        power_err=power_err,
-        df=df,
-        m=m,
-        poisson_power=poisson_power,
-    )
-    var10, var10_err = integrate_power_in_frequency_range(
-        frequency,
-        power,
-        freq_edges[1:3],
-        power_err=power_err,
-        df=df,
-        m=m,
-        poisson_power=poisson_power,
-    )
-    var11, var11_err = integrate_power_in_frequency_range(
-        frequency,
-        power,
-        freq_edges[3:5],
-        power_err=power_err,
-        df=df,
-        m=m,
-        poisson_power=poisson_power,
-    )
-    pc0 = var00 / var01
-    pc1 = var10 / var11
-    pc0_err = pc0 * (var00_err / var00 + var01_err / var01)
-    pc1_err = pc1 * (var10_err / var10 + var11_err / var11)
-    if return_log:
-        pc0_err = 1 / pc0 * pc0_err
-        pc1_err = 1 / pc1 * pc1_err
-        pc0 = np.log10(pc0)
-        pc1 = np.log10(pc1)
-    return pc0, pc0_err, pc1, pc1_err
-
-
-def hue_from_power_color(pc0, pc1, center=[4.51920, 0.453724]):
-    """Measure the angle of a point in the log-power color diagram wrt the center.
-
-    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
-    at -45 degrees wrt the horizontal axis.
-
-    See Heil et al. 2015, MNRAS, 448, 3348
-
-    Parameters
-    ----------
-    pc0 : float
-        The (linear, not log!) power color in the first frequency range
-    pc1 : float
-        The (linear, not log!) power color in the second frequency range
-
-    Other Parameters
-    ----------------
-    center : iterable, optional, default [4.51920, 0.453724]
-        The coordinates of the center of the power color diagram
-
-    Returns
-    -------
-    hue : float
-        The angle of the point wrt the center, in radians
-    """
-    pc0 = np.log10(pc0)
-    pc1 = np.log10(pc1)
-
-    center = np.log10(np.asarray(center))
-
-    return hue_from_logpower_color(pc0, pc1, center=center)
-
-
-def hue_from_logpower_color(log10pc0, log10pc1, center=(np.log10(4.51920), np.log10(0.453724))):
-    """Measure the angle of a point in the log-power color diagram wrt the center.
-
-    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
-    at -45 degrees wrt the horizontal axis.
-
-    See Heil et al. 2015, MNRAS, 448, 3348
-
-    Parameters
-    ----------
-    log10pc0 : float
-        The log10 power color in the first frequency range
-    log10pc1 : float
-        The log10 power color in the second frequency range
-
-    Other Parameters
-    ----------------
-    center : iterable, optional, default ``log10([4.51920, 0.453724])``
-        The coordinates of the center of the power color diagram
-
-    Returns
-    -------
-    hue : float
-        The angle of the point wrt the center, in radians
-    """
-    hue = 3 / 4 * np.pi - np.arctan2(log10pc1 - center[1], log10pc0 - center[0])
-    return hue
-
-
 def positive_fft_bins(n_bin, include_zero=False):
     """
     Give the range of positive frequencies of a complex FFT.
diff --git a/stingray/gti.py b/stingray/gti.py
index 99d7689c..4f633080 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -1,6 +1,5 @@
 import re
 import numpy as np
-import logging
 import warnings
 from collections.abc import Iterable
 import copy
diff --git a/stingray/io.py b/stingray/io.py
index 32a6b21d..e6a7d7e2 100644
--- a/stingray/io.py
+++ b/stingray/io.py
@@ -1,4 +1,3 @@
-import logging
 import math
 import copy
 import os
@@ -734,16 +733,12 @@ class EventReadOutput:
 
 
 def mkdir_p(path):  # pragma: no cover
-    """Safe ``mkdir`` function, found at [so-mkdir]_.
+    """Safe ``mkdir`` function
 
     Parameters
     ----------
     path : str
         The absolute path to the directory to be created
-
-    Notes
-    -----
-    .. [so-mkdir] http://stackoverflow.com/questions/600268/mkdir-p-functionality-in-python
     """
     import os
 
diff --git a/stingray/loggingconfig.py b/stingray/loggingconfig.py
index 72a93e71..579ebdc1 100644
--- a/stingray/loggingconfig.py
+++ b/stingray/loggingconfig.py
@@ -34,6 +34,6 @@ def setup_logger():
         formatter = CustomFormatter()
         handler.setFormatter(formatter)
         logger.addHandler(handler)
-        logger.setLevel(logging.DEBUG)
+        logger.setLevel(logging.WARNING)
 
     return logger
diff --git a/stingray/modeling/parameterestimation.py b/stingray/modeling/parameterestimation.py
index 4113c6cd..73441c01 100644
--- a/stingray/modeling/parameterestimation.py
+++ b/stingray/modeling/parameterestimation.py
@@ -1,7 +1,7 @@
 __all__ = ["OptimizationResults", "ParameterEstimation", "PSDParEst", "SamplingResults"]
 
 
-# check whether matplotlib is installed for easy plotting
+import logging
 import matplotlib.pyplot as plt
 from matplotlib.ticker import MaxNLocator
 
@@ -21,7 +21,6 @@ try:
 except ImportError:
     use_corner = False
 
-import logging
 from multiprocessing import Pool
 
 import numpy as np
@@ -46,7 +45,9 @@ from stingray.modeling.posterior import (
     logmin,
     fitter_to_model_params,
 )
-from stingray.loggingconfig import CustomFormatter
+from stingray.loggingconfig import CustomFormatter, setup_logger
+
+logger = setup_logger()
 
 
 class OptimizationResults(object):
@@ -428,7 +429,7 @@ class ParameterEstimation(object):
                 "TNC",
                 "SLSQP",
             ]:
-                logging.warning(
+                logger.warning(
                     "Fitting method %s " % self.fitmethod + "cannot incorporate the bounds you set!"
                 )
 
@@ -1213,7 +1214,7 @@ class SamplingResults(object):
 
                             ax.contour(xx, yy, zz, 7)
                         except ValueError:
-                            logging.info("Not making contours.")
+                            logger.info("Not making contours.")
 
         if save_plot:
             plt.savefig(filename, format="pdf")
@@ -1503,7 +1504,7 @@ class PSDParEst(ParameterEstimation):
                     sim_lpost1, t1, sim_lpost2, t2, neg=neg, max_post=max_post
                 )
             except RuntimeError:
-                logging.warning("Fitting was unsuccessful. " "Skipping this simulation!")
+                logger.warning("Fitting was unsuccessful. " "Skipping this simulation!")
                 continue
 
         return lrt_sim
@@ -1720,7 +1721,7 @@ class PSDParEst(ParameterEstimation):
                 max_y, maxfreq, maxind = self._compute_highest_outlier(sim_lpost, res, nmax=1)
                 max_y_all[i] = max_y[0]
             except RuntimeError:
-                logging.warning("Fitting unsuccessful! " "Skipping this simulation!")
+                logger.warning("Fitting unsuccessful! " "Skipping this simulation!")
                 continue
 
         return np.hstack(max_y_all)
diff --git a/stingray/power_colors.py b/stingray/power_colors.py
new file mode 100644
index 00000000..71b515e4
--- /dev/null
+++ b/stingray/power_colors.py
@@ -0,0 +1,495 @@
+import matplotlib.pyplot as plt
+from scipy.interpolate import interp1d
+import numpy as np
+from .fourier import integrate_power_in_frequency_range
+from collections.abc import Iterable
+
+DEFAULT_COLOR_CONFIGURATION = {
+    "center": [4.51920, 0.453724],
+    "ref_angle": 3 * np.pi / 4,
+    "state_definitions": {
+        "HSS": {"hue_limits": [300, 360], "color": "red"},
+        "LHS": {"hue_limits": [-20, 140], "color": "blue"},
+        "HIMS": {"hue_limits": [140, 220], "color": "green"},
+        "SIMS": {"hue_limits": [220, 300], "color": "yellow"},
+    },
+    "rms_spans": {
+        -20: [0.3, 0.7],
+        0: [0.3, 0.7],
+        10: [0.3, 0.6],
+        40: [0.25, 0.4],
+        100: [0.25, 0.35],
+        150: [0.2, 0.3],
+        170: [0.0, 0.3],
+        200: [0, 0.15],
+        370: [0, 0.15],
+    },
+}
+
+
+def _get_rms_span_functions(configuration=DEFAULT_COLOR_CONFIGURATION):
+    rms_spans = configuration["rms_spans"]
+
+    x = list(rms_spans.keys())
+
+    ymin = list([v[0] for v in rms_spans.values()])
+    ymax = list([v[1] for v in rms_spans.values()])
+
+    ymin_func = interp1d(x, ymin, kind="linear")
+    ymax_func = interp1d(x, ymax, kind="linear")
+    return ymin_func, ymax_func
+
+
+def _create_rms_hue_plot(polar=False, plot_spans=False, configuration=DEFAULT_COLOR_CONFIGURATION):
+    if polar:
+        fig, ax = plt.subplots(subplot_kw={"projection": "polar"})
+        ax.set_rmax(0.75)
+        ax.set_rticks([0, 0.25, 0.5, 0.75, 1])
+        ax.grid(True)
+        ax.set_rlim([0, 0.75])
+    else:
+        plt.figure()
+        plt.xlim(0, 360)
+        plt.ylim(0, 0.7)
+        plt.ylabel("Fractional rms")
+        plt.xlabel("Hue")
+        ax = plt.gca()
+
+    if not plot_spans:
+        return ax
+
+    ymin_func, ymax_func = _get_rms_span_functions(configuration)
+
+    for state in configuration["state_definitions"].keys():
+        color = configuration["state_definitions"][state]["color"]
+        xmin, xmax = configuration["state_definitions"][state]["hue_limits"]
+
+        def x_func(x):
+            return x
+
+        if polar:
+
+            def x_func(x):
+                return np.radians(x)
+
+        ax.fill_between(
+            np.linspace(x_func(xmin), x_func(xmax), 20),
+            ymin_func(np.linspace(xmin, xmax, 20)),
+            ymax_func(np.linspace(xmin, xmax, 20)),
+            color=color,
+            alpha=0.1,
+        )
+
+        if xmin < 0 and not polar:
+            ax.fill_between(
+                np.linspace(x_func(xmin + 360), x_func(360), 20),
+                ymin_func(np.linspace(xmin, 0, 20)),
+                ymax_func(np.linspace(xmin, 0, 20)),
+                color=color,
+                alpha=0.1,
+            )
+        if xmax > 360 and not polar:
+            ax.fill_between(
+                np.linspace(0, x_func(xmax - 360), 20),
+                ymin_func(np.linspace(360, xmax, 20)),
+                ymax_func(np.linspace(360, xmax, 20)),
+                color=color,
+                alpha=0.1,
+            )
+
+    return ax
+
+
+def _limit_angle_to_360(angle):
+    while angle >= 360:
+        angle = angle - 360
+    while angle < 0:
+        angle = angle + 360
+    return angle
+
+
+def _hue_line_data(center, angle, ref_angle=3 * np.pi / 4):
+    plot_angle = (-angle + ref_angle) % (np.pi * 2)
+
+    m = np.tan(plot_angle)
+    if np.isinf(m):
+        x = np.zeros_like(x) + center[0]
+        y = np.linspace(-4, 4, 20)
+    else:
+        x = np.linspace(0, 4, 20) * np.sign(np.cos(plot_angle)) + center[0]
+        y = center[1] + m * (x - center[0])
+    return x, y
+
+
+def _trace_states(ax, configuration=DEFAULT_COLOR_CONFIGURATION, **kwargs):
+    center = [np.log10(c) for c in configuration["center"]]
+    for state in configuration["state_definitions"].keys():
+        color = configuration["state_definitions"][state]["color"]
+        hue0, hue1 = configuration["state_definitions"][state]["hue_limits"]
+        hue_mean = (hue0 + hue1) / 2
+        hue_angle = (-np.radians(hue_mean) + 3 * np.pi / 4) % (np.pi * 2)
+
+        radius = 1.4
+        txt_x = radius * np.cos(hue_angle) + center[0]
+        txt_y = radius * np.sin(hue_angle) + center[1]
+        ax.text(txt_x, txt_y, state, color="k", ha="center", va="center")
+
+        hue0, hue1 = configuration["state_definitions"][state]["hue_limits"]
+        x0, y0 = _hue_line_data(center, np.radians(hue0), ref_angle=configuration["ref_angle"])
+
+        next_angle = hue0 + 5.0
+        x1, y1 = _hue_line_data(center, np.radians(hue0), ref_angle=configuration["ref_angle"])
+
+        while next_angle <= hue1:
+            x0, y0 = x1, y1
+            x1, y1 = _hue_line_data(
+                center, np.radians(next_angle), ref_angle=configuration["ref_angle"]
+            )
+            t1 = plt.Polygon(
+                [[x0[0], y0[0]], [x0[-1], y0[-1]], [x1[-1], y1[-1]]],
+                ls=None,
+                lw=0,
+                color=color,
+                **kwargs
+            )
+            ax.add_patch(t1)
+            next_angle += 5.0
+
+
+def _create_pc_plot(
+    xrange=[-2, 2],
+    yrange=[-2, 2],
+    plot_spans=False,
+    configuration=DEFAULT_COLOR_CONFIGURATION,
+):
+    """Creates an empty power color plot with labels in the right place."""
+    fig = plt.figure()
+    ax = plt.gca()
+
+    ax.set_aspect("equal")
+
+    ax.set_xlabel(r"log$_{10}$PC1")
+    ax.set_ylabel(r"log$_{10}$PC2")
+    ax.grid(False)
+
+    if not plot_spans:
+        ax.set_xlim(xrange)
+        ax.set_ylim(yrange)
+        return ax
+
+    center = np.log10(np.asarray(configuration["center"]))
+    ax.set_xlim(center[0] + np.asarray(xrange))
+    ax.set_ylim(center[1] + np.asarray(yrange))
+
+    for angle in range(0, 360, 20):
+        x, y = _hue_line_data(center, np.radians(angle), ref_angle=configuration["ref_angle"])
+
+        ax.plot(x, y, lw=0.2, ls=":", color="k", alpha=0.3, zorder=10)
+
+    ax.scatter(*center, marker="+", color="k")
+
+    limit_angles = set(
+        np.array(
+            [
+                configuration["state_definitions"][state]["hue_limits"]
+                for state in configuration["state_definitions"].keys()
+            ]
+        ).flatten()
+    )
+    limit_angles = [_limit_angle_to_360(angle) for angle in limit_angles]
+
+    for angle in limit_angles:
+        x, y = _hue_line_data(center, np.radians(angle), ref_angle=configuration["ref_angle"])
+
+        plt.plot(x, y, lw=1, ls=":", color="k", alpha=1, zorder=10)
+
+    _trace_states(ax, configuration=configuration, alpha=0.1)
+
+    return ax
+
+
+def plot_power_colors(
+    p1,
+    p1e,
+    p2,
+    p2e,
+    plot_spans=False,
+    configuration=DEFAULT_COLOR_CONFIGURATION,
+):
+    """
+    Plot power colors.
+
+    Parameters
+    ----------
+    p1 : float
+        The first power value.
+    p1e : float
+        The error in the first power value.
+    p2 : float
+        The second power value.
+    p2e : float
+        The error in the second power value.
+
+    Other Parameters
+    ----------------
+    center : (float, float)
+        The center coordinates of the plot. Default is (4.51920, 0.453724).
+    plot_spans : bool
+        Whether to plot the spans. Default is False.
+    configuration: bool
+        The color configuration to use. Default is DEFAULT_COLOR_CONFIGURATION.
+
+    Returns
+    -------
+    ax : :class:`matplotlib.AxesSubplot`
+        The matplotlib Axes object representing the power color plot.
+    """
+    p1e = 1 / p1 * p1e
+    p2e = 1 / p2 * p2e
+    p1 = np.log10(p1)
+    p2 = np.log10(p2)
+    # Create empty power color plot
+    ax = _create_pc_plot(plot_spans=plot_spans, configuration=configuration)
+    ax.errorbar(p1, p2, xerr=p1e, yerr=p2e, alpha=0.4, color="k")
+    ax.scatter(p1, p2, zorder=10, color="k")
+    return ax
+
+
+def plot_hues(
+    rms,
+    rmse,
+    pc1,
+    pc2,
+    polar=False,
+    plot_spans=False,
+    configuration=DEFAULT_COLOR_CONFIGURATION,
+):
+    hues = hue_from_power_color(pc1, pc2)
+
+    ax = _create_rms_hue_plot(polar=polar, plot_spans=plot_spans, configuration=configuration)
+    hues = hues % (np.pi * 2)
+    if not polar:
+        hues = np.degrees(hues)
+
+    ax.errorbar(hues, rms, yerr=rmse, fmt="o", alpha=0.5)
+    return ax
+
+
+def power_color(
+    frequency,
+    power,
+    power_err=None,
+    freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+    df=None,
+    m=1,
+    freqs_to_exclude=None,
+    poisson_power=0,
+    return_log=False,
+):
+    """
+    Calculate two power colors from a power spectrum.
+
+    Power colors are an alternative to spectral colors to understand the spectral state of an
+    accreting source. They are defined as the ratio of the power in two frequency ranges,
+    analogously to the colors calculated from electromagnetic spectra.
+
+    This function calculates two power colors, using the four frequency ranges contained
+    between the five frequency edges in ``freq_edges``. Given [f0, f1, f2, f3, f4], the
+    two power colors are calculated as the following ratios of the integrated power
+    (which are variances):
+
+    + PC0 = Var([f0, f1]) / Var([f2, f3])
+
+    + PC1 = Var([f1, f2]) / Var([f3, f4])
+
+    Errors are calculated using simple error propagation from the integrated power errors.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    frequency : iterable
+        The frequencies of the power spectrum
+    power : iterable
+        The power at each frequency
+
+    Other Parameters
+    ----------------
+    power_err : iterable
+        The power error bar at each frequency
+    freq_edges : iterable, optional, default ``[0.0039, 0.031, 0.25, 2.0, 16.0]``
+        The five edges defining the four frequency intervals to use to calculate the power color.
+        If empty, the power color is calculated using the frequencies from Heil et al. 2015.
+    df : float or float iterable, optional, default None
+        The frequency resolution of the input data. If None, it is calculated
+        from the median difference of input frequencies.
+    m : int, optional, default 1
+        The number of segments and/or contiguous frequency bins averaged to obtain power
+    freqs_to_exclude : 1-d or 2-d iterable, optional, default None
+        The ranges of frequencies to exclude from the calculation of the power color.
+        For example, the frequencies containing strong QPOs.
+        A 1-d iterable should contain two values for the edges of a single range. (E.g.
+        ``[0.1, 0.2]``). ``[[0.1, 0.2], [3, 4]]`` will exclude the ranges 0.1-0.2 Hz and 3-4 Hz.
+    poisson_power : float or iterable, optional, default 0
+        The Poisson noise level of the power spectrum. If iterable, it should have the same
+        length as ``frequency``. (This might apply to the case of a power spectrum with a
+        strong dead time distortion
+    return_log : bool, optional, default False
+        Return the base-10 logarithm of the power color and the errors
+
+    Returns
+    -------
+    PC0 : float
+        The first power color
+    PC0_err : float
+        The error on the first power color
+    PC1 : float
+        The second power color
+    PC1_err : float
+        The error on the second power color
+    """
+    freq_edges = np.asarray(freq_edges)
+    if len(freq_edges) != 5:
+        raise ValueError("freq_edges must have 5 elements")
+
+    frequency = np.asarray(frequency)
+    power = np.asarray(power)
+
+    if df is None:
+        df = np.median(np.diff(frequency))
+    input_frequency_low_edges = frequency - df / 2
+    input_frequency_high_edges = frequency + df / 2
+
+    if freq_edges.min() < input_frequency_low_edges[0]:
+        raise ValueError("The minimum frequency is larger than the first frequency edge")
+    if freq_edges.max() > input_frequency_high_edges[-1]:
+        raise ValueError("The maximum frequency is lower than the last frequency edge")
+
+    if power_err is None:
+        power_err = power / np.sqrt(m)
+    else:
+        power_err = np.asarray(power_err)
+
+    if freqs_to_exclude is not None:
+        if len(np.shape(freqs_to_exclude)) == 1:
+            freqs_to_exclude = [freqs_to_exclude]
+
+        if (
+            not isinstance(freqs_to_exclude, Iterable)
+            or len(np.shape(freqs_to_exclude)) != 2
+            or np.shape(freqs_to_exclude)[1] != 2
+        ):
+            raise ValueError("freqs_to_exclude must be of format [[f0, f1], [f2, f3], ...]")
+        for f0, f1 in freqs_to_exclude:
+            frequency_mask = (input_frequency_low_edges > f0) & (input_frequency_high_edges < f1)
+            idx0, idx1 = np.searchsorted(frequency, [f0, f1])
+            power[frequency_mask] = np.mean([power[idx0], power[idx1]])
+
+    var00, var00_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[:2],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var01, var01_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[2:4],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var10, var10_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[1:3],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    var11, var11_err = integrate_power_in_frequency_range(
+        frequency,
+        power,
+        freq_edges[3:5],
+        power_err=power_err,
+        df=df,
+        m=m,
+        poisson_power=poisson_power,
+    )
+    pc0 = var00 / var01
+    pc1 = var10 / var11
+    pc0_err = pc0 * (var00_err / var00 + var01_err / var01)
+    pc1_err = pc1 * (var10_err / var10 + var11_err / var11)
+    if return_log:
+        pc0_err = 1 / pc0 * pc0_err
+        pc1_err = 1 / pc1 * pc1_err
+        pc0 = np.log10(pc0)
+        pc1 = np.log10(pc1)
+    return pc0, pc0_err, pc1, pc1_err
+
+
+def hue_from_power_color(pc0, pc1, center=[4.51920, 0.453724]):
+    """Measure the angle of a point in the log-power color diagram wrt the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees wrt the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    pc0 : float
+        The (linear, not log!) power color in the first frequency range
+    pc1 : float
+        The (linear, not log!) power color in the second frequency range
+
+    Other Parameters
+    ----------------
+    center : iterable, optional, default [4.51920, 0.453724]
+        The coordinates of the center of the power color diagram
+
+    Returns
+    -------
+    hue : float
+        The angle of the point wrt the center, in radians
+    """
+    pc0 = np.log10(pc0)
+    pc1 = np.log10(pc1)
+
+    center = np.log10(np.asarray(center))
+
+    return hue_from_logpower_color(pc0, pc1, center=center)
+
+
+def hue_from_logpower_color(log10pc0, log10pc1, center=(np.log10(4.51920), np.log10(0.453724))):
+    """Measure the angle of a point in the log-power color diagram wrt the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees wrt the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348
+
+    Parameters
+    ----------
+    log10pc0 : float
+        The log10 power color in the first frequency range
+    log10pc1 : float
+        The log10 power color in the second frequency range
+
+    Other Parameters
+    ----------------
+    center : iterable, optional, default ``log10([4.51920, 0.453724])``
+        The coordinates of the center of the power color diagram
+
+    Returns
+    -------
+    hue : float
+        The angle of the point wrt the center, in radians
+    """
+    hue = 3 / 4 * np.pi - np.arctan2(log10pc1 - center[1], log10pc0 - center[0])
+    return hue
diff --git a/stingray/powerspectrum.py b/stingray/powerspectrum.py
index 4fd26243..d9d7dc17 100755
--- a/stingray/powerspectrum.py
+++ b/stingray/powerspectrum.py
@@ -77,6 +77,12 @@ class Powerspectrum(Crossspectrum):
         Note that for a single realization (``m=1``) the error is equal to the
         power.
 
+    unnorm_power: numpy.ndarray
+        The array of unnormalized powers
+
+    unnorm_power_err: numpy.ndarray
+        The uncertainties of ``unnorm_power``.
+
     df: float
         The frequency resolution.
 
@@ -458,7 +464,14 @@ class Powerspectrum(Crossspectrum):
 
     @staticmethod
     def from_events(
-        events, dt, segment_size=None, gti=None, norm="frac", silent=False, use_common_mean=True
+        events,
+        dt,
+        segment_size=None,
+        gti=None,
+        norm="frac",
+        silent=False,
+        use_common_mean=True,
+        save_all=False,
     ):
         """
         Calculate an average power spectrum from an event list.
@@ -494,6 +507,8 @@ class Powerspectrum(Crossspectrum):
             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.
         """
         if gti is None:
             gti = events.gti
@@ -505,6 +520,7 @@ class Powerspectrum(Crossspectrum):
             norm=norm,
             silent=silent,
             use_common_mean=use_common_mean,
+            save_all=save_all,
         )
 
     @staticmethod
@@ -517,6 +533,7 @@ class Powerspectrum(Crossspectrum):
         silent=False,
         use_common_mean=True,
         gti=None,
+        save_all=False,
     ):
         """Calculate AveragedPowerspectrum from a time series.
 
@@ -552,6 +569,8 @@ class Powerspectrum(Crossspectrum):
             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.
+        save_all : bool, default False
+            Save all intermediate PDSs used for the final average.
         """
         return powerspectrum_from_timeseries(
             ts,
@@ -562,11 +581,18 @@ class Powerspectrum(Crossspectrum):
             silent=silent,
             use_common_mean=use_common_mean,
             gti=gti,
+            save_all=save_all,
         )
 
     @staticmethod
     def from_lightcurve(
-        lc, segment_size=None, gti=None, norm="frac", silent=False, use_common_mean=True
+        lc,
+        segment_size=None,
+        gti=None,
+        norm="frac",
+        silent=False,
+        use_common_mean=True,
+        save_all=False,
     ):
         """
         Calculate a power spectrum from a light curve.
@@ -602,6 +628,8 @@ class Powerspectrum(Crossspectrum):
             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.
         """
         if gti is None:
             gti = lc.gti
@@ -612,11 +640,19 @@ class Powerspectrum(Crossspectrum):
             norm=norm,
             silent=silent,
             use_common_mean=use_common_mean,
+            save_all=save_all,
         )
 
     @staticmethod
     def from_lc_iterable(
-        iter_lc, dt, segment_size=None, gti=None, norm="frac", silent=False, use_common_mean=True
+        iter_lc,
+        dt,
+        segment_size=None,
+        gti=None,
+        norm="frac",
+        silent=False,
+        use_common_mean=True,
+        save_all=False,
     ):
         """
         Calculate the average power spectrum of an iterable collection of
@@ -653,6 +689,8 @@ class Powerspectrum(Crossspectrum):
             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.
         """
 
         return powerspectrum_from_lc_iterable(
@@ -663,6 +701,7 @@ class Powerspectrum(Crossspectrum):
             norm=norm,
             silent=silent,
             use_common_mean=use_common_mean,
+            save_all=save_all,
         )
 
     def _initialize_from_any_input(
@@ -801,8 +840,7 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         The array of mid-bin frequencies that the Fourier transform samples.
 
     power: numpy.ndarray
-        The array of normalized squared absolute values of Fourier
-        amplitudes.
+        The array of normalized powers
 
     power_err: numpy.ndarray
         The uncertainties of ``power``.
@@ -812,6 +850,12 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
         Note that for a single realization (``m=1``) the error is equal to the
         power.
 
+    unnorm_power: numpy.ndarray
+        The array of unnormalized powers
+
+    unnorm_power_err: numpy.ndarray
+        The uncertainties of ``unnorm_power``.
+
     df: float
         The frequency resolution.
 
@@ -824,6 +868,8 @@ class AveragedPowerspectrum(AveragedCrossspectrum, Powerspectrum):
     nphots: float
         The total number of photons in the light curve.
 
+    cs_all: list of :class:`Powerspectrum` objects
+        The list of all periodograms used to calculate the average periodogram.
     """
 
     def __init__(
diff --git a/stingray/utils.py b/stingray/utils.py
index d065d661..80090589 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -14,6 +14,9 @@ import scipy
 from numpy import histogram as histogram_np
 from numpy import histogram2d as histogram2d_np
 from numpy import histogramdd as histogramdd_np
+from .loggingconfig import setup_logger
+
+logger = setup_logger()
 
 try:
     import pyfftw
@@ -32,8 +35,8 @@ try:
 
     pyfftw.interfaces.cache.enable()
     HAS_PYFFTW = True
+    logger.info("Using PyFFTW")
 except ImportError:
-    warnings.warn("pyfftw not installed. Using standard scipy fft")
     from numpy.fft import ifft, fft, fftfreq, fftn, ifftn, fftshift, fft2, ifftshift, rfft, rfftfreq
 
     HAS_PYFFTW = False


Test patch:
diff --git a/stingray/modeling/tests/test_parameterestimation.py b/stingray/modeling/tests/test_parameterestimation.py
index 4e1f44ba..0fbbece5 100644
--- a/stingray/modeling/tests/test_parameterestimation.py
+++ b/stingray/modeling/tests/test_parameterestimation.py
@@ -11,7 +11,7 @@ from stingray import Powerspectrum, AveragedPowerspectrum
 from stingray.modeling import ParameterEstimation, PSDParEst, OptimizationResults, SamplingResults
 from stingray.modeling import PSDPosterior, set_logprior, PSDLogLikelihood, LogLikelihood
 from stingray.modeling.posterior import fitter_to_model_params
-from stingray.loggingconfig import CustomFormatter
+from stingray.loggingconfig import CustomFormatter, setup_logger
 
 try:
     from statsmodels.tools.numdiff import approx_hess
@@ -29,6 +29,8 @@ except ImportError:
 
 import matplotlib.pyplot as plt
 
+logger = setup_logger()
+
 pytestmark = pytest.mark.slow
 
 
diff --git a/stingray/tests/test_fourier.py b/stingray/tests/test_fourier.py
index 860a4152..beffcf1a 100644
--- a/stingray/tests/test_fourier.py
+++ b/stingray/tests/test_fourier.py
@@ -17,7 +17,7 @@ from stingray.fourier import normalize_periodograms, raw_coherence, estimate_int
 from stingray.fourier import bias_term, error_on_averaged_cross_spectrum, unnormalize_periodograms
 from stingray.fourier import impose_symmetry_lsft, lsft_slow, lsft_fast, rms_calculation
 from stingray.fourier import get_average_ctrate, normalize_leahy_from_variance
-from stingray.fourier import integrate_power_in_frequency_range, power_color, hue_from_power_color
+from stingray.fourier import integrate_power_in_frequency_range
 from stingray.fourier import get_rms_from_rms_norm_periodogram, get_rms_from_unnorm_periodogram
 
 from stingray.utils import check_allclose_and_print
@@ -768,82 +768,6 @@ class TestIntegration(object):
         assert np.allclose(powe, np.sqrt(2))
 
 
-class TestPowerColor(object):
-    @classmethod
-    def setup_class(cls):
-        cls.freq = np.arange(0.0001, 17, 0.00001)
-        cls.power = 1 / cls.freq
-
-    def test_power_color(self):
-        pc0, _, pc1, _ = power_color(self.freq, self.power)
-        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
-        assert np.isclose(pc0, 1)
-        assert np.isclose(pc1, 1)
-
-    def test_return_log(self):
-        pc0, _, pc1, _ = power_color(self.freq, self.power, return_log=True)
-        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
-        assert np.isclose(pc0, 0, atol=0.001)
-        assert np.isclose(pc1, 0, atol=0.001)
-
-    def test_bad_edges(self):
-        good = self.freq > 1 / 255  # the smallest frequency is 1/256
-        with pytest.raises(ValueError, match="The minimum frequency is larger "):
-            power_color(self.freq[good], self.power[good])
-
-        good = self.freq < 15  # the smallest frequency is 1/256
-        with pytest.raises(ValueError, match="The maximum frequency is lower "):
-            power_color(self.freq[good], self.power[good])
-
-        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
-            power_color(self.freq, self.power, freq_edges=[1])
-        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
-            power_color(self.freq, self.power, freq_edges=[1, 2, 3, 4, 5, 6])
-
-    def test_bad_excluded_interval(self):
-        for fte in ([1, 1.1, 3.0], [4], [[1, 1.1, 3.0]], 0, [[[1, 3]]]):
-            with pytest.raises(ValueError, match="freqs_to_exclude must be of "):
-                power_color(self.freq, self.power, freqs_to_exclude=fte)
-
-    def test_excluded_frequencies(self):
-        pc0, _, pc1, _ = power_color(self.freq, self.power, freqs_to_exclude=[1, 1.1])
-        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
-        # The excluded frequency interval is small enough that the approximation should work
-        assert np.isclose(pc0, 1, atol=0.001)
-        assert np.isclose(pc1, 1, atol=0.001)
-
-    def test_with_power_err(self):
-        pc0, pc0_err, pc1, pc1_err = power_color(
-            self.freq,
-            self.power,
-            power_err=self.power / 2,
-        )
-        pc0e, pc0e_err, pc1e, pc1e_err = power_color(
-            self.freq,
-            self.power,
-            power_err=self.power,
-        )
-        assert np.isclose(pc0, 1, atol=0.001)
-        assert np.isclose(pc1, 1, atol=0.001)
-        assert np.isclose(pc0e, 1, atol=0.001)
-        assert np.isclose(pc1e, 1, atol=0.001)
-        assert np.isclose(pc0e_err / pc0_err, 2, atol=0.001)
-        assert np.isclose(pc1e_err / pc1_err, 2, atol=0.001)
-
-    def test_hue(self):
-        center = (4.51920, 0.453724)
-        log_center = np.log10(np.asarray(center))
-        for angle in np.radians(np.arange(0, 380, 20)):
-            factor = rng.uniform(0.1, 10)
-            x = factor * np.cos(3 / 4 * np.pi - angle) + log_center[0]
-            y = factor * np.sin(3 / 4 * np.pi - angle) + log_center[1]
-            hue = hue_from_power_color(10**x, 10**y, center)
-            # Compare the angles in a safe way
-            c2 = (np.sin(hue) - np.sin(angle)) ** 2 + (np.cos(hue) - np.cos(angle)) ** 2
-            angle_diff = np.arccos((2.0 - c2) / 2.0)
-            assert np.isclose(angle_diff, 0, atol=0.001)
-
-
 class TestRMS(object):
     @classmethod
     def setup_class(cls):
diff --git a/stingray/tests/test_power_colors.py b/stingray/tests/test_power_colors.py
new file mode 100644
index 00000000..0f62e499
--- /dev/null
+++ b/stingray/tests/test_power_colors.py
@@ -0,0 +1,114 @@
+import pytest
+import numpy as np
+from stingray.power_colors import (
+    power_color,
+    hue_from_power_color,
+    plot_hues,
+    plot_power_colors,
+    DEFAULT_COLOR_CONFIGURATION,
+)
+
+rng = np.random.RandomState(1259723)
+
+
+class TestPowerColor(object):
+    @classmethod
+    def setup_class(cls):
+        cls.freq = np.arange(0.0001, 17, 0.00001)
+        cls.power = 1 / cls.freq
+        cls.pc0, cls.pc0e, cls.pc1, cls.pc1e = power_color(cls.freq, cls.power)
+        cls.lpc0, _, cls.lpc1, _ = power_color(cls.freq, cls.power, return_log=True)
+        cls.rms, cls.rmse = 0.1, 0.01
+        cls.configuration = DEFAULT_COLOR_CONFIGURATION
+        cls.configuration["state_definitions"]["HSS"]["hue_limits"] = [300, 370]
+
+    def test_power_color(self):
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        assert np.isclose(self.pc0, 1)
+        assert np.isclose(self.pc1, 1)
+
+    def test_return_log(self):
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        assert np.isclose(self.lpc0, 0, atol=0.001)
+        assert np.isclose(self.lpc1, 0, atol=0.001)
+
+    def test_bad_edges(self):
+        good = self.freq > 1 / 255  # the smallest frequency is 1/256
+        with pytest.raises(ValueError, match="The minimum frequency is larger "):
+            power_color(self.freq[good], self.power[good])
+
+        good = self.freq < 15  # the smallest frequency is 1/256
+        with pytest.raises(ValueError, match="The maximum frequency is lower "):
+            power_color(self.freq[good], self.power[good])
+
+        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
+            power_color(self.freq, self.power, freq_edges=[1])
+        with pytest.raises(ValueError, match="freq_edges must have 5 elements"):
+            power_color(self.freq, self.power, freq_edges=[1, 2, 3, 4, 5, 6])
+
+    def test_bad_excluded_interval(self):
+        for fte in ([1, 1.1, 3.0], [4], [[1, 1.1, 3.0]], 0, [[[1, 3]]]):
+            with pytest.raises(ValueError, match="freqs_to_exclude must be of "):
+                power_color(self.freq, self.power, freqs_to_exclude=fte)
+
+    def test_excluded_frequencies(self):
+        pc0, _, pc1, _ = power_color(self.freq, self.power, freqs_to_exclude=[1, 1.1])
+        # The colors calculated with these frequency edges on a 1/f spectrum should be 1
+        # The excluded frequency interval is small enough that the approximation should work
+        assert np.isclose(pc0, 1, atol=0.001)
+        assert np.isclose(pc1, 1, atol=0.001)
+
+    def test_with_power_err(self):
+        pc0, pc0_err, pc1, pc1_err = power_color(
+            self.freq,
+            self.power,
+            power_err=self.power / 2,
+        )
+        pc0e, pc0e_err, pc1e, pc1e_err = power_color(
+            self.freq,
+            self.power,
+            power_err=self.power,
+        )
+        assert np.isclose(pc0, 1, atol=0.001)
+        assert np.isclose(pc1, 1, atol=0.001)
+        assert np.isclose(pc0e, 1, atol=0.001)
+        assert np.isclose(pc1e, 1, atol=0.001)
+        assert np.isclose(pc0e_err / pc0_err, 2, atol=0.001)
+        assert np.isclose(pc1e_err / pc1_err, 2, atol=0.001)
+
+    def test_hue(self):
+        center = (4.51920, 0.453724)
+        log_center = np.log10(np.asarray(center))
+        for angle in np.radians(np.arange(0, 380, 20)):
+            factor = rng.uniform(0.1, 10)
+            x = factor * np.cos(3 / 4 * np.pi - angle) + log_center[0]
+            y = factor * np.sin(3 / 4 * np.pi - angle) + log_center[1]
+            hue = hue_from_power_color(10**x, 10**y, center)
+            # Compare the angles in a safe way
+            c2 = (np.sin(hue) - np.sin(angle)) ** 2 + (np.cos(hue) - np.cos(angle)) ** 2
+            angle_diff = np.arccos((2.0 - c2) / 2.0)
+            assert np.isclose(angle_diff, 0, atol=0.001)
+
+    @pytest.mark.parametrize("plot_spans", [True, False])
+    def test_plot_color(self, plot_spans):
+        plot_power_colors(
+            self.pc0,
+            self.pc0e,
+            self.pc1,
+            self.pc1e,
+            plot_spans=plot_spans,
+            configuration=self.configuration,
+        )
+
+    @pytest.mark.parametrize("plot_spans", [True, False])
+    @pytest.mark.parametrize("polar", [True, False])
+    def test_hues(self, plot_spans, polar):
+        plot_hues(
+            self.rms,
+            self.rmse,
+            self.pc0,
+            self.pc1,
+            plot_spans=plot_spans,
+            configuration=self.configuration,
+            polar=polar,
+        )