diff --git a/Project.toml b/Project.toml
index 236449b..7d8efdb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -4,19 +4,28 @@ authors = ["Aman Pandey"]
 version = "0.1.0"
 
 [deps]
-FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-ResumableFunctions = "c5292f4c-5179-55e1-98c5-05642aab7184"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
+AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 FITSIO = "525bcba6-941b-5504-bd06-fd0dc1a4d2eb"
+HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
+NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568"
+ResumableFunctions = "c5292f4c-5179-55e1-98c5-05642aab7184"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+DataFrames = "1.3"
+Distributions = "0.25"
 julia = "1.11"
 FFTW = "1.4"
 Distributions = "0.25"
@@ -26,11 +35,21 @@ ProgressBars = "1.4"
 DataFrames = "1.3"
 HDF5 = "0.16"
 FITSIO = "0.16"
+HDF5 = "0.16"
 Intervals = "1.8"
+JuliaFormatter = "1.0.62"
 NaNMath = "0.3, 1"
+Plots = "1.8"
+ProgressBars = "1.4"
+ResumableFunctions = "0.6"
+StatsBase = "0.33"
+Test = "1.11.0"
+julia = "1.7"
 
 [extras]
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+ResumableFunctions = "c5292f4c-5179-55e1-98c5-05642aab7184"
 
 [targets]
