import math
from math import atan, pi, tan
from typing import Any, Dict, List, Optional, Tuple

import numpy as np
from scipy.optimize import least_squares

from topaze.calc.calc_helpers import *
from topaze.toolbelt import PlgLogger
from topaze.toolbelt.i18n import tr

K_REFRACTION_DEFAULT = 0.13
EARTH_RADIUS_M_DEFAULT = 6371000.0  # meters


def _safe_hp(hp):
    return 0.0 if hp is None else float(hp)


def _elev_rad_from_av_gr(av_gr: float) -> float:
    """Convert zenith angle av (grads) to elevation (radians).
    Convention: av=0g zenith, 100g horizon, 200g nadir.
    """
    return (100.0 - av_gr) * math.pi / 200.0


def _apparent_level_adjustment(
    dh: float, k: float = 0.13, earth_radius: float = 6_370_000.0, verbose: bool = True
) -> float:
    """
    Apparent level correction (curvature + refraction):
    C = 0.5*(1 - k) * dh^2 / R   [meters]
    - dh: horizontal distance (m) (use dh already altered if you apply linear scale elsewhere)
    - k: refraction coefficient (default 0.13)
    - earth_radius: Earth radius in meters
    """
    coeff = 0.5 * (1.0 - float(k))  # 0.5*(1-k) ~ 0.435
    C = coeff * (float(dh) ** 2) / float(earth_radius)
    if verbose:
        print(
            tr("Apparent level correction {c} m (r={r:.3f} m, k={k:.3f})").format(
                c=C, r=dh, k=k
            )
        )
    return C


def _lsq_z_from_verticals(rows, z0, robust=True, loss="soft_l1", f_scale_gr=0.0005):
    """
    LSQ on vertical angles to estimate a single Z parameter.
    rows: list of dict with keys:
      - 'av_gr' (measured, grads)
      - 'r_m'   (horizontal range, m)
      - role-specific geometry:
        * INTERSECTION rows: 'zs','hi','hp','C' with ΔH = (z + hp) - (zs + hi)
        * RESECTION    rows: 'zr','hi','hp','C' with ΔH = (zr + hp) - (z + hi)
      - 'mode': 'intersection' or 'resection'
    z0: initial Z
    Returns: z_hat, stats dict (n, dof, rmse_gr, sigma0_gr, std_z_m, residuals[..])
    """
    g2r = pi / 200.0
    r2g = 200.0 / pi

    def residuals(z):
        res = []
        for rw in rows:
            r = rw["r_m"]
            C = rw.get("C", 0.0)
            if rw["mode"] == "intersection":
                deltaH = (z + rw["hp"]) - (rw["zs"] + rw["hi"])
            else:  # resection
                deltaH = (rw["zr"] + rw["hp"]) - (z + rw["hi"])
            u = (deltaH - C) / r
            av_hat = 100.0 - r2g * atan(u)  # grads
            res.append(rw["av_gr"] - av_hat)  # grads
        return np.array(res, dtype=float)

    # Solve
    if robust:
        res = least_squares(
            residuals,
            x0=np.array([z0], dtype=float),
            loss=loss,
            f_scale=f_scale_gr,
            jac="2-point",
        )
    else:
        res = least_squares(residuals, x0=np.array([z0], dtype=float), jac="2-point")

    z_hat = float(res.x[0])
    v = residuals(z_hat)
    n = v.size
    dof = max(n - 1, 0)
    rmse_gr = float(np.sqrt(np.mean(v**2))) if n else None
    sigma0_gr = None if dof == 0 else float(np.sqrt(np.sum(v**2) / dof))

    # Build Jacobian J (n x 1): dr_i/dz in grads per meter
    # av_hat = 100 - (200/pi)*atan(u), u = (ΔH - C)/r
    # d(av_hat)/dz = -(200/pi) * 1/(1+u^2) * d(u)/dz, with d(u)/dz = ±1/r
    # residual r_i = av_meas - av_hat → dr/dz = - d(av_hat)/dz
    J = []
    for rw in rows:
        r = rw["r_m"]
        C = rw.get("C", 0.0)
        if rw["mode"] == "intersection":
            deltaH = (z_hat + rw["hp"]) - (rw["zs"] + rw["hi"])
            du_dz = 1.0 / r
        else:
            deltaH = (rw["zr"] + rw["hp"]) - (z_hat + rw["hi"])
            du_dz = -1.0 / r
        u = (deltaH - C) / r
        d_avhat_dz = -(200.0 / pi) * (1.0 / (1.0 + u * u)) * du_dz
        dr_dz = -d_avhat_dz
        J.append([dr_dz])
    J = np.array(J, dtype=float)

    std_z = None
    if dof > 0:
        JTJ = float(J.T @ J)  # scalar
        if JTJ > 0 and sigma0_gr is not None:
            var_z = (sigma0_gr**2) / JTJ
            std_z = float(np.sqrt(var_z))

    # Build residual rows (keep what caller provided + residual fields)
    resid_rows = []
    for k, rw in enumerate(rows):
        resid_rows.append(
            {
                **rw,
                "av_comp_gr": rw["av_gr"] - float(v[k]),
                "residual_gr": float(v[k]),
            }
        )

    stats = {
        "n": n,
        "dof": dof,
        "rmse_gr": rmse_gr,
        "sigma0_gr": sigma0_gr,
        "std_z_m": std_z,
        "residuals": resid_rows,
    }
    return z_hat, stats


