diff --git a/CITATION.cff b/CITATION.cff
index 0df67a6..92b1531 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -10,8 +10,8 @@ authors:
 repository-code: https://github.com/JuliaAstro/SolarPosition.jl
 url: https://juliaastro.github.io/SolarPosition.jl/
 license: MIT
-version: 0.4.0
-date-released: 2026-01-05
+version: 0.4.1
+date-released: 2026-01-17
 abstract: >-
   SolarPosition.jl provides a simple, unified interface to a collection
   of validated solar position algorithms written in pure Julia.
diff --git a/Project.toml b/Project.toml
index 7450b97..3944d2b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "SolarPosition"
 uuid = "5b9d1343-a731-5a90-8730-7bf8d89bf3eb"
-version = "0.4.0"
+version = "0.4.1"
 authors = ["Stefan de Lange"]
 
 [workspace]
diff --git a/docs/make.jl b/docs/make.jl
index 0850fad..7fc11b1 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -51,6 +51,7 @@ makedocs(;
         "literature.md",
         "contributing.md",
     ],
+    draft = true,
 )
 
 deploydocs(;
diff --git a/docs/src/guides/plotting.md b/docs/src/guides/plotting.md
index 18570a6..a39e1b9 100644
--- a/docs/src/guides/plotting.md
+++ b/docs/src/guides/plotting.md
@@ -1,7 +1,13 @@
+```@meta
+Draft = false
+```
+
 # [Plotting with Makie.jl](@id plotting-examples)
 
 SolarPosition.jl provides a plotting extension for [Makie.jl](https://makie.juliaplots.org/stable/).
 
+The main plotting function is [`analemmas!`](@ref).
+
 To use it, simply import both the `SolarPosition` and `Makie` packages:
 
 ```@example plotting
@@ -11,62 +17,53 @@ using CairoMakie
 # supporting packages
 using Dates
 using TimeZones
-using DataFrames
 ```
 
 This example notebook is based on the [pvlib sun path example](https://pvlib-python.readthedocs.io/en/stable/gallery/solar-position/plot_sunpath_diagrams.html).
 
 ## Basic Sun Path Plotting
 
-Let's start by defining an observer location and calculating solar positions for a whole year:
+The plotting functions generate analemmas (figure-8 patterns showing the sun's position at
+each hour of the day throughout the year). You simply provide an observer location and
+the year you want to visualize:
 
 ```@example plotting
 # Define observer location (New Delhi, India)
 # Parameters: latitude, longitude, altitude in meters
-tz = tz"Asia/Kolkata"
 obs = Observer(28.6, 77.2, 0.0)
-
-# Generate hourly timestamps for a whole year
-times = collect(ZonedDateTime(DateTime(2019), tz):Hour(1):ZonedDateTime(DateTime(2020), tz))
-
-# This returns a StructVector with solar position data
-positions = solar_position(obs, times)
-
-# For plotting, we need to create a DataFrame that includes the timestamps
-df = DataFrame(positions)
-df.datetime = times
-
-# We can inspect the first few entries
-first(df, 5)
+tz = TimeZone("Asia/Kolkata")
+year = 2019
 ```
 
 ## Simple Sun Path Plot in Cartesian Coordinates
 
-We can visualize solar positions in cartesian coordinates using the `sunpathplot`
-function:
+We can visualize solar positions in cartesian coordinates using the `analemmas!`
+function. The function automatically generates analemmas for all 24 hours of the day:
 
 ```@example plotting
 fig = Figure(backgroundcolor = (:white, 0.0), textcolor= "#f5ab35")
 ax = Axis(fig[1, 1], backgroundcolor = (:white, 0.0))
-sunpathplot!(ax, df, hour_labels = false)
+analemmas!(ax, obs, year, hour_labels = false)
 fig
 ```
 
 ## Polar Coordinates with Hour Labels
 
-We can also work directly with a `DataFrame`. Note that for plotting we need to include
-the datetime information, so we add it to the DataFrame.
-
-Plotting in polar coordinates with `sunpathpolarplot` may yield a more intuitive
-representation of the solar path. Here, we also enable hourly labels for better
-readability:
+Plotting in polar coordinates with `analemmas!` may yield a more intuitive
+representation of the solar path. Here, we enable hourly labels for better readability:
 
 ```@example plotting
 fig2 = Figure(backgroundcolor = :transparent, textcolor= "#f5ab35", size = (800, 600))
 ax2 = PolarAxis(fig2[1, 1], backgroundcolor = "#1f2424")
-sunpathpolarplot!(ax2, df, hour_labels = true)
+analemmas!(ax2, obs, year, hour_labels = true)
+fig2
+```
 
-# Draw individual days
+Now let's manually plot the full solar path for specific dates March 21, June 21, and
+December 21. Also known as the vernal equinox, summer solstice, and winter solstice,
+respectively:
+
+```@example plotting
 line_objects = []
 for (date, label) in [(Date("2019-03-21"), "Mar 21"),
                       (Date("2019-06-21"), "Jun 21"),
@@ -74,9 +71,7 @@ for (date, label) in [(Date("2019-03-21"), "Mar 21"),
     times = collect(ZonedDateTime(DateTime(date), tz):Minute(5):ZonedDateTime(DateTime(date) + Day(1), tz))
     solpos = solar_position(obs, times)
     above_horizon = solpos.elevation .> 0
-    day_df = DataFrame(solpos)
-    day_df.datetime = times
-    day_filtered = day_df[above_horizon, :]
+    day_filtered = solpos[above_horizon]
     line_obj = lines!(ax2, deg2rad.(day_filtered.azimuth), day_filtered.zenith,
                       linewidth = 2, label = label)
     push!(line_objects, line_obj)
@@ -101,11 +96,23 @@ It tells us when the sun rises, reaches its highest point, and sets. And hence a
 length of the day. From the figure we can also read that in June the days are longest,
 while in December they are shortest.
 
-## Plotting without a custom axis
+## Custom Color Schemes
+
+You can customize the color scheme used for the analemmas by passing a `colorscheme` argument.
+Here's an example using the `:balance` colorscheme.
 
-Finally, we can also create plots without explicitly defining an axis beforehand.
-This is a more concise way to create plots, but it offers less customization:
+!!! info
+    More colorschemes are available in the [Makie documentation](https://docs.makie.org/dev/explanations/colors).
 
 ```@example plotting
-sunpathpolarplot(df, hour_labels = true, colorbar = true)
+fig = Figure(backgroundcolor = (:white, 0.0), textcolor= "#f5ab35")
+ax = Axis(fig[1, 1], backgroundcolor = (:white, 0.0))
+analemmas!(ax, obs, year, hour_labels = false, colorscheme = :balance)
+fig
+```
+
+## Docstrings
+
+```@docs
+analemmas!
 ```
diff --git a/examples/plotting.jl b/examples/plotting.jl
index 53934ea..cd578aa 100644
--- a/examples/plotting.jl
+++ b/examples/plotting.jl
@@ -1,35 +1,27 @@
 """Plot solar positions using SolarPosition.jl with hourly labels."""
 
 using Dates
-using TimeZones
-using DataFrames
 using CairoMakie
 using SolarPosition
 
 # define observer location (latitude, longitude, altitude in meters)
-tz = tz"Asia/Kolkata"
-obs = Observer(28.6, 77.2, 0.0)
+obs = Observer(28.6, 77.2, 0.0)  # New Delhi, India
+year = 2019
 
-# a whole year of hourly timestamps
-times = collect(ZonedDateTime(DateTime(2019), tz):Hour(1):ZonedDateTime(DateTime(2020), tz))
-positions = solar_position(obs, times)
-
-# plot positions from NamedTuple with hourly labels in polar coordinates
-df = DataFrame(positions)
-df.datetime = times
-sunpathplot(df)
-
-# plot DataFrame with hourly labels in cartesian coordinates
-sunpathpolarplot(df, hour_labels = true, colorbar = true)
-
-# plot DataFrame in polar coordinates with hourly labels
+# plot in cartesian coordinates with hourly labels
 fig = Figure()
-ax = PolarAxis(fig[1, 1], title = "Polar Coordinates with Hour Labels")
-sunpathpolarplot!(ax, df, hour_labels = true, colorbar = false)
+ax = Axis(fig[1, 1], title = "Cartesian Coordinates with Hour Labels")
+analemmas!(ax, obs, year, hour_labels = true)
 fig
 
-# example without hourly labels for comparison
+# plot in polar coordinates with hourly labels
 fig2 = Figure()
-ax2 = Axis(fig2[1, 1], title = "Cartesian Coordinates (No Labels)")
-sunpathplot!(ax2, df; hour_labels = false, colorbar = true)
+ax2 = PolarAxis(fig2[1, 1], title = "Polar Coordinates with Hour Labels")
+analemmas!(ax2, obs, year, hour_labels = true)
 fig2
+
+# example without hourly labels for comparison
+fig3 = Figure()
+ax3 = Axis(fig3[1, 1], title = "Cartesian Coordinates (No Labels)")
+analemmas!(ax3, obs, year, hour_labels = false)
+fig3
diff --git a/ext/SolarPositionMakieExt.jl b/ext/SolarPositionMakieExt.jl
index b746268..b06e4d6 100644
--- a/ext/SolarPositionMakieExt.jl
+++ b/ext/SolarPositionMakieExt.jl
@@ -1,205 +1,206 @@
 module SolarPositionMakieExt
 
-using Dates, Tables, Makie
+using Dates, Makie
 using SolarPosition
-import SolarPosition: sunpathpolarplot, sunpathpolarplot!, sunpathplot, sunpathplot!
-
-function _add_colorbar!(fig, plot_obj, position = fig[1, 2])
-    month_names =
-        ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
-    month_days = [15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349]
-    inverted_month_days = 365 .- reverse(month_days)
-    inverted_month_names = reverse(month_names)
-    return Colorbar(position, plot_obj, ticks = (inverted_month_days, inverted_month_names))
-end
-
-function _try_add_colorbar_to_recipe!(ax, plot_obj)
-    try
-        parent = ax.parent
-        while !isa(parent, Figure) &&
-                  !isnothing(parent) &&
-                  hasfield(typeof(parent), :parent)
-            parent = parent.parent
-        end
-        if isa(parent, Figure)
-            _add_colorbar!(parent, plot_obj, parent[1, 2])
-        end
-    catch
-        # silently fail if colorbar can't be added
-    end
-end
-
-function _add_hour_labels!(ax, dts, ze, el, az, coords)
-    hours = hour.(dts)
-    visible_mask = el .> 0
-    if !any(visible_mask)
-        return
+import SolarPosition: analemmas!
+
+"""
+    _generate_analemma_data(observer, year, hours; algorithm=PSA())
+
+Generate solar position data for analemmas at specified hours throughout the year.
+
+# Arguments
+- `observer::Observer`: Observer location
+- `year::Int`: Year for which to generate positions
+- `hours::AbstractVector{Int}`: Hours of the day to generate analemmas for (0-23)
+- `algorithm`: Solar position algorithm to use (default: PSA())
+
+# Returns
+A named tuple with fields: `datetime`, `zenith`, `azimuth`, `hour`, `dayofyear`
+"""
+function _generate_analemma_data(observer, year, hours; algorithm = PSA())
+    # Generate dates for every day of the year
+    datetimes = DateTime(year, 1, 1):Hour(1):DateTime(year, 12, 31)
+    n_points = length(datetimes)
+
+    # Pre-allocate arrays
+    zeniths = Vector{Float64}(undef, n_points)
+    azimuths = Vector{Float64}(undef, n_points)
+    hour_vals = Vector{Int}(undef, n_points)
+    doy_vals = Vector{Int}(undef, n_points)
+
+    # Calculate solar position
+    for (idx, dt) in enumerate(datetimes)
+        pos = solar_position(observer, dt, algorithm)
+        zeniths[idx] = pos.zenith
+        azimuths[idx] = pos.azimuth
+        hour_vals[idx] = Dates.hour(dt)
+        doy_vals[idx] = dayofyear(dt)
     end
 
-    ze_visible = ze[visible_mask]
-    el_visible = el[visible_mask]
-    az_visible = az[visible_mask]
-    hours_visible = hours[visible_mask]
+    return (
+        datetime = datetimes,
+        zenith = zeniths,
+        azimuth = azimuths,
+        hour = hour_vals,
+        dayofyear = doy_vals,
+    )
+end
 
-    for hour_val in unique(hours_visible)
-        hour_mask = hours_visible .== hour_val
-        if !any(hour_mask)
-            continue
+"""
+    _add_hour_labels!(ax, observer, year, hours, coords; algorithm=PSA())
+
+Add hour labels to a sun path plot at the position of maximum elevation for each hour.
+
+# Arguments
+- `ax`: The axis to add labels to
+- `observer::Observer`: Observer location
+- `year::Int`: Year for the solar positions
+- `hours::AbstractVector{Int}`: Hours to label (0-23)
+- `coords`: Coordinate system (`:polar` or `:cartesian`)
+- `algorithm`: Solar position algorithm to use (default: PSA())
+"""
+function _add_hour_labels!(ax, observer, year, hours, coords; algorithm = PSA())
+    # For each hour, find the day when the sun reaches maximum elevation at that hour
+    for hour in hours
+        # Sample a few days throughout the year to find max elevation
+        dates = [Date(year, m, 15) for m = 1:12]
+        max_el = -90.0
+        best_pos = nothing
+
+        for date in dates
+            dt = DateTime(date) + Hour(hour)
+            pos = solar_position(observer, dt, algorithm)
+            el = 90.0 - pos.zenith
+
+            if el > max_el
+                max_el = el
+                best_pos = pos
+            end
         end
 
-        if coords === :polar
-            idx = argmin(ze_visible[hour_mask])
-            label_az = az_visible[hour_mask][idx]
-            label_ze = ze_visible[hour_mask][idx]
-            x, y = deg2rad(label_az), label_ze
-        else  # cartesian coordinates
-            idx = argmax(el_visible[hour_mask])
-            label_az = az_visible[hour_mask][idx]
-            label_el = el_visible[hour_mask][idx]
-            offset = label_az < 180 ? -10 : 10
-            x, y = label_az + offset, label_el
+        # Only label if sun is above horizon
+        if max_el > 0 && !isnothing(best_pos)
+            if coords === :polar
+                x = deg2rad(best_pos.azimuth)
+                y = best_pos.zenith
+            else  # cartesian
+                offset = best_pos.azimuth < 180 ? -10 : 10
+                x = best_pos.azimuth + offset
+                y = max_el
+            end
+
+            text!(
+                ax,
+                x,
+                y,
+                text = lpad(string(hour), 2, '0'),
+                align = (:center, :bottom),
+                fontsize = 13,
+            )
         end
-
-        text!(
-            ax,
-            x,
-            y,
-            text = lpad(string(hour_val), 2, '0'),
-            align = (:center, :bottom),
-            fontsize = 13,
-        )
     end
 end
 
-
-@recipe(SunpathPlot) do scene
-    Theme(
-        colormap = :twilight,
-        backgroundcolor = :white,
-        markersize = 3,
-        hour_labels = true,
-        colorbar = true,
-    )
+"""
+    _configure_axis!(ax::Axis)
+
+Configure cartesian axis for analemma plot.
+"""
+function _configure_axis!(ax::Axis)
+    xlims!(ax, 0, 360)
+    ylims!(ax, 0, 90)
+    ax.xlabel = "Azimuth (°)"
+    ax.ylabel = "Elevation (°)"
+    ax.xticks = 0:30:360
+    ax.yticks = 0:10:90
 end
 
-@recipe(SunpathpolarPlot) do scene
-    Theme(colormap = :twilight, markersize = 3, hour_labels = true, colorbar = true)
+"""
+    _configure_axis!(ax::PolarAxis)
+
+Configure polar axis for analemma plot.
+"""
+function _configure_axis!(ax::PolarAxis)
+    ax.direction = -1
+    ax.theta_0 = -π / 2
+    ax.rlimits = (0, 90)
 end
 
-function Makie.plot!(sp::SunpathPlot{<:Tuple})
-    data = sp[1][]
-    dts = data.datetime
-    ze = data.zenith
-    az = data.azimuth
-    el = 90 .- ze
-    vals = 365 .- dayofyear.(dts)
+"""
+    _plot_analemmas!(sp, ax::Axis, data, observer, year, hours)
 
-    ax = current_axis()
-    if ax isa Axis
-        xlims!(ax, 0, 360)
-        ylims!(ax, 0, 90)
-        ax.xlabel = "Azimuth (°)"
-        ax.ylabel = "Elevation (°)"
-        ax.xticks = 0:30:360
-        ax.yticks = 0:10:90
-    end
+Plot analemmas on a cartesian axis.
+"""
+function _plot_analemmas!(sp, ax::Axis, data, observer, year, hours)
+    _configure_axis!(ax)
+
+    # Only plot points where sun is above horizon
+    el = 90 .- data.zenith
+    above_horizon = el .> 0
+    vals = 365 .- data.dayofyear[above_horizon]
 
     p = scatter!(
         sp,
-        az,
-        el;
+        data.azimuth[above_horizon],
+        el[above_horizon];
         color = vals,
-        colormap = sp.colormap[],
+        colormap = sp.colorscheme[],
         markersize = sp.markersize[],
     )
 
-    # add hourly labels if requested
+    # Add hour labels if requested
     if sp.hour_labels[]
-        _add_hour_labels!(ax, dts, ze, el, az, :cartesian)
-    end
-
-    # add colorbar if requested and possible
-    if sp.colorbar[]
-        _try_add_colorbar_to_recipe!(ax, p)
+        _add_hour_labels!(ax, observer, year, hours, :cartesian)
     end
 
     return p
 end
 
+"""
+    _plot_analemmas!(sp, ax::PolarAxis, data, observer, year, hours)
 
-function Makie.plot!(sp::SunpathpolarPlot{<:Tuple})
-    data = sp[1][]
-    dts = data.datetime
-    ze = data.zenith
-    az = data.azimuth
-    vals = 365 .- dayofyear.(dts)
+Plot analemmas on a polar axis.
+"""
+function _plot_analemmas!(sp, ax::PolarAxis, data, observer, year, hours)
+    _configure_axis!(ax)
 
-    ax = current_axis()
-    if ax isa PolarAxis
-        ax.direction = -1
-        ax.theta_0 = -π / 2
-        ax.rlimits = (0, 90)
-    end
+    # Only plot points where sun is above horizon
+    el = 90 .- data.zenith
+    above_horizon = el .> 0
+    vals = 365 .- data.dayofyear[above_horizon]
 
     p = scatter!(
         sp,
-        deg2rad.(az),
-        ze;
+        deg2rad.(data.azimuth[above_horizon]),
+        data.zenith[above_horizon];
         color = vals,
-        colormap = sp.colormap[],
+        colormap = sp.colorscheme[],
         markersize = sp.markersize[],
     )
 
-    # add hourly labels if requested
+    # Add hour labels if requested
     if sp.hour_labels[]
-        _add_hour_labels!(ax, dts, ze, 90 .- ze, az, :polar)
-    end
-
-    # add colorbar if requested and possible
-    if sp.colorbar[]
-        _try_add_colorbar_to_recipe!(ax, p)
+        _add_hour_labels!(ax, observer, year, hours, :polar)
     end
 
     return p
 end
 
-
-
-function Makie.convert_arguments(sp::Type{Union{<:SunpathPlot,<:SunpathpolarPlot}}, tbl)
-    cols = Tables.columns(tbl)
-    dts = Tables.getcolumn(cols, :datetime)
-    zenith = Tables.getcolumn(cols, :zenith)
-    azimuth = Tables.getcolumn(cols, :azimuth)
-    return ((datetime = dts, zenith = zenith, azimuth = azimuth),)
+@recipe(Analemmas) do scene
+    Theme(colorscheme = :twilight, markersize = 3, hour_labels = true)
 end
 
-function SolarPosition.sunpathpolarplot(
-    data;
-    hour_labels = true,
-    colorbar = true,
-    kwargs...,
-)
-    fig = Figure()
-    ax = PolarAxis(fig[1, 1])
-    p = sunpathpolarplot!(ax, data; hour_labels = hour_labels, kwargs...)
-
-    if colorbar
-        _add_colorbar!(fig, p)
-    end
-
-    return fig
-end
+function Makie.plot!(sp::Analemmas{<:Tuple{Observer,Int}})
+    observer = sp[1][]
+    year = sp[2][]
 
-function SolarPosition.sunpathplot(data; hour_labels = true, colorbar = true, kwargs...)
-    fig = Figure()
-    ax = Axis(fig[1, 1])
-    p = sunpathplot!(ax, data; hour_labels = hour_labels, kwargs...)
+    # Generate analemma data for all 24 hours
+    hours = collect(0:23)
+    data = _generate_analemma_data(observer, year, hours)
 
-    if colorbar
-        _add_colorbar!(fig, p)
-    end
-
-    return fig
+    ax = current_axis()
+    return _plot_analemmas!(sp, ax, data, observer, year, hours)
 end
 
 end # module
diff --git a/src/SolarPosition.jl b/src/SolarPosition.jl
index 37b1953..6716e2e 100644
--- a/src/SolarPosition.jl
+++ b/src/SolarPosition.jl
@@ -12,15 +12,65 @@ include("Utilities/Utilities.jl")
 @reexport using .Utilities
 
 # to make the makie extension work
-export sunpathplot
-export sunpathplot!
-export sunpathpolarplot
-export sunpathpolarplot!
-
-function sunpathplot end
-function sunpathplot! end
-function sunpathpolarplot end
-function sunpathpolarplot! end
+export analemmas!
+
+"""
+    $(TYPEDSIGNATURES)
+
+Plot analemmas (figure-8 patterns showing the sun's position at each hour throughout the
+year) for a given observer location and year.
+
+# Arguments
+- `ax`: A Makie [`Axis`](https://docs.makie.org/dev/reference/blocks/axis) or [`PolarAxis`](https://docs.makie.org/dev/reference/blocks/polaraxis) to plot on
+- `observer::Observer`: [`Observer`](@ref SolarPosition.Positioning.Observer) location (latitude, longitude, altitude)
+- `year::Int`: Year for which to generate the analemmas
+- `hour_labels::Bool=true`: Whether to add hour labels to the plot
+- `colorscheme::Symbol=:twilight`: Color scheme for the analemma points (any [Makie colormap](https://docs.makie.org/dev/explanations/colors))
+
+# Description
+This function automatically generates solar position data for all 24 hours of each day
+throughout the specified year and plots them as analemmas. The plot specializes to
+either [`PolarAxis`](https://docs.makie.org/dev/reference/blocks/polaraxis) or regular
+[`Axis`](https://docs.makie.org/dev/reference/blocks/axis) and adjusts the plot
+accordingly:
+
+- [`PolarAxis`](https://docs.makie.org/dev/reference/blocks/polaraxis): Plots in polar coordinates with azimuth as the angle and zenith as the radius
+- [`Axis`](https://docs.makie.org/dev/reference/blocks/axis): Plots in cartesian coordinates with azimuth on the x-axis and elevation on the y-axis
+
+The analemmas are colored by day of year. The default colorscheme is `:twilight`, but this
+can be customized using the `colorscheme` keyword argument.
+
+# Examples
+```julia
+using SolarPosition
+using CairoMakie
+
+# Define observer location (New Delhi, India)
+obs = Observer(28.6, 77.2, 0.0)
+year = 2019
+
+# Plot in cartesian coordinates
+fig = Figure()
+ax = Axis(fig[1, 1], title="Sun Path - Cartesian")
+analemmas!(ax, obs, year)
+fig
+
+# Plot in polar coordinates
+fig2 = Figure()
+ax2 = PolarAxis(fig2[1, 1], title="Sun Path - Polar")
+analemmas!(ax2, obs, year, hour_labels=true)
+fig2
+
+# Use a custom colorscheme
+fig3 = Figure()
+ax3 = Axis(fig3[1, 1], title="Sun Path - Custom Colors")
+analemmas!(ax3, obs, year, colorscheme=:viridis)
+fig3
+```
+
+See also: [`Observer`](@ref), [`solar_position`](@ref)
+"""
+function analemmas! end
 
 # to make the ModelingToolkit extension work
 export SolarPositionBlock