"""Small NumPy-only tools; positions use native zero-based FITS x,y centers."""

import math
import re
import numpy as np


def header_cards(raw):
    result = {}
    for offset in range(0, len(raw), 80):
        card = raw[offset : offset + 80].decode("ascii")
        if card[8:10] != "= ":
            continue
        key = card[:8].strip()
        value = card[10:]
        if value.lstrip().startswith("'"):
            match = re.match(r"\s*'((?:[^']|'')*)'", value)
            result[key] = match.group(1).replace("''", "'").strip() if match else value
        else:
            token = value.split("/")[0].strip()
            try:
                result[key] = float(token.replace("D", "E"))
            except ValueError:
                result[key] = token
    return result


def bilinear(image, x, y):
    x, y = np.broadcast_arrays(x, y)
    left, top = np.floor(x).astype(int), np.floor(y).astype(int)
    fraction_x, fraction_y = x - left, y - top
    result = np.zeros(x.shape, dtype=float)
    total_weight = np.zeros(x.shape, dtype=float)
    for offset_y in (0, 1):
        for offset_x in (0, 1):
            column, row = left + offset_x, top + offset_y
            weight = (fraction_x if offset_x else 1 - fraction_x) * (
                fraction_y if offset_y else 1 - fraction_y
            )
            inside = (
                (column >= 0)
                & (row >= 0)
                & (column < image.shape[1])
                & (row < image.shape[0])
            )
            values = image[
                np.clip(row, 0, image.shape[0] - 1),
                np.clip(column, 0, image.shape[1] - 1),
            ]
            valid = inside & np.isfinite(values)
            result += np.where(valid, values * weight, 0)
            total_weight += np.where(valid, weight, 0)
    return np.where(total_weight > 1 - 1e-8, result, np.nan)


def aperture(image, x, y, radius=3):
    margin = radius + 4
    left, right = max(0, math.floor(x - margin - 1)), min(image.shape[1], math.ceil(x + margin + 2))
    top, bottom = max(0, math.floor(y - margin - 1)), min(image.shape[0], math.ceil(y + margin + 2))
    image = image[top:bottom, left:right]
    rows, columns = np.indices(image.shape)
    rows, columns = rows + top, columns + left
    distance = np.hypot(columns - x, rows - y)
    background_mask = (
        (distance >= radius + 2) & (distance <= radius + 4) & np.isfinite(image)
    )
    if np.count_nonzero(background_mask) < 12:
        raise ValueError("Insufficient local sky samples")
    background = np.median(image[background_mask])
    weights = np.zeros(image.shape)
    for offset_y in (-0.375, -0.125, 0.125, 0.375):
        for offset_x in (-0.375, -0.125, 0.125, 0.375):
            weights += np.hypot(columns + offset_x - x, rows + offset_y - y) <= radius
    weights /= 16
    selected = weights > 0
    if not np.all(np.isfinite(image[selected])):
        raise ValueError("Aperture touches missing/bad pixels")
    return float(np.sum((image[selected] - background) * weights[selected]))


def remeasure_spot(image, original_x, original_y):
    """Local Gaussian + tilted-plane registration, not an Airy/photometric fit.

    A fixed five-pixel-radius window prevents the centroid from walking towards
    the rising stellar halo or a neighbouring speckle. Reported precision is not
    inferred from this unweighted fit; alternate templates test that limitation.
    """
    rows, columns = np.indices(image.shape)
    selected = (np.hypot(columns - original_x, rows - original_y) <= 4) & np.isfinite(
        image
    )
    xx, yy = columns[selected] - original_x, rows[selected] - original_y
    values = image[selected]
    if values.size < 30:
        raise ValueError("Too few pixels for local spot registration")
    best = None
    for offset_x in np.arange(-1, 1.001, 0.2):
        for offset_y in np.arange(-1, 1.001, 0.2):
            for sigma in (0.7, 0.9, 1.1, 1.3, 1.5):
                gaussian = np.exp(
                    -((xx - offset_x) ** 2 + (yy - offset_y) ** 2) / (2 * sigma**2)
                )
                design = np.column_stack((gaussian, np.ones(values.size), xx, yy))
                parameters = np.linalg.lstsq(design, values, rcond=None)[0]
                if parameters[0] <= 0:
                    continue
                score = float(np.sum((values - design @ parameters) ** 2))
                if best is None or score < best[0]:
                    best = (score, offset_x, offset_y, sigma)
    if best is None:
        raise ValueError("No positive local satellite-spot fit")
    if abs(best[1]) >= 0.99 or abs(best[2]) >= 0.99:
        raise ValueError(
            "Satellite-spot fit reaches the one-pixel local-search boundary"
        )
    return original_x + best[1], original_y + best[2]


def spot_template(image, x, y, half_width=6, aperture_radius=3):
    offsets = np.arange(-half_width, half_width + 1)
    xx, yy = np.meshgrid(offsets, offsets)
    stamp = bilinear(image, x + xx, y + yy)
    distance = np.hypot(xx, yy)
    background = np.nanmedian(stamp[distance >= half_width - 1])
    stamp = stamp - background
    if not np.all(np.isfinite(stamp)):
        raise ValueError("Missing satellite-template pixels")
    # Same finite aperture convention for the spot/star calibration and template.
    norm = aperture(stamp, half_width, half_width, aperture_radius)
    if norm <= 0:
        raise ValueError("Nonpositive integrated template normalization")
    return stamp / norm


def insert_template(image_shape, template, x, y):
    rows, columns = np.indices(image_shape)
    half = template.shape[0] // 2
    result = bilinear(template, columns - x + half, rows - y + half)
    return np.nan_to_num(result, nan=0)