"""Our analysis: trace the curves of Fig. 4 of Griv & Yahalom (2026), Astrophys. Space Sci. 371:118
(CC BY 4.0), from Springer's 1170 x 929 JPEG of the figure, and write them on a grid of x.

For each of the four panels (m = 1, 2, 3, 4) the script finds, in every pixel column inside the frame, the
orange pixels of the dashed curve (Im(omega*)/Omega) and the blue pixels of the solid curve (Re(omega*)/Omega),
takes their mean height, and converts pixels to values with the panel's gridlines: y = 0.5 and y = 1 at
rows 242 and 129 in the top panels and 728 and 615 in the bottom ones (226 px per unit, so y = 0 at 355 and
841), and x = 0 and x = 0.5 at columns 97 and 325.5 (left panels) or 699 and 927 (right panels), 457 px per
unit of x. The dashed curve has gaps; values on the grid are interpolated linearly between the nearest
columns that hold orange pixels, at most 40 px away on each side. Readings are good to about +-0.01-0.02.

Usage: python3 -I fig4_curves.py <Fig4.jpg> <out.csv>
The figure: https://media.springernature.com/lw1200/springer-static/image/art%3A10.1007%2Fs10509-026-04651-8/MediaObjects/10509_2026_4651_Fig4_HTML.jpg
Needs Pillow. Reads the image as pixels only.
"""
import csv
import sys

from PIL import Image

img = Image.open(sys.argv[1]).convert("RGB")
W, H = img.size
if (W, H) != (1170, 929):
    sys.exit(f"expected the 1170 x 929 image, got {W} x {H}")
px = img.load()


def is_orange(c):
    r, g, b = c
    return r > 170 and 50 < g < 160 and b < 110 and r - b > 100


def is_blue(c):
    r, g, b = c
    return b > 140 and r < 110 and b - r > 60


PANELS = {  # m: (x0 column, x0.5 column, y0.5 row, y1 row)
    1: (97.0, 325.5, 242.0, 129.0),
    2: (699.0, 927.0, 242.0, 129.0),
    3: (97.0, 325.5, 728.0, 615.0),
    4: (699.0, 927.0, 728.0, 615.0),
}
GRID = [round(i * 0.02, 2) for i in range(51)]


def trace(m):
    x0, xh, y05, y10 = PANELS[m]
    xunit = 2 * (xh - x0)
    unit = (y05 - y10) / 0.5
    yzero = y05 + 0.5 * unit
    im, re = {}, {}
    for col in range(int(x0) + 2, int(x0 + xunit) - 1):
        oy = [y for y in range(int(y10) - 140, int(yzero) + 2) if is_orange(px[col, y])]
        by = [y for y in range(int(y10) - 140, int(yzero) + 2) if is_blue(px[col, y])]
        xv = (col - x0) / xunit
        if oy:
            im[xv] = (yzero - sum(oy) / len(oy)) / unit
        if by:
            re[xv] = (yzero - sum(by) / len(by)) / unit
    return im, re, xunit


def on_grid(samples, xunit, x):
    if not samples:
        return None
    xs = sorted(samples)
    left = [s for s in xs if s <= x]
    right = [s for s in xs if s >= x]
    lim = 40 / xunit
    if left and right and x - left[-1] <= lim and right[0] - x <= lim:
        a, b = left[-1], right[0]
        if b == a:
            return samples[a]
        return samples[a] + (samples[b] - samples[a]) * (x - a) / (b - a)
    if x <= xs[0] and xs[0] - x <= lim:   # at the very left edge, the curve starts at zero
        return samples[xs[0]]
    if x >= xs[-1] and x - xs[-1] <= lim:
        return samples[xs[-1]]
    return None


rows = []
traces = {m: trace(m) for m in PANELS}
for x in GRID:
    row = {"x": f"{x:.2f}"}
    for m, (im, re, xunit) in traces.items():
        v = on_grid(im, xunit, x)
        row[f"im_m{m}"] = "" if v is None else f"{max(v, 0.0):.3f}"
    v = on_grid(traces[1][1], traces[1][2], x)
    row["re_m1"] = "" if v is None else f"{v:.3f}"
    rows.append(row)

with open(sys.argv[2], "w", newline="", encoding="utf-8") as f:
    f.write("# Values read by machine from Fig. 4 of Griv & Yahalom (2026), Astrophys. Space Sci. 371:118, "
            "https://doi.org/10.1007/s10509-026-04651-8, CC BY 4.0. Good to about +-0.01-0.02.\n")
    f.write("# x = k^2 c_r^2 / kappa^2 (the figure's axis); im_mN = Im(omega*)/Omega for azimuthal mode m = N "
            "(dashed curves); re_m1 = Re(omega*)/Omega for m = 1 (solid curve). Negative readings near zero are set to 0.\n")
    w = csv.DictWriter(f, fieldnames=["x", "im_m1", "im_m2", "im_m3", "im_m4", "re_m1"])
    w.writeheader()
    w.writerows(rows)
print("wrote", sys.argv[2], len(rows), "rows")
for r in rows[::5]:
    print(r)