Added algo calc, spar_error and plot scripts

src/algo1.py

@@ -3,39 +3,91 @@
 # Copyright (C) 2026 Association Exergie <association.exergie@gmail.com>
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-import time
+from __future__ import annotations
-import pandas as pd
+
-from pathlib import Path
 import argparse
-from tqdm import tqdm
+import sys
+from bisect import bisect_left
+from dataclasses import dataclass
+from pathlib import Path
 
-parser = argparse.ArgumentParser()
+import pandas as pd
-parser.add_argument("file", type=Path, help="Source data")
 
-args = parser.parse_args()
 
-run_data_path: Path = args.file
+DEFAULT_DATA_ROOTS = (
+    Path("~/projects/Exergie/Engine-Recorder/recordings-rpm").expanduser(),
+    Path("~/projects/Exergie/Engine-Recorder/Python/recordings_rpm").expanduser(),
+)
+GP_COLUMNS = ("gp0_falling", "gp1_falling")
+SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3")
+INTERPOLATED_COLUMNS = (
+    "algo1_continuous_angle_deg",
+    "algo1_revolution_index",
+    "algo1_crank_angle_deg",
+)
 
-df = pd.read_csv(run_data_path).set_index("time_us", drop=False)
+
+@dataclass(frozen=True)
+class EncoderPoint:
+    time_us: int
+    continuous_angle_deg: float
+
+
+@dataclass(frozen=True)
+class SparkPulse:
+    time_us: int
+    column: str
+
+
+def resolve_input_path(path: Path) -> Path:
+    expanded = path.expanduser()
+    if expanded.exists():
+        return expanded
+
+    for root in DEFAULT_DATA_ROOTS:
+        candidate = root / path
+        if candidate.exists():
+            return candidate
+
+    searched = ", ".join(str(root) for root in DEFAULT_DATA_ROOTS)
+    raise FileNotFoundError(f"{path} does not exist; also searched {searched}")
+
+
+def read_rpm(path: Path) -> pd.DataFrame:
+    df = pd.read_csv(path)
+    missing = {"time_us", *GP_COLUMNS} - set(df.columns)
+    if missing:
+        raise ValueError(f"missing columns: {', '.join(sorted(missing))}")
+
+    df["time_us"] = df["time_us"].astype("int64")
+    return df
+
+
+def gp_times(df: pd.DataFrame, column: str) -> list[int]:
+    pulse_rows = df.loc[df[column].fillna(0).astype(int) == 1]
+    return sorted(int(row.time_us) for row in pulse_rows.itertuples())
+
+
+def generate_sparks(
+    crank_times: list[int],
+    cam_times: list[int],
+    *,
+    spark_advance_deg: float,
+) -> list[SparkPulse]:
+    events: list[tuple[int, int, str]] = []
+    events.extend((time_us, 0, "cam") for time_us in cam_times)
+    events.extend((time_us, 1, "crank") for time_us in crank_times)
+    events.sort()
 
     sync_ok = False
     cycle_ctr = 0
     cam_ctr = 0
     exhaust_timestamp = 0
     intake_timestamp = 0
-spark_release_time_1 = 0
+    sparks: list[SparkPulse] = []
-spark_release_time_2 = 0
-spark_release_time_3 = 0
-spark_adv = 25
 
-rows = []
+    for time_us, _, event_type in events:
-
+        if event_type == "cam":
-for _, row in tqdm(df.iterrows(), total=len(df)):
-    time_us: int = row["time_us"]
-    crank: int = row["crank"]
-    cam: int = row["cam"]
-
-    if cam == 1:
             if cycle_ctr < 4:
                 sync_ok = False
             if not sync_ok and cycle_ctr == 4:
@@ -43,54 +95,224 @@ for _, row in tqdm(df.iterrows(), total=len(df)):
             if cam_ctr >= 2:
                 sync_ok = True
             cycle_ctr = 0
+            continue
 
-    if crank == 1:
         cycle_ctr += 1
         if cycle_ctr >= 5:
             sync_ok = False
 
         if cycle_ctr == 2:
             exhaust_timestamp = time_us
-
+        elif cycle_ctr == 3:
-        if cycle_ctr == 3:
             intake_timestamp = time_us
-
+        elif cycle_ctr == 4 and sync_ok:
-        if cycle_ctr == 4 and sync_ok:
             d1a = time_us - intake_timestamp
             d1b = intake_timestamp - exhaust_timestamp
-            spark_release_time_1 = int(
+            predicted_interval = 2 * d1a - d1b
-                time_us + (2 * d1a - d1b) * (45 - spark_adv) / 180
+            if d1a <= 0 or predicted_interval <= 0:
-            )
+                continue
+
+            spark_1_time = int(time_us + d1a * (45.0 - spark_advance_deg) / 180.0)
+            spark_2_time = int(time_us + d1b * (45.0 - spark_advance_deg) / 180.0)
 
             t_s = (
-                d1a * (135 + spark_adv) / 180 + (2 * d1a - d1b) * (45 - spark_adv) / 180
+                d1a * (135.0 + spark_advance_deg) / 180.0
+                + predicted_interval * (45.0 - spark_advance_deg) / 180.0
+            )
+            t_a = d1a / 2.0 + t_s / 2.0
+            spark_3_time = int(t_a * (45.0 - spark_advance_deg) / 180.0 + time_us)
+
+            sparks.extend(
+                (
+                    SparkPulse(spark_1_time, "spark_1"),
+                    SparkPulse(spark_2_time, "spark_2"),
+                    SparkPulse(spark_3_time, "spark_3"),
+                )
             )
 
-            t_a = d1a / 2 + t_s / 2
+    return sorted(set(sparks), key=lambda pulse: (pulse.time_us, pulse.column))
 
-            spark_release_time_2 = int(t_a * (45 - spark_adv) / 180 + time_us)
 
-            f_a = (315 + spark_adv) / (d1a) + (45 - spark_adv) / (2 * d1a - d1b)
+def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]:
+    required = {"time_us", "revolution_index", "crank_angle_deg"}
+    if not required <= set(df.columns):
+        return []
 
-            t_a = 1 / f_a
+    points: list[EncoderPoint] = []
+    angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"])
+    for row in angle_df.itertuples():
+        revolution_index = int(row.revolution_index)
+        if revolution_index < 0:
+            continue
+        angle = revolution_index * 360.0 + float(row.crank_angle_deg)
+        points.append(
+            EncoderPoint(time_us=int(row.time_us), continuous_angle_deg=angle)
+        )
 
-            spark_release_time_3 = int(2 * (t_a * (45 - spark_adv)) + time_us)
+    deduped: dict[int, EncoderPoint] = {}
+    for point in sorted(points, key=lambda item: item.time_us):
+        deduped[point.time_us] = point
+    return list(deduped.values())
 
