"""
/***************************************************************************
 Free Station
                                 A QGIS plugin
 Topaze
                             -------------------
        begin                : 2025-08-28
        git sha              : $Format:%H$
        copyright            : (C) 2025 by Jean-Marie ARSAC
        email                : jmarsac@arsac.wf
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
"""

import json
import math
import os
from typing import Any, Dict, List, Optional, Tuple

import numpy as np
from qgis.core import QgsProject
from scipy.optimize import least_squares

from topaze.calc.calc_helpers import angle_interne_gr
from topaze.calc.elevation_utils import (
    compute_free_station_z_from_verticals,
    compute_free_station_z_from_verticals_by_only_2_sights,
)
from topaze.file_utils import FileUtils
from topaze.report_helpers import fmt
from topaze.report_utils import ReportUtils
from topaze.toolbelt import PlgLogger
from topaze.toolbelt.i18n import tr
from topaze.topaze_calculator import TopazeCalculator

# --- LSQ helpers -------------------------------------------------------------


def _wrap_signed_gr(d: float) -> float:
    """Wrap grads to (-200, 200]."""
    d = float(d) % 400.0
    return d - 400.0 if d > 200.0 else d


# --- LSQ helpers: construire les paquets d’observations (angles / distances) ---
@staticmethod
def _collect_obs_for_lsq(
    station_id: str,
    obs_data: Dict[str, Any],
    stations_by_id: Dict[str, Dict[str, Any]],
    use_dirs: bool = True,
    use_dists: bool = True,
) -> Tuple[List[Dict[str, Any]], List[Dict[str, Any]]]:
    """
    Construit deux listes d’observations à partir de payload['obs'] :
      • dirs  : [{"to", "ah", "x", "y"}] (angles en grads)
      • dists : [{"to", "d", "x", "y"}]  (distances en mètres)
    en ne gardant que les cibles connues (présentes dans stations_by_id).
    """
    dirs, dists = [], []
    for o in obs_data.get("obs", []):
        if o.get("origine") != station_id:
            continue
        tgt = o.get("cible")
        if tgt not in stations_by_id:
            continue
        sx = stations_by_id[tgt].get("x")
        sy = stations_by_id[tgt].get("y")
        if sx is None or sy is None:
            continue
        if use_dirs:
            ah = o.get("ah")
            if ah is not None:
                dirs.append(
                    {"to": tgt, "ah": float(ah), "x": float(sx), "y": float(sy)}
                )
        if use_dists:
            di = o.get("di")
            if di is not None:
                dists.append(
                    {"to": tgt, "d": float(di), "x": float(sx), "y": float(sy)}
                )
    return dirs, dists