def compute_target_z_from_verticals_intersection(
    target_id: str,
    x_t: float,
    y_t: float,
    measures: list,
    stations_by_id: dict,
    verbose: bool = True,
    corrections: dict | None = None,
):
    """
    Compute target Z from vertical angles measured at stations with known Z.
    Supports apparent-level correction if corrections['niveau_apparent'] is True.
    """
    # Read corrections
    corr = corrections or {}
    apply_app_level = bool(corr.get("niveau_apparent", False))
    k_ref = float(corr.get("k_refraction", K_REFRACTION_DEFAULT))
    earth_R = float(corr.get("earth_radius_m", EARTH_RADIUS_M_DEFAULT))
    # NOTE: corr.get("k_alteration_lineaire") is intentionally ignored here.

    z_estimates = []
    rows = []

    for m in measures:
        if m.id_target != target_id:
            continue
        av = getattr(m, "av", None)
        if av is None:
            continue

        st = stations_by_id.get(m.id)
        if not st:
            continue
        zs = st["z"] if isinstance(st, dict) else st[2]
        if zs is None:
            continue

        hi = getattr(m, "hi", None)
        if hi is None:
            continue

        xs = st["x"] if isinstance(st, dict) else st[0]
        ys = st["y"] if isinstance(st, dict) else st[1]
        r = float(math.hypot(x_t - xs, y_t - ys))

        e_rad = _elev_rad_from_av_gr(float(av))
        hp = _safe_hp(getattr(m, "hp", None))

        # Apparent level correction term
        C = _apparent_level_adjustment(r, k_ref, earth_R) if apply_app_level else 0.0

        # Height difference ΔH = r*tan(e) (+ C if apparent level)
        z_i = (zs + float(hi)) + r * math.tan(e_rad) + C - hp
        z_estimates.append(z_i)

        rows.append(
            {
                "station": m.id,
                "r_m": r,
                "av_gr": float(av),
                "hi_m": float(hi),
                "hp_m": float(hp),
                "C_m": C,
                "z_est_m": z_i,
            }
        )

        if verbose:
            if apply_app_level:
                print(
                    tr(
                        "Z from {st}: r={r:.3f} m, av={av:.5f} gr, hi={hi:.3f} m, hp={hp:.3f} m, +C={C:.4f} m → Zt={zt:.4f} m"
                    ).format(
                        st=m.id,
                        r=r,
                        av=float(av),
                        hi=float(hi),
                        hp=float(hp),
                        C=C,
                        zt=z_i,
                    )
                )
            else:
                print(
                    tr(
                        "Z from {st}: r={r:.3f} m, av={av:.5f} gr, hi={hi:.3f} m, hp={hp:.3f} m → Zt={zt:.4f} m"
                    ).format(
                        st=m.id, r=r, av=float(av), hi=float(hi), hp=float(hp), zt=z_i
                    )
                )

    n = len(z_estimates)
    if n == 0:
        PlgLogger.log(
            message=tr("No usable vertical angles to compute target elevation."),
            log_level=1,
            push=True,
        )
        return None, None, None

    # 1) RAW: simple mean of direct Z_i (unchanged)
    z_raw = float(sum(z_estimates) / n)

    # 2) LSQ in angle space to get compensated Z
    rows_lsq = []
    for row in rows:  # rows already built per sight
        rows_lsq.append(
            {
                "mode": "intersection",
                "av_gr": row["av_gr"],
                "r_m": row["r_m"],
                "zs": (
                    stations_by_id[row["station"]]["z"]
                    if isinstance(stations_by_id[row["station"]], dict)
                    else stations_by_id[row["station"]][2]
                ),
                "hi": row["hi_m"],
                "hp": row["hp_m"],
                "C": row.get("C_m", 0.0),
                "station": row["station"],
            }
        )

    z_cmp, lsq_stats = _lsq_z_from_verticals(
        rows_lsq, z0=z_raw, robust=True, loss="soft_l1", f_scale_gr=0.0005
    )

    # 3) Metric residuals around compensated Z (keep per-sight diagnostics)
    v_z = [ri["z_est_m"] - z_cmp for ri in rows]
    dof = max(n - 1, 0)
    emq_z = None if dof == 0 else float(np.sqrt(np.sum(np.square(v_z)) / dof))

    # 4) Rebuild per-sight residual table with BOTH angular & metric fields
    res_rows = []
    for row in rows:
        st_rec = stations_by_id[row["station"]]
        zs = st_rec["z"] if isinstance(st_rec, dict) else st_rec[2]
        r = row["r_m"]
        hi = row["hi_m"]
        hp = row["hp_m"]
        C = row.get("C_m", 0.0)

        # From z_cmp, elevation back to AV_hat: tan(ê) = (ΔH - C) / r with ΔH = (z_cmp + hp) - (zs + hi)
        deltaH = (z_cmp + hp) - (zs + hi)
        elev_hat = math.atan((deltaH - C) / r)
        av_hat_gr = 100.0 - (elev_hat * 200.0 / math.pi)
        v_gr = wrap_gr_m200_200(row["av_gr"] - av_hat_gr)

        res_rows.append(
            {
                "station": row["station"],
                "r_m": r,
                "C_m": C,
                "av_meas_gr": row["av_gr"],  # << measured AV for the report
                "av_comp_gr": av_hat_gr,  # << computed AV at Z_comp
                "residual_gr": v_gr,
                "z_est_m": row["z_est_m"],  # << Z estimate from this single sight
                "residual_z_m": (
                    row["z_est_m"] - z_cmp
                ),  # << metric residual vs Z_comp
            }
        )

    # 5) Merge LSQ angular stats and our metric block
    stats = {
        **lsq_stats,
        "z_raw_m": z_raw,
        "z_cmp_m": z_cmp,
        "EMQ_Z_m": emq_z,
        "residuals": res_rows,  # << ensures the report has AV measured & Z estimates
    }
    return z_raw, z_cmp, stats


