import numpy as np
import plotly.graph_objects as go
from scipy.stats import f
import pandas as pd
import Berechnungen


class Genauigkeitsmaße:


    @staticmethod
    def berechne_s0apost(v: np.ndarray, P: np.ndarray, r: int) -> float:
        vTPv_matrix = v.T @ P @ v
        vTPv = float(vTPv_matrix.item())
        s0apost = np.sqrt(vTPv / r)
        return float(s0apost)


    @staticmethod
    def helmert_punktfehler(Qxx, s0_apost, unbekannten_liste, dim=3):
        diagQ = np.diag(Qxx)
        daten = []
        namen_str = [str(sym) for sym in unbekannten_liste]

        punkt_ids = []
        for n in namen_str:
            if n.upper().startswith('X'):
                punkt_ids.append(n[1:])

        for pid in punkt_ids:
            try:
                idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
                idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())

                qx = diagQ[idx_x]
                qy = diagQ[idx_y]
                qz = 0.0

                if dim == 3:
                    try:
                        idx_z = next(i for i, n in enumerate(namen_str) if n.upper() == f"Z{pid}".upper())
                        qz = diagQ[idx_z]
                    except StopIteration:
                        qz = 0.0

                sx = s0_apost * np.sqrt(qx)
                sy = s0_apost * np.sqrt(qy)
                sz = s0_apost * np.sqrt(qz) if dim == 3 else 0

                sP = s0_apost * np.sqrt(qx + qy + qz)

                daten.append([pid, float(sx), float(sy), float(sz), float(sP)])
            except:
                continue

        helmert_punktfehler = pd.DataFrame(daten, columns=["Punkt", "σx", "σy", "σz", f"σP_{dim}D"])
        return helmert_punktfehler



    @staticmethod
    def standardellipse(Qxx, s0_apost, unbekannten_liste):
        Qxx = np.asarray(Qxx, float)
        daten = []
        namen_str = [str(sym) for sym in unbekannten_liste]

        punkt_ids = []
        for n in namen_str:
            if n.upper().startswith('X'):
                punkt_ids.append(n[1:])

        for pid in punkt_ids:
            try:
                idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
                idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())

                qxx = Qxx[idx_x, idx_x]
                qyy = Qxx[idx_y, idx_y]
                qyx = Qxx[idx_y, idx_x]

                # Standardabweichungen
                sx = s0_apost * np.sqrt(qxx)
                sy = s0_apost * np.sqrt(qyy)
                sxy = (s0_apost ** 2) * qyx

                k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))

                # Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
                q_dmax = 0.5 * (qyy + qxx + k)
                q_dmin = 0.5 * (qyy + qxx - k)

                # Halbachsen
                s_max = s0_apost * np.sqrt(q_dmax)
                s_min = s0_apost * np.sqrt(q_dmin)

                # Richtungswinkel theta in gon:
                zaehler = 2 * qyx
                nenner = qxx - qyy
                t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)

                # Quadrantenabfrage
                if nenner > 0 and qyx > 0:  # Qxx - Qyy > 0 und Qyx > 0
                    t_gon = t_grund  # 0 - 50 gon
                elif nenner < 0 and qyx > 0:  # Qxx - Qyy < 0 und Qyx > 0
                    t_gon = 100 - t_grund  # 50 - 100 gon
                elif nenner < 0 and qyx < 0:  # Qxx - Qyy < 0 und Qyx < 0
                    t_gon = 100 + t_grund  # 100 - 150 gon
                elif nenner > 0 and qyx < 0:  # Qxx - Qyy > 0 und Qyx < 0
                    t_gon = 200 - t_grund  # 150 - 200 gon
                else:
                    t_gon = 0.0

                daten.append([
                    pid,
                    float(sx), float(sy), float(sxy),
                    float(s_max), float(s_min),
                    float(t_gon)
                ])

            except:
                continue

            standardellipse = pd.DataFrame(daten, columns=["Punkt", "σx", "σy", "σxy", "s_max", "s_min", "θ [gon]"])
        return standardellipse



    @staticmethod
    def konfidenzellipse(Qxx, s0_apost, unbekannten_liste, R, alpha):
        Qxx = np.asarray(Qxx, float)
        daten = []
        namen_str = [str(sym) for sym in unbekannten_liste]

        punkt_ids = [n[1:] for n in namen_str if n.upper().startswith('X')]

        # Faktor für Konfidenzellipse (F-Verteilung)
        kk = float(np.sqrt(2.0 * f.ppf(1.0 - alpha, 2, R)))

        for pid in punkt_ids:
            try:
                idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
                idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())

                qxx = Qxx[idx_x, idx_x]
                qyy = Qxx[idx_y, idx_y]
                qyx = Qxx[idx_y, idx_x]

                # Standardabweichungen
                sx = s0_apost * np.sqrt(qxx)
                sy = s0_apost * np.sqrt(qyy)
                sxy = (s0_apost ** 2) * qyx

                k = np.sqrt((qxx - qyy) ** 2 + 4 * (qyx ** 2))

                # Q_dmax/min = 0.5 * (Qyy + Qxx +/- k)
                q_dmax = 0.5 * (qyy + qxx + k)
                q_dmin = 0.5 * (qyy + qxx - k)

                # Halbachsen der Standardellipse
                s_max = s0_apost * np.sqrt(q_dmax)
                s_min = s0_apost * np.sqrt(q_dmin)

                # Halbachsen der Konfidenzellipse
                A_K = kk * s_max
                B_K = kk * s_min

                # Richtungswinkel theta in gon:
                zaehler = 2 * qyx
                nenner = qxx - qyy
                t_grund = 0.5 * np.arctan(abs(zaehler) / abs(nenner)) * (200 / np.pi)

                # Quadrantenabfrage
                if nenner > 0 and qyx > 0:
                    t_gon = t_grund  # 0 - 50 gon
                elif nenner < 0 and qyx > 0:
                    t_gon = 100 - t_grund  # 50 - 100 gon
                elif nenner < 0 and qyx < 0:
                    t_gon = 100 + t_grund  # 100 - 150 gon
                elif nenner > 0 and qyx < 0:
                    t_gon = 200 - t_grund  # 150 - 200 gon
                else:
                    t_gon = 0.0

                daten.append([
                    pid,
                    float(sx), float(sy), float(sxy),
                    float(A_K), float(B_K),
                    float(t_gon)
                ])

            except:
                continue

        konfidenzellipse = pd.DataFrame(daten, columns= ["Punkt", "σx", "σy", "σxy", "a_K", "b_K","θ [gon]"])

        return konfidenzellipse



    @staticmethod
    def konfidenzellipsoid(Qxx, s0_apost, unbekannten_liste, R, alpha, skala="f", return_2d_schnitte=True):

        Qxx = np.asarray(Qxx, float)
        namen_str = [str(sym) for sym in unbekannten_liste]
        punkt_ids = sorted({n[1:] for n in namen_str if n and n[0].upper() in ("X", "Y", "Z")})

        # Skalierungsfaktor für Konfidenzbereich
        if skala.lower() == "f":
            k2_3d = f.ppf(1.0 - alpha, df=3)
        elif skala.lower() == "f":
            k2_3d = 3.0 * f.ppf(1.0 - alpha, dfn=3, dfd=R)
        else:
            raise ValueError("skala muss 'chi2' oder 'f' sein.")

        daten = []

        for pid in punkt_ids:
            try:
                idx_x = next(i for i, n in enumerate(namen_str) if n.upper() == f"X{pid}".upper())
                idx_y = next(i for i, n in enumerate(namen_str) if n.upper() == f"Y{pid}".upper())
                idx_z = next(i for i, n in enumerate(namen_str) if n.upper() == f"Z{pid}".upper())
            except StopIteration:
                continue

            # 3x3-Block aus Qxx ziehen
            I = [idx_x, idx_y, idx_z]
            Qp = Qxx[np.ix_(I, I)]

            # Kovarianzmatrix (Sigma) des Punkts
            Sigma = (s0_apost ** 2) * Qp

            # Standardabweichungen
            sx = float(np.sqrt(Sigma[0, 0]))
            sy = float(np.sqrt(Sigma[1, 1]))
            sz = float(np.sqrt(Sigma[2, 2]))

            # Kovarianzen
            sxy = float(Sigma[0, 1])
            sxz = float(Sigma[0, 2])
            syz = float(Sigma[1, 2])

            # Eigenzerlegung (symmetrisch -> eigh)
            evals, evecs = np.linalg.eigh(Sigma)
            order = np.argsort(evals)[::-1]
            evals = evals[order]
            evecs = evecs[:, order]

            # Numerische Sicherheit: negative Mini-Eigenwerte durch Rundung abklemmen
            evals = np.clip(evals, 0.0, None)

            # Halbachsen des Konfidenzellipsoids:
            A, B, C = (np.sqrt(evals * k2_3d)).tolist()

            row = {
                "Punkt": pid,
                "σx": sx, "σy": sy, "σz": sz,
                "σxy": sxy, "σxz": sxz, "σyz": syz,
                "A_K": float(A), "B_K": float(B), "C_K": float(C),
                # Orientierung als Spaltenvektoren (Eigenvektoren)
                "evec_1": evecs[:, 0].tolist(),
                "evec_2": evecs[:, 1].tolist(),
                "evec_3": evecs[:, 2].tolist(),
                "skala_k2": float(k2_3d),
                "skala_typ": skala.lower()
            }

            # Optional: 2D-Schnitte (XY, XZ, YZ) als Ellipsenparameter
            if return_2d_schnitte:
                row.update(Genauigkeitsmaße.ellipsen_schnitt_2d(Sigma, alpha, R, skala))

            daten.append(row)

        return pd.DataFrame(daten)



    @staticmethod
    def ellipsen_schnitt_2d(Sigma3, alpha, R, skala):
        def ellipse_from_2x2(S2):
            # Skalierung für 2D
            if skala.lower() == "f":
                k2 = f.ppf(1.0 - alpha, df=2)
            else:
                k2 = 2.0 * f.ppf(1.0 - alpha, dfn=2, dfd=R)

            evals, evecs = np.linalg.eigh(S2)
            order = np.argsort(evals)[::-1]
            evals = np.clip(evals[order], 0.0, None)
            evecs = evecs[:, order]

            a, b = np.sqrt(evals * k2)

            # Winkel der Hauptachse (zu a) in der Ebene: atan2(vy, vx)
            vx, vy = evecs[0, 0], evecs[1, 0]
            theta_rad = np.arctan2(vy, vx)
            theta_gon = float(theta_rad * (200.0 / np.pi)) % 200.0

            return float(a), float(b), theta_gon

        # Submatrizen
        S_xy = Sigma3[np.ix_([0, 1], [0, 1])]
        S_xz = Sigma3[np.ix_([0, 2], [0, 2])]
        S_yz = Sigma3[np.ix_([1, 2], [1, 2])]

        axy, bxy, txy = ellipse_from_2x2(S_xy)
        axz, bxz, txz = ellipse_from_2x2(S_xz)
        ayz, byz, tyz = ellipse_from_2x2(S_yz)

        return {
            "aXY": axy, "bXY": bxy, "θXY [gon]": txy,
            "aXZ": axz, "bXZ": bxz, "θXZ [gon]": txz,
            "aYZ": ayz, "bYZ": byz, "θYZ [gon]": tyz,
        }



    @staticmethod
    def transform_q_with_your_functions(q_xyz, B, L):
        # East
        r11 = Berechnungen.E(L, 1, 0)
        r12 = Berechnungen.E(L, 0, 1)
        r13 = 0

        # North
        r21 = Berechnungen.N(B, L, 1, 0, 0)
        r22 = Berechnungen.N(B, L, 0, 1, 0)
        r23 = Berechnungen.N(B, L, 0, 0, 1)

        # Up
        r31 = Berechnungen.U(B, L, 1, 0, 0)
        r32 = Berechnungen.U(B, L, 0, 1, 0)
        r33 = Berechnungen.U(B, L, 0, 0, 1)

        R = np.array([
            [r11, r12, r13],
            [r21, r22, r23],
            [r31, r32, r33]
        ])

        q_enu = R @ q_xyz @ R.T
        return q_enu