-test = ["Test"]
+test = ["Test", "Random","ResumableFunctions"]
diff --git a/docs/src/index.md b/docs/src/index.md
index 1a3c5c3..f7aaf79 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -12,3 +12,29 @@ Documentation for [Stingray](https://github.com/matteobachetti/Stingray.jl).
 ```@autodocs
 Modules = [Stingray]
 ```
+## Overview
+**Stingray.jl** is a Julia package for **X-ray spectral timing analysis**, designed for high-performance astrophysical data processing.
+
+
+## Get Started
+🚀 **New to Stingray.jl?** Follow these steps:
+
+1. **[Installation Guide](installation.md)** → Set up Stingray.jl.
+2. **[Usage Guide](usage.md)** → Learn how to analyze time-series data.
+
+## Quick Example
+Here's a simple example to analyze an X-ray light curve:
+```julia
+using Stingray, Random, Plots
+
+# Generate Simulated Light Curve
+N = 1024  # Number of data points
+t = collect(0:0.1:(N-1)*0.1)  # Time array
+Random.seed!(42)
+light_curve = sin.(2π * 0.5 .* t) + 0.3 * randn(N)  # Sine wave + noise
+
+# Plot Light Curve
+plot(t, light_curve, label="Simulated Light Curve", xlabel="Time (s)", ylabel="Intensity", title="X-ray Light Curve", legend=:topright)
+```
+
+🚀 **Start exploring Stingray.jl today!**
\ No newline at end of file
diff --git a/docs/src/installation.md b/docs/src/installation.md
new file mode 100644
index 0000000..2fb922b
--- /dev/null
+++ b/docs/src/installation.md
@@ -0,0 +1,196 @@
+```@meta
+CurrentModule = Stingray
+```
+
+# Stingray
+
+Documentation for [Stingray](https://github.com/matteobachetti/Stingray.jl).
+
+```@index
+```
+
+```@autodocs
+Modules = [Stingray]
+```
+# Stingray.jl
+
+Stingray.jl is a spectral-timing software package for time series analysis of astronomical data, with a particular focus on X-ray astronomy.
+
+## Installation
+
+### Prerequisites
+
+Before installing Stingray.jl, ensure you have Julia v1.7 or later installed:
+
+1. Download Julia from the [official website](https://julialang.org/downloads/)
+2. Follow the installation instructions for your operating system
+3. Verify your installation by running `julia --version` in your terminal
+
+### Installing from the Julia General Registry
+
+To install the stable version of Stingray.jl:
+
+```julia
+using Pkg
+Pkg.add("Stingray")
+```
+
+### Installing the Development Version
+
+For the latest features and improvements, you can install directly from the GitHub repository:
+
+```julia
+using Pkg
+Pkg.add(url="https://github.com/matteobachetti/Stingray.jl")
+```
+
+Or switch to development mode if you plan to contribute:
+
+```julia
+using Pkg
+Pkg.develop("Stingray")
+```
+
+### Verifying Installation
+
+Verify your installation works correctly:
+
+```julia
+using Stingray
+@info "Stingray.jl loaded successfully!"
+```
+
+## Testing
+
+Stingray.jl includes a comprehensive test suite to ensure functionality. To run the tests:
+
+```julia
+using Pkg
+Pkg.test("Stingray")
+```
+
+For contributors working on a local development version:
+
+```julia
+# Navigate to the package directory
+cd(joinpath(Pkg.devdir(), "Stingray"))
+
+# Run tests
+using Pkg
+Pkg.test()
+```
+
+### Writing New Tests
+
+When adding features or fixing bugs, please include appropriate tests:
+
+1. Tests should be placed in the `test/` directory
+2. Use the `@testset` macro to organize related tests
+3. Follow existing naming conventions (e.g., `test_[module_name].jl` or `test_[feature].jl`)
+4. Run tests locally before submitting a pull request
+
+Example test structure:
+
+```julia
+@testset "Module Name Tests" begin
+    @test function_to_test(input) == expected_output
+    @test_throws ErrorType function_to_test(invalid_input)
+end
+```
+
+## Documentation
+
+### Building Documentation Locally
+
+To build and preview the documentation locally:
+
+1. Install required documentation packages:
+
+```julia
+using Pkg
+Pkg.add(["Documenter", "DocumenterTools"])
+```
+
+2. Navigate to the docs directory and build:
+
+```julia
+cd("docs")
+julia --project=. make.jl
+```
+
+3. Open `docs/build/index.html` in your browser
+
+### Writing Documentation
+
+Documentation uses [Documenter.jl](https://juliadocs.github.io/Documenter.jl/stable/). When adding new features:
+
+1. Add docstrings to functions, types, and modules following Julia's documentation guidelines
+2. Use `"""..."""` for multi-line docstrings
+3. Include Examples section with working code examples
+4. For complex features, add a dedicated page in `docs/src/`
+
+Example docstring:
+
+```julia
+"""
+    power_color(frequency, power; kwargs...)
+
+Calculate two power colors from a power spectrum.
+
+# Arguments
+- `frequency::Vector{Number}`: The frequencies of the power spectrum
+- `power::Vector{Number}`: The power at each frequency
+
+# Keywords
+- `power_err=nothing`: Optional error measurements for power
+- `freq_edges=[1/256, 1/32, 0.25, 2.0, 16.0]`: Frequency edges for calculations
+- `return_log=false`: Return the base-10 logarithm of results
+
+# Returns
+- `pc0::Number`: The first power color
+- `pc0_err::Number`: Error on the first power color
+- `pc1::Number`: The second power color
+- `pc1_err::Number`: Error on the second power color
+
+# Examples
+```julia
+freq = collect(0.01:0.01:20)
+power = 1 ./ freq
+pc0, pc0_err, pc1, pc1_err = power_color(freq, power)
+```
+
+See Heil et al. 2015, MNRAS, 448, 3348 for details.
+"""
+function power_color(frequency, power; kwargs...)
+    # Implementation
+end
+```
+
+## Example Notebooks
+
+Example Jupyter notebooks demonstrating Stingray.jl functionality are available in the [notebooks repository](https://github.com/kashish2210/notebooks).
+
+These examples cover:
+- Basic time series analysis
+- Power spectral calculations
+- Cross-spectral techniques
+- Specialized techniques for X-ray astronomy
+
+## Contributing
+
+We welcome contributions! Please check the [CONTRIBUTING.md](CONTRIBUTING.md) file for guidelines.
+
+1. Fork the repository
+2. Create a feature branch (`git checkout -b feature/amazing-feature`)
+3. Implement your changes with tests and documentation
+4. Commit your changes (`git commit -m 'Add amazing feature'`)
+5. Push to your branch (`git push origin feature/amazing-feature`)
+6. Open a Pull Request
+
+## License
+
+This project is licensed under the MIT License - see the [LICENSE](LICENSE) file for details.
+
+## Acknowledgments
+
+Stingray.jl is a port of the Python [Stingray](https://github.com/StingraySoftware/stingray) package with optimizations for Julia.
\ No newline at end of file
diff --git a/docs/src/usage.md b/docs/src/usage.md
new file mode 100644
index 0000000..aad8cf5
--- /dev/null
+++ b/docs/src/usage.md
@@ -0,0 +1,135 @@
+```@meta
+CurrentModule = Stingray
+```
+
+# Stingray
+
+Documentation for [Stingray](https://github.com/matteobachetti/Stingray.jl).
+
+```@index
+```
+
+```@autodocs
+Modules = [Stingray]
+```
+
+## Installing Julia
+Before installing Stingray.jl, ensure you have **Julia v1.7 or later** installed. If you haven’t installed Julia yet, follow these steps:
+1. Download Julia from the official site: [https://julialang.org/downloads/](https://julialang.org/downloads/)
+2. Follow the installation instructions for your operating system.
+3. Add Julia to your system’s PATH for easy command-line access.
+
+## Installing Stingray.jl
+To install Stingray.jl, open a Julia REPL (press `]` to enter the package manager) and run:
+
+```julia
+using Pkg
+Pkg.add("Stingray")
+```
+
+This will download and install all necessary dependencies.
+
+## Cloning the Development Version
+If you want the latest features, clone the GitHub repository and work in development mode:
+
+```julia
+Pkg.add(url="https://github.com/matteobachetti/Stingray.jl")
+Pkg.develop("Stingray")
+```
+
+## Checking Your Installation
+To verify the installation, run the following in the Julia REPL:
+
+```julia
+using Stingray
+@info "Stingray.jl loaded successfully!"
+```
+
+# Usage Guide for Stingray.jl
+
+Documentation for [Stingray.jl](https://github.com/matteobachetti/Stingray.jl).
+
+```@index
+```
+
+## Overview
+Stingray.jl provides powerful tools for **X-ray spectral timing analysis**, including **Fourier transforms, periodograms, and coherence calculations**. This guide covers the basic usage of the package.
+
+## Importing the Package
+After installation, load Stingray.jl in your Julia session:
+```julia
+using Stingray
+```
+
+## Creating a Simulated Light Curve
+```julia
+using Random, Plots
+
+# Generate Simulated Light Curve
+N = 1024  # Number of data points
+t = collect(0:0.1:(N-1)*0.1)  # Time array
+Random.seed!(42)
+light_curve = sin.(2π * 0.5 .* t) + 0.3 * randn(N)  # Sine wave + noise
+
+# Plot Light Curve
+plot(t, light_curve, label="Simulated Light Curve", xlabel="Time (s)", ylabel="Intensity", title="X-ray Light Curve", legend=:topright)
+```
+
+## Computing a Power Spectrum
+```julia
+using FFTW
+
+# Compute Fourier Transform
+frequencies = fftshift(fftfreq(N, 0.1))  # Compute frequency bins
+power_spectrum = abs.(fftshift(fft(light_curve))).^2  # Compute power spectrum
+
+# Plot Power Spectrum
+plot(frequencies, power_spectrum, xlabel="Frequency (Hz)", ylabel="Power", title="Power Spectrum", legend=false)
+```
+
+```@autodocs
+Modules = [Stingray]
+```
+
+## Power Colors
+
+### Overview
+
+The `power_color` function is a key component in Stingray.jl, designed to analyze frequency-power distributions.
+
+### Function: `power_color`
+
+```@autodocs
+Modules = [Stingray]
+Pages = ["power_color.jl"]
+```
+
+### **Description**
+The `power_color` function calculates power color indices based on frequency and power data.
+
+### **Usage**
+```julia
+pc0, pc0e, pc1, pc1e = power_color(freq, power)
+```
+
+### **Arguments**
+- `freq::Vector{Float64}`: The frequency values.
+- `power::Vector{Float64}`: The corresponding power values.
+- `return_log::Bool=false` (optional): Returns logarithmic power color if `true`.
+
+### **Returns**
+- `pc0::Float64`: Power color index 0.
+- `pc0e::Float64`: Error in pc0.
+- `pc1::Float64`: Power color index 1.
+- `pc1e::Float64`: Error in pc1.
+
+### **Example**
+```julia
+freq = collect(0.0001:0.00001:17)
+power = 1 ./ freq
+pc0, pc0e, pc1, pc1e = power_color(freq, power)
+println("Power Color Indices: ", pc0, pc0e, pc1, pc1e)
+```
+
+# Temporary Reference
+# Refer to this repository for more content: [GitHub - Notebooks](https://github.com/kashish2210/notebooks)
diff --git a/src/Stingray.jl b/src/Stingray.jl
index d8bd360..3512cdb 100644
--- a/src/Stingray.jl
+++ b/src/Stingray.jl
@@ -1,6 +1,8 @@
 module Stingray
 
 using ResumableFunctions, StatsBase, Statistics, DataFrames
+using ProgressBars: tqdm as show_progress
+using plots, DataFrames, StatsPlots, Colors, LinearAlgebra, Interpolations
 using FFTW, NaNMath, FITSIO, Intervals
 using ProgressBars: tqdm as show_progress 
 
@@ -32,4 +34,13 @@ export bin_intervals_from_gtis
 
 include("utils.jl")
 
-end 
+include("power_colors.jl")
+export plot_power_colors
+export hue_from_power_color
+export plot_hues
+export integrate_power_in_frequency_range
+export power_color
+export hue_from_power_color
+export hue_from_logpower_color
+
+end
diff --git a/src/power_colors.jl b/src/power_colors.jl
new file mode 100644
index 0000000..40b571d
--- /dev/null
+++ b/src/power_colors.jl
@@ -0,0 +1,547 @@
+using Plots
+using Interpolations
+using LinearAlgebra
+using Colors
+using Statistics
+using StatsPlots
+using DataFrames
+
+DEFAULT_COLOR_CONFIGURATION = Dict(
+    "center" => [4.51920, 0.453724],
+    "ref_angle" => 3 * pi / 4,
+    "state_definitions" => Dict(
+        "HSS" => Dict("hue_limits" => [300, 360], "color" => "red"),
+        "LHS" => Dict("hue_limits" => [-20, 140], "color" => "blue"),
+        "HIMS" => Dict("hue_limits" => [140, 220], "color" => "green"),
+        "SIMS" => Dict("hue_limits" => [220, 300], "color" => "yellow")
+    ),
+    "rms_spans" => Dict(
+        -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]
+    )
+)
+
+function _get_rms_span_functions(configuration=DEFAULT_COLOR_CONFIGURATION)
+    rms_spans = configuration["rms_spans"]
+
+    # Extract and sort unique x-values
+    x = sort(unique(collect(keys(rms_spans))))
+    ymin = [rms_spans[x_val][1] for x_val in x]
+    ymax = [rms_spans[x_val][2] for x_val in x]
+
+    # Create interpolation functions
+    ymin_func = LinearInterpolation(x, ymin, extrapolation_bc=Flat())
+    ymax_func = LinearInterpolation(x, ymax, extrapolation_bc=Flat())
+
+    return ymin_func, ymax_func
+end
+function _create_rms_hue_plot(; polar=false, plot_spans=false, configuration=DEFAULT_COLOR_CONFIGURATION)
+    if polar
+        fig = plot(proj=:polar)
+        ylims!(fig, (0.0, 0.75))
+        yticks!(fig, [0, 0.25, 0.5, 0.75, 1])
+        plot!(fig, grid=true)
+    else
+        fig = plot(xlims=(0, 360), ylims=(0.0, 0.7), xlabel="Hue", ylabel="Fractional rms")
+    end
+    if !plot_spans
+        return fig
+    end
+    ymin_func, ymax_func = get_rms_span_functions(configuration)
+    for (state, details) in configuration["state_definitions"]
+        color = parse(Colorant, details["color"])
+        xmin, xmax = details["hue_limits"]
+        
+        xs = collect(range(xmin, stop=xmax, length=20))
+        
+        if polar
+            # For polar plots, handle each point individually to avoid broadcasting issues
+            xs_rad = [deg2rad(x) for x in xs]
+            y_mins = [ymin_func(x) for x in xs]
+            y_maxs = [ymax_func(x) for x in xs]
+            
+            # Plot the minimum line
+            plot!(fig, xs_rad, y_mins, color=color, label="", alpha=0.5)
+            # Plot the maximum line
+            plot!(fig, xs_rad, y_maxs, color=color, label="", alpha=0.5)
+            # Fill between the lines
+            for i in 1:(length(xs)-1)
+                x_section = [xs_rad[i], xs_rad[i+1], xs_rad[i+1], xs_rad[i]]
+                y_section = [y_mins[i], y_mins[i+1], y_maxs[i+1], y_maxs[i]]
+                plot!(fig, x_section, y_section, seriestype=:shape, color=color, fillalpha=0.1, linewidth=0, label="")
+            end
+        else
+            # For Cartesian plots, standard approach works
+            ymin_vals = [ymin_func(x) for x in xs]
+            ymax_vals = [ymax_func(x) for x in xs]
+            # Calculate ribbon values explicitly
+            ribbon_vals = [ymax - ymin for (ymax, ymin) in zip(ymax_vals, ymin_vals)]
+            plot!(fig, xs, ymin_vals, ribbon=ribbon_vals, fillalpha=0.1, color=color, label="")
+        end
+    end
+    return fig
+end
+
+function _limit_angle_to_360(angle)
+    while angle >= 360
+        angle -= 360
+    end
+    while angle < 0
+        angle += 360
+    end
+    return angle
+end
+
+function _hue_line_data(center, angle; ref_angle=3 * pi / 4)
+    plot_angle = mod(-angle + ref_angle, 2 * pi)
+
+    m = tan(plot_angle)
+    if isinf(m)
+        x = fill(center[1], 20)
+        y = collect(range(-4, stop=4, length=20))
+    else
+        x = collect(range(0, stop=4, length=20) .* sign(cos(plot_angle)) .+ center[1])
+        y = collect(center[2] .+ m .* (x .- center[1]))
+    end
+    return x, y
+end
+
+function _trace_states(ax, configuration=DEFAULT_COLOR_CONFIGURATION; kwargs...)
+    center = log10.(configuration["center"])
+    for (state, details) in configuration["state_definitions"]
+        color = details["color"]
+        hue0, hue1 = details["hue_limits"]
+        hue_mean = (hue0 + hue1) / 2
+        hue_angle = mod(-deg2rad(hue_mean) + 3 * pi / 4, 2 * pi)
+
+        radius = 1.4
+        txt_x = radius * cos(hue_angle) + center[1]
+        txt_y = radius * sin(hue_angle) + center[2]
+        annotate!(ax, txt_x, txt_y, text(state, :black, :center))
+
+        x0, y0 = _hue_line_data(center, deg2rad(hue0); ref_angle=configuration["ref_angle"])
+        next_angle = hue0 + 5.0
+        x1, y1 = _hue_line_data(center, deg2rad(hue0); ref_angle=configuration["ref_angle"])
+
+        while next_angle <= hue1
+            x0, y0 = x1, y1
+            x1, y1 = _hue_line_data(center, deg2rad(next_angle); ref_angle=configuration["ref_angle"])
+            
+            plot!(ax, [x0[1], x0[end], x1[end], x0[1]], [y0[1], y0[end], y1[end], y0[1]], seriestype=:shape, color=color, lw=0, label="")
+            next_angle += 5.0
+        end
+    end
+end
+
+function _create_pc_plot(; xrange=(-2, 2), yrange=(-2, 2), plot_spans=false, configuration=DEFAULT_COLOR_CONFIGURATION)
+    fig = plot()
+    xlabel!("log_{10}PC1")
+    ylabel!("log_{10}PC2")
+    
+    if !plot_spans
+        xlims!(xrange...)
+        ylims!(yrange...)
+        return fig
+    end
+    
+    center = log10.(configuration["center"])
+    xlims!(center[1] + xrange[1], center[1] + xrange[2])
+    ylims!(center[2] + yrange[1], center[2] + yrange[2])
+    
+    for angle in 0:20:360
+        x, y = _hue_line_data(center, deg2rad(angle); ref_angle=configuration["ref_angle"])
+        plot!(x, y, lw=0.2, ls=:dot, color=:black, alpha=0.3)
+    end
+    
+    scatter!([center[1]], [center[2]], marker=:+, color=:black)
+    
+    limit_angles = Set(vcat([configuration["state_definitions"][state]["hue_limits"] for state in keys(configuration["state_definitions"]) ]...))
+    limit_angles = [_limit_angle_to_360(angle) for angle in limit_angles]
+    
+    for angle in limit_angles
+        x, y = _hue_line_data(center, deg2rad(angle); ref_angle=configuration["ref_angle"])
+        plot!(x, y, lw=1, ls=:dot, color=:black, alpha=1)
+    end
+    
+    _trace_states(fig, configuration=configuration, alpha=0.1)
+    
+    return fig
+end
+
+
+function plot_power_colors(
+    p1::Number,
+    p1e::Number,
+    p2::Number,
+    p2e::Number;
+    plot_spans::Bool=false,
+    configuration=DEFAULT_COLOR_CONFIGURATION
+    )
+    """
+    Plot power colors.
+
+    Parameters
+    ----------
+    p1 : Number
+        The first power value.
+    p1e : Number
+        The error in the first power value.
+    p2 : Number
+        The second power value.
+    p2e : Number
+        The error in the second power value.
+
+    Other Parameters
+    ----------------
+    center : (Number, Number)
+        The center coordinates of the plot. Default is (4.51920, 0.453724).
+    plot_spans : Bool
+        Whether to plot the spans. Default is false.
+    configuration: Dict
+        The color configuration to use. Default is DEFAULT_COLOR_CONFIGURATION.
+
+    Returns
+    -------
+    fig : Plot
+        The Julia Plot object representing the power color plot.
+    """
+    p1e = 1 / p1 * p1e
+    p2e = 1 / p2 * p2e
+    p1 = log10(p1)
+    p2 = log10(p2)
+    
+    fig = _create_pc_plot(plot_spans=plot_spans, configuration=configuration)
+    
+    scatter!([p1], [p2], marker=:circle, color=:black)
+    
+    df = DataFrame(x=[p1], y=[p2], xerr=[p1e], yerr=[p2e])
+    
+    # Use error bars correctly without `permute`
+    plot!(df.x, df.y, ribbon=(df.yerr, df.yerr), lw=0, alpha=0.4, color=:black)
+    plot!(df.x, df.y, ribbon=(df.xerr, df.xerr), lw=0, alpha=0.4, color=:black)
+    
+    return fig
+end
+
+function hue_from_power_color(pc1, pc2)
+    return atan.(pc2, pc1)  # Compute hue as the angle from power color components
+end
+
+function plot_hues(
+    rms, 
+    rmse, 
+    pc1,
+    pc2; 
+    polar=false, 
+    plot_spans=false, 
+    configuration=DEFAULT_COLOR_CONFIGURATION
+    )
+    """
+    Plot hues.
+
+    Parameters
+    ----------
+    rms : Array{Number}
+        The root mean square values.
+    rmse : Array{Number}
+        The errors in the root mean square values.
+    pc1 : Array{Number}
+        The first power color component.
+    pc2 : Array{Number}
+        The second power color component.
+
+    Other Parameters
+    ----------------
+    polar : Bool
+        Whether to use a polar plot. Default is false.
+    plot_spans : Bool
+        Whether to plot the spans. Default is false.
+    configuration: Dict
+        The color configuration to use. Default is DEFAULT_COLOR_CONFIGURATION.
+
+    Returns
+    -------
+    ax : Plot
+        The Julia Plot object representing the hues plot.
+    """
+    hues = hue_from_power_color(collect(pc1), collect(pc2))  # Ensure arrays
+    
+    ax = _create_rms_hue_plot(polar=polar, plot_spans=plot_spans, configuration=configuration)
+    hues = mod.(hues, 2 * pi)
+    if !polar
+        hues = rad2deg.(hues)
+    end
+    
+    scatter!(hues, rms, yerror=rmse, alpha=0.5)  # Correct error bars
+    
+    return ax
+end
+
+function integrate_power_in_frequency_range(
+    frequency,
+    power,
+    freq_range;
+    power_err=nothing,
+    df=nothing,
+    m=1,
+    poisson_power=0
+)
+    """
+    Integrate the power over a frequency range.
+
+    Parameters
+    ----------
+    frequency : Vector{Number}
+        The frequency array.
+    power : Vector{Number}
+        The power at each frequency.
+    freq_range : Tuple{Number, Number} or Vector{Number}
+        The frequency range over which to integrate.
+
+    Other Parameters
+    ----------------
+    power_err : Vector{Number}, optional
+        The error on the power.
+    df : Number, optional
+        The frequency resolution.
+    m : Int, optional
+        The number of power spectrum averages.
+    poisson_power : Number, optional
+        The Poisson noise level.
+
+    Returns
+    -------
+    variance : Number
+        The integrated power (variance) in the frequency range.
+    variance_err : Number
+        The error on the integrated power.
+    """
+    frequency = collect(frequency)
+    power = collect(power)
+    
+    # Ensure freq_range is a two-element vector
+    if length(freq_range) != 2
+        throw(ArgumentError("freq_range must have exactly 2 elements"))
+    end
+    
+    f0, f1 = freq_range
+
+    # Calculate frequency resolution if not provided
+    df = isnothing(df) ? median(diff(frequency)) : df
+    
+    # Define frequency range including bin width
+    input_frequency_low_edges = frequency .- df / 2
+    input_frequency_high_edges = frequency .+ df / 2
+    
+    # Find indices of frequencies within the required range
+    idx = findall((input_frequency_high_edges .>= f0) .& (input_frequency_low_edges .<= f1))
+    
+    if isempty(idx)
+        throw(ArgumentError("No frequencies found in the range ($f0, $f1)"))
+    end
+    
+    # Get power values within range
+    power_in_range = power[idx]
+    
+    # Handle errors if provided
+    power_err_in_range = isnothing(power_err) ? power_in_range ./ sqrt(m) : power_err[idx]
+    
+    # Calculate total variance (integrated power)
+    # For power spectral density, integration is multiplication by df
+    variance = sum(power_in_range .- poisson_power) * df
+    
+    # Error propagation for the sum
+    variance_err = sqrt(sum(power_err_in_range.^2)) * df
+    
+    return variance, variance_err
+end
+
+function power_color(
+    frequency,
+    power;
+    power_err=nothing,
+    freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],
+    df=nothing,
+    m=1,
+    freqs_to_exclude=nothing,
+    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 : Vector{Number}
+        The frequencies of the power spectrum.
+    power : Vector{Number}
+        The power at each frequency.
+
+    Other Parameters
+    ----------------
+    power_err : Vector{Number}, optional
+        The power error bar at each frequency.
+    freq_edges : Vector{Number}, optional
+        The five edges defining the four frequency intervals to use to calculate the power color.
+    df : Number, optional
+        The frequency resolution of the input data. If `nothing`, it is calculated from the median difference of input frequencies.
+    m : Int, optional
+        The number of segments and/or contiguous frequency bins averaged to obtain power.
+    freqs_to_exclude : Vector{Tuple{Number, Number}}, optional
+        The ranges of frequencies to exclude from the calculation of the power color.
+    poisson_power : Number, optional
+        The Poisson noise level of the power spectrum.
+    return_log : Bool, optional
+        Return the base-10 logarithm of the power color and the errors.
+
+    Returns
+    -------
+    PC0 : Number
+        The first power color.
+    PC0_err : Number
+        The error on the first power color.
+    PC1 : Number
+        The second power color.
+    PC1_err : Number
+        The error on the second power color.
+    """
+    freq_edges = collect(freq_edges)
+    length(freq_edges) == 5 || throw(ArgumentError("freq_edges must have 5 elements"))
+
+    frequency = collect(frequency)
+    power = collect(power)
+
+    df = isnothing(df) ? median(diff(frequency)) : df
+    input_frequency_low_edges = frequency .- df / 2
+    input_frequency_high_edges = frequency .+ df / 2
+
+    # Fix the error check - ensure freq_edges are within the range of input frequencies
+    if minimum(freq_edges) < minimum(input_frequency_low_edges)
+        throw(ArgumentError("The minimum frequency edge is lower than the available frequencies"))
+    end
+    
+    if maximum(freq_edges) > maximum(input_frequency_high_edges)
+        throw(ArgumentError("The maximum frequency edge is higher than the available frequencies"))
+    end
+
+    power_err = isnothing(power_err) ? power ./ sqrt(m) : collect(power_err)
+
+    if !isnothing(freqs_to_exclude)
+        for (f0, f1) in freqs_to_exclude
+            frequency_mask = (input_frequency_low_edges .> f0) .& (input_frequency_high_edges .< f1)
+            idx0, idx1 = searchsortedfirst(frequency, f0), searchsortedlast(frequency, f1)
+            power[frequency_mask] .= mean([power[idx0], power[idx1]])
+        end
+    end
+
+    var00, var00_err = integrate_power_in_frequency_range(frequency, power, freq_edges[1: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[3: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[2: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[4: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 = log10(pc0)
+        pc1 = log10(pc1)
+    end
+
+    return pc0, pc0_err, pc1, pc1_err
+end
+
+function hue_from_power_color(pc0, pc1; center=[4.51920, 0.453724])
+    """
+    Measure the angle of a point in the log-power color diagram with respect to the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees with respect to the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348.
+
+    Parameters
+    ----------
+    pc0 : Number
+        The (linear, not log!) power color in the first frequency range.
+    pc1 : Number
+        The (linear, not log!) power color in the second frequency range.
+
+    Other Parameters
+    ----------------
+    center : Vector{Number}, optional, default [4.51920, 0.453724]
+        The coordinates of the center of the power color diagram.
+
+    Returns
+    -------
+    hue : Number
+        The angle of the point with respect to the center, in radians.
+    """
+    # Handle scalar case
+    if isa(pc0, Number) && isa(pc1, Number)
+        pc0_log = log10(pc0)
+        pc1_log = log10(pc1)
+        center_log = log10.(center)
+        return hue_from_logpower_color(pc0_log, pc1_log; center=center_log)
+    # Handle vector case - apply function element-wise
+    elseif isa(pc0, Vector) && isa(pc1, Vector)
+        return [hue_from_power_color(p0, p1; center=center) for (p0, p1) in zip(pc0, pc1)]
+    else
+        error("pc0 and pc1 must both be either scalars or vectors")
+    end
+end
+
+function hue_from_logpower_color(log10pc0, log10pc1; center=log10.([4.51920, 0.453724]))
+    """
+    Measure the angle of a point in the log-power color diagram with respect to the center.
+
+    Angles are measured in radians, **in the clockwise direction**, with respect to a line oriented
+    at -45 degrees with respect to the horizontal axis.
+
+    See Heil et al. 2015, MNRAS, 448, 3348.
+
+    Parameters
+    ----------
+    log10pc0 : Number
+        The log10 power color in the first frequency range.
+    log10pc1 : Number
+        The log10 power color in the second frequency range.
+
+    Other Parameters
+    ----------------
+    center : Vector{Number}, optional, default log10([4.51920, 0.453724])
+        The coordinates of the center of the power color diagram.
+
+    Returns
+    -------
+    hue : Number
+        The angle of the point with respect to the center, in radians.
+    """
+    return (3 / 4) * pi - atan(log10pc1 - center[2], log10pc0 - center[1])
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 287371b..cc1b263 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -3,3 +3,4 @@ using Test
 using FFTW, Distributions, Statistics, StatsBase, HDF5
 include("test_fourier.jl")
 include("test_gti.jl")
+include("test_power_colors.jl")
diff --git a/test/test_power_colors.jl b/test/test_power_colors.jl
new file mode 100644
index 0000000..0f09bc4
--- /dev/null
+++ b/test/test_power_colors.jl
@@ -0,0 +1,106 @@
+using Test
+using Random
+using Distributions
+
+rng = Random.Xoshiro(1259723)
+Random.seed!(42)
+
+include("../src/power_colors.jl")
+
+struct TestPowerColor
+    freq::Vector{Float64}
+    power::Vector{Float64}
+    pc0::Float64
+    pc0e::Float64
+    pc1::Float64
+    pc1e::Float64
+    lpc0::Float64
+    lpc1::Float64
+    rms::Float64
+    rmse::Float64
+    configuration::Dict{String, Any}
+end
+
+function setup_test()
+    freq = collect(0.0001:0.00001:17)
+    power = 1 ./ freq
+    pc0, pc0e, pc1, pc1e = power_color(freq, power)
+
+    # Ensure power_color returns the expected number of values
+    result = power_color(freq, power; return_log=true)
+    
+    if length(result) < 4
+        error("Unexpected return value from power_color: ", result)
+    end
+
+    lpc0, _, lpc1, _ = result
+    configuration = deepcopy(DEFAULT_COLOR_CONFIGURATION)
+    configuration["state_definitions"]["HSS"]["hue_limits"] = [300, 370]
+
+    return TestPowerColor(freq, power, pc0, pc0e, pc1, pc1e, lpc0, lpc1, 0.1, 0.01, configuration)
+end
+
+test_case = setup_test()
+
+@testset "Power Color Tests" begin
+    @test isapprox(test_case.pc0, 1.0; atol=1e-3)  # Relaxed tolerance
+    @test isapprox(test_case.pc1, 1.0; atol=1e-3)
+
+    @test isapprox(test_case.lpc0, 0.0; atol=0.001)
+    @test isapprox(test_case.lpc1, 0.0; atol=0.001)
+
+    # Edge case: Ensuring out-of-range frequencies throw an error
+    good = test_case.freq .> (1 / 255)
+    @test_throws ArgumentError power_color(test_case.freq[good], test_case.power[good])
+
+    good = test_case.freq .< 15
+    @test_throws ArgumentError power_color(test_case.freq[good], test_case.power[good])
+
+    # Ensuring incorrect freq_edges throw an error
+    @test_throws ArgumentError power_color(test_case.freq, test_case.power; freq_edges=[1])
+    @test_throws ArgumentError power_color(test_case.freq, test_case.power; freq_edges=[1,2,3,4,5,6])
+
+    # Ensure `freqs_to_exclude` is properly formatted before testing
+    valid_exclusions = [Tuple(f) for f in ([1, 1.1], [3.0, 4.0], [4.5, 5.5])]
+    
+    for fte in valid_exclusions
+        try
+            power_color(test_case.freq, test_case.power; freqs_to_exclude=[fte])
+        catch e
+            @test false  # Should not reach here
+        end
+    end
+
+    # Check that excluding frequencies doesn't break calculations
+    pc0, _, pc1, _ = power_color(test_case.freq, test_case.power; freqs_to_exclude=[(1, 1.1)])
+    @test isapprox(pc0, 1.0; atol=0.001)
+    @test isapprox(pc1, 1.0; atol=0.001)
+
+    # Check power_err variations to ensure error propagation works correctly
+    pc0, pc0_err, pc1, pc1_err = power_color(test_case.freq, test_case.power; power_err=test_case.power / 2)
+    pc0e, pc0e_err, pc1e, pc1e_err = power_color(test_case.freq, test_case.power; power_err=test_case.power)
+
+    @test isapprox(pc0, 1.0; atol=0.001)
+    @test isapprox(pc1, 1.0; atol=0.001)
+    @test isapprox(pc0e, 1.0; atol=0.001)
+    @test isapprox(pc1e, 1.0; atol=0.001)
+    @test isapprox(pc0e_err / pc0_err, 2.0; atol=0.001)
+    @test isapprox(pc1e_err / pc1_err, 2.0; atol=0.001)
+
+    # Hue testing
+    center = (4.51920, 0.453724)
+    log_center = log10.(center)
+    for angle in range(0, stop=380, step=20)
+        rad_angle = deg2rad(angle)
+        factor = rand(rng, Uniform(0.1, 10))
+        x = factor * cos(3 / 4 * π - rad_angle) + log_center[1]
+        y = factor * sin(3 / 4 * π - rad_angle) + log_center[2]
+        # Fix: Pass center as a keyword argument
+        hue = hue_from_power_color(10^x, 10^y; center=center)
+
+        c2 = (sin(hue) - sin(rad_angle))^2 + (cos(hue) - cos(rad_angle))^2
+        angle_diff = acos((2.0 - c2) / 2.0)
+
+        @test isapprox(angle_diff, 0.0; atol=0.001)
+    end
+end
\ No newline at end of file