"""
    # Residuals (vertical) wrt z_cmp
    v_z = [zi - z_cmp for zi in z_estimates]
    dof = max(n - 1, 0)
    emq_z = None if dof == 0 else float(math.sqrt(sum((dz) ** 2 for dz in v_z) / dof))

    # Angular residuals (recompute av̂ from z_cmp: tan(ê) = (ΔH - C)/r)
    res_rows = []
    v_ang = []
    for row in rows:
        st = stations_by_id[row["station"]]
        zs = st["z"] if isinstance(st, dict) else st[2]
        xs = st["x"] if isinstance(st, dict) else st[0]
        ys = st["y"] if isinstance(st, dict) else st[1]
        r = row["r_m"]
        hi = row["hi_m"]
        hp = row["hp_m"]
        C = row["C_m"]

        deltaH = (z_cmp + hp) - (zs + hi)
        elev_hat = math.atan((deltaH - C) / r)
        av_hat_gr = 100.0 - (elev_hat * 200.0 / math.pi)
        v_gr = wrap_gr_m200_200(row["av_gr"] - av_hat_gr)
        v_ang.append(v_gr)

        res_rows.append(
            {
                "station": row["station"],
                "r_m": r,
                "C_m": C,
                "av_meas_gr": row["av_gr"],
                "av_comp_gr": av_hat_gr,
                "residual_gr": v_gr,
                "residual_z_m": (row["z_est_m"] - z_cmp),
            }
        )

    rmse_gr = float(np.sqrt(np.mean(np.square(v_ang)))) if v_ang else None
    s0_gr = (
        None
        if dof == 0 or not v_ang
        else float(np.sqrt(np.sum(np.square(v_ang)) / dof))
    )

    stats = {
        "n": n,
        "dof": dof,
        "z_raw_m": z_raw,
        "z_cmp_m": z_cmp,
        "EMQ_Z_m": emq_z,
        "rmse_gr": rmse_gr,
        "sigma0_gr": s0_gr,
        "apparent_level": {
            "enabled": apply_app_level,
            "k_refraction": k_ref,
            "earth_radius_m": earth_R,
        },
        "residuals": res_rows,
    }
    return z_raw, z_cmp, stats
"""


