Example - a multi-star spectrum viewer¶

A complete blueprint for a small tool that does this:

"I have 10 FITS spectra of different stars. I want to load them, normalise the continuum, compute a quality SNR, identify the Balmer lines and measure H-alpha for each - and see the result in one interactive page."

You will build it twice, with the same workflow, just two different doors into spectro-kernel:

Mode A - library. A Python script that imports spectro_kernel and does everything in-process. Simplest, fastest, fully offline once the kernel is installed.
Mode B - MCP server. The same workflow, but every step is a call to a deployed spectro-kernel MCP server. Useful when the analysis runs outside Python (a web backend in Node/Go, an AI agent) or you want one shared kernel for several apps.

Both modes produce the same numbers and the same HTML page. That parity is the whole point of the catalogue.

0 · The scenario¶

You have a folder of FITS spectra, one per star:

data/
├── vega.fits
├── deneb.fits
├── altair.fits
├── arcturus.fits
├── betelgeuse.fits
├── rigel.fits
├── sirius.fits
├── procyon.fits
├── capella.fits
└── pollux.fits

(Don't have them? Use any of your own; or generate ten synthetic ones - see Your first spectrum for the helper snippet.)

Goal - for each spectrum:

Read the FITS file.
Normalise the continuum (polynomial, order 3).
Compute the signal-to-noise ratio (DER_SNR).
Detect the Balmer lines, identify them.
Fit H-alpha with a Gaussian - record its centre, FWHM and equivalent width.

…then render one interactive Plotly page that overlays the ten normalised spectra and shows a summary table.

1 · Install once¶

uv venv --python 3.12
uv pip install "spectro-kernel[viz]"      # library + Plotly viz
uv pip install fastmcp                    # only needed for Mode B

(Or pip install if you prefer - see Getting started.)

2 · Mode A - library¶

A single self-contained script. viewer_lib.py:

"""Multi-star spectrum viewer - library mode."""
from pathlib import Path
import plotly.graph_objects as go

from spectro_kernel import WorkContext, run_algorithm
from spectro_kernel.io import read_fits

DATA = Path("data")
HALPHA_REST = 6562.79          # Å
HTML_OUT   = Path("viewer.html")


def analyse_one(fits_path: Path) -> dict:
    """Run the same pipeline on one FITS file; return a tidy summary."""
    ctx = WorkContext(spectrum=read_fits(fits_path))

    # 1. normalise so every spectrum sits around flux = 1
    run_algorithm("normalize_polynomial", ctx, {"order": 3})

    # 2. data-quality metric
    run_algorithm("snr_der", ctx)

    # 3. find Balmer lines (matches against the bundled catalogue)
    run_algorithm("detect_lines", ctx, {"catalog": "balmer"})

    # 4. fit H-alpha - may fail (e.g. line not in range); swallow gracefully
    fit_out = run_algorithm(
        "fit_gaussian_line", ctx,
        {"line_center_angstrom": HALPHA_REST,
         "window_angstrom": 40.0, "label": "H-alpha"},
        raise_on_error=False,
    )

    fit = ctx.line_fits.get("H-alpha")
    return {
        "star": fits_path.stem,
        "spectrum": ctx.spectrum,                                  # normalised
        "snr": ctx.metrics["snr_der"],
        "n_lines_matched": int(ctx.metrics.get("n_lines_matched", 0)),
        "halpha_fitted": fit_out.success and fit is not None,
        "halpha_fwhm_a": fit.fwhm_angstrom if fit else None,
        "halpha_ew_a":   fit.equivalent_width_angstrom if fit else None,
        "history": ctx.history,                                    # full audit trail
    }


def main() -> None:
    fits_files = sorted(DATA.glob("*.fits"))
    if not fits_files:
        raise SystemExit(f"No FITS files in {DATA}/")

    results = [analyse_one(p) for p in fits_files]

    # ---- Plotly overlay --------------------------------------------------
    fig = go.Figure()
    for r in results:
        s = r["spectrum"]
        fig.add_scatter(x=s.wavelength, y=s.flux, mode="lines", name=r["star"])
    fig.update_layout(
        title="Ten stars - continuum-normalised spectra",
        xaxis_title="Wavelength (Å)",
        yaxis_title="Normalised flux",
        template="plotly_dark",
        legend_title_text="star",
    )
    fig.write_html(HTML_OUT, include_plotlyjs="cdn")

    # ---- Summary table ---------------------------------------------------
    print(f"\n{'Star':<14} {'SNR':>6}  {'lines':>5}  "
          f"{'Hα FWHM (Å)':>12}  {'Hα EW (Å)':>10}")
    print("-" * 56)
    for r in results:
        fwhm = f"{r['halpha_fwhm_a']:>12.2f}" if r["halpha_fwhm_a"] else f"{'—':>12}"
        ew   = f"{r['halpha_ew_a']:>10.3f}"  if r["halpha_ew_a"]   else f"{'—':>10}"
        print(f"{r['star']:<14} {r['snr']:>6.0f}  "
              f"{r['n_lines_matched']:>5d}  {fwhm}  {ew}")
    print(f"\nOpen {HTML_OUT} to inspect the overlay.")


if __name__ == "__main__":
    main()

Run it:

python viewer_lib.py
open viewer.html

You get a viewer.html with the 10 spectra overlaid, plus the summary table on stdout. Each call's parameters, version and hashes are recorded in results[i]["history"] if you want to log them.

3 · Mode B - through the MCP server¶

Same goal, same numbers. The only thing that changes is where the kernel runs: now it sits on a hosted server (any container host - DigitalOcean App Platform, Fly.io, Render, …) and we talk to it over the Model Context Protocol.

viewer_mcp.py:

"""Multi-star spectrum viewer - MCP client mode.

Same analysis as viewer_lib.py, driven through a hosted spectro-kernel MCP server.
"""
from __future__ import annotations
import asyncio
import json
from pathlib import Path

import plotly.graph_objects as go
from fastmcp import Client

MCP_URL    = "https://your-app.ondigitalocean.app/mcp"
DATA       = Path("data")
HALPHA_REST = 6562.79
HTML_OUT   = Path("viewer-mcp.html")


def _unwrap(result):
    """Pull the JSON payload out of a FastMCP CallToolResult, tolerantly."""
    if getattr(result, "structured_content", None):
        return result.structured_content
    if getattr(result, "data", None) is not None:
        return result.data
    return json.loads(result.content[0].text)  # last-resort text fallback


async def analyse_one(client: Client, fits_path: Path) -> dict:
    star = fits_path.stem

    # Each star gets its own server-side session; results don't leak between stars.
    sid = _unwrap(await client.call_tool("create_session", {}))["session_id"]

    try:
        await client.call_tool(
            "load_spectrum",
            {"session_id": sid, "path": str(fits_path.resolve())},
        )
        await client.call_tool(
            "normalize_polynomial",
            {"session_id": sid, "params": {"order": 3}},
        )
        await client.call_tool("snr_der", {"session_id": sid})
        await client.call_tool(
            "detect_lines",
            {"session_id": sid, "params": {"catalog": "balmer"}},
        )
        fit_result = _unwrap(await client.call_tool(
            "fit_gaussian_line",
            {"session_id": sid,
             "params": {"line_center_angstrom": HALPHA_REST,
                        "window_angstrom": 40.0, "label": "H-alpha"}},
        ))

        # Pull the consolidated state and a downsampled spectrum for plotting.
        state = _unwrap(await client.call_tool(
            "get_session_state", {"session_id": sid}))
        preview = _unwrap(await client.call_tool(
            "get_spectrum_preview", {"session_id": sid}))

        return {
            "star": star,
            "snr": state["metrics"].get("snr_der"),
            "n_lines_matched": int(state["metrics"].get("n_lines_matched", 0)),
            "halpha_fwhm_a": fit_result.get("metrics", {}).get("fwhm_angstrom"),
            "halpha_ew_a":   fit_result.get("artifacts", {})
                                       .get("line_fit", {})
                                       .get("equivalent_width_angstrom"),
            "preview": preview,
        }
    finally:
        await client.call_tool("end_session", {"session_id": sid})


async def main() -> None:
    fits_files = sorted(DATA.glob("*.fits"))
    if not fits_files:
        raise SystemExit(f"No FITS files in {DATA}/")

    async with Client(MCP_URL) as client:
        results = await asyncio.gather(
            *(analyse_one(client, p) for p in fits_files)
        )

    # Plotly overlay - using the previews returned by the server.
    fig = go.Figure()
    for r in results:
        p = r["preview"]
        fig.add_scatter(x=p["wavelength_sample"], y=p["flux_sample"],
                        mode="lines", name=r["star"])
    fig.update_layout(
        title="Ten stars - via the MCP server",
        xaxis_title="Wavelength (Å)",
        yaxis_title="Normalised flux",
        template="plotly_dark",
    )
    fig.write_html(HTML_OUT, include_plotlyjs="cdn")

    print(f"\n{'Star':<14} {'SNR':>6}  {'lines':>5}  {'Hα FWHM (Å)':>12}")
    print("-" * 44)
    for r in results:
        fwhm = f"{r['halpha_fwhm_a']:>12.2f}" if r["halpha_fwhm_a"] else f"{'—':>12}"
        print(f"{r['star']:<14} {r['snr']:>6.0f}  "
              f"{r['n_lines_matched']:>5d}  {fwhm}")


if __name__ == "__main__":
    asyncio.run(main())

Run it:

python viewer_mcp.py
open viewer-mcp.html

Where do the FITS files live?

The example passes a local path to load_spectrum. The MCP server loads that path, not the client - so the file must be reachable from the server. Either deploy with the data baked in, or pass an https:// URL (FITS over HTTP is supported by the kernel's reader), or add an upload tool that takes the bytes.

4 · Same numbers, every door¶

viewer_lib.py and viewer_mcp.py use the same algorithms with the same parameters; they should print identical SNRs and identical H-alpha FWHM values (to floating-point precision). That parity is what makes the kernel useful: a measurement quoted from your dashboard is exactly the same as one quoted from a notebook, or from a Claude session.

Want to assert it programmatically? Add a tiny pytest that runs both scripts and diffs the metrics - that is how spectro-kernel itself enforces parity in CI.

5 · Bonus - let an agent do it¶

With the MCP server deployed and configured in Claude Desktop, the same workflow can be driven in natural language:

"You have access to the spectro-kernel MCP server. For each FITS file in data/, create a session, normalise with normalize_polynomial, then snr_der, then fit_gaussian_line on H-alpha at 6562.79 Å with a 40 Å window. Return a markdown table of the SNR and the H-alpha FWHM per star."

Claude discovers the tools (list_algorithms, describe_algorithm), composes the calls - identical to what the script does, including the same session-per-star isolation - and writes the table. The kernel never knows whether the orchestrator is a Python script or an AI agent.

6 · Next steps for a real project¶

When this becomes a proper application (a "starsight" mini-app, say):

Wrap it in a Streamlit or Dash UI - three lines around the analysis.
Replace the file glob with a database query (you already have ctx.metrics and ctx.line_fits to persist).
Add a few more lines: fit_lorentzian_line for the wings, measure_radial_velocity if your spectra are calibrated, lomb_scargle if you also have a time series.
Deploy the UI to Netlify / Vercel, point it at the kernel deployed on DO - and you've got a viewer running entirely in the cloud.

Use this document as the spec when you spin up that repo.