from __future__ import annotations

import json
import math
from dataclasses import dataclass
from pathlib import Path
from typing import Any, Dict, List, Optional

import numpy as np

from topaze.report_utils import ReportUtils
from topaze.toolbelt.i18n import tr

# ---------------------------------------------------------------------------
# Dataclass de résultat
# ---------------------------------------------------------------------------


@dataclass
class TrilaterationResult:
    """
    Result of a trilateration computation for a single unknown point.
    """

    point_id: str

    x0: float
    y0: float
    z0: Optional[float]

    x: float
    y: float
    z: Optional[float]

    sigma_x: float
    sigma_y: float
    sigma_z: Optional[float]

    sigma0_xy: float
    sigma0_z: Optional[float]

    cov_xy: float

    n_obs_xy: int
    n_obs_z: int
    dof_xy: int
    dof_z: int

    distance_residuals: List[Dict[str, Any]]
    vertical_residuals: List[Dict[str, Any]]

    meta: Dict[str, Any]
    config: Dict[str, Any]


# ---------------------------------------------------------------------------
# Lecture JSON + calcul
# ---------------------------------------------------------------------------


def _load_trilateration_json(json_path: Path) -> Dict[str, Any]:
    with json_path.open("r", encoding="utf-8") as f:
        return json.load(f)


def _initial_guess_by_circle_intersections(coords_xy: np.ndarray, H_obs: np.ndarray):
    n = len(H_obs)
    if n < 2:
        raise ValueError(
            tr("Initial guess by circle intersections needs at least 2 observations.")
        )

    intersection_points = []
    weights = []

    for i in range(n):
        for j in range(i + 1, n):
            x1, y1 = coords_xy[i]
            x2, y2 = coords_xy[j]
            r1, r2 = H_obs[i], H_obs[j]
            d = np.hypot(x2 - x1, y2 - y1)

            # Vérifier si les cercles s'intersectent
            if d > r1 + r2 or d < abs(r1 - r2):
                continue

            # Calcul des points d'intersection
            a = (r1**2 - r2**2 + d**2) / (2 * d)
            h = np.sqrt(r1**2 - a**2)
            xm = x1 + a * (x2 - x1) / d
            ym = y1 + a * (y2 - y1) / d

            # Deux points d'intersection
            points = [
                (xm + h * (y2 - y1) / d, ym - h * (x2 - x1) / d),
                (xm - h * (y2 - y1) / d, ym + h * (x2 - x1) / d),
            ]

            # Filtrer : garder le point le plus proche des autres stations
            for point in points:
                dx = point[0] - coords_xy[:, 0]
                dy = point[1] - coords_xy[:, 1]
                d_calc = np.hypot(dx, dy)
                # Score = inverse de la distance moyenne aux autres observations
                score = 1.0 / (1.0 + np.mean(np.abs(d_calc - H_obs)))
                intersection_points.append(point)
                weights.append(score)

    if not intersection_points:
        return float(coords_xy[:, 0].mean()), float(coords_xy[:, 1].mean())

    # Moyenne pondérée
    x0 = np.average([p[0] for p in intersection_points], weights=weights)
    y0 = np.average([p[1] for p in intersection_points], weights=weights)

    return x0, y0


