from dataclasses import dataclass
from typing import Sequence, List, Dict
import sympy as sp
import numpy as np
import math
from decimal import Decimal
import matplotlib.pyplot as plt


@dataclass
class Genauigkeitsmaße:
    def __init__(self):
        pass

    @staticmethod
    def s0apost(v, P, r):
        vv = (v.T * P * v)[0, 0]
        vv = float(vv)  # Sympy -> float

        if r <= 0:
            raise ValueError(f"Redundanz r muss > 0 sein, ist {r}.")
        if not math.isfinite(vv):
            raise ValueError(f"vv ist nicht endlich (NaN/Inf). vv={vv}")
        if vv < 0:
            raise ValueError(f"vv ist negativ. vv={vv}")

        return math.sqrt(vv / float(r))



    def helmertscher_punktfehler_3D(self, sigma_x: float, sigma_y: float, sigma_z: float) -> float:
        sx = sp.sympify(sigma_x)
        sy = sp.sympify(sigma_y)
        sz = sp.sympify(sigma_z)

        helmert_pf_3D = sp.sqrt(sx**2 + sy**2 + sz**2)
        return float(helmert_pf_3D)


    def helmertscher_punktfehler_3D_alle(
        self,
        sigma_x_list: Sequence[float],
        sigma_y_list: Sequence[float],
        sigma_z_list: Sequence[float],
    ) -> List[float]:
        if not (
            len(sigma_x_list) == len(sigma_y_list) == len(sigma_z_list)
        ):
            raise ValueError("Listen sigma_x, sigma_y, sigma_z müssen gleich lang sein.")

        return [
            self.helmertscher_punktfehler_3D(sx, sy, sz)
            for sx, sy, sz in zip(sigma_x_list, sigma_y_list, sigma_z_list)
        ]


    def fehlerellipse_standard_2D(
        self,
        Q_xx: float,
        Q_yy: float,
        Q_xy: float,
        sigma_0: float,
    ) -> Dict[str, float]:

        Q_xx_s = sp.sympify(Q_xx)
        Q_yy_s = sp.sympify(Q_yy)
        Q_xy_s = sp.sympify(Q_xy)
        sigma0_s = sp.sympify(sigma_0)

        k = sp.sqrt((Q_xx_s - Q_yy_s)**2 + 4 * Q_xy_s**2)

        Q_dmax = sp.Rational(1, 2) * (Q_xx_s + Q_yy_s + k)
        Q_dmin = sp.Rational(1, 2) * (Q_xx_s + Q_yy_s - k)

        s_max = sigma0_s * sp.sqrt(Q_dmax)
        s_min = sigma0_s * sp.sqrt(Q_dmin)

        theta_rad = sp.Rational(1, 2) * sp.atan2(2 * Q_xy_s, Q_xx_s - Q_yy_s)
        theta_gon = theta_rad * 200 / sp.pi

        return {
            "Q_dmax": float(Q_dmax),
            "Q_dmin": float(Q_dmin),
            "s_max": float(s_max),
            "s_min": float(s_min),
            "theta_rad": float(theta_rad),
            "theta_gon": float(theta_gon),
        }



    def fehlerellipse_konfidenz_2D(
        self,
        s_max: float,
        s_min: float,
        scale_value: float,
    ) -> Dict[str, float]:

        s_max_s = sp.sympify(s_max)
        s_min_s = sp.sympify(s_min)
        scale_s = sp.sympify(scale_value)

        faktor = sp.sqrt(scale_s)

        A_K = faktor * s_max_s
        B_K = faktor * s_min_s

        return {
            "A_K": float(A_K),
            "B_K": float(B_K),
        }



    def plot_ellipsen_alle(
        self,
        x_list: Sequence[float],
        y_list: Sequence[float],
        Q_xx_list: Sequence[float],
        Q_yy_list: Sequence[float],
        Q_xy_list: Sequence[float],
        sigma_0: float,
        scale_value: float | None = None,  # None = Standardellipse, sonst Konfidenzellipse
    ) -> plt.Axes:

        n = len(x_list)
        if not (
            n == len(y_list) == len(Q_xx_list) == len(Q_yy_list) == len(Q_xy_list)
        ):
            raise ValueError("Alle Listen müssen gleich lang sein (ein Eintrag pro Punkt).")

        fig, ax = plt.subplots()

        for i in range(n):
            std = self.fehlerellipse_standard_2D(
                Q_xx_list[i], Q_yy_list[i], Q_xy_list[i], sigma_0
            )
            s_max = std["s_max"]
            s_min = std["s_min"]
            theta = std["theta_rad"]

            # ggf. Konfidenzellipse statt Standardellipse
            if scale_value is not None:
                konf = self.fehlerellipse_konfidenz_2D(s_max, s_min, scale_value)
                A = konf["A_K"]
                B = konf["B_K"]
            else:
                A = s_max
                B = s_min

            t = np.linspace(0, 2 * np.pi, 200)
            ct = np.cos(theta)
            st = np.sin(theta)

            x_local = A * np.cos(t)
            y_local = B * np.sin(t)

            x_ell = x_list[i] + ct * x_local - st * y_local
            y_ell = y_list[i] + st * x_local + ct * y_local

            ax.plot(x_ell, y_ell, linewidth=0.8)
            ax.plot(x_list[i], y_list[i], marker=".", color="k", markersize=3)

        ax.set_aspect("equal", "box")
        ax.grid(True)
        ax.set_xlabel("x")
        ax.set_ylabel("y")
        return ax