def compute_station_z_from_verticals_resection(
    station_id: str,
    x_s: float,
    y_s: float,
    measures: list,
    refs_by_id: dict,
    verbose: bool = True,
    corrections: dict | None = None,
):
    """
    Compute station Z from vertical angles to references with known Z.
    Supports apparent-level correction if corrections['niveau_apparent'] is True.
    """
    # Read corrections
    corr = corrections or {}
    apply_app_level = bool(corr.get("niveau_apparent", False))
    k_ref = float(corr.get("k_refraction", K_REFRACTION_DEFAULT))
    earth_R = float(corr.get("earth_radius_m", EARTH_RADIUS_M_DEFAULT))
    # NOTE: corr.get("k_alteration_lineaire") is intentionally ignored here.

    z_estimates = []
    rows = []

    for m in measures:
        if m.id != station_id:
            continue
        av = getattr(m, "av", None)
        if av is None:
            continue

        ref = refs_by_id.get(m.id_target)
        if not ref:
            continue
        zr = ref["z"] if isinstance(ref, dict) else ref[2]
        if zr is None:
            continue

        hi = getattr(m, "hi", None)
        if hi is None:
            continue

        xr = ref["x"] if isinstance(ref, dict) else ref[0]
        yr = ref["y"] if isinstance(ref, dict) else ref[1]
        r = float(math.hypot(xr - x_s, yr - y_s))
        # (facultatif) pareil pour l'altération linéaire :
        # k_lin_ppm = float(corrections.get("k_alteration_lineaire", 0.0))
        # r_eff = r * (1.0 + k_lin_ppm * 1e-6)
        r_eff = r

        e_rad = _elev_rad_from_av_gr(float(av))
        hp = _safe_hp(getattr(m, "hp", None))

        C = (
            _apparent_level_adjustment(r_eff, k_ref, earth_R)
            if apply_app_level
            else 0.0
        )

        # ΔH = r*tan(e) (+ C if apparent level), so:
        # Zs = (Zr + hp) - [r*tan(e) + C] - hi
        z_i = (zr + hp) - (r_eff * math.tan(e_rad) - C) - float(hi)
        z_estimates.append(z_i)

        rows.append(
            {
                "ref": m.id_target,
                "r_m": r,
                "av_gr": float(av),
                "hi_m": float(hi),
                "hp_m": float(hp),
                "C_m": C,
                "z_est_m": z_i,
            }
        )

        if verbose:
            if apply_app_level:
                print(
                    tr(
                        "Z of station from ref {ref}: r={r:.3f} m, av={av:.5f} gr, hi={hi:.3f} m, hp={hp:.3f} m, +C={C:.4f} m → Zs={zs:.4f} m"
                    ).format(
                        ref=m.id_target,
                        r=r,
                        av=float(av),
                        hi=float(hi),
                        hp=float(hp),
                        C=C,
                        zs=z_i,
                    )
                )
            else:
                print(
                    tr(
                        "Z of station from ref {ref}: r={r:.3f} m, av={av:.5f} gr, hi={hi:.3f} m, hp={hp:.3f} m → Zs={zs:.4f} m"
                    ).format(
                        ref=m.id_target,
                        r=r,
                        av=float(av),
                        hi=float(hi),
                        hp=float(hp),
                        zs=z_i,
                    )
                )

    n = len(z_estimates)
    if n == 0:
        PlgLogger.log(
            message=tr("No usable vertical angles to compute station elevation."),
            log_level=1,
            push=True,
        )
        return None, None, None

    # 1) RAW
    z_raw = float(sum(z_estimates) / n)

    # 2) LSQ on angles
    rows_lsq = []
    for row in rows:
        rows_lsq.append(
            {
                "mode": "resection",
                "av_gr": row["av_gr"],
                "r_m": row["r_m"],
                "zr": (
                    refs_by_id[row["ref"]]["z"]
                    if isinstance(refs_by_id[row["ref"]], dict)
                    else refs_by_id[row["ref"]][2]
                ),
                "hi": row["hi_m"],
                "hp": row["hp_m"],
                "C": row.get("C_m", 0.0),
                "ref": row["ref"],
            }
        )

    z_cmp, lsq_stats = _lsq_z_from_verticals(
        rows_lsq, z0=z_raw, robust=True, loss="soft_l1", f_scale_gr=0.0005
    )

    # 3) Metric residuals around compensated Z
    v_z = [ri["z_est_m"] - z_cmp for ri in rows]
    dof = max(n - 1, 0)
    emq_z = None if dof == 0 else float(np.sqrt(np.sum(np.square(v_z)) / dof))

    # 4) Rebuild per-sight residual table
    res_rows = []
    for row in rows:
        ref_rec = refs_by_id[row["ref"]]
        zr = ref_rec["z"] if isinstance(ref_rec, dict) else ref_rec[2]
        r = row["r_m"]
        hi = row["hi_m"]
        hp = row["hp_m"]
        C = row.get("C_m", 0.0)

        # From z_cmp: ΔH = (zr + hp) - (z_cmp + hi) ; tan(ê) = (ΔH - C) / r
        deltaH = (zr + hp) - (z_cmp + hi)
        elev_hat = math.atan((deltaH - C) / r)
        av_hat_gr = 100.0 - (elev_hat * 200.0 / math.pi)
        v_gr = wrap_gr_m200_200(row["av_gr"] - av_hat_gr)

        res_rows.append(
            {
                "ref": row["ref"],
                "r_m": r,
                "C_m": C,
                "av_meas_gr": row["av_gr"],  # measured AV
                "av_comp_gr": av_hat_gr,  # computed AV
                "residual_gr": v_gr,
                "z_est_m": row["z_est_m"],  # per-sight Z estimate
                "residual_z_m": (row["z_est_m"] - z_cmp),
            }
        )

    stats = {
        **lsq_stats,
        "z_raw_m": z_raw,
        "z_cmp_m": z_cmp,
        "EMQ_Z_m": emq_z,
        "residuals": res_rows,  # << contains full fields expected by the report
    }
    return z_raw, z_cmp, stats