def _adjust_xy(
    coords_xy: np.ndarray,
    H_obs: np.ndarray,
    max_iter: int = 20,
    tol: float = 1e-6,
):
    """
    Non-linear LSQ adjustment for XY using horizontal distances.

    Minimise Σ (d_calc - H_obs)² par Gauss-Newton.
    """
    n = H_obs.size
    if n < 2:
        raise ValueError(tr("XY adjustment needs at least 2 horizontal distances."))

    # Initial guess: centroid of known stations
    x0, y0 = _initial_guess_by_circle_intersections(coords_xy, H_obs)
    x = float(x0)
    y = float(y0)

    # Equal weights
    P = np.eye(n, dtype=float)
    N = None

    for _ in range(max_iter):
        dx = x - coords_xy[:, 0]
        dy = y - coords_xy[:, 1]
        d_calc = np.hypot(dx, dy)
        d_calc = np.where(d_calc < 1e-9, 1e-9, d_calc)

        # residuals r = d_calc - H_obs
        r = d_calc - H_obs

        # Jacobian J = ∂d_calc/∂(X,Y)
        J = np.column_stack((dx / d_calc, dy / d_calc))  # shape (n,2)

        N = J.T @ P @ J
        u = J.T @ P @ r.reshape(-1, 1)

        try:
            delta = -np.linalg.solve(N, u)
        except np.linalg.LinAlgError:
            delta = -np.linalg.pinv(N) @ u

        dx_corr = float(delta[0, 0])
        dy_corr = float(delta[1, 0])

        x += dx_corr
        y += dy_corr

        if math.hypot(dx_corr, dy_corr) < tol:
            break

    # Final residuals & precision estimates
    dx = x - coords_xy[:, 0]
    dy = y - coords_xy[:, 1]
    d_calc = np.hypot(dx, dy)
    d_calc = np.where(d_calc < 1e-9, 1e-9, d_calc)

    v = H_obs - d_calc  # obs - computed
    dof = max(1, n - 2)
    sigma0_sq = float((v @ v) / dof)
    sigma0 = math.sqrt(sigma0_sq)

    if N is None:
        raise RuntimeError(tr("Normal matrix not computed in XY adjustment."))

    try:
        Qxx = np.linalg.inv(N)
    except np.linalg.LinAlgError:
        Qxx = np.linalg.pinv(N)

    sigma_x = sigma0 * math.sqrt(float(Qxx[0, 0]))
    sigma_y = sigma0 * math.sqrt(float(Qxx[1, 1]))
    cov_xy = sigma0**2 * float(Qxx[0, 1])

    return x0, y0, x, y, sigma0, sigma_x, sigma_y, cov_xy, v, d_calc