-    spark_1 = 1 if spark_release_time_1 == time_us and sync_ok else 0
+
-    spark_2 = 1 if spark_release_time_2 == time_us and sync_ok else 0
+def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None:
-    spark_3 = 1 if spark_release_time_3 == time_us and sync_ok else 0
+    if not points:
-    rows.append(
+        return None
-        {
+
-            "time_us": time_us,
+    times = [point.time_us for point in points]
-            "crank": crank,
+    index = bisect_left(times, time_us)
-            "cam": cam,
+
-            "spark_1": spark_1,
+    if index < len(points) and points[index].time_us == time_us:
-            "spark_2": spark_2,
+        return points[index].continuous_angle_deg
-            "spark_3": spark_3,
+    if index == 0 or index == len(points):
+        return None
+
+    before = points[index - 1]
+    after = points[index]
+    span = after.time_us - before.time_us
+    if span <= 0:
+        return before.continuous_angle_deg
+
+    fraction = (time_us - before.time_us) / span
+    return before.continuous_angle_deg + fraction * (
+        after.continuous_angle_deg - before.continuous_angle_deg
+    )
+
+
+def add_spark_rows(df: pd.DataFrame, sparks: list[SparkPulse]) -> pd.DataFrame:
+    for column in SPARK_COLUMNS:
+        df[column] = 0
+    for column in INTERPOLATED_COLUMNS:
+        if column not in df.columns:
+            df[column] = pd.NA
+
+    rows_by_time = {
+        int(row.time_us): index for index, row in enumerate(df.itertuples())
     }
+    new_rows_by_time: dict[int, dict[str, object]] = {}
+
+    for pulse in sparks:
+        if pulse.time_us in rows_by_time:
+            df.at[rows_by_time[pulse.time_us], pulse.column] = 1
+            continue
+
+        if pulse.time_us not in new_rows_by_time:
+            row = {column: pd.NA for column in df.columns}
+            row["time_us"] = pulse.time_us
+            for column in SPARK_COLUMNS:
+                row[column] = 0
+            new_rows_by_time[pulse.time_us] = row
+        new_rows_by_time[pulse.time_us][pulse.column] = 1
+
+    if new_rows_by_time:
+        new_rows = list(new_rows_by_time.values())
+        df = pd.concat(
+            [df, pd.DataFrame(new_rows, columns=df.columns)], ignore_index=True
         )
 
-out_df = pd.DataFrame(rows)
+    return df.sort_values("time_us", kind="stable").reset_index(drop=True)
-base_name, _ = run_data_path.stem.rsplit("_", 1)
+
-output = run_data_path.parent / f"{base_name}_output.csv"
+
-out_df.to_csv(output)
+def add_interpolated_angles(
+    df: pd.DataFrame, points: list[EncoderPoint]
+) -> pd.DataFrame:
+    for index, row in df.iterrows():
+        needs_angle = bool(sum(int(row[column] or 0) for column in SPARK_COLUMNS))
+        if not needs_angle:
+            continue
+
+        continuous_angle = interpolate_angle(points, int(row["time_us"]))
+        if continuous_angle is None:
+            continue
+
+        revolution_index = int(continuous_angle // 360.0)
+        crank_angle = continuous_angle - revolution_index * 360.0
+        df.at[index, "algo1_continuous_angle_deg"] = round(continuous_angle, 6)
+        df.at[index, "algo1_revolution_index"] = revolution_index
+        df.at[index, "algo1_crank_angle_deg"] = round(crank_angle, 6)
+
+    return df
+
+
+def output_path_for(input_path: Path, output_path: Path | None) -> Path:
+    if output_path is not None:
+        return output_path.expanduser()
+
+    stem = input_path.stem
+    if stem.endswith("_rpm"):
+        stem = stem[: -len("_rpm")]
+    return input_path.with_name(f"{stem}_algo1.csv")
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Run algo1 crank/cam angle prediction against one RPM dataset."
+    )
+    parser.add_argument(
+        "file",
+        type=Path,
+        help="RPM CSV file. Bare names are searched in recorder RPM folders.",
+    )
+    parser.add_argument("-o", "--output", type=Path, help="Output CSV path.")
+    parser.add_argument(
+        "--crank-gp",
+        choices=GP_COLUMNS,
+        default="gp0_falling",
+        help="GP falling column to treat as crank pulses.",
+    )
+    parser.add_argument(
+        "--cam-gp",
+        choices=GP_COLUMNS,
+        default="gp1_falling",
+        help="GP falling column to treat as cam pulses.",
+    )
+    parser.add_argument(
+        "--spark-advance",
+        type=float,
+        default=25.0,
+        help="Spark advance in crank degrees, matching the original algo1 constant.",
+    )
+    return parser.parse_args()
+
+
+def main() -> int:
+    args = parse_args()
+    if args.crank_gp == args.cam_gp:
+        print(
+            "--crank-gp and --cam-gp must select different GP columns.", file=sys.stderr
+        )
+        return 1
+
+    try:
+        input_path = resolve_input_path(args.file)
+        df = read_rpm(input_path)
+    except (OSError, ValueError) as exc:
+        print(f"Could not read input: {exc}", file=sys.stderr)
+        return 1
+
+    crank_times = gp_times(df, args.crank_gp)
+    cam_times = gp_times(df, args.cam_gp)
+    sparks = generate_sparks(
+        crank_times,
+        cam_times,
+        spark_advance_deg=args.spark_advance,
+    )
+
+    output_df = add_spark_rows(df.copy(), sparks)
+    output_df = add_interpolated_angles(output_df, encoder_points(df))
+
+    output_path = output_path_for(input_path, args.output)
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+    output_df.to_csv(output_path, index=False)
+
+    print(
+        f"Wrote {output_path} "
+        f"crank={args.crank_gp}({len(crank_times)}) "
+        f"cam={args.cam_gp}({len(cam_times)}) "
+        f"sparks={len(sparks)}"
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

src/plot.py

@@ -2,73 +2,380 @@
 # Copyright (C) 2026 Association Exergie <association.exergie@gmail.com>
 # SPDX-License-Identifier: GPL-3.0-or-later
 
-import pandas as pd
+from __future__ import annotations
-from pathlib import Path
+
 import argparse
-from tqdm import tqdm
+import sys
+from bisect import bisect_left
+from dataclasses import dataclass
+from pathlib import Path
+
 import matplotlib.pyplot as plt
+import pandas as pd
+from matplotlib.ticker import MultipleLocator
 
-parser = argparse.ArgumentParser()
-parser.add_argument("directory", type=Path, help="Source data directory")
 
-args = parser.parse_args()
+GP_COLUMNS = ("gp0_falling", "gp1_falling")
+SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3")
+PULSE_STYLES = {
+    "gp0_falling": ("gp0 falling", "tab:green", 5.0, 0.18),
+    "gp1_falling": ("gp1 falling", "tab:red", 5.0, 0.18),
+    "spark_1": ("spark 1", "tab:blue", 2.4, 0.65),
+    "spark_2": ("spark 2", "tab:orange", 2.4, 0.65),
+    "spark_3": ("spark 3", "tab:purple", 2.4, 0.65),
+}
 
-directory: Path = args.directory
 
-if not directory.is_dir():
+@dataclass(frozen=True)
-    parser.error(f"{directory} is not a valid directory")
+class EncoderPoint:
+    time_us: int
+    continuous_angle_deg: float
 
-print(f"Processing data in: {directory}")
 
-files: list[Path] = []
+def read_csv(path: Path) -> pd.DataFrame:
+    df = pd.read_csv(path)
+    if "time_us" not in df.columns:
+        raise ValueError("missing column: time_us")
+    df["time_us"] = df["time_us"].astype("int64")
+    return df
 
-for path in directory.glob("*.csv"):
-    stem = path.stem
 
+def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]:
+    required = {"time_us", "revolution_index", "crank_angle_deg"}
+    if not required <= set(df.columns):
+        return []
+
+    points: list[EncoderPoint] = []
+    angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"])
+    for row in angle_df.itertuples():
+        revolution_index = int(row.revolution_index)
+        if revolution_index < 0:
+            continue
+        continuous_angle = revolution_index * 360.0 + float(row.crank_angle_deg)
+        points.append(EncoderPoint(int(row.time_us), continuous_angle))
+
+    deduped: dict[int, EncoderPoint] = {}
+    for point in sorted(points, key=lambda item: item.time_us):
+        deduped[point.time_us] = point
+    return list(deduped.values())
+
+
+def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None:
+    if not points:
+        return None
+
+    times = [point.time_us for point in points]
+    index = bisect_left(times, time_us)
+    if index < len(points) and points[index].time_us == time_us:
+        return points[index].continuous_angle_deg
+    if index == 0 or index == len(points):
+        return None
+
+    before = points[index - 1]
+    after = points[index]
+    span = after.time_us - before.time_us
+    if span <= 0:
+        return before.continuous_angle_deg
+
+    fraction = (time_us - before.time_us) / span
+    return before.continuous_angle_deg + fraction * (
+        after.continuous_angle_deg - before.continuous_angle_deg
+    )
+
+
+def iter_csv_paths(path: Path) -> list[Path]:
+    if path.is_dir():
+        return sorted(path.rglob("*.csv"))
+    return [path]
+
+
+def has_angle_data(points: list[EncoderPoint]) -> bool:
+    return len(points) >= 2
+
+
+def rpm_rows(df: pd.DataFrame) -> pd.DataFrame:
+    if "rpm_raw" not in df.columns or "rpm_lpf" not in df.columns:
+        return pd.DataFrame()
+    return df.dropna(subset=["rpm_raw", "rpm_lpf"])
+
+
+def plot_time(df: pd.DataFrame, ax: plt.Axes) -> None:
+    rpm_df = rpm_rows(df)
+    if not rpm_df.empty:
+        time_s = rpm_df["time_us"] / 1_000_000.0
+        ax.plot(time_s, rpm_df["rpm_raw"], label="raw", linewidth=0.75, alpha=0.35)
+        ax.plot(time_s, rpm_df["rpm_lpf"], label="lpf", linewidth=1.4)
+    ax.set_xlabel("time (s)")
+
+
+def plot_angle_overlay(df: pd.DataFrame, ax: plt.Axes) -> None:
+    rpm_df = rpm_rows(df)
+    if rpm_df.empty:
+        return
+
+    raw_label_used = False
+    lpf_label_used = False
+    for revolution in sorted(rpm_df["revolution_index"].dropna().astype(int).unique()):
+        if revolution < 0:
+            continue
+        rev_df = rpm_df[rpm_df["revolution_index"].astype(int) == revolution]
+        ax.plot(
+            rev_df["crank_angle_deg"],
+            rev_df["rpm_raw"],
+            color="0.55",
+            linewidth=0.5,
+            alpha=0.12,
+            label="raw" if not raw_label_used else None,
+        )
+        ax.plot(
+            rev_df["crank_angle_deg"],
+            rev_df["rpm_lpf"],
+            linewidth=0.8,
+            alpha=0.28,
+            label="lpf per rev" if not lpf_label_used else None,
+        )
+        raw_label_used = True
+        lpf_label_used = True
+
+    ax.set_xlim(0, 360)
+    ax.set_xlabel("crank angle (deg)")
+
+
+def plot_angle_continuous(df: pd.DataFrame, ax: plt.Axes) -> None:
+    rpm_df = rpm_rows(df)
+    if rpm_df.empty:
+        return
+
+    angle_df = rpm_df.dropna(subset=["revolution_index", "crank_angle_deg"])
+    angle_df = angle_df[angle_df["revolution_index"].astype(int) >= 0]
+    if angle_df.empty:
+        return
+
+    continuous_angle = (
+        angle_df["revolution_index"].astype(float) * 360.0
+        + angle_df["crank_angle_deg"].astype(float)
+    )
+    ax.plot(continuous_angle, angle_df["rpm_raw"], label="raw", linewidth=0.75, alpha=0.35)
+    ax.plot(continuous_angle, angle_df["rpm_lpf"], label="lpf", linewidth=1.4)
+
+    max_angle = float(continuous_angle.max())
+    turn_angle = 0.0
+    while turn_angle <= max_angle:
+        ax.axvline(turn_angle, color="0.7", linewidth=0.5, alpha=0.35)
+        turn_angle += 360.0
+
+    ax.set_xlim(0, max_angle)
+    ax.xaxis.set_major_locator(MultipleLocator(360.0))
+    ax.xaxis.set_minor_locator(MultipleLocator(180.0))
+    ax.set_xlabel("continuous crank angle (deg)")
+
+
+def selected_pulse_columns(
+    df: pd.DataFrame,
+    gp_mode: str,
+    spark_columns: list[str],
+) -> list[str]:
+    columns: list[str] = []
+    if gp_mode in {"gp0", "both"} and "gp0_falling" in df.columns:
+        columns.append("gp0_falling")
+    if gp_mode in {"gp1", "both"} and "gp1_falling" in df.columns:
+        columns.append("gp1_falling")
+    columns.extend(column for column in spark_columns if column in df.columns)
+    return columns
+
+
+def pulse_x_value(
+    row: pd.Series,
+    *,
+    x_axis: str,
+    points: list[EncoderPoint],
+) -> float | None:
+    time_us = int(row["time_us"])
+    if x_axis == "time":
+        return time_us / 1_000_000.0
+
+    continuous_angle = interpolate_angle(points, time_us)
+    if continuous_angle is None:
+        return None
+    if x_axis in {"angle", "angle-overlay"}:
+        return continuous_angle % 360.0
+    return continuous_angle
+
+
+def plot_pulse_bars(
+    df: pd.DataFrame,
+    ax: plt.Axes,
+    *,
+    x_axis: str,
+    gp_mode: str,
+    spark_columns: list[str],
+    points: list[EncoderPoint],
+) -> None:
+    for column in selected_pulse_columns(df, gp_mode, spark_columns):
+        label, color, linewidth, alpha = PULSE_STYLES[column]
+        pulse_df = df.loc[df[column].fillna(0).astype(int) == 1]
+        plotted_label = False
+        for _, row in pulse_df.iterrows():
+            x_value = pulse_x_value(row, x_axis=x_axis, points=points)
+            if x_value is None:
+                continue
+            ax.axvline(
+                x_value,
+                color=color,
+                alpha=alpha,
+                linewidth=linewidth,
+                label=label if not plotted_label else None,
+            )
+            plotted_label = True
+
+
+def plot_file(
+    csv_path: Path,
+    output: Path | None,
+    *,
+    x_axis: str,
+    gp_mode: str,
+    spark_columns: list[str],
+) -> bool:
     try:
-        base_name, channel = stem.rsplit("_", 1)
+        df = read_csv(csv_path)
-    except ValueError:
+    except (OSError, ValueError) as exc:
-        print(f"Skipping badly named file: {path}")
+        print(f"Could not read {csv_path}: {exc}", file=sys.stderr)
-        continue
+        return False
 
-    if channel != "output":
+    points = encoder_points(df)
-        print(f"Skipping unknown file: {path}")
+    if x_axis.startswith("angle") and not has_angle_data(points):
-        continue
+        print(f"Could not plot crank angle for {csv_path}: missing angle columns", file=sys.stderr)
+        return False
 
-    files.append(path)
+    fig, ax = plt.subplots(figsize=(12, 5))
+    if x_axis in {"angle", "angle-overlay"}:
+        plot_angle_overlay(df, ax)
+    elif x_axis == "angle-continuous":
+        plot_angle_continuous(df, ax)
+    else:
+        plot_time(df, ax)
 
-
+    plot_pulse_bars(
-for file in files:
+        df,
-    output_df = pd.read_csv(file).set_index("time_us", drop=False)
+        ax,
-    fig, ax = plt.subplots(figsize=(12, 6))
+        x_axis=x_axis,
-
+        gp_mode=gp_mode,
-    # RPM
+        spark_columns=spark_columns,
-    ax.plot(output_df["time_us"] / 1_000_000, output_df["crank"], label="crank")
+        points=points,
-    ax.plot(
-        output_df["time_us"] / 1_000_000, output_df["cam"], label="cam", color="red"
     )
-    ax.plot(
-        output_df["time_us"] / 1_000_000,
-        output_df["spark_1"],
-        label="spark_1",
-        color="green",
-    )
-    ax.plot(
-        output_df["time_us"] / 1_000_000,
-        output_df["spark_2"],
-        label="spark_2",
-        color="yellow",
-    )
-    ax.plot(
-        output_df["time_us"] / 1_000_000,
-        output_df["spark_3"],
-        label="spark_3",
-        color="orange",
-    )
-    ax.set_xlabel("Time (s)")
-    ax.set_ylim(0, 1)
-    ax.grid(True)
 
+    ax.set_ylabel("RPM")
+    ax.set_title(str(csv_path))
+    ax.grid(alpha=0.25)
+    if x_axis == "angle-continuous":
+        ax.grid(which="minor", axis="x", alpha=0.15, linewidth=0.6)
+    ax.legend()
     fig.tight_layout()
+
+    if output is not None:
+        output.parent.mkdir(parents=True, exist_ok=True)
+        fig.savefig(output, dpi=160)
+        print(f"Wrote {output}")
+        plt.close(fig)
+    else:
         plt.show()
         plt.close(fig)
+
+    return True
+
+
+def output_path_for(csv_path: Path, output: Path | None, multiple: bool) -> Path | None:
+    if output is None:
+        return None
+    if not multiple:
+        return output
+    return output / f"{csv_path.stem}.png"
+
+
+def parse_spark_columns(value: str) -> list[str]:
+    if value == "none":
+        return []
+    if value == "all":
+        return list(SPARK_COLUMNS)
+
+    columns = [column.strip() for column in value.split(",") if column.strip()]
+    if not columns:
+        raise argparse.ArgumentTypeError("must specify at least one spark column")
+
+    invalid = sorted(set(columns) - set(SPARK_COLUMNS))
+    if invalid:
+        valid = ", ".join(("none", "all", *SPARK_COLUMNS))
+        raise argparse.ArgumentTypeError(
+            f"invalid spark column(s): {', '.join(invalid)}; choose from {valid}"
+        )
+
+    deduped: list[str] = []
+    for column in columns:
+        if column not in deduped:
+            deduped.append(column)
+    return deduped
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Plot algo1 RPM output CSV files.")
+    parser.add_argument("csv_path", type=Path, help="Algo1 CSV file or directory of CSV files.")
+    parser.add_argument("-o", "--output", type=Path, help="Optional image output path or directory.")
+    parser.add_argument(
+        "--x",
+        choices=("time", "angle", "angle-overlay", "angle-continuous"),
+        default="time",
+        help="Plot RPM against time, overlaid 0-360 crank angle, or continuous crank angle.",
+    )
+    parser.add_argument(
+        "--gp",
+        choices=("none", "gp0", "gp1", "both"),
+        default="both",
+        help="Overlay GP falling events as vertical bars.",
+    )
+    parser.add_argument(
+        "--sparks",
+        type=parse_spark_columns,
+        default="all",
+        help="Overlay spark events: none, all, one spark, or a comma-separated subset.",
+    )
+    return parser.parse_args()
+
+
+def main() -> int:
+    args = parse_args()
+
+    if not args.csv_path.exists():
+        print(f"Path does not exist: {args.csv_path}", file=sys.stderr)
+        return 1
+
+    csv_paths = iter_csv_paths(args.csv_path)
+    if not csv_paths:
+        print(f"No CSV files found under {args.csv_path}", file=sys.stderr)
+        return 1
+
+    multiple = len(csv_paths) > 1
+    if multiple and args.output is not None and args.output.suffix:
+        print("When plotting a directory, --output must be a directory.", file=sys.stderr)
+        return 1
+
+    spark_columns = (
+        args.sparks if isinstance(args.sparks, list) else parse_spark_columns(args.sparks)
+    )
+    failures = 0
+    for csv_path in csv_paths:
+        output = output_path_for(csv_path, args.output, multiple)
+        if not plot_file(
+            csv_path,
+            output,
+            x_axis=args.x,
+            gp_mode=args.gp,
+            spark_columns=spark_columns,
+        ):
+            failures += 1
+
+    return 1 if failures else 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

src/spark_error.py

@@ -0,0 +1,837 @@
+# Copyright (C) 2026 Hector van der Aa <hector@h3cx.dev>
+# Copyright (C) 2026 Association Exergie <association.exergie@gmail.com>
+# SPDX-License-Identifier: GPL-3.0-or-later
+
+from __future__ import annotations
+
+import argparse
+import math
+import statistics
+import sys
+from bisect import bisect_left, bisect_right
+from dataclasses import dataclass
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import pandas as pd
+
+
+GP_COLUMNS = ("gp0_falling", "gp1_falling")
+SPARK_COLUMNS = ("spark_1", "spark_2", "spark_3")
+DEFAULT_PPR = 256
+CURVE_STYLES = {
+    "spark_1": {"linestyle": (0, (1.0, 1.2)), "linewidth": 1.7},
+    "spark_2": {"linestyle": (0, (1.0, 2.2)), "linewidth": 2.3},
+    "spark_3": {"linestyle": (0, (1.0, 3.2)), "linewidth": 2.9},
+}
+
+
+@dataclass(frozen=True)
+class EncoderPoint:
+    time_us: int
+    continuous_angle_deg: float
+
+
+@dataclass(frozen=True)
+class SparkError:
+    spark_column: str
+    spark_time_us: int
+    reference_crank_time_us: int
+    expected_delta_deg: float
+    actual_delta_deg: float
+    error_deg: float
+
+
+@dataclass(frozen=True)
+class AngleEvent:
+    column: str
+    time_us: int
+    continuous_angle_deg: float
+
+
+@dataclass(frozen=True)
+class ErrorStats:
+    count: int
+    minimum: float
+    maximum: float
+    mean: float
+    median: float
+    stdev: float
+    rms: float
+    abs_mean: float
+    abs_max: float
+
+
+def read_algo1(path: Path) -> pd.DataFrame:
+    df = pd.read_csv(path)
+    missing = {"time_us", "revolution_index", "crank_angle_deg"} - set(df.columns)
+    if missing:
+        raise ValueError(f"missing columns: {', '.join(sorted(missing))}")
+
+    spark_missing = set(SPARK_COLUMNS) - set(df.columns)
+    if spark_missing:
+        raise ValueError(f"missing spark columns: {', '.join(sorted(spark_missing))}")
+
+    df["time_us"] = df["time_us"].astype("int64")
+    return df.sort_values("time_us", kind="stable").reset_index(drop=True)
+
+
+def encoder_points(df: pd.DataFrame) -> list[EncoderPoint]:
+    angle_df = df.dropna(subset=["revolution_index", "crank_angle_deg"])
+    points: list[EncoderPoint] = []
+
+    for row in angle_df.itertuples():
+        revolution_index = int(row.revolution_index)
+        if revolution_index < 0:
+            continue
+        continuous_angle = revolution_index * 360.0 + float(row.crank_angle_deg)
+        points.append(EncoderPoint(int(row.time_us), continuous_angle))
+
+    deduped: dict[int, EncoderPoint] = {}
+    for point in sorted(points, key=lambda item: item.time_us):
+        deduped[point.time_us] = point
+    return list(deduped.values())
+
+
+def interpolate_angle(points: list[EncoderPoint], time_us: int) -> float | None:
+    if not points:
+        return None
+
+    times = [point.time_us for point in points]
+    index = bisect_left(times, time_us)
+    if index < len(points) and points[index].time_us == time_us:
+        return points[index].continuous_angle_deg
+    if index == 0 or index == len(points):
+        return None
+
+    before = points[index - 1]
+    after = points[index]
+    span_us = after.time_us - before.time_us
+    if span_us <= 0:
+        return before.continuous_angle_deg
+
+    fraction = (time_us - before.time_us) / span_us
+    return before.continuous_angle_deg + fraction * (
+        after.continuous_angle_deg - before.continuous_angle_deg
+    )
+
+
+def pulse_times(df: pd.DataFrame, column: str) -> list[int]:
+    pulse_rows = df.loc[df[column].fillna(0).astype(int) == 1]
+    return sorted(int(row.time_us) for row in pulse_rows.itertuples())
+
+
+def spark_events(df: pd.DataFrame, spark_columns: list[str]) -> list[tuple[int, str]]:
+    events: list[tuple[int, str]] = []
+    for column in spark_columns:
+        for time_us in pulse_times(df, column):
+            events.append((time_us, column))
+    return sorted(events)
+
+
+def previous_time(times: list[int], time_us: int) -> int | None:
+    index = bisect_right(times, time_us) - 1
+    if index < 0:
+        return None
+    return times[index]
+
+
+def wrap_signed_180(angle_deg: float) -> float:
+    return (angle_deg + 180.0) % 360.0 - 180.0
+
+
+def calculate_errors(
+    df: pd.DataFrame,
+    *,
+    crank_gp: str,
+    spark_columns: list[str],
+    spark_advance_deg: float,
+) -> list[SparkError]:
+    points = encoder_points(df)
+    crank_times = pulse_times(df, crank_gp)
+    expected_delta_deg = 45.0 - spark_advance_deg
+    errors: list[SparkError] = []
+
+    for spark_time_us, spark_column in spark_events(df, spark_columns):
+        reference_time_us = previous_time(crank_times, spark_time_us)
+        if reference_time_us is None:
+            continue
+
+        reference_angle = interpolate_angle(points, reference_time_us)
+        spark_angle = interpolate_angle(points, spark_time_us)
+        if reference_angle is None or spark_angle is None:
+            continue
+
+        actual_delta_deg = spark_angle - reference_angle
+        error_deg = wrap_signed_180(actual_delta_deg - expected_delta_deg)
+        errors.append(
+            SparkError(
+                spark_column=spark_column,
+                spark_time_us=spark_time_us,
+                reference_crank_time_us=reference_time_us,
+                expected_delta_deg=expected_delta_deg,
+                actual_delta_deg=actual_delta_deg,
+                error_deg=error_deg,
+            )
+        )
+
+    return errors
+
+
+def stats_for(values: list[float]) -> ErrorStats | None:
+    if not values:
+        return None
+
+    return ErrorStats(
+        count=len(values),
+        minimum=min(values),
+        maximum=max(values),
+        mean=statistics.fmean(values),
+        median=statistics.median(values),
+        stdev=statistics.stdev(values) if len(values) >= 2 else 0.0,
+        rms=math.sqrt(statistics.fmean(value * value for value in values)),
+        abs_mean=statistics.fmean(abs(value) for value in values),
+        abs_max=max(abs(value) for value in values),
+    )
+
+
+def format_stats(label: str, stats: ErrorStats | None) -> list[str]:
+    if stats is None:
+        return [f"{label}: no valid spark samples"]
+
+    return [
+        f"{label}:",
+        f"  count: {stats.count}",
+        f"  min error deg: {stats.minimum:.6f}",
+        f"  max error deg: {stats.maximum:.6f}",
+        f"  mean error deg: {stats.mean:.6f}",
+        f"  median error deg: {stats.median:.6f}",
+        f"  stdev error deg: {stats.stdev:.6f}",
+        f"  rms error deg: {stats.rms:.6f}",
+        f"  mean abs error deg: {stats.abs_mean:.6f}",
+        f"  max abs error deg: {stats.abs_max:.6f}",
+    ]
+
+
+def build_report(
+    input_path: Path,
+    errors: list[SparkError],
+    *,
+    crank_gp: str,
+    spark_columns: list[str],
+    spark_advance_deg: float,
+    ppr: int,
+    band_size_deg: float,
+) -> str:
+    expected_delta_deg = 45.0 - spark_advance_deg
+    encoder_step_deg = 360.0 / ppr
+    lines = [
+        f"input: {input_path}",
+        f"crank reference GP: {crank_gp}",
+        f"spark columns: {', '.join(spark_columns)}",
+        f"spark advance deg: {spark_advance_deg:.6f}",
+        f"expected angle after 45 BTDC crank pulse deg: {expected_delta_deg:.6f}",
+        f"encoder ppr: {ppr}",
+        f"encoder step deg: {encoder_step_deg:.6f}",
+        "",
+    ]
+
+    all_errors = [error.error_deg for error in errors]
+    lines.extend(format_stats("all sparks", stats_for(all_errors)))
+
+    for column in spark_columns:
+        column_errors = [
+            error.error_deg for error in errors if error.spark_column == column
+        ]
+        lines.append("")
+        lines.extend(format_stats(column, stats_for(column_errors)))
+
+    lines.append("")
+    lines.append("samples:")
+    lines.append(
+        "spark,time_us,reference_crank_time_us,error_band,expected_delta_deg,"
+        "actual_delta_deg,error_deg"
+    )
+    for error in errors:
+        lines.append(
+            f"{error.spark_column},{error.spark_time_us},"
+            f"{error.reference_crank_time_us},"
+            f"{band_label(error.error_deg, band_size_deg)},"
+            f"{error.expected_delta_deg:.6f},"
+            f"{error.actual_delta_deg:.6f},"
+            f"{error.error_deg:.6f}"
+        )
+
+    return "\n".join(lines) + "\n"
+
+
+def plot_distribution(
+    errors: list[SparkError],
+    output_path: Path,
+    *,
+    spark_columns: list[str],
+    band_size_deg: float,
+) -> None:
+    fig, axes = plt.subplots(
+        len(spark_columns),
+        1,
+        figsize=(10, max(3.0, 2.5 * len(spark_columns))),
+        sharex=True,
+    )
+    if len(spark_columns) == 1:
+        axes = [axes]
+
+    all_values = [error.error_deg for error in errors]
+    if all_values:
+        minimum = math.floor(min(all_values) / band_size_deg) * band_size_deg
+        maximum = math.ceil(max(all_values) / band_size_deg) * band_size_deg
+        bins = []
+        edge = minimum
+        while edge <= maximum + band_size_deg * 0.5:
+            bins.append(edge)
+            edge += band_size_deg
+        if len(bins) < 2:
+            bins.append(minimum + band_size_deg)
+    else:
+        bins = 10
+
+    for ax, column in zip(axes, spark_columns):
+        values = [error.error_deg for error in errors if error.spark_column == column]
+        if values:
+            ax.hist(
+                values,
+                bins=bins,
+                density=True,
+                alpha=0.75,
+                color="tab:blue",
+                edgecolor="black",
+            )
+        else:
+            ax.text(
+                0.5,
+                0.5,
+                "No valid samples",
+                ha="center",
+                va="center",
+                transform=ax.transAxes,
+            )
+
+        ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65)
+        ax.text(
+            0.0,
+            0.92,
+            "target",
+            ha="center",
+            va="top",
+            transform=ax.get_xaxis_transform(),
+            fontsize=8,
+        )
+        ax.set_ylabel("density")
+        ax.set_title(column)
+        ax.grid(alpha=0.25)
+
+    axes[-1].set_xlabel("angular error (deg)")
+    fig.suptitle(f"Spark angular error distribution ({band_size_deg:g} deg bins)")
+    fig.tight_layout()
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(output_path, dpi=160)
+    plt.close(fig)
+
+
+def density_curve(values: list[float], points: int = 240) -> tuple[list[float], list[float]]:
+    minimum = min(values)
+    maximum = max(values)
+    count = len(values)
+
+    if count >= 2:
+        stdev = statistics.stdev(values)
+    else:
+        stdev = 0.0
+
+    span = maximum - minimum
+    bandwidth = 1.06 * stdev * count ** (-1.0 / 5.0) if stdev > 0.0 else 0.0
+    if bandwidth <= 0.0:
+        bandwidth = max(span / 8.0, 1.0)
+
+    margin = max(span * 0.15, bandwidth * 3.0, 1.0)
+    start = minimum - margin
+    stop = maximum + margin
+    step = (stop - start) / (points - 1)
+
+    x_values = [start + index * step for index in range(points)]
+    normalizer = count * bandwidth * math.sqrt(2.0 * math.pi)
+    y_values = [
+        sum(
+            math.exp(-0.5 * ((x_value - value) / bandwidth) ** 2)
+            for value in values
+        )
+        / normalizer
+        for x_value in x_values
+    ]
+    return x_values, y_values
+
+
+def band_edges(error_deg: float, band_size_deg: float = 1.0) -> tuple[float, float]:
+    lower = math.floor(error_deg / band_size_deg) * band_size_deg
+    return lower, lower + band_size_deg
+
+
+def format_band_value(value: float) -> str:
+    if value == 0:
+        return "0"
+    if value.is_integer():
+        return f"{int(value):+d}"
+    return f"{value:+.3f}".rstrip("0").rstrip(".")
+
+
+def band_label(error_deg: float, band_size_deg: float = 1.0) -> str:
+    lower, upper = band_edges(error_deg, band_size_deg)
+    return f"{format_band_value(lower)}_to_{format_band_value(upper)}_deg"
+
+
+def band_directory_name(error_deg: float, band_size_deg: float = 1.0) -> str:
+    return band_label(error_deg, band_size_deg).replace("+", "pos").replace("-", "neg")
+
+
+def angle_events(
+    df: pd.DataFrame,
+    points: list[EncoderPoint],
+    columns: list[str],
+) -> list[AngleEvent]:
+    events: list[AngleEvent] = []
+    for column in columns:
+        if column not in df.columns:
+            continue
+        for time_us in pulse_times(df, column):
+            angle = interpolate_angle(points, time_us)
+            if angle is None:
+                continue
+            events.append(AngleEvent(column, time_us, angle))
+    return sorted(events, key=lambda event: (event.continuous_angle_deg, event.column))
+
+
+def cycle_rows(
+    df: pd.DataFrame,
+    center_angle_deg: float,
+    before_deg: float,
+    after_deg: float,
+) -> pd.DataFrame:
+    if "rpm_raw" not in df.columns or "rpm_lpf" not in df.columns:
+        return pd.DataFrame()
+
+    rows = df.dropna(subset=["revolution_index", "crank_angle_deg", "rpm_raw", "rpm_lpf"]).copy()
+    if rows.empty:
+        return rows
+
+    rows = rows[rows["revolution_index"].astype(int) >= 0]
+    rows["continuous_angle_deg"] = (
+        rows["revolution_index"].astype(float) * 360.0
+        + rows["crank_angle_deg"].astype(float)
+    )
+    rows["relative_angle_deg"] = rows["continuous_angle_deg"] - center_angle_deg
+    return rows[
+        (rows["relative_angle_deg"] >= -before_deg)
+        & (rows["relative_angle_deg"] <= after_deg)
+    ]
+
+
+def plot_cycle(
+    df: pd.DataFrame,
+    points: list[EncoderPoint],
+    error: SparkError,
+    output_path: Path,
+    *,
+    before_deg: float,
+    after_deg: float,
+    band_size_deg: float,
+) -> bool:
+    spark_angle = interpolate_angle(points, error.spark_time_us)
+    if spark_angle is None:
+        return False
+
+    rows = cycle_rows(df, spark_angle, before_deg, after_deg)
+    events = angle_events(df, points, [*GP_COLUMNS, *SPARK_COLUMNS])
+    window_events = [
+        event
+        for event in events
+        if -before_deg <= event.continuous_angle_deg - spark_angle <= after_deg
+    ]
+
+    fig, (rpm_ax, angle_ax, pulse_ax) = plt.subplots(
+        3,
+        1,
+        figsize=(12, 8),
+        sharex=True,
+        gridspec_kw={"height_ratios": [2.0, 1.2, 1.0]},
+    )
+
+    if not rows.empty:
+        rpm_ax.plot(
+            rows["relative_angle_deg"],
+            rows["rpm_raw"],
+            label="rpm raw",
+            linewidth=0.75,
+            alpha=0.35,
+        )
+        rpm_ax.plot(
+            rows["relative_angle_deg"],
+            rows["rpm_lpf"],
+            label="rpm lpf",
+            linewidth=1.4,
+        )
+        angle_ax.plot(
+            rows["relative_angle_deg"],
+            rows["crank_angle_deg"],
+            label="crank angle",
+            color="tab:orange",
+            linewidth=1.2,
+        )
+
+    for event in window_events:
+        relative_angle = event.continuous_angle_deg - spark_angle
+        if event.column == "gp0_falling":
+            pulse_ax.vlines(relative_angle, 0.0, 0.8, color="tab:green", linewidth=1.6)
+        elif event.column == "gp1_falling":
+            pulse_ax.vlines(relative_angle, 1.0, 1.8, color="tab:red", linewidth=1.6)
+        elif event.time_us == error.spark_time_us and event.column == error.spark_column:
+            pulse_ax.vlines(relative_angle, 2.0, 3.0, color="black", linewidth=2.8)
+        else:
+            pulse_ax.vlines(relative_angle, 2.0, 2.8, color="tab:blue", linewidth=1.4, alpha=0.6)
+
+    for ax in (rpm_ax, angle_ax, pulse_ax):
+        ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65)
+        ax.grid(alpha=0.25)
+
+    reference_angle = interpolate_angle(points, error.reference_crank_time_us)
+    if reference_angle is not None:
+        reference_relative = reference_angle - spark_angle
+        target_relative = reference_relative + error.expected_delta_deg
+        for ax in (rpm_ax, angle_ax, pulse_ax):
+            ax.axvline(reference_relative, color="0.35", linestyle="--", linewidth=1.0, alpha=0.65)
+            ax.axvline(target_relative, color="tab:purple", linestyle=":", linewidth=2.0, alpha=0.85)
+        pulse_ax.text(
+            target_relative,
+            3.05,
+            "target",
+            color="tab:purple",
+            ha="center",
+            va="bottom",
+            fontsize=8,
+        )
+
+    rpm_ax.set_ylabel("RPM")
+    rpm_ax.legend(loc="upper right")
+    angle_ax.set_ylabel("crank angle (deg)")
+    angle_ax.set_ylim(0, 360)
+    angle_ax.legend(loc="upper right")
+    pulse_ax.set_ylabel("pulses")
+    pulse_ax.set_yticks([0.4, 1.4, 2.5])
+    pulse_ax.set_yticklabels(["gp0", "gp1", "spark"])
+    pulse_ax.set_ylim(-0.2, 3.2)
+    pulse_ax.set_xlabel("crank angle relative to analyzed spark (deg)")
+    pulse_ax.set_xlim(-before_deg, after_deg)
+    fig.suptitle(
+        f"{error.spark_column} at {error.spark_time_us} us "
+        f"error={error.error_deg:.3f} deg "
+        f"band={band_label(error.error_deg, band_size_deg)}"
+    )
+    fig.tight_layout()
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(output_path, dpi=160)
+    plt.close(fig)
+    return True
+
+
+def write_cycle_plots(
+    df: pd.DataFrame,
+    errors: list[SparkError],
+    output_root: Path,
+    *,
+    before_deg: float,
+    after_deg: float,
+    band_size_deg: float,
+) -> int:
+    points = encoder_points(df)
+    written = 0
+    errors_by_band: dict[str, list[SparkError]] = {}
+
+    for index, error in enumerate(errors, start=1):
+        band_name = band_directory_name(error.error_deg, band_size_deg)
+        errors_by_band.setdefault(band_name, []).append(error)
+        directory = output_root / band_name
+        output_path = directory / (
+            f"{index:04d}_{error.spark_column}_{error.spark_time_us}_"
+            f"err_{error.error_deg:+.3f}.png"
+        )
+        if not plot_cycle(
+            df,
+            points,
+            error,
+            output_path,
+            before_deg=before_deg,
+            after_deg=after_deg,
+            band_size_deg=band_size_deg,
+        ):
+            continue
+        written += 1
+
+    for band_name, band_errors in errors_by_band.items():
+        output_path = output_root / band_name / "_stacked_cycles.png"
+        if plot_stacked_cycles(
+            df,
+            points,
+            band_errors,
+            output_path,
+            before_deg=before_deg,
+            after_deg=after_deg,
+            band_name=band_name,
+        ):
+            written += 1
+
+    return written
+
+
+def plot_stacked_cycles(
+    df: pd.DataFrame,
+    points: list[EncoderPoint],
+    errors: list[SparkError],
+    output_path: Path,
+    *,
+    before_deg: float,
+    after_deg: float,
+    band_name: str,
+) -> bool:
+    if not errors:
+        return False
+
+    fig, (rpm_ax, angle_ax) = plt.subplots(
+        2,
+        1,
+        figsize=(12, 7),
+        sharex=True,
+        gridspec_kw={"height_ratios": [2.0, 1.2]},
+    )
+    plotted = False
+
+    for error in errors:
+        spark_angle = interpolate_angle(points, error.spark_time_us)
+        reference_angle = interpolate_angle(points, error.reference_crank_time_us)
+        if spark_angle is None or reference_angle is None:
+            continue
+
+        rows = cycle_rows(df, spark_angle, before_deg, after_deg)
+        if rows.empty:
+            continue
+
+        rpm_ax.plot(
+            rows["relative_angle_deg"],
+            rows["rpm_raw"],
+            color="0.45",
+            linewidth=0.55,
+            alpha=0.18,
+        )
+        rpm_ax.plot(
+            rows["relative_angle_deg"],
+            rows["rpm_lpf"],
+            color="tab:blue",
+            linewidth=0.8,
+            alpha=0.28,
+        )
+        angle_ax.plot(
+            rows["relative_angle_deg"],
+            rows["crank_angle_deg"],
+            color="tab:orange",
+            linewidth=0.7,
+            alpha=0.22,
+        )
+        target_relative = reference_angle + error.expected_delta_deg - spark_angle
+        for ax in (rpm_ax, angle_ax):
+            ax.axvline(target_relative, color="tab:purple", linestyle=":", linewidth=0.8, alpha=0.18)
+        plotted = True
+
+    if not plotted:
+        plt.close(fig)
+        return False
+
+    for ax in (rpm_ax, angle_ax):
+        ax.axvline(0.0, color="black", linewidth=1.0, alpha=0.65)
+        ax.grid(alpha=0.25)
+
+    rpm_ax.plot([], [], color="0.45", linewidth=0.8, alpha=0.45, label="rpm raw")
+    rpm_ax.plot([], [], color="tab:blue", linewidth=1.0, alpha=0.55, label="rpm lpf")
+    rpm_ax.axvline(float("nan"), color="tab:purple", linestyle=":", linewidth=1.2, label="target")
+    rpm_ax.set_ylabel("RPM")
+    rpm_ax.legend(loc="upper right")
+    angle_ax.set_ylabel("crank angle (deg)")
+    angle_ax.set_ylim(0, 360)
+    angle_ax.set_xlabel("crank angle relative to analyzed spark (deg)")
+    angle_ax.set_xlim(-before_deg, after_deg)
+    fig.suptitle(f"Stacked cycles for {band_name} ({len(errors)} samples)")
+    fig.tight_layout()
+    output_path.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(output_path, dpi=160)
+    plt.close(fig)
+    return True
+
+
+def default_output_path(input_path: Path, suffix: str) -> Path:
+    stem = input_path.stem
+    return input_path.with_name(f"{stem}_{suffix}")
+
+
+def parse_spark_columns(value: str) -> list[str]:
+    if value == "all":
+        return list(SPARK_COLUMNS)
+
+    columns = [column.strip() for column in value.split(",") if column.strip()]
+    if not columns:
+        raise argparse.ArgumentTypeError("must specify at least one spark column")
+
+    invalid = sorted(set(columns) - set(SPARK_COLUMNS))
+    if invalid:
+        valid = ", ".join((*SPARK_COLUMNS, "all"))
+        raise argparse.ArgumentTypeError(
+            f"invalid spark column(s): {', '.join(invalid)}; choose from {valid}"
+        )
+
+    deduped: list[str] = []
+    for column in columns:
+        if column not in deduped:
+            deduped.append(column)
+    return deduped
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(
+        description="Measure algo1 spark angular error from an algo1 output CSV."
+    )
+    parser.add_argument("file", type=Path, help="CSV generated by src/algo1.py.")
+    parser.add_argument(
+        "--crank-gp",
+        choices=GP_COLUMNS,
+        default="gp0_falling",
+        help="GP falling column used as the crank pulse reference.",
+    )
+    parser.add_argument(
+        "--sparks",
+        type=parse_spark_columns,
+        default="all",
+        help="Spark columns to analyze: all, one spark, or a comma-separated subset.",
+    )
+    parser.add_argument(
+        "--spark-advance",
+        type=float,
+        default=25.0,
+        help="Test spark advance in crank degrees before TDC.",
+    )
+    parser.add_argument(
+        "--ppr",
+        type=int,
+        default=DEFAULT_PPR,
+        help="Rotary encoder pulses per revolution, used for report context.",
+    )
+    parser.add_argument("-o", "--output", type=Path, help="Distribution plot path.")
+    parser.add_argument("--log", type=Path, help="Text report path.")
+    parser.add_argument(
+        "--cycle-plots",
+        type=Path,
+        help="Optional root folder for per-spark engine-cycle plots grouped by error band.",
+    )
+    parser.add_argument(
+        "--cycle-before-deg",
+        type=float,
+        default=540.0,
+        help="Crank-angle range before the analyzed spark for each cycle plot.",
+    )
+    parser.add_argument(
+        "--cycle-after-deg",
+        type=float,
+        default=180.0,
+        help="Crank-angle range after the analyzed spark for each cycle plot.",
+    )
+    parser.add_argument(
+        "--band-size",
+        type=float,
+        default=1.0,
+        help="Angular error band size in degrees for categorized cycle plot folders.",
+    )
+    return parser.parse_args()
+
+
+def main() -> int:
+    args = parse_args()
+    input_path = args.file.expanduser()
+    if not input_path.exists():
+        print(f"Path does not exist: {input_path}", file=sys.stderr)
+        return 1
+    if args.band_size <= 0.0:
+        print("--band-size must be greater than zero.", file=sys.stderr)
+        return 1
+    if args.cycle_before_deg <= 0.0:
+        print("--cycle-before-deg must be greater than zero.", file=sys.stderr)
+        return 1
+    if args.cycle_after_deg <= 0.0:
+        print("--cycle-after-deg must be greater than zero.", file=sys.stderr)
+        return 1
+
+    spark_columns = (
+        args.sparks if isinstance(args.sparks, list) else parse_spark_columns(args.sparks)
+    )
+
+    try:
+        df = read_algo1(input_path)
+    except (OSError, ValueError) as exc:
+        print(f"Could not read input: {exc}", file=sys.stderr)
+        return 1
+
+    if args.crank_gp not in df.columns:
+        print(f"Missing crank GP column: {args.crank_gp}", file=sys.stderr)
+        return 1
+
+    errors = calculate_errors(
+        df,
+        crank_gp=args.crank_gp,
+        spark_columns=spark_columns,
+        spark_advance_deg=args.spark_advance,
+    )
+    report = build_report(
+        input_path,
+        errors,
+        crank_gp=args.crank_gp,
+        spark_columns=spark_columns,
+        spark_advance_deg=args.spark_advance,
+        ppr=args.ppr,
+        band_size_deg=args.band_size,
+    )
+
+    log_path = args.log.expanduser() if args.log else default_output_path(input_path, "spark_error.txt")
+    plot_path = args.output.expanduser() if args.output else default_output_path(input_path, "spark_error.png")
+
+    log_path.parent.mkdir(parents=True, exist_ok=True)
+    log_path.write_text(report)
+    plot_distribution(
+        errors,
+        plot_path,
+        spark_columns=spark_columns,
+        band_size_deg=args.band_size,
+    )
+
+    print(report, end="")
+    print(f"Wrote {log_path}")
+    print(f"Wrote {plot_path}")
+    if args.cycle_plots is not None:
+        cycle_root = args.cycle_plots.expanduser()
+        written = write_cycle_plots(
+            df,
+            errors,
+            cycle_root,
+            before_deg=args.cycle_before_deg,
+            after_deg=args.cycle_after_deg,
+            band_size_deg=args.band_size,
+        )
+        print(f"Wrote {written} categorized cycle plots under {cycle_root}")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())