@staticmethod
def _free_station_lsq_xytheta(
    dirs: List[Dict[str, Any]],
    dists: List[Dict[str, Any]],
    x0: float,
    y0: float,
    robust: bool = True,
    loss: str = "soft_l1",
    sigma_dir_gr: float = 0.0003,  # ~3 cc, à adapter à l’instrument
    sigma_dist_m: float = 0.005,  # ~5 mm, à adapter
) -> Tuple[float, float, float, Dict[str, Any]]:
    """
    LSQ mixte sur (X,Y,θ) avec observations d’angles (gr) et distances (m).
    Modèle :
      • direction i :   ah_i ≈ θ + az(X,Y→P_i) + v_i        [grads]
      • distance  j :   d_j ≈ ||(X,Y) − P_j|| + e_j          [m]

    Retourne: (X̂, Ŷ, Θ̂_gr, stats) où
      stats = { "n","dof","emqx_z","rmse_gr","sigma0_gr","std_x_m","std_y_m","std_theta_gr",
                "residuals_dir":[{to, ah_meas_gr, ah_comp_gr, res_gr},...],
                "residuals_dist":[{to, d_meas_m, d_comp_m, res_m},...] }
    """
    use_dirs = len(dirs) > 0
    use_dists = len(dists) > 0
    if not (use_dirs or use_dists):
        raise ValueError(
            tr("No observations available for LSQ (neither directions nor distances).")
        )

    # Matrices / vecteurs pour empiler les résidus (angles en g, distances en m)
    xd = np.array([o["x"] for o in dirs], dtype=float) if use_dirs else np.empty((0,))
    yd = np.array([o["y"] for o in dirs], dtype=float) if use_dirs else np.empty((0,))
    zd = np.array([o["ah"] for o in dirs], dtype=float) if use_dirs else np.empty((0,))

    xt = np.array([o["x"] for o in dists], dtype=float) if use_dists else np.empty((0,))
    yt = np.array([o["y"] for o in dists], dtype=float) if use_dists else np.empty((0,))
    dt = np.array([o["d"] for o in dists], dtype=float) if use_dists else np.empty((0,))

    inv_sig_dir = 1.0 / float(sigma_dir_gr)
    inv_sig_dst = 1.0 / float(sigma_dist_m)

    def residuals(u):
        X, Y, TH = float(u[0]), float(u[1]), float(u[2])  # TH en grads
        res_list = []

        if use_dirs:
            az = np.degrees(np.arctan2(xd - X, yd - Y)) * (200.0 / 180.0)
            az = np.mod(az, 400.0)
            v = zd - (TH + az)
            v = (v + 200.0) % 400.0 - 200.0  # wrap (-200,200]
            res_list.append(inv_glob * v)  # “unitless” pour σ0

        if use_dists:
            r = np.sqrt((xt - X) ** 2 + (yt - Y) ** 2)
            e = dt - r
            res_list.append(inv_sig_dst * e)

        if not res_list:
            return np.zeros((0,))
        return np.concatenate(res_list)

    # initial TH depuis la première direction si dispo, sinon 0
    if use_dirs and zd.size:
        az0 = (np.degrees(np.arctan2(xd[0] - x0, yd[0] - y0)) * (200.0 / 180.0)) % 400.0
        th0 = (zd[0] - az0) % 400.0
    else:
        th = 0.0

    inv_glob = 1.0 / sigma_dir_gr if use_dirs else 1.0
    inv_dl = 1.0 / sigma_dist_m if use_dists else 1.0

    def residuals(u):
        X, Y, TH = float(u[0]), float(u[1]), float(u[2] % 400.0)
        res = []

        if use_dirs and zd.size:
            az = np.degrees(np.arctan2(xd - X, yd - Y)) * (200.0 / 180.0)
            az = np.mod(az, 400.0)
            v = zd - (TH + az)
            v = (v + 200.0) % 400.0 - 200.0
            res.append(v * inv_glob)

        if use_dists and dt.size:
            r = np.sqrt((xt - X) ** 2 + (yt - Y) ** 2)
            e = dt - r
            res.append(e * inv_dl)

        return np.concatenate(res) if res else np.zeros((0,))

    u0 = np.array([x0, y0, th0 if use_dirs and zd.size else 0.0], dtype=float)

    # LSQ robuste par défaut
    if robust:
        out = least_squares(
            residuals,
            u0,
            method="trf",
            loss="soft_l1",
            f_scale=1.0,
            jac="2-point",
            xtol=1e-12,
            ftol=1e-12,
            gtol=1e-12,
            max_nfev=300,
        )
    else:
        out = least_squares(
            residuals,
            u0,
            method="lm",
            jac="2-point",
            xtol=1e-12,
            ftol=1e-12,
            gtol=1e-12,
        )

    Xh, Yh, THh = float(out.x[0]), float(out.x[1]), float(out.x[2] % 400.0)

    # Résidus “non scalés”
    v = residuals(out.x)
    n_obs = len(dirs) + len(dists)
    dof = max(n_obs - 3, 0)
    # σ0 et RMSE (sur résidus scalés)
    rmse = float(np.sqrt(np.mean(v**2))) if n_obs else None
    sigma0 = None if dof == 0 else float(np.sqrt(np.sum(v**2) / dof))

    # Empilement des résidus séparés (non scalés) pour le tableau
    residuals_dir = []
    residuals_dst = []
    if use_dirs and zd.size:
        az = np.degrees(np.arctan2(xd - Xh, yd - Yh)) * (200.0 / 180.0)
        az = np.mod(az, 400.0)
        vdir = zd - (THh + az)
        vdir = (vdir + 200.0) % 400.0 - 200.0
        for o, az_i, vi in zip(dirs, az, vdir):
            residuals_dir.append(
                {
                    "to": o["to"],
                    "ah_meas_gr": float(o["ah"]),
                    "ah_comp_gr": float((THh + az_i) % 400.0),
                    "res_gr": float(vi),
                }
            )
    if use_dists and dt.size:
        rhat = np.sqrt((xt - Xh) ** 2 + (yt - Yh) ** 2)
        vdst = dt - rhat
        for o, rh, ve in zip(dists, rhat, vdst):
            residuals_int = {
                "to": o["to"],
                "d_meas_m": float(o["d"]),
                "d_comp_m": float(rh),
                "res_m": float(ve),
            }
            residuals_dst.append(residuals_int)

    # NB: EMQ X/Y au sens “écart type des inconnues” = racine des diag de (σ0^2 * (JᵀJ)⁻¹)
    std_x = std_y = std_th = None
    if dof > 0 and hasattr(out, "jac") and out.jac is not None:
        JTJ = out.jac.T @ out.jac
        try:
            Cov = sigma0**2 * np.linalg.inv(JTJ)
            std_x = float(np.sqrt(Cov[0, 0]))
            std_y = float(np.sqrt(Cov[1, 1]))
            std_th = float(np.sqrt(Cov[2, 2]))
        except np.linalg.LinAlgError:
            pass

    return (
        Xh,
        Yh,
        THh,
        {
            "n": int(n_obs),
            "dof": int(dof),
            "emqx_z": None,  # rempli plus tard si EMQ Z
            "rmse_gr": rmse,
            "sigma0_gr": sigma0,
            "std_x_m": std_x,
            "std_y_m": float(std_y) if std_y is not None else None,
            "std_theta_gr": std_th,
            "residuals_dir": residuals_dir,
            "residuals_dist": residuals_dst,
        },
    )


def _auto_choose_candidate_by_direction_order(
    A: Tuple[float, float],
    B: Tuple[float, float],
    ah_A: Optional[float],
    ah_B: Optional[float],
    candidates: List[Tuple[float, float]],
    tol_gr: float = 0.01,  # tolérance absolue sur le résidu retenu
    amb_margin_gr: float = 0.001,  # marge d’indécision entre les deux résidus
    amb_zone_gr: float = 0.002,  # zone autour de 0g (et ±200g) rendant le signe peu fiable
    max_tol_gr: float = 0.10,  # si les deux résidus > max_tol → données incohérentes
) -> Optional[Tuple[Tuple[float, float], int, str]]:
    """
    Choose candidate whose signed delta dir(M→A)→dir(M→B) matches measured signed delta.
    Returns (chosen_xy, index, reason) or None if undecidable.

    Exemple concret où elle renvoie None

    ah_A = 123.4560 g, ah_B = 123.4561 g → Δ_meas ≈ 0.0001 g (≈ 0.3″)
    =>  amb_zone_gr le considère comme trop proche de 0g ⇒ None (ambigu).

    Ou encore : deux résidus r1=0.0042 g, r2=0.0043 g
    => |r1−r2|=0.0001 g < amb_margin_gr => None.

    Ou enfin : r1=0.12 g, r2=0.15 g (> max_tol_gr=0.10 g)
    => None (données incohérentes).
    """
    if not candidates or ah_A is None or ah_B is None:
        return None
    if len(candidates) == 1:
        return (tuple(candidates[0]), 0, tr("Only one candidate available."))

    d_meas = _wrap_signed_gr(ah_B - ah_A)

    # Zone d’ambiguïté intrinsèque autour de 0g / ±200g (perte d’information de sens)
    if abs(d_meas) < amb_zone_gr or abs(abs(d_meas) - 200.0) < amb_zone_gr:
        return None

    xa, ya = A
    xb, yb = B

    # Résidus pour les deux candidats
    items = []
    for i, (xm, ym) in enumerate(candidates):
        dir_MA = TopazeCalculator.bearing_gon(xm, ym, xa, ya)
        dir_MB = TopazeCalculator.bearing_gon(xm, ym, xb, yb)
        d_calc = _wrap_signed_gr(dir_MB - dir_MA)
        res = abs(_wrap_signed_gr(d_calc - d_meas))
        items.append((res, i, (xm, ym)))

    # Trie par résidu croissant
    items.sort(key=lambda t: t[0])
    r1, i1, pt1 = items[0]
    r2, i2, pt2 = items[1]

    # Si les deux trop grands → incohérent
    if min(r1, r2) > max_tol_gr:
        return None

    # Si quasi identiques → indécidable
    if abs(r1 - r2) < amb_margin_gr:
        return None

    # Sinon on choisit le meilleur
    reason = tr("Chosen by direction order: signed delta residual = {r:.5f} gr").format(
        r=r1
    )
    return (pt1, i1, reason)