def _adjust_z_mean_and_weighted(
    obs_info: List[Dict[str, Any]],
    config: Dict[str, Any],
) -> Dict[str, Any]:
    """
    Compute both unweighted and weighted Z estimates.

    Returns:
        A dictionary containing unweighted mean Z, weighted Z,
        their standard deviations and residuals.
    """

    z_cfg = config.get("z_adjustment", {})
    if not isinstance(z_cfg, dict):
        z_cfg = {}

    sigma_station_z_default = float(
        z_cfg.get("sigma_station_z", config.get("station_z_stddev", 0.005))
    )
    sigma_hi_default = float(
        z_cfg.get("sigma_hi", config.get("instrument_height_stddev", 0.002))
    )
    sigma_hp_default = float(
        z_cfg.get("sigma_hp", config.get("prism_height_stddev", 0.002))
    )
    sigma_di_const = float(
        z_cfg.get("sigma_di_const", config.get("slope_distance_stddev", 0.003))
    )
    sigma_di_ppm = float(
        z_cfg.get("sigma_di_ppm", config.get("slope_distance_stddev_ppm", 2.0))
    )
    sigma_av_gon_default = float(
        z_cfg.get("sigma_av_gon", config.get("vertical_angle_stddev_gon", 0.0005))
    )

    min_variance = float(z_cfg.get("min_variance", 1e-12))

    z_values: List[float] = []
    variances: List[float] = []

    for info in obs_info:
        di = float(info["di"])
        av = float(info["av"])
        hi = float(info["hi"])
        hp = float(info["hp"])
        station_z = float(info["station_z"])

        inc_rad = (100.0 - av) * math.pi / 200.0

        # Same formula as your current code:
        # Z_unknown = Z_station + hp - hi - dZ_geom
        Zi = station_z + hp - hi - float(info["dZ_geom"])
        z_values.append(Zi)

        sigma_station_z = float(info.get("station_sigma_z", sigma_station_z_default))
        sigma_hi = float(info.get("sigma_hi", sigma_hi_default))
        sigma_hp = float(info.get("sigma_hp", sigma_hp_default))
        sigma_av_gon = float(info.get("sigma_av_gon", sigma_av_gon_default))

        sigma_di = math.sqrt(sigma_di_const**2 + (sigma_di_ppm * di * 1e-6) ** 2)

        # Zi = Z_station + hp - hi - di * sin(inc)
        dZi_ddi = -math.sin(inc_rad)
        dZi_dav_gon = di * math.cos(inc_rad) * math.pi / 200.0

        var_i = (
            sigma_station_z**2
            + sigma_hp**2
            + sigma_hi**2
            + (dZi_ddi * sigma_di) ** 2
            + (dZi_dav_gon * sigma_av_gon) ** 2
        )

        variances.append(max(float(var_i), min_variance))

    Z_arr = np.array(z_values, dtype=float)
    var_arr = np.array(variances, dtype=float)

    n_obs_z = len(Z_arr)
    dof_z = max(1, n_obs_z - 1)

    # ------------------------------------------------------------------
    # 1. Unweighted mean = ordinary least squares with equal weights
    # ------------------------------------------------------------------
    z_mean = float(Z_arr.mean())
    v_mean = Z_arr - z_mean

    sigma0_z_mean = float(math.sqrt(float((v_mean @ v_mean) / dof_z)))
    sigma_z_mean = float(sigma0_z_mean / math.sqrt(n_obs_z))

    vertical_residuals_mean: List[Dict[str, Any]] = []

    for info, Zi, vZi in zip(obs_info, Z_arr, v_mean):
        vertical_residuals_mean.append(
            {
                "obs_index": info["index"],
                "target_id": info["target_id"],
                "measured_vertical_angle": info["av"],
                "di": info["di"],
                "geometric_dZ": info["dZ_geom"],
                "z_from_observation": float(Zi),
                "residual": float(vZi),
                "abs_residual": float(abs(vZi)),
                "weight": 1.0,
                "sigma_z_observation": None,
                "method": "unweighted_mean",
            }
        )

    # ------------------------------------------------------------------
    # 2. Weighted mean = weighted least squares
    # ------------------------------------------------------------------
    w_abs = 1.0 / var_arr

    # Relative weights: easier to interpret in reports.
    weights = w_abs / float(np.mean(w_abs))

    sum_w = float(np.sum(weights))

    z_weighted = float(np.sum(weights * Z_arr) / sum_w)
    v_weighted = Z_arr - z_weighted

    sigma0_z_weighted = float(
        math.sqrt(float(np.sum(weights * v_weighted * v_weighted)) / dof_z)
    )
    sigma_z_weighted = float(sigma0_z_weighted * math.sqrt(1.0 / sum_w))

    vertical_residuals_weighted: List[Dict[str, Any]] = []

    for info, Zi, vZi, wi, var_i in zip(
        obs_info,
        Z_arr,
        v_weighted,
        weights,
        var_arr,
    ):
        q_v_i = max(0.0, 1.0 / float(wi) - 1.0 / sum_w)

        if sigma0_z_weighted > 0.0 and q_v_i > 0.0:
            normalized_residual = float(vZi / (sigma0_z_weighted * math.sqrt(q_v_i)))
        else:
            normalized_residual = None

        vertical_residuals_weighted.append(
            {
                "obs_index": info["index"],
                "target_id": info["target_id"],
                "measured_vertical_angle": info["av"],
                "di": info["di"],
                "geometric_dZ": info["dZ_geom"],
                "z_from_observation": float(Zi),
                "residual": float(vZi),
                "abs_residual": float(abs(vZi)),
                "weight": float(wi),
                "sigma_z_observation": float(math.sqrt(var_i)),
                "normalized_residual": normalized_residual,
                "method": "weighted_mean",
            }
        )

    return {
        "z_mean": z_mean,
        "sigma_z_mean": sigma_z_mean,
        "sigma0_z_mean": sigma0_z_mean,
        "vertical_residuals_mean": vertical_residuals_mean,
        "z_weighted": z_weighted,
        "sigma_z_weighted": sigma_z_weighted,
        "sigma0_z_weighted": sigma0_z_weighted,
        "vertical_residuals_weighted": vertical_residuals_weighted,
        "n_obs_z": n_obs_z,
        "dof_z": dof_z,
    }