def plot_netz_komplett_final(x_vektor, unbekannten_labels, beobachtungs_labels, Qxx, sigma0_apost,
                             k_faktor=2.447, v_faktor=1000):
    """
    Optimierter Plot für Jupyter Notebook:
    - k_faktor: Statistischer Sicherheitsfaktor (2.447 entspricht 95% für 2D)
    - v_faktor: Optische Überhöhung der Ellipsen (z.B. 1000 = mm werden als m dargestellt)
    """

    x_vektor = np.asarray(x_vektor, float).reshape(-1)
    Qxx = np.asarray(Qxx, float)

    # 1. Datenaufbereitung
    coords = {}
    punkt_ids = sorted(set(str(l)[1:] for l in unbekannten_labels if str(l).startswith(('X', 'Y', 'Z'))))
    pts_data = []

    for pid in punkt_ids:
        try:
            ix = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"X{pid}")
            iy = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"Y{pid}")
            x, y = float(x_vektor[ix]), float(x_vektor[iy])
            coords[pid] = (x, y)

            # Kovarianzmatrix extrahieren und mit s0^2 skalieren
            q_idx = [ix, iy]
            Q_sub = Qxx[np.ix_(q_idx, q_idx)] * (sigma0_apost ** 2)
            pts_data.append({'id': pid, 'x': x, 'y': y, 'Q': Q_sub})
        except:
            continue

    if len(pts_data) == 0:
        raise ValueError(
            "Keine Netzpunkte extrahiert. Prüfe: x_vektor Form (u,) und Qxx Form (u,u) sowie Labels 'X<id>'/'Y<id>'.")

    fig = go.Figure()

    # 2. Beobachtungen (Gruppiert)
    beob_typen = {
        'GNSS-Basislinien': {'pattern': 'gnss', 'color': 'rgba(255, 100, 0, 0.4)'},
        'Nivellement': {'pattern': 'niv', 'color': 'rgba(0, 200, 100, 0.4)'},
        'Tachymeter': {'pattern': '', 'color': 'rgba(100, 100, 100, 0.3)'}
    }

    for typ, info in beob_typen.items():
        x_l, y_l = [], []
        for bl in beobachtungs_labels:
            bl_str = str(bl).lower()
            if (info['pattern'] in bl_str and info['pattern'] != '') or (
                    info['pattern'] == '' and 'gnss' not in bl_str and 'niv' not in bl_str):
                pts = [pid for pid in coords if f"_{pid}" in str(bl) or str(bl).startswith(f"{pid}_")]
                if len(pts) >= 2:
                    x_l.extend([coords[pts[0]][0], coords[pts[1]][0], None])
                    y_l.extend([coords[pts[0]][1], coords[pts[1]][1], None])

        if x_l:
            fig.add_trace(go.Scatter(x=x_l, y=y_l, mode='lines', name=typ, line=dict(color=info['color'], width=1)))

    # 3. Konfidenzellipsen mit v_faktor
    for pt in pts_data:
        vals, vecs = np.linalg.eigh(pt['Q'])
        order = vals.argsort()[::-1]
        vals, vecs = vals[order], vecs[:, order]

        theta = np.degrees(np.arctan2(vecs[1, 0], vecs[0, 0]))
        # Skalierung: k_faktor (Statistik) * v_faktor (Optik)
        a = k_faktor * np.sqrt(vals[0]) * v_faktor
        b = k_faktor * np.sqrt(vals[1]) * v_faktor

        t = np.linspace(0, 2 * np.pi, 40)
        e_x = a * np.cos(t)
        e_y = b * np.sin(t)
        R = np.array([[np.cos(np.radians(theta)), -np.sin(np.radians(theta))],
                      [np.sin(np.radians(theta)), np.cos(np.radians(theta))]])
        rot = np.dot(R, np.array([e_x, e_y]))

        fig.add_trace(go.Scatter(
            x=rot[0, :] + pt['x'], y=rot[1, :] + pt['y'],
            mode='lines', line=dict(color='red', width=1.5),
            name=f"Ellipsen (Vergrößert {v_faktor}x)",
            legendgroup="Ellipsen",
            showlegend=(pt == pts_data[0]),  # Nur einmal in der Legende zeigen
            hoverinfo='skip'
        ))

    # 4. Punkte
    df_pts = pd.DataFrame(pts_data)
    fig.add_trace(go.Scatter(
        x=df_pts['x'], y=df_pts['y'], mode='markers+text',
        text=df_pts['id'], textposition="top center",
        marker=dict(size=8, color='black'), name="Netzpunkte"
    ))

    # 5. Layout & Notebook-Größe
    fig.update_layout(
        title=f"Netzausgleichung: Ellipsen {v_faktor}-fach vergrößert (k={k_faktor})",
        xaxis=dict(title="X [m]", tickformat="f", separatethousands=True, scaleanchor="y", scaleratio=1, showgrid=True,
                   gridcolor='lightgrey'),
        yaxis=dict(title="Y [m]", tickformat="f", separatethousands=True, showgrid=True, gridcolor='lightgrey'),
        width=1100,  # Breite angepasst
        height=900,  # Höhe deutlich vergrößert für Jupiter Notebook
        plot_bgcolor='white',
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01, bgcolor="rgba(255,255,255,0.8)")
    )

    # Info-Annotation als Ersatz für einen physischen Maßstabstab
    fig.add_annotation(
        text=f"<b>Maßstab Ellipsen:</b><br>Dargestellte Größe = Wahre Ellipse × {v_faktor}",
        align='left', showarrow=False, xref='paper', yref='paper', x=0.02, y=0.05,
        bgcolor="white", bordercolor="black", borderwidth=1)

    fig.show(config={'scrollZoom': True})