"""
    v_z = [zi - z_cmp for zi in z_estimates]
    dof = max(n - 1, 0)
    emq_z = None if dof == 0 else float(math.sqrt(sum((dz) ** 2 for dz in v_z) / dof))

    # Angular residuals (from z_cmp: tan(ê) = (ΔH - C)/r, here ΔH = (Zr + hp) - (z_cmp + hi))
    res_rows = []
    v_ang = []
    for row in rows:
        ref = refs_by_id[row["ref"]]
        xr = ref["x"] if isinstance(ref, dict) else ref[0]
        yr = ref["y"] if isinstance(ref, dict) else ref[1]
        zr = ref["z"] if isinstance(ref, dict) else ref[2]

        r = row["r_m"]
        hi = row["hi_m"]
        hp = row["hp_m"]
        C = row["C_m"]
        deltaH = (zr + hp) - (z_cmp + hi)
        elev_hat = math.atan((deltaH - C) / r)
        av_hat_gr = 100.0 - (elev_hat * 200.0 / math.pi)
        v_gr = wrap_gr_m200_200(row["av_gr"] - av_hat_gr)
        v_ang.append(v_gr)

        res_rows.append(
            {
                "ref": row["ref"],
                "r_m": r,
                "C_m": C,
                "av_meas_gr": row["av_gr"],
                "av_comp_gr": av_hat_gr,
                "residual_gr": v_gr,
                "residual_z_m": (row["z_est_m"] - z_cmp),
            }
        )

    rmse_gr = float(np.sqrt(np.mean(np.square(v_ang)))) if v_ang else None
    s0_gr = (
        None
        if dof == 0 or not v_ang
        else float(np.sqrt(np.sum(np.square(v_ang)) / dof))
    )

    stats = {
        "n": n,
        "dof": dof,
        "z_raw_m": z_raw,
        "z_cmp_m": z_cmp,
        "EMQ_Z_m": emq_z,
        "rmse_gr": rmse_gr,
        "sigma0_gr": s0_gr,
        "apparent_level": {
            "enabled": apply_app_level,
            "k_refraction": k_ref,
            "earth_radius_m": earth_R,
        },
        "residuals": res_rows,
    }
    return z_raw, z_cmp, stats
"""