def _compute_solution(data: Dict[str, Any]) -> TrilaterationResult:
    """
    Utilise la structure trilateration.json :

      - meta: infos projet
      - config: compute_z, linear_scale_factor, tolerances...
      - calcul: ["S35", ...] -> point inconnu
      - stations: [{matricule, x, y, z}, ...] -> points connus
      - obs: liste d'observations avec :
          * type == "reference"
          * origine == point inconnu
          * cible == station connue
          * di = distance de pente (slope)
          * av = angle vertical en gon (100g = horizontale)
          * hi / hp
    """
    meta = data.get("meta", {})
    config = data.get("config", {})

    # Seul config.compute_z est pris en compte (meta.compute_z ignoré)
    compute_z = bool(config.get("compute_z", True))
    scale = float(config.get("linear_scale_factor", 1.0))

    # Unknown point id
    calcul_list = data.get("calcul") or []
    if not calcul_list:
        raise ValueError(tr("Trilateration JSON file has no `calcul` block."))
    unknown_id = calcul_list[0]

    # Known stations
    stations = data.get("stations", [])
    stations_by_id = {s["matricule"]: s for s in stations}
    if not stations_by_id:
        raise ValueError(tr("Trilateration JSON has no `stations` block."))

    # Observations: from unknown point to known stations, type == "reference"
    raw_obs: List[tuple[int, Dict[str, Any]]] = []
    for idx, obs in enumerate(data.get("obs", []), start=1):
        if obs.get("type") != "reference":
            continue
        if obs.get("origine") != unknown_id:
            continue

        target_id = obs.get("cible")
        if target_id not in stations_by_id:
            raise ValueError(
                tr("Observation #{idx} references unknown station `{pid}`.").format(
                    idx=idx, pid=target_id
                )
            )
        raw_obs.append((idx, obs))

    if len(raw_obs) < 3:
        raise ValueError(
            tr(
                "Trilateration needs at least 3 distance observations from "
                "point `{pid}` to known stations, got {n}."
            ).format(pid=unknown_id, n=len(raw_obs))
        )

    # Build arrays for XY adjustment
    coords_xy: List[tuple[float, float]] = []
    H_obs: List[float] = []
    obs_info: List[Dict[str, Any]] = []

    for idx, obs in raw_obs:
        target_id = obs["cible"]
        st = stations_by_id[target_id]

        xk = float(st["x"])
        yk = float(st["y"])
        zk = float(st.get("z", 0.0))

        # Slope distance + vertical angle
        di = float(obs["di"]) * scale
        av = float(obs["av"])
        hi = float(obs.get("hi", 0.0))
        hp = float(obs.get("hp", 0.0))

        # Convention verticale :
        #   av en gons, 0g = zénith, 100g = horizontale, 200g = nadir
        # On utilise l'angle d'inclinaison i = 100g - av
        inc_gon = 100.0 - av
        inc_rad = inc_gon * math.pi / 200.0

        # Horizontal distance (from slope)
        dH = di * math.cos(inc_rad)
        # Vertical geometric component (axe instrument -> prisme)
        dZ = di * math.sin(inc_rad)

        coords_xy.append((xk, yk))
        H_obs.append(dH)

        obs_info.append(
            {
                "index": idx,
                "target_id": target_id,
                "station_x": xk,
                "station_y": yk,
                "station_z": zk,
                "di": di,
                "av": av,
                "hi": hi,
                "hp": hp,
                "dH_obs": dH,
                "dZ_geom": dZ,
            }
        )

    coords_xy_arr = np.array(coords_xy, dtype=float)
    H_obs_arr = np.array(H_obs, dtype=float)
    n_obs_xy = len(H_obs_arr)

    # XY non-linear LSQ
    X0, Y0, X, Y, sigma0_xy, sigma_x, sigma_y, cov_xy, v_H, H_calc = _adjust_xy(
        coords_xy_arr, H_obs_arr
    )

    # Tolerance for horizontal distance residuals
    dist_res_tol = float(
        config.get(
            "distance_residual_tolerance",
            config.get("distance_back_forward_tolerance", 0.01),
        )
    )

    distance_residuals: List[Dict[str, Any]] = []
    for info, v_i, Hc in zip(obs_info, v_H, H_calc):
        distance_residuals.append(
            {
                "obs_index": info["index"],
                "target_id": info["target_id"],
                "measured_horizontal": info["dH_obs"],
                "computed_horizontal": float(Hc),
                "residual": float(v_i),
                "abs_residual": float(abs(v_i)),
                "within_tolerance": abs(v_i) <= dist_res_tol,
            }
        )

    # Z adjustment (optional)
    z_est: Optional[float] = None
    sigma_z: Optional[float] = None
    sigma0_z: Optional[float] = None
    vertical_residuals: List[Dict[str, Any]] = []
    n_obs_z = 0
    dof_z = 0
    ## if weighted
    z_mean: Optional[float] = None
    z_weighted: Optional[float] = None

    sigma_z_mean: Optional[float] = None
    sigma0_z_mean: Optional[float] = None

    sigma_z_weighted: Optional[float] = None
    sigma0_z_weighted: Optional[float] = None

    vertical_residuals_mean: List[Dict[str, Any]] = []
    vertical_residuals_weighted: List[Dict[str, Any]] = []
    if compute_z:
        z_results = _adjust_z_mean_and_weighted(obs_info, config)

        z_mean = z_results["z_mean"]
        sigma_z_mean = z_results["sigma_z_mean"]
        sigma0_z_mean = z_results["sigma0_z_mean"]
        vertical_residuals_mean = z_results["vertical_residuals_mean"]

        z_weighted = z_results["z_weighted"]
        sigma_z_weighted = z_results["sigma_z_weighted"]
        sigma0_z_weighted = z_results["sigma0_z_weighted"]
        vertical_residuals_weighted = z_results["vertical_residuals_weighted"]

        n_obs_z = z_results["n_obs_z"]
        dof_z = z_results["dof_z"]

        # Valeur principale retenue pour le point calculé.
        # Ici je propose de retenir la valeur pondérée.
        z_est = z_weighted
        sigma_z = sigma_z_weighted
        sigma0_z = sigma0_z_weighted
        vertical_residuals = vertical_residuals_weighted

    dof_xy = max(1, n_obs_xy - 2)

    return TrilaterationResult(
        point_id=unknown_id,
        x0=float(X0),
        y0=float(Y0),
        z0=z_mean,
        x=float(X),
        y=float(Y),
        z=z_weighted,
        sigma_x=float(sigma_x),
        sigma_y=float(sigma_y),
        sigma_z=float(sigma_z) if sigma_z is not None else None,
        sigma0_xy=float(sigma0_xy),
        sigma0_z=float(sigma0_z) if sigma0_z is not None else None,
        cov_xy=float(cov_xy),
        n_obs_xy=n_obs_xy,
        n_obs_z=n_obs_z,
        dof_xy=dof_xy,
        dof_z=dof_z,
        distance_residuals=distance_residuals,
        vertical_residuals=vertical_residuals,
        meta=meta,
        config=config,
    )


