#!/usr/bin/env python3
"""Analyze radial profiles of the current TN ring model."""

from __future__ import annotations

import json
import math
import sys
from pathlib import Path

import bpy
import numpy as np


def main() -> None:
    source = Path(sys.argv[sys.argv.index("--") + 1]).resolve()
    output_dir = Path(sys.argv[sys.argv.index("--") + 2]).resolve()
    output_dir.mkdir(parents=True, exist_ok=True)
    bpy.ops.import_scene.gltf(filepath=str(source))
    obj = max((candidate for candidate in bpy.context.scene.objects if candidate.type == "MESH"), key=lambda item: len(item.data.vertices))
    coordinates = np.array([obj.matrix_world @ vertex.co for vertex in obj.data.vertices], dtype=np.float64)
    center_xy = (coordinates[:, :2].min(axis=0) + coordinates[:, :2].max(axis=0)) * 0.5
    local = coordinates.copy()
    local[:, :2] -= center_xy
    radii = np.hypot(local[:, 0], local[:, 1])
    angles = np.arctan2(local[:, 1], local[:, 0])

    bins = 180
    profiles = []
    for index in range(bins):
        minimum = -math.pi + (2 * math.pi * index / bins)
        maximum = -math.pi + (2 * math.pi * (index + 1) / bins)
        mask = (angles >= minimum) & (angles < maximum)
        values_r = radii[mask]
        values_z = local[mask, 2]
        profiles.append({
            "angle_deg": round(math.degrees((minimum + maximum) * 0.5), 3),
            "count": int(mask.sum()),
            "r_min": float(np.quantile(values_r, 0.01)) if len(values_r) else None,
            "r_median": float(np.quantile(values_r, 0.5)) if len(values_r) else None,
            "r_max": float(np.quantile(values_r, 0.99)) if len(values_r) else None,
            "z_min": float(np.quantile(values_z, 0.01)) if len(values_z) else None,
            "z_median": float(np.quantile(values_z, 0.5)) if len(values_z) else None,
            "z_max": float(np.quantile(values_z, 0.99)) if len(values_z) else None,
        })

    stable = [p for p in profiles if p["count"] >= 100]
    result = {
        "source": str(source),
        "vertices": len(obj.data.vertices),
        "center_xy": center_xy.tolist(),
        "bounds": {"min": coordinates.min(axis=0).tolist(), "max": coordinates.max(axis=0).tolist()},
        "global": {
            "r_q01": float(np.quantile(radii, 0.01)),
            "r_q25": float(np.quantile(radii, 0.25)),
            "r_q50": float(np.quantile(radii, 0.50)),
            "r_q75": float(np.quantile(radii, 0.75)),
            "r_q99": float(np.quantile(radii, 0.99)),
            "z_q01": float(np.quantile(local[:, 2], 0.01)),
            "z_q50": float(np.quantile(local[:, 2], 0.50)),
            "z_q99": float(np.quantile(local[:, 2], 0.99)),
        },
        "profiles": profiles,
    }
    output = output_dir / "analysis.json"
    output.write_text(json.dumps(result, indent=2), encoding="utf-8")
    print("TN_ANALYSIS", json.dumps(result["global"]))
    print("TN_CENTER", json.dumps(result["center_xy"]))
    print("TN_PROFILE_EXTREMES", json.dumps({
        "widest": sorted(stable, key=lambda p: (p["r_max"] - p["r_min"]), reverse=True)[:12],
        "tallest": sorted(stable, key=lambda p: (p["z_max"] - p["z_min"]), reverse=True)[:12],
    }))


if __name__ == "__main__":
    main()