class FreeStation:
    """
    Free-station (AFT method) for Topaze
    - Reads /tmp Topaze temp file "free_station.json" (via FileUtils)
    - Computes M from A,B, AM distance and angle at M (grads) between MA and MB
    - Produces a Markdown report + HTML/PDF (like resection/intersection)

    Example of free_station.json:
    {
      "calcul": ["M"],
      "stations": [
        { "matricule": "A1", "x": 905477.302, "y": 6602977.797 },
        { "matricule": "B1", "x": 905509.941, "y": 6603428.027 }
      ],
      "free_station": {
        "A": "A1",
        "B": "B1",
        "choose": "left"
      },
      "obs": [
        {
          "origine": "M",
          "hi": 1.61,
          "cible": "A1",
          "ah": 100.123,
          "av": null,
          "di": 123.456,
          "hp": null
        },
        {
          "origine": "M",
          "hi": 1.61,
          "cible": "B1",
          "ah": 112.4686,
          "av": null,
          "di": null,
          "hp": null
        }
      ]
    }

    """

    @staticmethod
    @staticmethod
    def refine_free_station_with_lsq_and_report(
        payload: Dict[str, Any],
        obs_data: Dict[str, Any],
        *,
        use_dirs: bool = True,
        use_dists: bool = True,
        robust: bool = True,
        sigma_dir_gr: float = 0.0003,
        sigma_dist_m: float = 0.005,
        out_basename: str = "free_station_report",
    ) -> Dict[str, Any]:
        """
        Exploite toutes les visées (angles et/ou distances) pour affiner (X,Y,θ) par LSQ
        à partir du candidat choisi, puis régénère le Markdown/HTML/PDF.
        """
        res = payload.get("result") or {}
        if "chosen" not in res:
            PlgLogger.log(
                tr("No chosen candidate; LSQ refinement skipped."),
                log_level=1,
                push=True,
            )
            return payload

        X0, Y0 = res["chosen"]
        stations = {s["matricule"]: s for s in obs_data.get("stations", [])}
        dirs, dists = _collect_obs_for_lsq(
            payload["station_id"],
            obs_data,
            stations,
            use_dirs=use_dirs,
            use_dists=use_dists,
        )

        if not dirs and not dists:
            PlgLogger.log(
                tr("No usable observations for LSQ refinement."),
                log_level=1,
                push=True,
            )
            return payload

        try:
            Xh, Yh, THh, stats = _free_station_lsq_xytheta(
                dirs,
                dists,
                x0=X0,
                y0=Y0,
                robust=robust,
                sigma_dir_gr=sigma_dir_gr,
                sigma_dist_m=sigma_dist_m,
            )
        except Exception as ex:
            PlgLogger.log(
                tr("LSQ refinement failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )
            return payload

        payload.setdefault("lsq_block", {})
        payload["lsq_block"].update({"x": Xh, "y": Yh, "theta_gr": THh, "stats": stats})

        # Régénérer le rapport avec section LSQ
        md_report = FreeStation.export_markdown_report(payload)
        out_dir = QgsProject.instance().homePath() + "/rapport/"
        os.makedirs(out_dir, exist_ok=True)
        md_name = out_basename + ".md"
        FileUtils.save_temp_file(md_name, md_report, out_dir)
        md_path = os.path.join(out_dir, md_name)

        try:
            ReportUtils.markdown_to_html(md_path, title=tr("Free Station Report"))
        except Exception as ex:
            PlgLogger.log(
                tr("HTML export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )
        try:
            ReportUtils.markdown_to_pdf(
                md_path, title=tr("Free Station Report"), margins_mm=(7, 15, 7, 15)
            )
        except Exception as ex:
            PlgLogger.log(
                tr("PDF export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )

        PlgLogger.log(
            tr("LSQ refinement done: X={x:.4f}, Y={y:.4f}, θ={t:.4f} g").format(
                x=Xh, y=Yh, t=THh
            ),
            log_level=0,
            push=True,
        )
        return payload

    @staticmethod
    def compute_free_station():
        """
        Root function to compute free-station from free_station.json.
        Similar spirit to Resection.compute_resection().

        Returns:
            payload dict (see exporter); or None-equivalent payload on error.
        """
        data = FileUtils.load_temp_file("free_station.json")
        if not data:
            PlgLogger.log(
                message=tr("No data from free_station.json"),
                log_level=2,
                push=True,
            )
            out_dir = QgsProject.instance().homePath() + "/rapport/"
            md_name = "free_station_report.md"
            FileUtils.remove_temp_file(md_name, out_dir)
            return None

        try:
            obs_data = json.loads(data)
        except Exception as ex:
            PlgLogger.log(
                message=tr("Invalid JSON: {err}").format(err=str(ex)),
                log_level=2,
                push=True,
            )
            return None

        station_id = obs_data["calcul"][0]
        corrections = obs_data.get("corrections", {})
        stations = {s["matricule"]: s for s in obs_data.get("stations", [])}
        fs = obs_data.get("free_station") or {}

        A_id = fs.get("A")
        B_id = fs.get("B")
        choose = (fs.get("choose") or "auto").lower()
        use_apparent_level = bool(corrections.get("niveau_apparent", True))

        if A_id not in stations or B_id not in stations:
            PlgLogger.log(
                tr("Unknown A or B id in 'free_station' block."),
                log_level=2,
                push=True,
            )
            return None

        A = stations[A_id]
        B = stations[B_id]

        # Find the two sights from M to A and M to B
        ah_A = ah_B = None
        di_A = None  # AM distance
        for o in obs_data.get("obs", []):
            if o.get("origine") != station_id:
                continue
            if o.get("cible") == A_id:
                ah_A = o.get("ah")
                di_A = o.get("di")
            elif o.get("cible") == B_id:
                ah_B = o.get("ah")

        if ah_A is None or ah_B is None or di_A is None:
            PlgLogger.log(
                tr("Missing AHs or AM distance: require ah(M→A), ah(M→B), di(M→A)."),
                log_level=2,
                push=True,
            )
            return None

        if A["x"] is None or A["y"] is None or B["x"] is None or B["y"] is None:
            PlgLogger.log(
                tr("Missing coordinates for A or B station."),
                log_level=2,
                push=True,
            )
            return None
        # Extract coordinates
        xa = float(A["x"])
        ya = float(A["y"])
        xb = float(B["x"])
        yb = float(B["y"])

        # Build AM and M_gr
        AM = float(di_A)
        M_gr = angle_interne_gr(float(ah_B) - float(ah_A))

        if abs(M_gr) < 1e-4 or abs(200.0 - M_gr) < 1e-4:
            PlgLogger.log(
                tr("Degenerate triangle: angle at M too small or straight."),
                log_level=2,
                push=True,
            )
            return None

        # Compute free station solution
        try:
            res = FreeStation.free_station_aft_single(
                (xa, ya),
                (xb, yb),
                AM,
                M_gr,
                choose=choose,
                ah_A=float(ah_A),
                ah_B=float(ah_B),
            )
        except Exception as ex:
            PlgLogger.log(
                message=tr("Free station computation failed: {err}").format(
                    err=str(ex)
                ),
                log_level=2,
                push=True,
            )
            return None

        # Choose XY for Z computation
        chosen_xy = None
        if "chosen" in res:
            chosen_xy = tuple(res["chosen"])
        elif res.get("candidates"):
            chosen_xy = tuple(res["candidates"][0])
        else:
            PlgLogger.log(
                tr("No candidate to compute station elevation."),
                log_level=2,
                push=True,
            )
            return None

        x_result, y_result = chosen_xy

        # Compute Z from verticals if available
        z_raw = z_cmp = None
        z_stats = None
        z_block = compute_free_station_z_from_verticals_by_only_2_sights(
            station_id=station_id,
            chosen_xy=chosen_xy,
            stations_by_id=stations,
            obs_list=obs_data.get("obs", []),
            use_apparent_level=use_apparent_level,
            default_hp=0.0,
            weight_mode="inv_dH",
        )

        # Build payload for reporting (déjà dans ton code)
        payload = {
            "station_id": station_id,
            "A": {"id": A.get("matricule", "A"), "x": xa, "y": ya},
            "B": {"id": B.get("matricule", "B"), "x": xb, "y": yb},
            "result": res,
            "z_block": None,
            "use_apparent_level": use_apparent_level,
            "stations": obs_data.get("stations", []),  # utile pour le LSQ helper
        }
        payload["z_block"] = z_block

        # --- LSQ (XY, θ) optionnel si un candidat est déjà fixé automatiquement ---
        lsq_cfg = obs_data.get(
            "lsq", {}
        )  # ex. {"enabled": true, "use_dirs": true, "use_dists": true, "robust": true}
        lsq_enabled = bool(lsq_cfg.get("enabled", False))
        use_dirs = bool(lsq_cfg.get("use_dirs", True))
        use_dists = bool(lsq_cfg.get("use_dists", True))
        robust = bool(lsq_cfg.get("robust", True))
        sigma_dir_gr = float(lsq_cfg.get("sigma_dir_gr", 0.0003))  # ~3cc
        sigma_dist_m = float(lsq_cfg.get("sigma_dist_m", 0.005))  # 5 mm

        if lsq_enabled and ("chosen" in res):
            # On affine (X,Y,θ) en exploitant toutes les visées sélectionnées
            try:
                dirs, dists = _collect_obs_for_lsq(
                    station_id,
                    obs_data,
                    {s["matricule"]: s for s in obs_data.get("stations", [])},
                    use_dirs=use_dirs,
                    use_dists=use_dists,
                )
                if dirs or dists:
                    X0, Y0 = res["chosen"]
                    Xh, Yh, THh, stats = _free_station_lsq_xytheta(
                        dirs,
                        dists,
                        x0=X0,
                        y0=Y0,
                        robust=robust,
                        sigma_dir_gr=sigma_dir_gr,
                        sigma_dist_m=sigma_dist_m,
                    )
                    payload["lsq_block"] = {
                        "x": Xh,
                        "y": Yh,
                        "theta_gr": THh,
                        "stats": stats,
                    }
                else:
                    PlgLogger.log(
                        tr("No usable observations for LSQ refinement."),
                        log_level=1,
                    )
            except Exception as ex:
                PlgLogger.log(
                    tr("LSQ refinement failed: {e}").format(e=str(ex)),
                    log_level=2,
                    push=True,
                )

            # Export report
            md_report = FreeStation.export_markdown_report(payload)

            out_dir = QgsProject.instance().homePath() + "/rapport/"
            os.makedirs(out_dir, exist_ok=True)
            md_name = "free_station_report.md"
            FileUtils.save_temp_file(md_name, md_report, out_dir)
            md_path = os.path.join(out_dir, md_name)

            try:
                ReportUtils.markdown_to_html(md_path, title=tr("Free Station Report"))
            except Exception as ex:
                PlgLogger.log(
                    message=tr("HTML export failed: {e}").format(e=str(ex)),
                    log_level=2,
                    push=True,
                )
            try:
                ReportUtils.markdown_to_pdf(
                    md_path,
                    title=tr("Free Station Report"),
                    margins_mm=(7, 15, 7, 15),
                )
            except Exception as ex:
                PlgLogger.log(
                    message=tr("PDF export failed: {e}").format(e=str(ex)),
                    log_level=2,
                    push=True,
                )

            PlgLogger.log(
                tr("Free station computed: X={x:.4f}, Y={y:.4f}").format(
                    x=x_result, y=y_result
                ),
                log_level=0,
            )
            return payload

    @staticmethod
    def choose_free_station_candidate(
        payload: Dict[str, Any],
        side: Optional[str] = None,  # "left" | "right" | None
        index: Optional[int] = None,  # 0-based index if preferred
    ) -> Dict[str, Any]:
        """
        Fix the chosen XY candidate in the payload, based on side or index.
        Mutates and returns the payload (adds 'result.chosen', 'result.side', 'result.chosen_index').
        """
        res = payload.get("result") or {}
        cands = res.get("candidates") or []
        if not cands:
            raise ValueError(tr("No candidates available to choose from."))

        # By index (highest priority if provided)
        if index is not None:
            if index < 0 or index >= len(cands):
                raise IndexError(tr("Candidate index out of range."))
            chosen = tuple(cands[index])
            res["chosen"] = chosen
            res["chosen_index"] = int(index)
            res["side"] = None
            PlgLogger.log(tr("Candidate chosen by index: {i}").format(i=index))
            return payload

        # By side (left/right) in AB local frame
        if side:
            side = side.lower()
            if side not in ("left", "right"):
                raise ValueError(tr("Invalid side; expected 'left' or 'right'."))

            A = payload["A"]
            res_meta = res
            xa, ya = float(A["x"]), float(A["y"])
            gab = float(res_meta["gab_gr"])

            best = None
            for i, pt in enumerate(cands):
                xloc, yloc = FreeStation.world_to_local(xa, ya, gab, pt[0], pt[1])
                this_side = "left" if yloc >= 0 else "right"
                if this_side == side:
                    best = (pt, i)
                    break

            if best is None:
                # fallback: take first candidate (rare corner case)
                best = (tuple(cands[0]), 0)

            chosen, idx = best
            res["chosen"] = chosen
            res["chosen_index"] = int(idx)
            res["side"] = side
            PlgLogger.log(
                tr("Candidate chosen by side: {s} (index {i})").format(s=side, i=idx)
            )
            return payload

        # If nothing specified, keep as-is (no choice)
        PlgLogger.log(tr("No side/index provided; no candidate fixed."), log_level=1)
        return payload

    @staticmethod
    def finalize_free_station_with_z_and_report(
        payload: Dict[str, Any],
        obs_data: Dict[str, Any],
        out_basename: str = "free_station_report",
    ) -> Dict[str, Any]:
        """
        After the user has chosen a candidate, compute station Z from verticals (if any),
        attach it to payload['z_block'], and re-export the report (MD + HTML/PDF).
        """
        res = payload.get("result") or {}
        if "chosen" not in res:
            PlgLogger.log(
                tr(
                    "No chosen candidate. Please choose a candidate before computing elevation."
                ),
                log_level=2,
                push=True,
            )
            return payload

        chosen_xy = tuple(res["chosen"])
        stations = {s["matricule"]: s for s in obs_data.get("stations", [])}

        # Read corrections block (optional)
        corr = obs_data.get("corrections") or {}
        use_apparent_level = bool(corr.get("niveau_apparent", False))

        # Compute Z from vertical observations
        # station id to use for filtering 'obs' (fallbacks if missing)
        sid = payload.get("station_id") or (obs_data.get("calcul") or ["M"])[0]

        z_block = compute_free_station_z_from_verticals(
            station_id=sid,
            chosen_xy=chosen_xy,
            stations_by_id=stations,
            obs_list=obs_data.get("obs", []),
            use_apparent_level=use_apparent_level,
            default_hp=0.0,
            weight_mode="inv_dH",
        )
        payload["z_block"] = z_block

        # Re-export report with Z section
        md_report = FreeStation.export_markdown_report(payload)
        out_dir = QgsProject.instance().homePath() + "/rapport/"
        os.makedirs(out_dir, exist_ok=True)
        md_name = out_basename + ".md"
        FileUtils.save_temp_file(md_name, md_report, out_dir)
        md_path = os.path.join(out_dir, md_name)

        # HTML/PDF
        try:
            ReportUtils.markdown_to_html(md_path, title=tr("Free station Report"))
        except Exception as ex:
            PlgLogger.log(
                message=tr("HTML export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )
        try:
            ReportUtils.markdown_to_pdf(
                md_path, title=tr("Free station Report"), margins_mm=(7, 15, 7, 15)
            )
        except Exception as ex:
            PlgLogger.log(
                message=tr("PDF export failed: {e}").format(e=str(ex)),
                log_level=2,
                push=True,
            )

        PlgLogger.log(tr("Report updated with elevation (Z)."))
        return payload

    @staticmethod
    def gr2rad(g: float) -> float:
        """Convert grads to radians"""
        return float(g) * math.pi / 200.0

    @staticmethod
    def rad2gr(r: float) -> float:
        """Convert radians to grads"""
        return float(r) * 200.0 / math.pi

    @staticmethod
    def rotate_local_to_world(
        xa: float, ya: float, gab_gr: float, x_local: float, y_local: float
    ) -> Tuple[float, float]:
        """
        Map local (x along AB, y left of AB) to world XY using AB bearing.
        gab_gr: bearing AB in grads.
        """
        t = FreeStation.gr2rad(gab_gr)
        ct, st = math.cos(t), math.sin(t)
        X = xa + x_local * st + y_local * ct
        Y = ya + x_local * ct - y_local * st
        return X, Y

    @staticmethod
    def world_to_local(
        xa: float, ya: float, gab_gr: float, X: float, Y: float
    ) -> Tuple[float, float]:
        """
        Map world XY back to local (x along AB, y left of AB).
        """
        t = FreeStation.gr2rad(gab_gr)
        ct, st = math.cos(t), math.sin(t)
        dx, dy = (X - xa), (Y - ya)
        xloc = st * dx + ct * dy
        yloc = ct * dx - st * dy
        return xloc, yloc

    @staticmethod
    def free_station_aft_single(
        A: Tuple[float, float],
        B: Tuple[float, float],
        AM: float,
        M_gr: float,
        choose: str = "auto",
        ah_A: Optional[float] = None,
        ah_B: Optional[float] = None,
    ) -> Dict[str, Any]:
        """
        Compute station M from AFT 'free-station' with one pair (A,B).

        Inputs:
          A = (xA, yA), B = (xB, yB)
          AM = distance A->M (m)
          M_gr = angle at M (grads) between MA and MB
          choose: "auto" | "left" | "right"
            - "auto": return both candidates; caller decides later.
            - "left"/"right": pick the solution with y>0 (left) or y<0 (right) w.r.t AB.

        Returns:
          dict with keys:
            "candidates": [ (XM1, YM1), (XM2, YM2) ]
            "chosen": (XM, YM) if choose != "auto"
            "side": "left"/"right"
            "triangles": [ { "A_gr", "B_gr", "M_gr", "AB", "BM", "candidate_pair" }, ... ]
            "gab_gr": bearing AB
            "AB": distance AB
            "AM": distance AM
            "M_gr": angle at M
        """
        xa, ya = A
        xb, yb = B
        AB = math.hypot(xb - xa, yb - ya)

        if AB <= 0.0:
            raise ValueError(tr("A and B must be distinct."))
        if AM <= 0.0:
            raise ValueError(tr("AM must be > 0."))

        M_rad = FreeStation.gr2rad(M_gr)
        s = (AM / AB) * math.sin(M_rad)

        # Numeric guard
        if abs(s) > 1.0:
            s = max(-1.0, min(1.0, s))

        # Two possible angles at B (principal and supplementary)
        B1_rad = math.asin(max(-1.0, min(1.0, s)))
        B2_rad = math.pi - B1_rad

        gab = TopazeCalculator.bearing_gon(xa, ya, xb, yb)
        triangles = []
        uniq = {}

        sinM = math.sin(M_rad)
        if abs(sinM) < 1e-12:
            raise ValueError(tr("Degenerate triangle: M angle too small."))

        for B_rad in (B1_rad,):
            B_gr = FreeStation.rad2gr(B_rad)
            A_gr = 200.0 - (B_gr + float(M_gr))
            A_rad = FreeStation.gr2rad(A_gr)

            # BM from sine law: AB/sin(M) = BM/sin(A) → BM = AB * sin(A) / sin(M)
            BM = AB * math.sin(A_rad) / sinM

            # Intersection of circles in local AB frame
            x_local = (AM * AM - BM * BM + AB * AB) / (2.0 * AB)
            y_sq = AM * AM - x_local * x_local
            if y_sq < 0.0:
                y_sq = 0.0  # numeric guard
            y_local = math.sqrt(y_sq)

            M_left = FreeStation.rotate_local_to_world(xa, ya, gab, x_local, +y_local)
            M_right = FreeStation.rotate_local_to_world(xa, ya, gab, x_local, -y_local)

            triangles.append(
                {
                    "A_gr": A_gr,
                    "B_gr": B_gr,
                    "M_gr": float(M_gr),
                    "AB": AB,
                    "BM": BM,
                    "candidate_pair": [M_left, M_right],
                }
            )

            # Keep unique candidates (round to mm to merge duplicates if any)
            for pt in (M_left, M_right):
                key = (round(pt[0] * 1000), round(pt[1] * 1000))
                uniq[key] = pt

            # Filter valid triangles (angles strictly positive)
            valid_tris = []
            for tri in triangles:
                A_gr = tri["A_gr"]
                B_gr = tri["B_gr"]
                if 0.0 < A_gr < 200.0 and 0.0 < B_gr < 200.0:
                    valid_tris.append(tri)

            # Fallback if nothing passes
            if not valid_tris:
                valid_tris = triangles[:]

            # Choose the most stable triangle (max |y|)
            def y_abs_max_for_triangle(tri):
                (xL, yL), (xR, yR) = tri["candidate_pair"]
                _, yLloc = FreeStation.world_to_local(xa, ya, gab, xL, yL)
                _, yRloc = FreeStation.world_to_local(xa, ya, gab, xR, yR)
                return max(abs(yLloc), abs(yRloc))

            best_tri = max(valid_tris, key=y_abs_max_for_triangle)

            # Build candidates from best triangle only
            candidates = list(best_tri["candidate_pair"])
            triangles = [best_tri]

            # auto selection when 'choose' is "auto"
            auto_choice = None
            if (choose is None) or (choose == "auto"):
                if len(candidates) >= 2 and ah_A is not None and ah_B is not None:
                    auto = _auto_choose_candidate_by_direction_order(
                        A=(xa, ya),
                        B=(xb, yb),
                        ah_A=float(ah_A),
                        ah_B=float(ah_B),
                        candidates=candidates,
                        tol_gr=0.01,
                        amb_margin_gr=0.001,
                        amb_zone_gr=0.002,
                        max_tol_gr=0.10,
                    )
                    if auto is not None:
                        chosen_xy, idx, why = auto
                        _, yloc = FreeStation.world_to_local(
                            xa, ya, gab, chosen_xy[0], chosen_xy[1]
                        )
                        side = "left" if yloc >= 0 else "right"
                        PlgLogger.log(tr("Auto-selected candidate: {w}").format(w=why))
                        auto_choice = {
                            "chosen": chosen_xy,
                            "chosen_index": idx,
                            "side": side,
                        }

        out = {
            "candidates": candidates,
            "triangles": triangles,
            "gab_gr": gab,
            "AB": AB,
            "AM": AM,
            "M_gr": float(M_gr),
        }

        # If we auto-decided, store and skip manual side selection
        if auto_choice:
            out.update(auto_choice)
        # Otherwise allow explicit side selection from caller (as before)
        elif choose in ("left", "right") and candidates:
            t = FreeStation.gr2rad(gab)
            ct, st = math.cos(t), math.sin(t)
            chosen = None
            chosen_side = None
            for pt in candidates:
                dx, dy = (pt[0] - xa), (pt[1] - ya)
                xloc = st * dx + ct * dy
                yloc = ct * dx - st * dy
                side = "left" if yloc >= 0 else "right"
                if side == choose:
                    chosen = pt
                    chosen_side = side
                    break
            if chosen is None and candidates:
                chosen = candidates[0]
                chosen_side = "left" if choose == "left" else "right"
            out["chosen"] = chosen
            out["side"] = chosen_side
            out["chosen_index"] = next(
                (i for i, pt in enumerate(candidates) if pt == chosen), None
            )
        return out

    @staticmethod
    def export_markdown_report(payload: Dict[str, Any]) -> str:
        """
        Export a Markdown report for free-station computation.
        Similar to Resection.export_markdown_report().

        Args:
            payload: Dict containing A, B, result, z_block, etc.

        Returns:
            Markdown string
        """
        L: List[str] = []
        L.append("# " + tr("Free-station (AFT)"))
        L.append("")

        A = payload["A"]
        B = payload["B"]
        res = payload["result"]

        L.append("## " + tr("Inputs"))
        L.append("")
        L.append("- {}: X={}, Y={}".format(tr("Point A"), fmt(A["x"]), fmt(A["y"])))
        L.append("- {}: X={}, Y={}".format(tr("Point A"), fmt(A["x"]), fmt(A["y"])))
        L.append("- {}: X={}, Y={}".format(tr("Point B"), fmt(B["x"]), fmt(B["y"])))
        L.append("- {}: {}".format(tr("Distance AM (m)"), fmt(res["AM"])))
        L.append("- {}: {}".format(tr("Angle at M (gr)"), fmt(res["M_gr"])))
        L.append(
            "- {}: {}  —  {}: {}".format(
                tr("AB (m)"),
                fmt(res["AB"]),
                tr("Bearing AB (gr)"),
                fmt(res["gab_gr"]),
            )
        )
        L.append("")

        rho = res["AM"] / res["AB"] if res.get("AB") else float("nan")

        # Signed distances to AB (local y) for each candidate
        y_list = []
        for pt in res.get("candidates", []):
            _, yloc = FreeStation.world_to_local(
                A["x"], A["y"], res["gab_gr"], pt[0], pt[1]
            )
            y_list.append(yloc)

        y_abs = [abs(v) for v in y_list]
        y_min = min(y_abs) if y_abs else None
        y_min_pct = (
            (100.0 * y_min / res["AB"]) if (y_min is not None and res["AB"]) else None
        )

        L.append("## " + tr("Geometry check"))
        L.append("")
        L.append("- {}: {}".format(tr("AM/AB ratio"), fmt(rho)))
        L.append("- {}: {}".format(tr("M (gr) at station"), fmt(res["M_gr"])))
        if y_list:
            y_str = "; ".join(["{} {}".format(fmt(v), tr("m")) for v in y_abs])
            L.append(
                "- {}: {}".format(tr("Distance to AB (|y|) for candidates"), y_str)
            )
        if y_min_pct is not None:
            flag = ""
            if y_min_pct < 1.0:
                flag = " " + tr("(warning: near-tangent configuration)")
            L.append(
                "- {}: {}%{}".format(
                    tr("Smallest |y| / AB"),
                    fmt(y_min_pct, "{:.2f}"),
                    flag,
                )
            )
        L.append("")

        L.append("## " + tr("Triangle solutions"))
        L.append("")
        L.append(
            "| "
            + tr("A (gr)")
            + " | "
            + tr("B (gr)")
            + " | "
            + tr("M (gr)")
            + " | "
            + tr("BM (m)")
            + " | "
            + tr("Candidates (X,Y)")
            + " |"
        )
        L.append("|---:|---:|---:|---:|---|")
        for tri in res["triangles"]:
            pair = tri["candidate_pair"]
            cand_str = "; ".join([f"({fmt(p[0])}, {fmt(p[1])})" for p in pair])
            L.append(
                f"| {fmt(tri['A_gr'])} | {fmt(tri['B_gr'])} | {fmt(tri['M_gr'])} | {fmt(tri['BM'])} | {cand_str} |"
            )
        L.append("")

        L.append("## " + tr("Results"))

        L.append("")
        if "chosen" in res:
            L.append(
                "- {} ({}): X={}, Y={}".format(
                    tr("Chosen solution"),
                    res.get("side", "?"),
                    fmt(res["chosen"][0]),
                    fmt(res["chosen"][1]),
                )
            )
        else:
            if res["candidates"]:
                for i, pt in enumerate(res["candidates"], 1):
                    L.append(
                        "- {} {}: X={}, Y={}".format(
                            tr("Candidate"),
                            i,
                            fmt(pt[0]),
                            fmt(pt[1]),
                        )
                    )
            else:
                L.append("- " + tr("No valid candidate."))
        L.append("")

        # --- LSQ (XY + orientation) si dispo ---
        lsq = payload.get("lsq_block")
        if lsq:
            S = lsq.get("stats", {})
            L.append("## " + tr("Least Squares (XY, orientation)"))
            L.append("")
            L.append("- X / Y: {} / {}".format(fmt(lsq.get("x")), fmt(lsq.get("y"))))
            L.append(
                "- {}: {}".format(tr("Orientation θ (gr)"), fmt(lsq.get("theta_gr")))
            )
            L.append(
                "- {}: {}  —  {}: {} ".format(
                    tr("Observations"),
                    fmt(S.get("n"), "{:d}"),
                    tr("DOF"),
                    fmt(S.get("dof"), "{:d}"),
                )
            )
            L.append(
                "- σ₀ (unitless): {}  —  RMSE (gr): {}".format(
                    fmt(S.get("sigma0_gr")),
                    fmt(S.get("rmse_gr")),
                )
            )
            L.append(
                "- {} (m): {}  —  {} (m): {}  —  {} (gr): {}".format(
                    tr("EMQ X (1σ)"),
                    fmt(S.get("std_x_m"), "{:.4f}"),
                    tr("EMQ Y (1σ)"),
                    fmt(S.get("std_y_m"), "{:.4f}"),
                    tr("EMQ θ (1σ)"),
                    fmt(S.get("std_theta_gr")),
                )
            )
            L.append("")

            if S.get("residuals_dir"):
                L.append("### " + tr("Per-sight direction residuals"))
                L.append(
                    "| "
                    + tr("To")
                    + " | "
                    + tr("Measured (gr)")
                    + " | "
                    + tr("Computed (gr)")
                    + " | "
                    + tr("Residual (gr)")
                    + " |"
                )
                L.append("|---|---:|---:|---:|")
                for r in S["residuals_dir"]:
                    L.append(
                        f"| {r['to']} | {fmt(r['ah_meas_gr'])} | {fmt(r['ah_comp_gr'])} | {fmt(r['res_gr'])} |"
                    )
                L.append("")

            if S.get("residuals_dist"):
                L.append("### " + tr("Per-sight distance residuals"))
                L.append(
                    "| "
                    + tr("To")
                    + " | "
                    + tr("Measured (m)")
                    + " | "
                    + tr("Computed (m)")
                    + " | "
                    + tr("Residual (m)")
                    + " |"
                )
                L.append("|---|---:|---:|---:|")
                for r in S["residuals_dist"]:
                    L.append(
                        f"| {r['to']} | {fmt(r['d_meas_m'])} | {fmt(r['d_comp_m'])} | {fmt(r['res_m'])} |"
                    )
                L.append("")

        # Elevation (Z) block
        FreeStation._add_z_block(L, payload)

        return "\n".join(L)

    @staticmethod
    def _add_z_block(L: List[str], payload: Dict[str, Any]):
        """Append the 'Elevation (Z)' block for a free-station."""
        L.append("## " + tr("Elevation (Z)"))
        L.append("")

        zb = payload.get("z_block")
        if not zb:
            L.append("- " + tr("No vertical observations available for elevation."))
            L.append("")
            return

        L.append(
            "- {}".format(tr("Station elevation (raw / compensated)"))
            + " {} / {}".format(fmt(zb.get("z_raw")), fmt(zb.get("z_cmp")))
            + tr("m")
        )

        use_apparent_level = payload.get("use_apparent_level", False)
        L.append(
            "- {}: {}".format(
                tr("Apparent level correction"),
                tr("ON") if use_apparent_level else tr("OFF"),
            )
        )

        L.append(
            "- {}: {}  —  {}: {}".format(
                tr("Observations"),
                fmt(zb.get("n"), "{:d}"),
                tr("DOF"),
                fmt(zb.get("dof"), "{:d}"),
            )
        )
        emqx = zb.get("emqx_z")  # dispersion métrique des z_est (m)
        stdz = zb.get("std_z_m")  # EMQ Z (1σ) théorique via covariance LSQ (m)
        rmse = zb.get("rmse_gr")  # RMSE sur résidus angulaires (gr)
        s0 = zb.get("sigma0_gr")  # σ₀ a posteriori (gr)
        L.append(
            "- {} (m): {}  —  {} (m): {}  —  RMSE (gr): {}  —  σ₀ (gr): {}".format(
                tr("EMQ Z (1σ)"),
                fmt(emqx, "{:.4f}"),
                tr("Std Z (1σ, LSQ)"),
                fmt(stdz, "{:.4f}"),
                fmt(rmse),
                fmt(s0),
            )
        )
        L.append("")

        # --- Warning thresholds
        STDZ_WARN = 0.010  # 1 cm
        SIGMA0_WARN_GR = 0.020  # 0.02 g = 2 cc

        # Show warnings only if LSQ exists (dof > 0)
        dof_z = zb.get("dof", 0)
        if dof_z and (stdz is not None or s0 is not None):
            warn_msgs = []
            if (stdz is not None) and (stdz > STDZ_WARN):
                warn_msgs.append(
                    tr("Std Z exceeds {thr} m").format(thr=fmt(STDZ_WARN, "{:.3f}"))
                )
            if (s0 is not None) and (s0 > SIGMA0_WARN_GR):
                warn_msgs.append(
                    tr("σ₀ exceeds {thr} gr").format(thr=fmt(SIGMA0_WARN_GR))
                )
            if warn_msgs:
                # Markdown quote block to make it stand out
                L.append("> " + tr("Warning") + ": " + "; ".join(warn_msgs))
                L.append("")

        # Per-sight table
        per_sight = zb.get("per_sight", [])
        if per_sight:
            L.append("### " + tr("Per-sight verticals"))
            L.append("")
            L.append(
                "| "
                + tr("To")
                + " | "
                + tr("Distance (m)")
                + " | "
                + tr("AV measured (gr)")
                + " | "
                + tr("AV computed (gr)")
                + " | "
                + tr("Residual (gr)")
                + " | "
                + tr("Z estimate (m)")
                + " | "
                + tr("Z residual (m)")
                + " |"
            )
            L.append("|---|---:|---:|---:|---:|---:|---:|")
            for p in per_sight:
                L.append(
                    f"| {p.get('to','')} | {fmt(p.get('dH'))} | {fmt(p.get('av_gr'))} | "
                    f"{fmt(p.get('av_calc_gr', float('nan')))} | "
                    f"{fmt(p.get('res_gr', float('nan')))} | "
                    f"{fmt(p.get('z_est'))} | {fmt(p.get('res_m', 0.0))} |"
                )
            L.append("")


# Test function
def _smoke_test_free_station():
    # Cas synthétique simple
    A = (0.0, 0.0)
    B = (100.0, 0.0)
    AM = 80.0

    M_gr = 5e-3
    try:
        res = FreeStation.free_station_aft_single(A, B, AM, M_gr, choose="left")
        assert False, "Should have raised an error for degenerate triangle"
    except ValueError:
        pass  # Expected

    M_gr = 50.0  # ~45° en rad, triangle non dégénéré
    res = FreeStation.free_station_aft_single(A, B, AM, M_gr, choose="left")
    assert res and "candidates" in res and len(res["candidates"]) >= 1
    xm, ym = res.get("chosen", res["candidates"][0])
    assert -1000.0 < xm < 1100.0 and -1000.0 < ym < 1100.0

    AM = 50.0
    M_gr = 100.0
    res = FreeStation.free_station_aft_single(A, B, AM, M_gr, choose="right")
    assert res and "candidates" in res and len(res["candidates"]) >= 1
    xm, ym = res.get("chosen", res["candidates"][0])
    assert -1000.0 < xm < 1100.0 and -1000.0 < ym < 1100.0

    print(tr("Smoke test passed!"))


if __name__ == "__main__":
    _smoke_test_free_station()