# ---------------------------------------------------------------------------
# Génération du rapport Markdown
# ---------------------------------------------------------------------------


def _format_float(value: Optional[float], ndigits: int = 4) -> str:
    if value is None:
        return ""
    return f"{value:.{ndigits}f}"


def _build_markdown_report(data: Dict[str, Any], result: TrilaterationResult) -> str:
    """
    Construit le rapport Markdown, dans l'esprit de _build_markdown_report()
    de closed_traverse.py, mais adapté à une trilatération ponctuelle.
    """
    meta = data.get("meta", {})
    config = data.get("config", {})

    project_name = meta.get("project_name", "")
    observer = meta.get("observer", "")
    instrument = meta.get("instrument", "")
    date = meta.get("date", "")
    angle_unit = meta.get("angle_unit", "gon")
    distance_unit = meta.get("distance_unit", "m")

    # Tolérance utilisée pour interpréter les résidus de distances horizontales
    dist_tol = config.get("distance_residual_tolerance")
    if dist_tol is None:
        dist_tol = config.get("distance_back_forward_tolerance")

    md_parts: List[str] = []

    md_parts.append("# {}".format(tr("Trilateration report")))
    md_parts.append("")
    md_parts.append("{}- **{}**: {}".format("", tr("Project"), project_name))
    md_parts.append("- **{}**: {}".format(tr("Observer"), observer))
    md_parts.append("- **{}**: {}".format(tr("Instrument"), instrument))
    md_parts.append("- **{}**: {}".format(tr("Date"), date))
    md_parts.append("- **{}**: {}".format(tr("Unknown point"), result.point_id))
    md_parts.append(
        "- **{}**: {}".format(
            tr("Number of horizontal distance observations"), result.n_obs_xy
        )
    )
    md_parts.append("")

    # Bloc tolérances : on ne garde que la tolérance sur les résidus de distance
    if dist_tol is not None:
        md_parts.append("### {}".format(tr("Check tolerances")))
        md_parts.append("")
        md_parts.append(
            "- **{}**: ".format(tr("Horizontal distance residual tolerance"))
            + f"{dist_tol} {distance_unit}"
        )
        md_parts.append("")

    # Coordonnées
    md_parts.append("## {}".format(tr("Coordinates of unknown point")))
    md_parts.append("")
    md_parts.append(tr("| Coordinate| Raw (m)| Adjusted (m) | Std. dev. (m) |"))
    md_parts.append("|-------|------------|-----------|-----------|")
    md_parts.append(
        "| X | {} | {} | {} |".format(
            _format_float(result.x0, 4),
            _format_float(result.x, 4),
            _format_float(result.sigma_x, 4),
        )
    )
    md_parts.append(
        "| Y | {} | {} | {} |".format(
            _format_float(result.y0, 4),
            _format_float(result.y, 4),
            _format_float(result.sigma_y, 4),
        )
    )
    if result.z is not None and result.sigma_z is not None:
        md_parts.append(
            "| Z | {} | {} | {} |".format(
                _format_float(result.z0, 4),
                _format_float(result.z, 4),
                _format_float(result.sigma_z, 4),
            )
        )
    md_parts.append("")

    # Résidus horizontaux
    md_parts.append("## {}".format(tr("Horizontal distance residuals")))
    md_parts.append("")
    md_parts.append(
        tr("| Target | Measured Hd (m) | Computed Hd (m) | Residual (m) | Status |")
    )
    md_parts.append("|---|---:|---:|:---:|:---:|")

    for r in result.distance_residuals:
        ok_str = ""
        if r.get("within_tolerance") is np.True_:
            ok_str = "ok"
        elif r.get("within_tolerance") is np.False_:
            ok_str = "**KO**"

        md_parts.append(
            "| {pid} | {mH} | {cH} | {ares} | {ok} |".format(
                pid=r["target_id"],
                mH=_format_float(r["measured_horizontal"], 4),
                cH=_format_float(r["computed_horizontal"], 4),
                ares=_format_float(r["abs_residual"], 4),
                ok=ok_str,
            )
        )

    md_parts.append("")
    md_parts.append(
        "_{}_".format(
            tr(
                "Horizontal distances are derived from raw slope distances and vertical angles."
            )
        )
    )
    md_parts.append("")

    # Partie verticale (si calcul Z activé)
    if result.n_obs_z > 0:
        md_parts.append("## {}".format(tr("Vertical information")))
        md_parts.append("")
        md_parts.append(
            "- **{}**: {}".format(tr("Number of vertical observations"), result.n_obs_z)
        )
        md_parts.append(
            "- **{}**: {}".format(tr("Degrees of freedom (Z)"), result.dof_z)
        )
        md_parts.append("")
        md_parts.append(
            tr(
                "| Target | va ({ang}) | sd ({dist}) | "
                "dZ geom. (m) | Z from obs (m) | Residual (m) |"
            ).format(ang=angle_unit, dist=distance_unit)
        )
        md_parts.append("|---|---:|---:|---:|---:|---:|")

        for r in result.vertical_residuals:
            md_parts.append(
                "| {pid} | {av} | {di} | {dz} | {zobs} | {ares} |".format(
                    pid=r["target_id"],
                    av=_format_float(r["measured_vertical_angle"], 4),
                    di=_format_float(r["di"], 4),
                    dz=_format_float(r["geometric_dZ"], 4),
                    zobs=_format_float(r["z_from_observation"], 4),
                    ares=_format_float(r["abs_residual"], 4),
                )
            )

        md_parts.append("")
        md_parts.append(
            "_{}_".format(
                tr(
                    "Each vertical observation provides an individual estimate of the unknown point elevation."
                )
            )
        )
        md_parts.append(
            "_{}_".format(
                tr(
                    "The final Z is obtained by a least-squares average of these individual estimates."
                )
            )
        )
        md_parts.append("")

    return "\n".join(md_parts)