def plot_netz_final_mit_df_ellipsen(x_vektor, unbekannten_labels, beobachtungs_labels, df_ellipsen, v_faktor=1000):
    # 1. Punkte extrahieren
    coords = {}
    # Wir nehmen an, dass die Reihenfolge im x_vektor X, Y, Z pro Punkt ist
    punkt_ids = sorted(set(str(l)[1:] for l in unbekannten_labels if str(l).startswith(('X', 'Y', 'Z'))))

    for pid in punkt_ids:
        try:
            ix = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"X{pid}")
            iy = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"Y{pid}")
            coords[pid] = (float(x_vektor[ix]), float(x_vektor[iy]))
        except:
            continue

    fig = go.Figure()

    # 2. Beobachtungslinien (Gruppiert)
    beob_typen = {
        'GNSS-Basislinien': {'pattern': 'gnss', 'color': 'rgba(255, 100, 0, 0.4)'},
        'Nivellement': {'pattern': 'niv', 'color': 'rgba(0, 200, 100, 0.4)'},
        'Tachymeter': {'pattern': '', 'color': 'rgba(100, 100, 100, 0.3)'}
    }

    for typ, info in beob_typen.items():
        x_l, y_l = [], []
        for bl in beobachtungs_labels:
            bl_str = str(bl).lower()
            # Einfache Logik zur Typtrennung
            if (info['pattern'] in bl_str and info['pattern'] != '') or \
                    (info['pattern'] == '' and 'gnss' not in bl_str and 'niv' not in bl_str):
                pts = [pid for pid in coords if f"_{pid}" in str(bl) or str(bl).startswith(f"{pid}_")]
                if len(pts) >= 2:
                    x_l.extend([coords[pts[0]][0], coords[pts[1]][0], None])
                    y_l.extend([coords[pts[0]][1], coords[pts[1]][1], None])

        if x_l:
            fig.add_trace(go.Scatter(x=x_l, y=y_l, mode='lines', name=typ, line=dict(color=info['color'], width=1)))

    # 3. Ellipsen aus dem DataFrame zeichnen
    for _, row in df_ellipsen.iterrows():
        pid = str(row['Punkt'])
        if pid in coords:
            x0, y0 = coords[pid]

            # Werte aus DF (mit v_faktor skalieren)
            a = row['a_K'] * v_faktor
            b = row['b_K'] * v_faktor
            theta_gon = row['θ [gon]']

            # Umrechnung: gon -> rad für die Rotation
            # Da im Plot X horizontal und Y vertikal ist, entspricht theta_gon dem Winkel zur X-Achse
            theta_rad = theta_gon * (np.pi / 200.0)

            # Ellipsen berechnen
            t = np.linspace(0, 2 * np.pi, 50)
            e_x = a * np.cos(t)
            e_y = b * np.sin(t)

            # Ausrichtung der Ellipsen
            R = np.array([[np.cos(theta_rad), -np.sin(theta_rad)],
                          [np.sin(theta_rad), np.cos(theta_rad)]])

            rot = np.dot(R, np.array([e_x, e_y]))

            fig.add_trace(go.Scatter(
                x=rot[0, :] + x0, y=rot[1, :] + y0,
                mode='lines', line=dict(color='red', width=1.5),
                name='Konfidenzellipsen',
                legendgroup='Ellipsen',
                showlegend=(pid == df_ellipsen.iloc[0]['Punkt']),
                hoverinfo='text',
                text=f"Punkt {pid}<br>a_K: {row['a_K']:.4f}m<br>b_K: {row['b_K']:.4f}m"
            ))

    # Punkte plotten
    df_pts = pd.DataFrame([(pid, c[0], c[1]) for pid, c in coords.items()], columns=['ID', 'X', 'Y'])
    fig.add_trace(go.Scatter(
        x=df_pts['X'], y=df_pts['Y'], mode='markers+text',
        text=df_pts['ID'], textposition="top center",
        marker=dict(size=8, color='black'), name="Netzpunkte"))

    # Layout
    fig.update_layout(
        title=f"Netzplot (Ellipsen {v_faktor}x überhöht)",
        xaxis=dict(title="X [m]", tickformat="f", separatethousands=True, scaleanchor="y", scaleratio=1,
                   showgrid=True, gridcolor='lightgrey'),
        yaxis=dict(title="Y [m]", tickformat="f", separatethousands=True, showgrid=True, gridcolor='lightgrey'),
        width=1100, height=900,
        plot_bgcolor='white')

    # Maßstabsangabe
    fig.add_annotation(
        text=f"<b>Skalierung:</b><br>Ellipsengröße im Plot = {v_faktor} × Realität",
        align='left', showarrow=False, xref='paper', yref='paper', x=0.02, y=0.02,
        bgcolor="rgba(255,255,255,0.8)", bordercolor="black", borderwidth=1)

    fig.show(config={'scrollZoom': True})

    import plotly.graph_objects as go
    import numpy as np



