"""
/***************************************************************************
 Resection
                                 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 os
from math import atan, atan2, cos, hypot, pi, sin, tan

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

from topaze.calc.calc_helpers import (
    az_gr_north_cw,
    line_intersection,
    wrap_gr_0_400,
    wrap_gr_m200_200,
)
from topaze.calc.elevation_utils import compute_station_z_from_verticals_resection
from topaze.calc_utils import *
from topaze.file_utils import FileUtils
from topaze.ptopo import Ptopo
from topaze.report_helpers import fmt
from topaze.report_utils import ReportUtils
from topaze.shared import shared_ptopo_array
from topaze.toolbelt import PlgLogger
from topaze.toolbelt.i18n import tr
from topaze.topaze_utils import TopazeUtils
from topaze.topo_sight import TopoSight


class Resection:
    @staticmethod
    def compute_resection():
        x_raw = y_raw = x_lsq = y_lsq = z_raw = z_lsq = None
        data = FileUtils.load_temp_file("relevement.json")
        if not data:
            PlgLogger.log(
                message=tr("No data from relevement.json"), log_level=2, push=True
            )
            out_dir = QgsProject.instance().homePath() + "/rapport/"
            md_name = "resection_report.md"
            FileUtils.remove_temp_file(md_name, out_dir)
            return (None, None, None), (None, None, None), None

        obs_data = json.loads(data)
        corrections = obs_data.get("corrections", {})  # may be empty

        station_id = obs_data["calcul"][0]
        measures = []
        for o in obs_data["obs"]:
            m = None
            if (
                o.get("ah", None) is not None
                or o.get("av", None) is not None
                or o.get("di", None) is not None
            ):
                m = TopoSight(
                    id=o["origine"],
                    id_target=o["cible"],
                    hi=o["hi"],
                    ah=o["ah"],
                    av=o["av"],
                    di=o["di"],
                    hp=o["hp"],
                )
                measures.append(m)
        for m in measures:
            for station in obs_data["stations"]:
                if station["matricule"] == m.id:
                    m.x = station["x"]
                    m.y = station["y"]
                    m.z = station["z"]
                if station["matricule"] == m.id_target:
                    m.x_target = station["x"]
                    m.y_target = station["y"]
                    m.z_target = station["z"]
            """
            print(
                f"obs is id: {m.id} x: {m.x} y: {m.y} hi: {m.hi} "
                f"id_target: {m.id_target} x_target: {m.x_target} "
                f"y_target: {m.y_target} ah: {m.ah} hp: {m.hp}"
            )
            """
        meas = []
        for m in measures:
            if m.ah is not None:
                meas.append(m)
        # Normalize directions and sort once
        for mm in meas:
            mm.ah = (mm.ah % 400.0 + 400.0) % 400.0
        h_mea = sorted(meas, key=lambda t: t.ah)
        try:
            # Raw (unadjusted) coordinates using Tienstra
            x_raw, y_raw = Resection.resection_raw_tienstra(h_mea)
            # print(
            #    tr("Raw resection (Tienstra): ({x}, {y})").format(x=x_raw, y=y_raw)
            # )
        except Exception as ex:
            PlgLogger.log(
                message=tr("Raw resection failed: {err}").format(err=str(ex)),
                log_level=2,
                push=True,
            )
            x_raw = y_raw = None
            return (None, None, None), (None, None, None), None

        # Compensated (barycentric method)
        try:
            x_cmp, y_cmp = Resection.barycentric_resection(obs_data["calcul"][0], h_mea)
            # print(
            #    tr("Compensated resection (barycentric): ({x}, {y})").format(
            #        x=x_cmp, y=y_cmp
            #    )
            # )
        except Exception as ex:
            PlgLogger.log(
                message=tr("Compensated resection failed: {err}").format(err=str(ex)),
                log_level=2,
                push=True,
            )
            return (None, None, None), (None, None, None), None

        # print("raw : ", (x_raw, y_raw), "cmp : ", (x_cmp, y_cmp))

        # LSQ refinement (absolute directions) using compensated as initial guess if available, else raw
        x_start, y_start = (
            (x_cmp, y_cmp)
            if (x_cmp is not None and y_cmp is not None)
            else (x_raw, y_raw)
        )
        # print("start : ", (x_start, y_start))
        lsq_dict = None
        try:
            (
                x_lsq,
                y_lsq,
                theta_lsq,
                resids_lsq,
                stats_lsq,
            ) = Resection.resection_lsq_orientation(h_mea, x_start, y_start)
            # print(
            #    tr(
            #        "LSQ resection (X, Y, orientation): ({x}, {y}), theta={t} gr"
            #    ).format(x=x_lsq, y=y_lsq, t=theta_lsq)
            # )
            lsq_dict = {
                "x": x_lsq,
                "y": y_lsq,
                "theta_gr": theta_lsq,
                "residuals": resids_lsq,
                "stats": stats_lsq,
            }

        except Exception as ex:
            PlgLogger.log(
                message=tr("LSQ resection failed: {err}").format(err=str(ex)),
                log_level=2,
                push=True,
            )
            return (None, None, None), (None, None, None), None

        # Build references dict with Z (known points)
        refs_by_id = {
            st["matricule"]: {"x": st["x"], "y": st["y"], "z": st.get("z")}
            for st in obs_data["stations"]
        }

        v_mea = []
        for m in measures:
            if m.av is not None:
                v_mea.append(m)
        z_raw = z_cmp = None
        z_stats = None
        if x_cmp is not None and y_cmp is not None:
            station_id = obs_data["calcul"][0]  # already have this above in your code
            z_raw, z_cmp, z_stats = compute_station_z_from_verticals_resection(
                station_id,
                x_cmp,
                y_cmp,
                measures,
                refs_by_id,
                verbose=True,
                corrections=corrections,
            )
            if z_cmp is not None:
                print(
                    tr("Station elevation (raw/comp): {zr:.4f} / {zc:.4f} m").format(
                        zr=z_raw, zc=z_cmp
                    )
                )

        # Build Markdown report with LSQ block if available
        md_report = Resection.export_markdown_report(
            station_id,
            measures,
            x_raw,
            y_raw,
            x_cmp,
            y_cmp,
            lsq=lsq_dict,
            z_raw=z_raw,
            z_cmp=z_cmp,
            z_stats=z_stats,
        )

        # Save to temp file for later viewing
        out_dir = QgsProject.instance().homePath() + "/rapport/"
        os.makedirs(out_dir, exist_ok=True)
        md_name = "resection_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("Resection 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("Resection 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,
            )

        lsq_dict["zstats"] = z_stats
        return (x_raw, y_raw, z_raw), (x_lsq, y_lsq, z_cmp), lsq_dict

    @staticmethod
    def _cot(x):
        s = sin(x)
        if abs(s) < 1e-14:
            return None
        return cos(x) / s

    @staticmethod
    def _tienstra_triplet(A, B, C, alpha_g, beta_g, gamma_g):
        """
        Tienstra on one triplet:
          - A,B,C are known points (x,y)
          - alpha,beta,gamma are measured angles at the UNKNOWN P that subtend
            sides (BC), (CA), (AB), in grads, with alpha+beta+gamma = 400 g.
        Uses K = 1/(cot(VertexAngle) - cot(MeasuredAtP)) weights.
        """
        # Vertex angles of triangle ABC (radians)
        Aang = Resection._angle_vertex_rad(B, A, C)
        Bang = Resection._angle_vertex_rad(C, B, A)
        Cang = Resection._angle_vertex_rad(A, C, B)

        # Measured angles at P (radians)
        ar = alpha_g * pi / 200.0
        br = beta_g * pi / 200.0
        gr = gamma_g * pi / 200.0

        cA, cB, cC = Resection._cot(Aang), Resection._cot(Bang), Resection._cot(Cang)
        ca, cb, cg = Resection._cot(ar), Resection._cot(br), Resection._cot(gr)
        if None in (cA, cB, cC, ca, cb, cg):
            raise ZeroDivisionError(
                tr("Degenerate triangle or angle leads to undefined cotangent.")
            )

        dA = cA - ca
        dB = cB - cb
        dC = cC - cg
        if abs(dA) < 1e-14 or abs(dB) < 1e-14 or abs(dC) < 1e-14:
            raise ZeroDivisionError(tr("Near-singular Tienstra denominator."))

        KA = 1.0 / dA
        KB = 1.0 / dB
        KC = 1.0 / dC
        K = KA + KB + KC

        x = (KA * A[0] + KB * B[0] + KC * C[0]) / K
        y = (KA * A[1] + KB * B[1] + KC * C[1]) / K
        return x, y

    @staticmethod
    def barycentric_resection(station_id, measures):
        """
        Compensated (barycentric) resection:
        - Build Tienstra solutions on all consecutive triplets along the 0..400 g circle.
        - Return the simple average of those triplet solutions (robust and unbiased in practice).
        With N=3, this is exactly classical Tienstra. With N>=4, it matches LSQ at the cm/mm level.
        """
        n = len(measures)
        if n < 3:
            raise ValueError(
                tr("At least 3 directions are required for barycentric resection.")
            )

        # Normalize absolute directions and circularly sort by direction (0..400 g)
        for m in measures:
            m.ah = (m.ah % 400.0 + 400.0) % 400.0
        measures = sorted(measures, key=lambda mm: mm.ah)

        # Convenience arrays
        X = [m.x_target for m in measures]
        Y = [m.y_target for m in measures]
        D = [m.ah for m in measures]  # absolute directions in grads (0..400)

        sols = []
        for i in range(n):
            # Define consecutive triplet (A=Ti, B=T{i+1}, C=T{i+2})
            iA = i
            iB = (i + 1) % n
            iC = (i + 2) % n

            A = (X[iA], Y[iA])
            B = (X[iB], Y[iB])
            C = (X[iC], Y[iC])

            # Forward oriented sectors in grads so that alpha+beta+gamma == 400 g
            gamma_g = Resection._wrap_gr_0_400(D[iB] - D[iA])  # P: angle subtending AB
            alpha_g = Resection._wrap_gr_0_400(D[iC] - D[iB])  # P: angle subtending BC
            beta_g = 400.0 - (gamma_g + alpha_g)  # P: angle subtending CA

            if beta_g <= 0.0:
                # Not a valid forward sweep for this triplet; skip
                continue

            try:
                x, y = Resection._tienstra_triplet(A, B, C, alpha_g, beta_g, gamma_g)
                sols.append((x, y))
            except ZeroDivisionError:
                # Skip near-degenerate triplet
                continue

        if not sols:
            raise RuntimeError(
                tr("Barycentric resection failed: no valid triplet solution.")
            )

        # Mean of triplet solutions (we can switch to median if we want extra robustness)
        x_mean = sum(p[0] for p in sols) / len(sols)
        y_mean = sum(p[1] for p in sols) / len(sols)
        return x_mean, y_mean

    @staticmethod
    def _tienstra(A, B, C, alpha, beta, gamma):
        """
        Tienstra formula for angular resection (raw, unadjusted).
        A, B, C are known points as (x, y).
        alpha, beta, gamma are the angles at the unknown point P
        subtending sides (BC), (CA), (AB) respectively (in radians).
        The angles must sum to pi (180°).
        """
        (xA, yA), (xB, yB), (xC, yC) = A, B, C

        # Triangle side lengths opposite to A, B, C
        a = hypot(xB - xC, yB - yC)
        b = hypot(xC - xA, yC - yA)
        c = hypot(xA - xB, yA - yB)

        # Weights per Tienstra
        wA = a / sin(alpha)
        wB = b / sin(beta)
        wC = c / sin(gamma)
        s = wA + wB + wC

        x = (wA * xA + wB * xB + wC * xC) / s
        y = (wA * yA + wB * yB + wC * yC) / s
        return x, y

    @staticmethod
    def _angle_vertex_rad(P, Q, R):
        """
        Angle at vertex Q of triangle P-Q-R, in [0, π].
        """
        v1x, v1y = P[0] - Q[0], P[1] - Q[1]
        v2x, v2y = R[0] - Q[0], R[1] - Q[1]
        cross = v1x * v2y - v1y * v2x
        dot = v1x * v2x + v1y * v2y
        return abs(atan2(cross, dot))  # [0, π]

    @staticmethod
    def _forward_sectors_gr(measures):
        """
        Forward *oriented* sector angles S_i between consecutive directions, in grads.
        Assumes measures are sorted by .ah on the 0..400 circle.
        """
        n = len(measures)
        S = []
        for i in range(n):
            j = (i + 1) % n
            S.append(Resection._wrap_gr_0_400(measures[j].ah - measures[i].ah))
        return S  # sum should be ~400g

    @staticmethod
    def resection_raw_tienstra(measures):
        """
        RAW (unadjusted) position via **correct Tienstra** on consecutive triplets.
        α,β,γ are the *measured* angles at the unknown point (from directions),
        A,B,C are the **triangle vertex angles** at the known stations.
        """
        n = len(measures)
        if n < 3:
            raise ValueError(
                tr("At least 3 directions are required for a raw resection.")
            )

        # Normalize directions and sort by direction value (0..400 g)
        for m in measures:
            m.ah = (m.ah % 400.0 + 400.0) % 400.0
        measures = sorted(measures, key=lambda mm: mm.ah)

        # Build forward sectors; check they wrap to one full turn
        S = Resection._forward_sectors_gr(measures)
        ssum = sum(S)
        if not (390.0 <= ssum <= 410.0):
            PlgLogger.log(
                message=tr("Warning: sector sum is {s} g (expected ~400 g).").format(
                    s=ssum
                ),
                log_level=2,
                push=True,
            )

        # Convenience arrays
        pts = [(m.x_target, m.y_target) for m in measures]

        sols = []
        for i in range(n):
            # Triplet (A=Ti, B=T{i+1}, C=T{i+2}) in the *same* forward sweep
            A = pts[i]
            B = pts[(i + 1) % n]
            C = pts[(i + 2) % n]

            # Measured angles at P between lines to known points:
            # γ = sector Ti -> T{i+1}, α = sector T{i+1} -> T{i+2}, β = 400 - (γ+α)
            gamma_g = S[i]
            alpha_g = S[(i + 1) % n]
            beta_g = 400.0 - (gamma_g + alpha_g)
            if beta_g <= 0.0:
                # Not a valid forward triplet; skip
                continue

            # Triangle vertex angles at A,B,C (from coordinates), in radians
            Aang = Resection._angle_vertex_rad(B, A, C)
            Bang = Resection._angle_vertex_rad(C, B, A)
            Cang = Resection._angle_vertex_rad(A, C, B)

            # Convert measured α,β,γ to radians
            ar = alpha_g * pi / 200.0
            br = beta_g * pi / 200.0
            gr = gamma_g * pi / 200.0

            # cot(x) = cos/sin ; guard against tiny sines
            def cot(x):
                s = sin(x)
                if abs(s) < 1e-14:
                    return None
                return cos(x) / s

            cA, cB, cC = cot(Aang), cot(Bang), cot(Cang)
            ca, cb, cg = cot(ar), cot(br), cot(gr)
            if None in (cA, cB, cC, ca, cb, cg):
                continue

            # Tienstra weights: K = 1 / (cot(VertexAngle) - cot(MeasuredAngleAtP))
            denA = cA - ca
            denB = cB - cb
            denC = cC - cg
            # Guard against near-singular denominators
            if abs(denA) < 1e-14 or abs(denB) < 1e-14 or abs(denC) < 1e-14:
                continue

            KA = 1.0 / denA
            KB = 1.0 / denB
            KC = 1.0 / denC
            Ksum = KA + KB + KC

            x = (KA * A[0] + KB * B[0] + KC * C[0]) / Ksum
            y = (KA * A[1] + KB * B[1] + KC * C[1]) / Ksum
            sols.append((x, y))

        if not sols:
            raise RuntimeError(
                tr(
                    "Raw resection failed: no valid triplet (degenerate geometry or angle coincidences)."
                )
            )

        # Robust aggregation of triplet solutions (median keeps it 'raw' and stable)
        xs = sorted(p[0] for p in sols)
        ys = sorted(p[1] for p in sols)
        k = len(xs) // 2
        if len(xs) % 2:
            x_raw, y_raw = xs[k], ys[k]
        else:
            x_raw = 0.5 * (xs[k - 1] + xs[k])
            y_raw = 0.5 * (ys[k - 1] + ys[k])

        return x_raw, y_raw

    # --- Angle helpers (radians ↔ grads) ---
    @staticmethod
    def _azimuth_gr(dx, dy):
        # Azimuth from dX,dY in grads, wrapped into [0,200)
        # atan2 returns angle from X-axis; we want from +X (east) counterclockwise -> same convention
        ang_rad = atan2(dy, dx)  # [-pi, pi]
        ang_gr = ang_rad * 200.0 / pi  # grads
        ang_gr = ang_gr % 200.0
        return ang_gr

    @staticmethod
    def _rad_to_gr(rad):
        # Convert radians to grads (200g = π rad)
        return rad * 200.0 / pi

    @staticmethod
    def _wrap_gr_0_200(g):
        # Wrap grads into [0,200)
        g = g % 200.0
        return g if g >= 0.0 else g + 200.0

    @staticmethod
    def _wrap_gr_0_400(g):
        return g % 400.0

    @staticmethod
    def _wrap_gr_minus200_200(g):
        g = (g + 200.0) % 400.0 - 200.0
        return 200.0 if abs(g + 200.0) < 1e-12 else g

    @staticmethod
    def _azimuth_gr_north_cw(dx, dy):
        """
        Geodetic surveying convention: 0g = North, increasing clockwise.
        For vector (dx,dy) = (Xt-Xp, Yt-Yp), az = atan2(dx, dy) in radians.
        """
        from math import atan2, pi

        ang_gr = atan2(dx, dy) * 200.0 / pi
        return ang_gr % 400.0

    @staticmethod
    def _wrap_gr_minus100_100(g):
        # Wrap grads into (-100, 100]
        g = (g + 100.0) % 200.0 - 100.0
        # Map -100 to +100 for nicer symmetry (optional)
        return 100.0 if abs(g + 100.0) < 1e-12 else g

    @staticmethod
    def _angle_at_P_radians(xp, yp, xi, yi, xj, yj):
        """
        Smallest positive angle at P between vectors P->Ti and P->Tj, in radians, in [0, π].
        """
        v1x, v1y = xi - xp, yi - yp
        v2x, v2y = xj - xp, yj - yp
        cross = v1x * v2y - v1y * v2x
        dot = v1x * v2x + v1y * v2y
        ang = abs(atan2(cross, dot))  # returns in [-π, π], take absolute -> [0, π]
        return ang

    @staticmethod
    def _sectors_gr(measures):
        """
        Build consecutive sector angles S_i between directions (in grads), each in (0, 200].
        Assumes measures are circularly ordered by increasing 'ah' on [0, 400).
        """
        n = len(measures)
        S = []
        for i in range(n):
            j = (i + 1) % n
            # forward difference on 0..400
            a = Resection._wrap_gr_0_400(measures[j].ah - measures[i].ah)
            # use the smaller arc (<= 200g)
            if a > 200.0:
                a = 400.0 - a
            S.append(a)
        return S  # length n, each in (0,200]

    @staticmethod
    def _measured_sector_angles_gr(measures):
        """
        Build measured sector angles between consecutive targets, in grads (0..200).
        Uses the same logic as your barycentric method (differences of measured angles).
        """
        n = len(measures)
        sectors = []
        for i in range(n):
            j = (i + 1) % n
            a = measures[j].ah - measures[i].ah  # grads
            if a < 0.0:
                a += 200.0
            sectors.append(a)  # already in grads
        # Rotate (-1) to align with your compensated method
        from collections import deque

        d = deque(sectors)
        d.rotate(-1)
        return list(d)

    @staticmethod
    def _computed_sector_angles_gr(measures, xP, yP):
        """
        Compute theoretical sector angles (in grads) at solution P between consecutive targets.
        """
        n = len(measures)
        angles_gr = []
        for i in range(n):
            j = (i + 1) % n
            ang_rad = Resection._angle_at_P_radians(
                xP,
                yP,
                measures[i].x_target,
                measures[i].y_target,
                measures[j].x_target,
                measures[j].y_target,
            )
            angles_gr.append(Resection._rad_to_gr(ang_rad))
        # Same rotation (-1) to keep indexing consistent with measured sectors
        from collections import deque

        d = deque(angles_gr)
        d.rotate(-1)
        return list(d)

    @staticmethod
    def _residuals_by_sector(measures, xP, yP):
        """
        Build residuals per observation sector: measured - computed (in grads, wrapped to (-100,100]).
        Returns a list of dicts with sector labels and residuals.
        Sector i is defined by target Ti -> T{i+1}.
        """
        n = len(measures)
        meas_gr = Resection._measured_sector_angles_gr(measures)
        comp_gr = Resection._computed_sector_angles_gr(measures, xP, yP)
        items = []
        for i in range(n):
            tgt_i = measures[i].id_target
            tgt_j = measures[(i + 1) % n].id_target
            resid = Resection._wrap_gr_minus100_100(meas_gr[i] - comp_gr[i])
            items.append(
                {
                    "sector": f"{tgt_i} → {tgt_j}",
                    "measured_gr": meas_gr[i],
                    "computed_gr": comp_gr[i],
                    "residual_gr": resid,
                }
            )
        return items

    @staticmethod
    def resection_lsq_orientation(measures, x0, y0):
        """
        LSQ on absolute directions: measured_i ≈ theta + Az_NorthCW(P→Ti).
        Unknowns: X, Y, theta (grads on 0..400).
        """
        import numpy as np
        from scipy.optimize import least_squares

        xs = np.array([m.x_target for m in measures], dtype=float)
        ys = np.array([m.y_target for m in measures], dtype=float)
        z = np.array(
            [m.ah for m in measures], dtype=float
        )  # measured directions (0..400)
        n = len(measures)
        if n < 2:
            raise ValueError(
                tr("At least 2 directions are required for LSQ orientation.")
            )

        # Initial theta from first sight
        az0 = Resection._azimuth_gr_north_cw(xs[0] - x0, ys[0] - y0)
        theta0 = Resection._wrap_gr_0_400(z[0] - az0)

        def fun(u):
            x, y, theta = u[0], u[1], u[2]
            az = np.arctan2(xs - x, ys - y) * (200.0 / np.pi)  # North-CW
            az = np.mod(az, 400.0)
            v = z - (theta + az)  # grads
            v = (v + 200.0) % 400.0 - 200.0  # wrap to (-200,200]
            return v

        u0 = np.array([x0, y0, theta0], dtype=float)
        res = least_squares(fun, u0, method="lm")

        x_hat, y_hat, theta_hat = res.x[0], res.x[1], Resection._wrap_gr_0_400(res.x[2])
        v = fun(np.array([x_hat, y_hat, theta_hat], dtype=float))
        dof = max(n - 3, 0)
        rmse = float(np.sqrt(np.mean(v**2)))
        sigma0 = None if dof == 0 else float(np.sqrt(np.sum(v**2) / dof))

        # --- parameter covariance and EMQs (only if DOF>0 and Jacobian OK) ---
        std_x = std_y = std_theta = None
        cov = None
        if sigma0 is not None and hasattr(res, "jac") and res.jac is not None:
            J = (
                res.jac
            )  # shape (n, 3) — derivatives of residuals (gr) w.r.t [X,Y,theta]
            JTJ = J.T @ J
            try:
                JTJ_inv = np.linalg.inv(JTJ)
                cov = (sigma0**2) * JTJ_inv  # units: [X,Y] in m^2, theta in gr^2
                std_x = float(np.sqrt(cov[0, 0]))  # EMQ X (m)
                std_y = float(np.sqrt(cov[1, 1]))  # EMQ Y (m)
                std_theta = float(np.sqrt(cov[2, 2]))  # EMQ theta (gr)
            except np.linalg.LinAlgError:
                pass  # leave std_* as None

        # Per-sight residuals
        residuals = []
        for m in measures:
            # Azimuth (North-clockwise) and wrap to 0..400 g
            az_gr = Resection._azimuth_gr_north_cw(
                m.x_target - x_hat, m.y_target - y_hat
            )
            comp = Resection._wrap_gr_0_400(theta_hat + az_gr)

            # Angular residual (measured - computed), wrapped to (-200,200]
            resi_gr = Resection._wrap_gr_minus200_200(m.ah - comp)

            # Range from solution to target (meters in EPSG:2154)
            dx = m.x_target - x_hat
            dy = m.y_target - y_hat
            R = math.hypot(dx, dy)

            # Signed transverse linear residual at the target:
            # positive residual => clockwise => to the right of the ray
            eps_rad = resi_gr * math.pi / 200.0
            trans_m = R * math.tan(eps_rad)  # exact (for small angles, ~ R*eps_rad)

            # Components of the right-hand normal (East, North)
            az_rad = az_gr * math.pi / 200.0
            nE = math.cos(az_rad)  # right normal unit vector (E component)
            nN = -math.sin(az_rad)  # right normal unit vector (N component)

            dE = trans_m * nE  # signed transverse offset, East
            dN = trans_m * nN  # signed transverse offset, North

            residuals.append(
                {
                    "id_from": m.id,
                    "id_to": m.id_target,
                    "measured_gr": m.ah,
                    "computed_gr": comp,
                    "residual_gr": resi_gr,
                    # metric equivalents:
                    "range_m": R,
                    "transverse_m": trans_m,
                    "transverse_mm": trans_m * 1000.0,
                    "dE_m": dE,
                    "dN_m": dN,
                }
            )

        stats = {
            "rmse_gr": rmse,
            "sigma0_gr": sigma0,
            "n": n,
            "dof": dof,
            "std_x_m": std_x,
            "std_y_m": std_y,
            "std_theta_gr": std_theta,
            "cov": cov.tolist() if cov is not None else None,
        }
        return x_hat, y_hat, theta_hat, residuals, stats

    @staticmethod
    def export_markdown_report(
        station_id,
        measures,
        x_raw,
        y_raw,
        x_cmp,
        y_cmp,
        lsq=None,
        z_raw=None,
        z_cmp=None,
        z_stats=None,
    ):
        """
        Export a Markdown report with raw (Tienstra) and compensated (barycentric) resection results.

                'lsq' optional dict:
          {
            "x": x_hat, "y": y_hat, "theta_gr": theta_hat,
            "residuals": [ {id_from, id_to, measured_gr, computed_gr, residual_gr}, ... ],
            "stats": { "rmse_gr": ..., "sigma0_gr": ..., "n": ..., "dof": ... }
          }

        """
        lines = []
        lines.append("# {}".format(tr("Resection computing")))
        lines.append("")
        lines.append("**{}:** {}".format(tr("Station ID"), station_id))
        lines.append("")
        lines.append("## " + tr("Observations"))
        lines.append("")
        lines.append(
            "| "
            + tr("From")
            + " | "
            + tr("To")
            + " | "
            + tr("X")
            + " | "
            + tr("Y")
            + " | "
            + tr("Horizontal Angle (grads)")
            + " |"
        )
        lines.append("|---|---|---:|---:|---:|")
        for m in measures:
            lines.append(
                f"| {m.id} | {m.id_target} | {m.x_target:.3f} | {m.y_target:.3f} | {m.ah:.4f} |"
            )

        lines.append("")
        lines.append("## " + tr("Results"))
        lines.append("")
        lines.append("| " + tr("Method") + " | X | Y |")
        lines.append("|---|---:|---:|")
        if x_raw is not None and y_raw is not None:
            lines.append(
                "| {} | {:.4f} | {:.4f} |".format(
                    tr("Tienstra"),
                    x_raw,
                    y_raw,
                )
            )
        if x_cmp is not None and y_cmp is not None:
            lines.append(
                "| {} | {:.4f} | {:.4f} |".format(
                    tr("Barycentric"),
                    x_cmp,
                    y_cmp,
                )
            )

        lines.append("")
        lines.append("")
        # LSQ (absolute directions) section
        if lsq is not None:
            stats = lsq.get("stats", {}) if isinstance(lsq, dict) else {}
            n = stats.get("n")
            dof = stats.get("dof")
            rmse = stats.get("rmse_gr")
            s0 = stats.get("sigma0_gr")
            sx = stats.get("std_x_m")
            sy = stats.get("std_y_m")
            st = stats.get("std_theta_gr")  # if present

            lines.append("## " + tr("Least Squares summary"))
            lines.append("")
            lines.append(
                "- {}: {}  —  {}: {}".format(
                    tr("Observations"),
                    fmt(n, "{:d}"),
                    tr("DOF"),
                    fmt(dof, "{:d}"),
                )
            )
            lines.append("- RMSE (gr): {}  —  σ₀ (gr): {}".format(fmt(rmse), fmt(s0)))

            # EMQs per coordinate (show dashes if unavailable)
            lines.append(
                "- {} (m): {}  —  {} (m): {}".format(
                    tr("EMQ X (1σ)"),
                    fmt(sx, "{:.4f}"),
                    tr("EMQ Y (1σ)"),
                    fmt(sy, "{:.4f}"),
                )
            )

            # Optional orientation EMQ if your solver provides it
            if st is not None:
                lines.append("- {} (gr): {}".format(tr("EMQ θ (1σ)"), fmt(st)))

            lines.append("")

            # Per-sight residuals (guard against missing key)
            residuals = stats.get("residuals", [])
            if residuals:
                lines.append("### " + tr("Per-sight residuals"))
                lines.append("")
                lines.append(
                    "| "
                    + tr("From")
                    + " | "
                    + tr("To")
                    + " | "
                    + tr("Measured (gr)")
                    + " | "
                    + tr("Computed (gr)")
                    + " | "
                    + tr("Residual (gr)")
                    + " |"
                )
                lines.append("|---|---|---:|---:|---:|")
                for r in residuals:
                    lines.append(
                        f"| {r.get('id_from','')} | {r.get('id_to','')} | "
                        f"{fmt(r.get('measured_gr'))} | {fmt(r.get('computed_gr'))} | {fmt(r.get('residual_gr'))} |"
                    )
                lines.append("")

                # (Optional) metric equivalents table if we compute them
                if "range_m" in residuals[0] and "transverse_mm" in residuals[0]:
                    lines.append("### " + tr("Per-sight residuals — metric equivalent"))
                    lines.append("")
                    lines.append(
                        "| "
                        + tr("From")
                        + " | "
                        + tr("To")
                        + " | "
                        + tr("Distance (m)")
                        + " | "
                        + tr("Transverse (m)")
                        + " | "
                        + tr("Transverse (mm)")
                        + " |"
                    )
                    lines.append("|---|---|---:|---:|---:|")
                    for r in residuals:
                        lines.append(
                            f"| {r.get('id_from','')} | {r.get('id_to','')} | "
                            f"{fmt(r.get('range_m'), '{:.3f}')} | {fmt(r.get('transverse_m'), '{:.4f}')} | "
                            f"{fmt(r.get('transverse_mm'), '{:.1f}')} |"
                        )
                    lines.append("")

        # Elevation (Z) block for resection
        add_z_block_resection(lines, z_raw, z_cmp, z_stats)

        return "\n".join(lines)


def add_z_block_resection(L, z_raw, z_cmp, zstats):
    """Append the 'Elevation (Z)' block for a resection station."""
    L.append("## " + tr("Elevation (Z)"))
    L.append("")

    if z_raw is None and z_cmp is None:
        L.append(tr("No vertical angles available to compute station elevation."))
        L.append("")
        return

    # Summary line
    L.append(
        "- {}: {} / {} {}".format(
            tr("Station elevation (raw / compensated)"),
            fmt(z_raw, "{:.4f}"),
            fmt(z_cmp, "{:.4f}"),
            tr("m"),
        )
    )
    if zstats:
        n = zstats.get("n")
        dof = zstats.get("dof")
        emq = zstats.get("EMQ_Z_m")
        rm = zstats.get("rmse_gr")
        s0 = zstats.get("sigma0_gr")

        # NEW: apparent level line
        ap = zstats.get("apparent_level", {}) or {}
        enabled = bool(ap.get("enabled", False))
        k_ref = ap.get("k_refraction", None)
        earth_R = ap.get("earth_radius_m", None)
        onoff = tr("ON") if enabled else tr("OFF")
        line = "- {}: {}".format(tr("Apparent level correction"), onoff)
        if k_ref is not None and earth_R is not None:
            line += "- {}: {}  —  {}: {}".format(
                tr("Apparent level correction"),
                onoff,
                tr("k"),
                fmt(k_ref, "{:.3f}"),
                tr("R"),
                fmt(earth_R, "{:.0f}"),
                tr("m"),
            )
        L.append(line)

        L.append(
            "- {}: {}  —  {}: {}".format(
                tr("Observations"),
                fmt(n, "{:d}"),
                tr("DOF"),
                fmt(dof, "{:d}"),
            )
        )
        L.append(
            "- {}: {}  —  {}: {}".format(
                tr("EMQ Z (1σ)"), fmt(emq, "{:.4f}"), tr("RMSE"), fmt(rm)
            )
        )
        L.append("")

        residuals = zstats.get("residuals") or []
        if residuals:
            # NEW: add C (m) column if present
            show_C = any("C_m" in r for r in residuals)
            L.append("### " + tr("Per-sight verticals"))
            L.append("")
            if show_C:
                L.append(
                    "| "
                    + tr("Reference")
                    + " | "
                    + tr("Distance (m)")
                    + " | C (m) | "
                    + tr("AV measured (gr)")
                    + " | "
                    + tr("AV computed (gr)")
                    + " | "
                    + tr("Residual (gr)")
                    + " | "
                    + tr("Z estimate (m)")
                    + " | "
                    + tr("Z residual (m)")
                    + " |"
                )
                L.append("|---|---:|---:|---:|---:|---:|---:|---:|")
            else:
                L.append(
                    "| "
                    + tr("Reference")
                    + " | "
                    + 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 r in residuals:
                if show_C:
                    L.append(
                        f"| {r.get('ref','')} | "
                        f"{fmt(r.get('r_m'), '{:.3f}')} | {fmt(r.get('C_m'), '{:.4f}')} | "
                        f"{fmt(r.get('av_meas_gr'))} | {fmt(r.get('av_comp_gr'))} | "
                        f"{fmt(r.get('residual_gr'))} | "
                        f"{fmt(r.get('z_est_m'), '{:.4f}')} | {fmt(r.get('residual_z_m'), '{:.4f}')} |"
                    )
                else:
                    L.append(
                        f"| {r.get('ref','')} | "
                        f"{fmt(r.get('r_m'), '{:.3f}')} | "
                        f"{fmt(r.get('av_meas_gr'))} | {fmt(r.get('av_comp_gr'))} | "
                        f"{fmt(r.get('residual_gr'))} | "
                        f"{fmt(r.get('z_est_m'), '{:.4f}')} | {fmt(r.get('residual_z_m'), '{:.4f}')} |"
                    )
            L.append("")