# ---------------------------------------------------------------------------
# Classe Trilateration avec méthode de classe (même style que ClosedTraverse)
# ---------------------------------------------------------------------------


class Trilateration:
    """
    Helper class for computing an isolated point by trilateration.

    Main entry point:
        Trilateration.compute_trilateration_from_json(json_path, output_dir)
    """

    @classmethod
    def compute_trilateration_from_json(
        cls,
        json_path: Path,
        output_path: Path,
    ) -> Dict[str, Any]:
        """
        Calcule un point isolé par trilatération à partir d'un fichier JSON
        (par ex. 'trilateration.json').

        Étapes:
          - lecture du JSON,
          - conversion distances de pente + angles verticaux → distances horizontales,
          - ajustement non-linéaire XY sur les distances horizontales,
          - calcul optionnel de Z depuis les composantes verticales,
          - calcul des résidus (horizontaux + verticaux),
          - génération d'un rapport Markdown `trilateration.md`,
          - conversion éventuelle en HTML, PDF et ODT si les fonctions
            markdown_to_html / markdown_to_pdf / markdown_to_odt sont disponibles.

        Retourne un dict contenant le résultat et les chemins de rapport.
        """
        output_path.mkdir(parents=True, exist_ok=True)

        data = _load_trilateration_json(json_path)
        result = _compute_solution(data)

        # Rapport Markdown
        md_content = _build_markdown_report(data, result)
        md_path = output_path / "trilateration.md"
        md_path.write_text(md_content, encoding="utf-8")

        # Chemins pour conversions éventuelles
        html_path = output_path / "trilateration.html"
        pdf_path = output_path / "trilateration.pdf"
        odt_path = output_path / "trilateration.odt"

        # Conversions HTML / PDF / ODT (même pattern que closed_traverse)
        ReportUtils.markdown_to_html(md_path, html_path.as_posix())
        ReportUtils.markdown_to_pdf(md_path, pdf_path.as_posix())
        ReportUtils.markdown_to_odt(md_path, odt_path.as_posix())

        return {
            "result": result,
            "report_paths": {
                "md": md_path,
                "html": html_path if html_path.exists() else None,
                "pdf": pdf_path if pdf_path.exists() else None,
                "odt": odt_path if odt_path.exists() else None,
            },
        }


if __name__ == "__main__":
    # Exemple d'appel direct pour tests:
    Trilateration.compute_trilateration_from_json(
        Path("c:/temp/trilateration.json"),
        Path("c:/temp"),
    )
    pass