def compute_free_station_z_from_verticals_by_only_2_sights(
    station_id: str,
    chosen_xy: Tuple[float, float],
    stations_by_id: Dict[str, Dict[str, Any]],
    obs_list: List[Dict[str, Any]],
    use_apparent_level: bool = False,
    default_hp: float = 0.0,
    weight_mode: str = "inv_dH",  # "equal" | "inv_dH"
) -> Optional[Dict[str, Any]]:
    """
    Estimate station elevation Z_M from vertical angles measured at station_id
    towards known-Z targets.

    Inputs:
      - station_id: ID of the station in 'obs' (filters origine == station_id)
      - chosen_xy: (X_M, Y_M) chosen free-station solution.
      - stations_by_id: {matricule: {"x","y","z",...}} (Z must be present for used targets).
      - obs_list: observations [{'origine','cible','av','hi','hp',...}] in grads/meters.
      - use_apparent_level: if True, apply apparent-level correction (curvature+refraction) to ΔH.
      - default_hp: fallback prism height if missing (meters).
      - weight_mode: "equal" or "inv_dH" (weights 1/dH).

    Returns:
      dict with keys:
        {
          "z_raw": float,
          "z_cmp": float,
          "n": int, "dof": int,
          "emqx_z": float,         # 1σ (m)
          "rmse_gr": float, "sigma0_gr": float,
          "per_sight": [
              {
                "to": "ID",
                "dH": float, "av_gr": float,
                "z_est": float, "res_m": float,
                "av_calc_gr": float, "res_gr": float
              }, ...
          ]
        }
      or None if not enough observations.
    """
    XM, YM = chosen_xy
    per = []

    if not station_id:
        PlgLogger.log(
            tr("Missing station_id for elevation computation."),
            log_level=2,
            push=True,
        )
        return None

    # Collect per-sight Z estimates
    for o in obs_list:
        if o.get("origine") != station_id:
            continue
        tid = o.get("cible")
        if not tid or tid not in stations_by_id:
            continue
        target = stations_by_id[tid]
        if "z" not in target or target["z"] is None:
            continue  # need known Z on target

        av = o.get("av")
        if av is None:
            continue  # need vertical angle
        hi = float(o.get("hi") or 0.0)
        hp = float(o.get("hp") if o.get("hp") is not None else default_hp)

        XP, YP = float(target["x"]), float(target["y"])
        ZP = float(target["z"])

        # horizontal distance from chosen M to target P
        dH = math.hypot(XP - XM, YP - YM)
        if dH <= 0.0:
            continue  # degenerate

        e_rad = _elev_rad_from_av_gr(float(av))
        dZ_line = dH * math.tan(
            e_rad
        )  # vertical diff along level from M-axis to target (sign follows tan)

        # Apparent-level correction (curvature+refraction) if requested
        if use_apparent_level:
            # Convention: subtract the correction from the geometric vertical difference
            # so that Δ is referred to the level line (same signe que dZ_line).
            try:
                from topaze.topaze_calculator import TopazeCalculator as TC

                corr = TC.apparent_level_adjustment(dH)
            except Exception:
                # Fallback: 0.42 * d^2 / R (R ~ 6370000 m);
                corr = 0.42 * (dH**2) / 6370000.0
            dZ_line = dZ_line - corr

        # Z_M estimate from this sight:
        # (ZP + hp) - (ZM + hi) = dZ_line  =>  ZM = ZP + hp - hi - dZ_line
        z_est = ZP + hp - hi - dZ_line

        per.append({"to": str(tid), "dH": dH, "av_gr": float(av), "z_est": z_est})

    n = len(per)
    if n == 0:
        PlgLogger.log(
            tr("No usable vertical observations to compute station elevation."),
            log_level=1,
        )
        return None
    if n == 1:
        # With one sight, raw == compensated; residuals = 0
        z_raw = z_cmp = per[0]["z_est"]
        for p in per:
            p.update({"av_calc_gr": per[0]["av_gr"], "res_gr": 0.0, "res_m": 0.0})
        return {
            "z_raw": z_raw,
            "z_cmp": z_cmp,
            "n": 1,
            "dof": 0,
            "emqx_z": 0.0,
            "rmse_gr": 0.0,
            "sigma0_gr": 0.0,
            "per_sight": per,
        }

    # Raw = plain (unweighted) mean
    z_raw = sum(p["z_est"] for p in per) / n

    # Weighted mean (simple model): w_i = 1/dH or = 1
    if weight_mode == "inv_dH":
        weights = [1.0 / p["dH"] for p in per]
    else:
        weights = [1.0] * n
    wsum = sum(weights)
    z_cmp = sum(w * p["z_est"] for w, p in zip(weights, per)) / wsum

    # Residuals (meters) and angle residuals (grads) from compensated Z
    sum_r2 = 0.0
    sum_v2 = 0.0
    for p in per:
        r_m = p["z_est"] - z_cmp  # meters
        p["res_m"] = r_m
        # Compute computed angle from compensated Z:
        # Δ = (ZP + hp) - (ZM + hi)  => e_calc = atan(Δ / dH)  => av_calc = 100g - e_calc
        tid = p["to"]
        t = stations_by_id[tid]
        ZP = float(t["z"])
        hi = 0.0  # will try to read per-observation hi/hp again
        hp = 0.0
        # find the original obs to get its hi/hp again
        for o in obs_list:
            if o.get("origine") == station_id and o.get("cible") == tid:
                hi = float(o.get("hi") or 0.0)
                hp = float(o.get("hp") if o.get("hp") is not None else default_hp)
                break
        dH = p["dH"]
        Delta_calc = (ZP + hp) - (z_cmp + hi)
        e_calc = math.atan2(Delta_calc, dH)
        av_calc_gr = 100.0 - (e_calc * 200.0 / math.pi)

        v_gr = float(p["av_gr"]) - av_calc_gr  # angle residual in grads
        p["av_calc_gr"] = av_calc_gr
        p["res_gr"] = v_gr

        sum_r2 += r_m * r_m
        sum_v2 += v_gr * v_gr

    dof = n - 1  # 1 unknown (Z_M)
    emqx_z = math.sqrt(sum_r2 / dof) if dof > 0 else 0.0
    rmse_gr = math.sqrt(sum_v2 / dof) if dof > 0 else 0.0
    sigma0_gr = rmse_gr  # here identical (single variance component)

    return {
        "z_raw": z_raw,
        "z_cmp": z_cmp,
        "n": n,
        "dof": dof,
        "emqx_z": emqx_z,
        "rmse_gr": rmse_gr,
        "sigma0_gr": sigma0_gr,
        "per_sight": per,
    }