def plot_netz_3D(x_vektor, unbekannten_labels, beobachtungs_labels, df_ellipsen, v_faktor=1000):
    """
    Erzeugt einen interaktiven 3D-Plot des Netzes.
    - v_faktor: Vergrößerung der Genauigkeits-Achsen (z.B. 1000 für mm -> m)
    """
    # 1. Punkte extrahieren
    pts = {}
    punkt_ids = sorted(set(str(l)[1:] for l in unbekannten_labels if str(l).startswith(('X', 'Y', 'Z'))))

    for pid in punkt_ids:
        try:
            ix = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"X{pid}")
            iy = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"Y{pid}")
            iz = next(i for i, s in enumerate(unbekannten_labels) if str(s) == f"Z{pid}")
            pts[pid] = (float(x_vektor[ix]), float(x_vektor[iy]), float(x_vektor[iz]))
        except:
            continue

    fig = go.Figure()

    # 2. Beobachtungen (Linien im Raum)
    # Wir zeichnen hier einfach alle Verbindungen
    x_line, y_line, z_line = [], [], []
    for bl in beobachtungs_labels:
        p_in_l = [pid for pid in pts if f"_{pid}" in str(bl) or str(bl).startswith(f"{pid}_")]
        if len(p_in_l) >= 2:
            p1, p2 = pts[p_in_l[0]], pts[p_in_l[1]]
            x_line.extend([p1[0], p2[0], None])
            y_line.extend([p1[1], p2[1], None])
            z_line.extend([p1[2], p2[2], None])

    fig.add_trace(go.Scatter3d(
        x=x_line, y=y_line, z=z_line,
        mode='lines', line=dict(color='gray', width=2),
        name='Beobachtungen'
    ))

    # 3. Punkte & "Fehler-Kreuze" (als Ersatz für Ellipsoide)
    # Ein echtes 3D-Ellipsoid ist grafisch schwer, daher zeichnen wir 3 Achsen
    for pid, coord in pts.items():
        # Hier könnten wir die echten Halbachsen aus der 3D-Eigenwertanalyse nutzen
        # Für den Anfang plotten wir die Standardabweichungen sX, sY, sZ als Kreuz
        fig.add_trace(go.Scatter3d(
            x=[coord[0]], y=[coord[1]], z=[coord[2]],
            mode='markers+text', text=[pid],
            marker=dict(size=4, color='black'), name=f'Punkt {pid}'
        ))

    # 4. Layout
    fig.update_layout(
        scene=dict(
            xaxis_title='X [m]',
            yaxis_title='Y [m]',
            zaxis_title='Z [m]',
            aspectmode='data'  # WICHTIG: Verhältnisse 1:1:1 bewahren
        ),
        width=1000, height=800,
        title="Geozentrisches Netz in 3D"
    )

    fig.show()

    # Aufruf