def compute_free_station_z_from_verticals(
    station_id: str,
    chosen_xy: Tuple[float, float],
    stations_by_id: Dict[str, Dict[str, Any]],
    obs_list: List[Dict[str, Any]],
    use_apparent_level: bool = False,
    default_hp: float = 0.0,
    weight_mode: str = "inv_dH",  # "equal" | "inv_dH"  (pour z_raw uniquement)
    lsq: bool = True,  # activer la compensation LSQ (angles)
    robust: bool = True,  # LSQ robuste (soft_l1)
    sigma_av_gr: float = 0.0003,  # écart-type angulaire (gr) ~ 3 cc
) -> Optional[Dict[str, Any]]:
    """
    Estimate station elevation Z_M from vertical angles measured at station_id
    towards known-Z targets (supports N ≥ 2).

    Returns dict with:
      - z_raw, z_cmp, n, dof, emqx_z (m),
      - rmse_gr, sigma0_gr (grads),
      - std_z_m (m)  # EMQ Z (1σ) théorique via covariance LSQ,
      - per_sight: [{to, dH, av_gr, z_est, av_calc_gr, res_gr, res_m}, ...]
    """

    def _wrap_signed_gr(d: float) -> float:
        d = float(d) % 400.0
        return d - 400.0 if d > 200.0 else d

    def _elev_rad_from_av_gr(av_gr: float) -> float:
        # av: zénithal en grads (0g zenith, 100g horizon)
        return (100.0 - float(av_gr)) * math.pi / 200.0

    XM, YM = chosen_xy
    per: List[Dict[str, Any]] = []

    if not station_id:
        PlgLogger.log(
            tr("Missing station_id for elevation computation."),
            log_level=2,
            push=True,
        )
        return None

    # --- Collecte des visées verticales vers des cibles de Z connu
    for o in obs_list:
        if o.get("origine") != station_id:
            continue
        tid = o.get("cible")
        if not tid or tid not in stations_by_id:
            continue
        target = stations_by_id[tid]
        if (
            target.get("z") is None
            or target.get("x") is None
            or target.get("y") is None
        ):
            continue  # besoin de X,Y,Z cible

        av = o.get("av")
        if av is None:
            continue
        hi = float(o.get("hi") or 0.0)
        hp = float(o.get("hp") if o.get("hp") is not None else default_hp)

        XP, YP = float(target["x"]), float(target["y"])
        ZP = float(target["z"])

        dH = math.hypot(XP - XM, YP - YM)
        if dH <= 0.0:
            continue

        e_rad = _elev_rad_from_av_gr(float(av))
        dZ_line = dH * math.tan(e_rad)

        # Correction de niveau apparent (courbure + réfraction), si demandée
        if use_apparent_level:
            try:
                from topaze.topaze_calculator import TopazeCalculator as TC

                C = TC.apparent_level_adjustment(dH)  # >0
            except Exception:
                C = 0.42 * (dH**2) / 6_370_000.0
            dZ_line = dZ_line - C  # référencé à la ligne de niveau

        # Z_M estimé par cette visée: (ZP + hp) - (ZM + hi) = dZ_line -> ZM = ZP + hp - hi - dZ_line
        z_est = ZP + hp - hi - dZ_line

        per.append(
            {
                "to": str(tid),
                "dH": dH,
                "av_gr": float(av),
                "hi": hi,
                "hp": hp,
                "ZP": ZP,
                "z_est": z_est,
            }
        )

    n = len(per)
    if n == 0:
        PlgLogger.log(
            tr("No usable vertical observations to compute station elevation."),
            log_level=1,
        )
        return None
    if n == 1:
        z_raw = z_cmp = per[0]["z_est"]
        per[0].update({"av_calc_gr": per[0]["av_gr"], "res_gr": 0.0, "res_m": 0.0})
        return {
            "z_raw": z_raw,
            "z_cmp": z_cmp,
            "n": 1,
            "dof": 0,
            "emqx_z": 0.0,
            "rmse_gr": 0.0,
            "sigma0_gr": 0.0,
            "std_z_m": 0.0,
            "per_sight": per,
        }

    # --- z_raw : moyenne (pondérée optionnelle)
    if weight_mode == "inv_dH":
        w = [1.0 / p["dH"] for p in per]
    else:
        w = [1.0] * n
    wsum = sum(w)
    z_raw = sum(wi * p["z_est"] for wi, p in zip(w, per)) / wsum

    # --- LSQ sur les angles (si activé et n≥2)
    # Modèle: av_meas ≈ 100g - (200/pi)*atan(u), u = ((ZP+hp) - (Z + hi) - C)/dH
    z_cmp = z_raw
    rmse_gr = 0.0
    sigma0_gr = 0.0
    std_z_m = None

    if lsq and n >= 2:
        # Prépare les arrays pour la fonction de résidus (en grads non scalés)
        dH_arr = np.array([p["dH"] for p in per], dtype=float)
        ZP_arr = np.array([p["ZP"] for p in per], dtype=float)
        hi_arr = np.array([p["hi"] for p in per], dtype=float)
        hp_arr = np.array([p["hp"] for p in per], dtype=float)
        av_arr = np.array([p["av_gr"] for p in per], dtype=float)

        # Recompute C_i pour chaque visée si apparent level actif
        if use_apparent_level:
            try:
                from topaze.topaze_calculator import TopazeCalculator as TC

                C_arr = np.array(
                    [TC.apparent_level_adjustment(d) for d in dH_arr], dtype=float
                )
            except Exception:
                C_arr = 0.42 * (dH_arr**2) / 6_370_000.0
        else:
            C_arr = np.zeros_like(dH_arr)

        g2r = math.pi / 200.0
        r2g = 200.0 / math.pi

        def res_grads(z):
            Z = float(z[0])
            u = ((ZP_arr + hp_arr) - (Z + hi_arr) - C_arr) / dH_arr
            av_hat = 100.0 - r2g * np.arctan(u)  # grads
            v = av_arr - av_hat  # grads
            v = (v + 200.0) % 400.0 - 200.0  # wrap (-200,200]
            # Échelle (optionnelle) par sigma_av_gr → résidus “unitless” pour σ0
            return v / float(sigma_av_gr)

        z0 = np.array([z_raw], dtype=float)
        if robust:
            sol = least_squares(
                res_grads,
                z0,
                method="trf",
                loss="soft_l1",
                f_scale=1.0,
                jac="2-point",
                xtol=1e-12,
                ftol=1e-12,
                gtol=1e-12,
                max_nfev=200,
            )
        else:
            sol = least_squares(
                res_grads,
                z0,
                method="lm",
                jac="2-point",
                xtol=1e-12,
                ftol=1e-12,
                gtol=1e-12,
            )

        z_cmp = float(sol.x[0])

        # Résidus non scalés (en grads) à la solution
        # (reprendre res_grads et “déscaler”)
        def res_grads_unscaled(Z):
            u = ((ZP_arr + hp_arr) - (Z + hi_arr) - C_arr) / dH_arr
            av_hat = 100.0 - r2g * np.arctan(u)
            v = av_arr - av_hat
            return (v + 200.0) % 400.0 - 200.0

        v_gr = res_grads_unscaled(z_cmp)
        dof = max(n - 1, 0)  # 1 inconnue: Z
        rmse_gr = float(np.sqrt(np.mean(v_gr**2)))
        sigma0_gr = None if dof == 0 else float(np.sqrt(np.sum(v_gr**2) / dof))

        # Jacobienne (non scalée) J = ∂(res_grads_unscaled)/∂Z
        # av_hat = 100 - (200/pi)*atan(u) ; u = ((ZP+hp)-(Z+hi)-C)/dH => d u/dZ = -1/dH
        # d(av_hat)/dZ = -(200/pi) * 1/(1+u^2) * (-1/dH) = + (200/pi)/(dH*(1+u^2))
        # r = av_meas - av_hat => dr/dZ = - d(av_hat)/dZ = - (200/pi)/(dH*(1+u^2))
        u = ((ZP_arr + hp_arr) - (z_cmp + hi_arr) - C_arr) / dH_arr
        J = -(r2g) / (dH_arr * (1.0 + u * u))  # (n x 1) en gr/m

        # Covariance sur Z (m^2): Σ_Z = σ0^2 * (J^T J)^-1   (J en gr/m, σ0 en gr)
        JTJ = float(np.sum(J * J))
        if dof > 0 and JTJ > 0.0 and sigma0_gr is not None:
            var_z = (sigma0_gr**2) / JTJ
            std_z_m = float(np.sqrt(var_z))
        else:
            std_z_m = None
    else:
        dof = n - 1
        std_z_m = None

    # --- Remplir le tableau per_sight (angles & métrique)
    sum_r2 = 0.0
    sum_v2 = 0.0
    for p in per:
        r_m = p["z_est"] - z_cmp
        p["res_m"] = r_m

        # AV calculé à partir de z_cmp
        dH = p["dH"]
        ZP = p["ZP"]
        hi = p["hi"]
        hp = p["hp"]
        if use_apparent_level:
            try:
                from topaze.topaze_calculator import TopazeCalculator as TC

                C = TC.apparent_level_adjustment(dH)
            except Exception:
                C = 0.42 * (dH**2) / 6_370_000.0
        else:
            C = 0.0

        Delta_calc = (ZP + hp) - (z_cmp + hi)  # m
        e_calc = math.atan2((Delta_calc - C), dH)
        av_calc_gr = 100.0 - (e_calc * 200.0 / math.pi)
        p["av_calc_gr"] = av_calc_gr

        v_gr = float(p["av_gr"]) - av_calc_gr
        v_gr = _wrap_signed_gr(v_gr)
        p["res_gr"] = v_gr

        sum_r2 += r_m * r_m
        sum_v2 += v_gr * v_gr

    # EMQ (métrique) des estimations autour de z_cmp
    emqx_z = math.sqrt(sum_r2 / max(n - 1, 1))
    if sigma0_gr is None:
        sigma0_gr = math.sqrt(sum_v2 / max(n - 1, 1))

    return {
        "z_raw": z_raw,
        "z_cmp": z_cmp,
        "n": n,
        "dof": max(n - 1, 0),
        "emqx_z": emqx_z,  # dispersion métrique des z_est
        "rmse_gr": rmse_gr,  # RMSE des résidus angulaires (gr)
        "sigma0_gr": sigma0_gr,  # σ₀ a posteriori (gr)
        "std_z_m": std_z_m,  # EMQ Z (1σ) via covariance LSQ (m) si calculable
        "per_sight": per,
    }
