Ellipsoid-Auswahl angepasst

revert Delete Tests/gha_resultsKarney.pkl

ES/Hansen_ES_CMA.py

@@ -0,0 +1,157 @@
+import numpy as np
+from numpy.typing import NDArray
+
+
+def felli(x: NDArray) -> float:
+    N = x.shape[0]
+    if N < 2:
+        raise ValueError("dimension must be greater than one")
+    exponents = np.arange(N) / (N - 1)
+    return float(np.sum((1e6 ** exponents) * (x ** 2)))
+
+
+def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000,
+          func_args=(), func_kwargs=None, seed=0,
+          bestEver=np.inf, noImproveGen=0, absTolImprove=1e-12, maxNoImproveGen=100, sigmaImprove=1e-12):
+    if func_kwargs is None:
+        func_kwargs = {}
+
+    if seed is not None:
+        np.random.seed(seed)
+
+    # Initialization (aus Parametern statt hart verdrahtet)
+    if xmean is None:
+        xmean = np.random.rand(N)
+    else:
+        xmean = np.asarray(xmean, dtype=float)
+        N = xmean.shape[0]
+
+    if stopeval is None:
+        stopeval = int(1e3 * N ** 2)
+
+    # Strategy parameter setting: Selection
+    lambda_ = 4 + int(np.floor(3 * np.log(N)))
+    mu = lambda_ / 2.0
+
+    # muXone recombination weights
+    weights = np.log(mu + 0.5) - np.log(np.arange(1, int(mu) + 1))
+    mu = int(np.floor(mu))
+    weights = weights / np.sum(weights)
+    mueff = np.sum(weights) ** 2 / np.sum(weights ** 2)
+
+    # Strategy parameter setting: Adaptation
+    cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N)
+    cs = (mueff + 2) / (N + mueff + 5)
+    c1 = 2 / ((N + 1.3) ** 2 + mueff)
+    cmu = min(1 - c1,
+              2 * (mueff - 2 + 1 / mueff) / ((N + 2) ** 2 + 2 * mueff))
+    damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs
+
+    # Initialize dynamic (internal) strategy parameters and constants
+    pc = np.zeros(N)
+    ps = np.zeros(N)
+    B = np.eye(N)
+    D = np.eye(N)
+    C = B @ D @ (B @ D).T
+    eigeneval = 0
+    chiN = np.sqrt(N) * (1 - 1 / (4 * N) + 1 / (21 * N ** 2))
+
+    # Generation Loop
+    counteval = 0
+    arx = np.zeros((N, lambda_))
+    arz = np.zeros((N, lambda_))
+    arfitness = np.zeros(lambda_)
+
+    gen = 0
+
+    # print(f'    [CMA-ES] Start: lambda = {lambda_}, sigma ={round(sigma, 6)}, stopeval = {stopeval}')
+
+    while counteval < stopeval:
+        gen += 1
+
+        # Generate and evaluate lambda offspring
+        for k in range(lambda_):
+            arz[:, k] = np.random.randn(N)
+            arx[:, k] = xmean + sigma * (B @ D @ arz[:, k])
+            arfitness[k] = float(func(arx[:, k], *func_args, **func_kwargs))  # <-- allgemein
+            counteval += 1
+
+        # Sort by fitness and compute weighted mean into xmean
+        idx = np.argsort(arfitness)
+        arfitness = arfitness[idx]
+        arindex = idx
+
+        xold = xmean.copy()
+        xmean = arx[:, arindex[:mu]] @ weights
+        zmean = arz[:, arindex[:mu]] @ weights
+
+        # Stagnation check
+        fbest = arfitness[0]
+        if bestEver - fbest > absTolImprove:
+            bestEver = fbest
+            noImproveGen = 0
+        else:
+            noImproveGen += 1
+
+        if gen == 1 or gen % 50 == 0:
+            # print(f'    [CMA-ES] Gen {gen}, best = {round(fbest, 6)}, sigma = {sigma:.3g}')
+            pass
+
+        if noImproveGen >= maxNoImproveGen:
+            # print(f'    [CMA-ES] Abbruch: keine Verbesserung > {round(absTolImprove, 3)} in {maxNoImproveGen} Generationen.')
+            break
+
+        if sigma < sigmaImprove:
+            # print(f'    [CMA-ES] Abbruch: sigma zu klein {sigma:.3g}')
+            break
+
+        # Cumulation: Update evolution paths
+        ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * (B @ zmean)
+        norm_ps = np.linalg.norm(ps)
+        hsig = norm_ps / np.sqrt(1 - (1 - cs) ** (2 * counteval / lambda_)) / chiN < (1.4 + 2 / (N + 1))
+        hsig = 1.0 if hsig else 0.0
+
+        pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (B @ D @ zmean)
+
+        # Adapt covariance matrix C
+        BDz = B @ D @ arz[:, arindex[:mu]]
+        C = (1 - c1 - cmu) * C \
+            + c1 * (np.outer(pc, pc) + (1 - hsig) * cc * (2 - cc) * C) \
+            + cmu * BDz @ np.diag(weights) @ BDz.T
+
+        # Adapt step-size sigma
+        sigma = sigma * np.exp((cs / damps) * (norm_ps / chiN - 1))
+
+        # Update B and D from C (Eigenzerlegung, O(N^2))
+        if counteval - eigeneval > lambda_ / ((c1 + cmu) * N * 10):
+            eigeneval = counteval
+            # enforce symmetry
+            C = (C + C.T) / 2.0
+            eigvals, B = np.linalg.eigh(C)
+            D = np.diag(np.sqrt(eigvals))
+
+        # Break, if fitness is good enough
+        if arfitness[0] <= stopfitness:
+            break
+
+        # Escape flat fitness, or better terminate?
+        if arfitness[0] == arfitness[int(np.ceil(0.7 * lambda_)) - 1]:
+            sigma = sigma * np.exp(0.2 + cs / damps)
+            # print('    [CMA-ES] stopfitness erreicht.')
+            # print("warning: flat fitness, consider reformulating the objective")
+            break
+
+        # print(f"{counteval}: {arfitness[0]}")
+
+    # Final Message
+    # print(f"{counteval}: {arfitness[0]}")
+    xmin = arx[:, arindex[0]]
+    bestValue = arfitness[0]
+    # print(f'    [CMA-ES] Ende: Gen = {gen}, best = {round(bestValue, 6)}')
+    return xmin
+
+
+if __name__ == "__main__":
+    xmin = escma(felli, N=10)  # <-- Zielfunktion wird übergeben
+    print("Bestes gefundenes x:", xmin)
+    print("f(xmin) =", felli(xmin))

ES/gha1_ES.py

@@ -0,0 +1,247 @@
+from __future__ import annotations
+
+from typing import List, Tuple
+
+import numpy as np
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+from ES.Hansen_ES_CMA import escma
+from GHA_triaxial.gha1_ana import gha1_ana
+from GHA_triaxial.gha1_approx import gha1_approx
+from GHA_triaxial.utils import jacobi_konstante
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_mpi_pi
+
+
+def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
+        -> Tuple[NDArray, NDArray, NDArray, float, float, NDArray]:
+    """
+    Analytische ENU-Basis in ellipsoidische Koordinaten (β, ω) nach Karney (2025), S. 2
+    :param beta: Beta Koordinate
+    :param omega: Omega Koordinate
+    :param ell: Ellipsoid
+    :return:    E_hat = Einheitsrichtung entlang wachsendem ω (East)
+                N_hat = Einheitsrichtung entlang wachsendem β (North)
+                U_hat = Einheitsnormale (Up)
+                En & Nn = Längen der unnormierten Ableitungen
+                R (XYZ) = Punkt in XYZ
+    """
+    # Berechnungshilfen
+    omega = wrap_mpi_pi(omega)
+    cb = np.cos(beta)
+    sb = np.sin(beta)
+    co = np.cos(omega)
+    so = np.sin(omega)
+    # D = sqrt(a^2 - c^2)
+    D = np.sqrt(ell.ax*ell.ax - ell.b*ell.b)
+    # Sx = sqrt(a^2 - b^2 sin^2β - c^2 cos^2β)
+    Sx = np.sqrt(ell.ax*ell.ax - ell.ay*ell.ay*(sb*sb) - ell.b*ell.b*(cb*cb))
+    # Sz = sqrt(a^2 sin^2ω + b^2 cos^2ω - c^2)
+    Sz = np.sqrt(ell.ax*ell.ax*(so*so) + ell.ay*ell.ay*(co*co) - ell.b*ell.b)
+
+    # Karney Gl. (4)
+    X = ell.ax * co * Sx / D
+    Y = ell.ay * cb * so
+    Z = ell.b * sb * Sz / D
+    R = np.array([X, Y, Z], dtype=float)
+
+    # --- Ableitungen - Karney Gl. (5a,b,c)---
+    # E = ∂R/∂ω
+    dX_dw = -ell.ax * so * Sx / D
+    dY_dw = ell.ay * cb * co
+    dZ_dw = ell.b * sb * (so * co * (ell.ax*ell.ax - ell.ay*ell.ay) / Sz) / D
+    E = np.array([dX_dw, dY_dw, dZ_dw], dtype=float)
+
+    # N = ∂R/∂β
+    dX_db = ell.ax * co * (sb * cb * (ell.b*ell.b - ell.ay*ell.ay) / Sx) / D
+    dY_db = -ell.ay * sb * so
+    dZ_db = ell.b * cb * Sz / D
+    N = np.array([dX_db, dY_db, dZ_db], dtype=float)
+
+    # U = Grad(x^2/a^2 + y^2/b^2 + z^2/c^2 - 1)
+    U = np.array([X/(ell.ax*ell.ax), Y/(ell.ay*ell.ay), Z/(ell.b*ell.b)], dtype=float)
+
+    En = float(np.linalg.norm(E))
+    Nn = float(np.linalg.norm(N))
+    Un = float(np.linalg.norm(U))
+
+    N_hat = N / Nn
+    E_hat = E / En
+    U_hat = U / Un
+    E_hat -= float(np.dot(E_hat, N_hat)) * N_hat
+    E_hat = E_hat / max(np.linalg.norm(E_hat), 1e-18)
+
+    return E_hat, N_hat, U_hat, En, Nn, R
+
+
+def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_hat_curr: NDArray, U_hat_curr: NDArray) -> float:
+    """
+    Berechnet das Azimut in der lokalen Tangentialebene am aktuellen Punkt P_curr, gemessen
+    an der Bewegungsrichtung vom vorherigen Punkt P_prev nach P_curr.
+    :param P_prev: vorheriger Punkt
+    :param P_curr: aktueller Punkt
+    :param E_hat_curr: Einheitsvektor der lokalen Tangentialrichtung am Punkt P_curr
+    :param N_hat_curr: Einheitsvektor der lokalen Tangentialrichtung am Punkt P_curr
+    :param U_hat_curr: Einheitsnormalenvektor am Punkt P_curr
+    :return: Azimut in Radiant
+    """
+    v = (P_curr - P_prev).astype(float)
+    vT = v - float(np.dot(v, U_hat_curr)) * U_hat_curr
+    vTn = max(np.linalg.norm(vT), 1e-18)
+    vT_hat = vT / vTn
+    sE = float(np.dot(vT_hat, E_hat_curr))
+    sN = float(np.dot(vT_hat, N_hat_curr))
+
+    return wrap_mpi_pi(float(np.arctan2(sE, sN)))
+
+
+def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float, gamma0: float,
+                        ell: EllipsoidTriaxial, maxSegLen: float = 1000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]:
+    """
+    Berechnung der 1. GHA mithilfe der CMA-ES.
+    Die CMA-ES optimiert sukzessive einen Punkt, der maxSegLen vom vorherigen Punkt entfernt und zusätzlich auf der
+    geodätischen Linien liegt. Somit entsteht ein Geodäten ähnlicher Polygonzug auf der Oberfläche des dreiachsigen Ellipsoids.
+    :param beta_i: Beta Koordinate am Punkt i
+    :param omega_i: Omega Koordinate am Punkt i
+    :param alpha_i: Azimut am Punkt i
+    :param ds: Gesamtlänge
+    :param gamma0: Jacobi-Konstante am Startpunkt
+    :param ell: Ellipsoid
+    :param maxSegLen: maximale Segmentlänge
+    :param sigma0:
+    :return:
+    """
+    # Startbasis
+    E_i, N_i, U_i, En_i, Nn_i, P_i = ENU_beta_omega(beta_i, omega_i, ell)
+    # Prediktor: dβ ≈ ds cosα / |N|, dω ≈ ds sinα / |E|
+
+    En_eff = max(En_i, 1e-9)
+    Nn_eff = max(Nn_i, 1e-9)
+
+    d_beta = ds * np.cos(alpha_i) / Nn_eff
+    d_omega = ds * np.sin(alpha_i) / En_eff
+
+    # optional: harte Schritt-Clamps (verhindert wrap-chaos)
+    d_beta = float(np.clip(d_beta, -0.2, 0.2))  # rad
+    d_omega = float(np.clip(d_omega, -0.2, 0.2))  # rad
+
+    # d_beta = ds * float(np.cos(alpha_i)) / Nn_i
+    # d_omega = ds * float(np.sin(alpha_i)) / En_i
+    beta_pred = beta_i + d_beta
+    omega_pred = wrap_mpi_pi(omega_i + d_omega)
+
+    xmean = np.array([beta_pred, omega_pred], dtype=float)
+
+    if sigma0 is None:
+        R0 = (ell.ax + ell.ay + ell.b) / 3
+        sigma0 = 1e-5 * (ds / R0)
+
+    def fitness(x: NDArray) -> float:
+        """
+        Fitnessfunktion: Fitnesscheck erfolgt anhand der Segmentlänge und der Jacobi-Konstante.
+        Die Segmentlänge muss möglichst gut zum Sollwert passen. Die Jacobi-Konstante am Punkt x muss zur
+        Jacobi-Konstanten am Startpunkt passen, damit der Polygonzug auf derselben geodätischen Linie bleibt.
+        :param x: Koordinate in beta, lambda aus der CMA-ES
+        :return: Fitnesswert (f)
+        """
+        beta = x[0]
+        omega = wrap_mpi_pi(x[1])
+
+        P = ell.ell2cart_karney(beta, omega)  # in kartesischer Koordinaten
+        d = float(np.linalg.norm(P - P_i))  # Distanz zwischen
+
+        # maxSegLen einhalten
+        J_len = ((d - ds) / ds) ** 2
+        w_len = 1.0
+
+        # Azimut für Jacobi-Konstante
+        E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta, omega, ell)
+        alpha_end = azimuth_at_ESpoint(P_i, P, E_j, N_j, U_j)
+
+        # Jacobi-Konstante
+        g_end = jacobi_konstante(beta, omega, alpha_end, ell)
+        J_gamma = (g_end - gamma0) ** 2
+        w_gamma = 10
+
+        f = float(w_len * J_len + w_gamma * J_gamma)
+
+        return f
+
+    xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0)  # Aufruf CMA-ES
+
+    beta_best = xb[0]
+    omega_best = wrap_mpi_pi(xb[1])
+    P_best = ell.ell2cart_karney(beta_best, omega_best)
+    E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta_best, omega_best, ell)
+    alpha_end = azimuth_at_ESpoint(P_i, P_best, E_j, N_j, U_j)
+
+    return beta_best, omega_best, P_best, alpha_end
+
+
+def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000, all_points: bool = False)\
+        -> Tuple[NDArray, float, NDArray] | Tuple[NDArray, float]:
+    """
+    Aufruf der 1. GHA mittels CMA-ES
+    :param ell: Ellipsoid
+    :param beta0: Beta Startkoordinate
+    :param omega0: Omega Startkoordinate
+    :param alpha0: Azimut Startkoordinate
+    :param s_total: Gesamtstrecke
+    :param maxSegLen: maximale Segmentlänge
+    :param all_points: Alle Punkte ausgeben?
+    :return: Zielpunkt Pk, Azimut am Zielpunkt und Punktliste
+    """
+    beta = float(beta0)
+    omega = wrap_mpi_pi(float(omega0))
+    alpha = wrap_mpi_pi(float(alpha0))
+
+    gamma0 = jacobi_konstante(beta, omega, alpha, ell)  # Referenz-γ0
+
+    P_all: List[NDArray] = [ell.ell2cart_karney(beta, omega)]
+    alpha_end: List[float] = [alpha]
+
+    s_acc = 0.0
+    step = 0
+    nsteps_est = int(np.ceil(s_total / maxSegLen))
+    while s_acc < s_total - 1e-9:
+        step += 1
+        ds = min(maxSegLen, s_total - s_acc)
+        # print(f"[GHA1-ES] Step {step}/{nsteps_est}  ds={ds:.3f} m  s_acc={s_acc:.3f} m  beta={beta:.6f}  omega={omega:.6f}  alpha={alpha:.6f}")
+
+        beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_i=alpha, ds=ds, gamma0=gamma0,
+                                                    ell=ell, maxSegLen=maxSegLen)
+        s_acc += ds
+        P_all.append(P)
+        alpha_end.append(wrap_mpi_pi(alpha))
+        if step > nsteps_est + 50:
+            raise RuntimeError("GHA1_ES: Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
+    Pk = P_all[-1]
+    alpha1 = float(alpha_end[-1])
+
+    if all_points:
+        return Pk, alpha1, np.array(P_all)
+    else:
+        return Pk, alpha1
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    s = 180000
+    # alpha0 = 3
+    alpha0 = wu.gms2rad([5, 0, 0])
+    beta = 0
+    omega = 0
+    P0 = ell.ell2cart(beta, omega)
+    point1, alpha1 = gha1_ana(ell, P0, alpha0=alpha0, s=s, maxM=100, maxPartCircum=32)
+    point1app, alpha1app = gha1_approx(ell, P0, alpha0=alpha0, s=s, ds=1000)
+
+    res, alpha, points = gha1_ES(ell, beta0=beta, omega0=-omega, alpha0=alpha0, s_total=s, maxSegLen=1000)
+
+    print(point1)
+    print(res)
+    print(alpha)
+    print(points)
+    # print("alpha1 (am Endpunkt):", res.alpha1)
+    print(res - point1)
+    print(point1app - point1, "approx")

ES/gha2_ES.py

@@ -0,0 +1,182 @@
+from typing import Tuple
+
+import numpy as np
+import plotly.graph_objects as go
+from numpy.typing import NDArray
+
+from ES.Hansen_ES_CMA import escma
+from GHA_triaxial.gha2_num import gha2_num
+from GHA_triaxial.utils import sigma2alpha
+from ellipsoid_triaxial import EllipsoidTriaxial
+
+
+def Sehne(P1: NDArray, P2: NDArray) -> float:
+    """
+    Berechnung der 3D-Distanz zwischen zwei kartesischen Punkten
+    :param P1: kartesische Koordinate Punkt 1
+    :param P2: kartesische Koordinate Punkt 2
+    :return: Bogenlänge s
+    """
+    R12 = P2-P1
+    s = float(np.linalg.norm(R12))
+
+    return s
+
+
+def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, all_points: bool = False) -> Tuple[float, float, float, NDArray] | Tuple[float, float, float]:
+    """
+    Berechnen der 2. GHA mithilfe der CMA-ES.
+    Die CMA-ES optimiert sukzessive den Mittelpunkt zwischen Start- und Zielpunkt. Der Abbruch der Berechnung erfolgt, wenn alle Segmentlängen <= maxSegLen sind.
+    Die Distanzen zwischen den einzelnen Punkten werden als direkte 3D-Distanzen berechnet und aufaddiert.
+    :param ell: Ellipsoid
+    :param P0: Startpunkt
+    :param Pk: Zielpunkt
+    :param maxSegLen: maximale Segmentlänge
+    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird
+    :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge
+    """
+    P_left: NDArray = None
+    P_right: NDArray = None
+
+    def midpoint_fitness(x: tuple) -> float:
+        """
+        Fitness für einen Mittelpunkt P_middle zwischen P_left und P_right auf dem triaxialen Ellipsoid:
+        - Minimiert d(P_left, P_middle) + d(P_middle, P_right)
+        - Erzwingt d(P_left,P_middle) ≈ d(P_middle,P_right)  (echter Mittelpunkt im Sinne der Polygonkette)
+        :param x: enthält die Startwerte von u und v
+        :return: Fitnesswert (f)
+        """
+        nonlocal P_left, P_right, ell
+
+        u, v = x
+        P_middle = ell.para2cart(u, v)
+        d1 = Sehne(P_left, P_middle)
+        d2 = Sehne(P_middle, P_right)
+        base = d1 + d2
+
+        # midpoint penalty (dimensionslos)
+        # relative Differenz, skaliert über verschiedene Segmentlängen
+        denom = max(base, 1e-9)
+        pen_equal = ((d1 - d2) / denom) ** 2
+        w_equal = 10.0
+
+        f = base + denom * w_equal * pen_equal
+
+        return f
+
+    R0 = (ell.ax + ell.ay + ell.b) / 3
+    if maxSegLen is None:
+        maxSegLen = R0 * 1 / (637.4*2)  # 10km Segment bei mittleren Erdradius
+
+    sigma_uv_nom = 1e-3 * (maxSegLen / R0)  # ~1e-5
+    points: list[NDArray] = [P0, Pk]
+    startIter = 0
+    level = 0
+
+    while True:
+        seg_lens = [Sehne(points[i], points[i+1]) for i in range(len(points)-1)]
+        max_len = max(seg_lens)
+        if max_len <= maxSegLen:
+            break
+
+        level += 1
+        new_points: list[NDArray] = [points[0]]
+        for i in range(len(points) - 1):
+            A = points[i]
+            B = points[i+1]
+            dAB = Sehne(A, B)
+            # print(dAB)
+
+            if dAB > maxSegLen:
+                # global P_left, P_right
+                P_left, P_right = A, B
+                Au, Av = ell.cart2para(A)
+                Bu, Bv = ell.cart2para(B)
+                u0 = (Au + Bu) / 2
+                v0 = Av + 0.5 * np.arctan2(np.sin(Bv - Av), np.cos(Bv - Av))
+                xmean = [u0, v0]
+
+                sigmaStep = sigma_uv_nom * (Sehne(A, B) / maxSegLen)
+
+                u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep)  # Aufruf CMA-ES
+
+                P_next = ell.para2cart(u, v)
+                new_points.append(P_next)
+                startIter += 1
+                maxIter = 10000
+                if startIter > maxIter:
+                    raise RuntimeError("GHA2_ES: maximale Iterationen überschritten")
+            new_points.append(B)
+
+        points = new_points
+        # print(f"[Level {level}] Punkte: {len(points)} | max Segment: {max_len:.3f} m")
+
+    P_all = np.vstack(points)
+    totalLen = float(np.sum(np.linalg.norm(P_all[1:] - P_all[:-1], axis=1)))
+
+    if len(points) >= 3:
+        p0i = ell.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0))
+        sigma0 = (p0i - P0) / np.linalg.norm(p0i - P0)
+        alpha0 = sigma2alpha(ell, sigma0, P0)
+
+        p1i = ell.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2]))
+        sigma1 = (Pk - p1i) / np.linalg.norm(Pk - p1i)
+        alpha1 = sigma2alpha(ell, sigma1, Pk)
+    else:
+        alpha0 = None
+        alpha1 = None
+
+    if all_points:
+        return alpha0, alpha1, totalLen, P_all
+    return alpha0, alpha1, totalLen
+
+
+def show_points(points: NDArray, pointsES: NDArray, p0: NDArray, p1: NDArray):
+    """
+    Anzeigen der Punkte
+    :param points: wahre Punkte der Linie
+    :param pointsES: Punkte der Linie aus ES
+    :param p0: wahrer Startpunkt
+    :param p1: wahrer Endpunkt
+    """
+    fig = go.Figure()
+
+    fig.add_scatter3d(x=pointsES[:, 0], y=pointsES[:, 1], z=pointsES[:, 2],
+                      mode='lines', line=dict(color="green", width=3), name="Numerisch")
+    fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2],
+                      mode='lines', line=dict(color="red", width=3), name="ES")
+    fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
+                      mode='markers', marker=dict(color="black"), name="P0")
+    fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
+                      mode='markers', marker=dict(color="black"), name="P1")
+
+    fig.update_layout(
+        scene=dict(xaxis_title='X [km]',
+                   yaxis_title='Y [km]',
+                   zaxis_title='Z [km]',
+                   aspectmode='data'))
+
+    fig.show()
+
+
+if __name__ == '__main__':
+    ell = EllipsoidTriaxial.init_name("Bursa1970")
+
+    beta0, lamb0 = (0.2, 0.1)
+    P0 = ell.ell2cart(beta0, lamb0)
+    beta1, lamb1 = (0.3, 0.2)
+    P1 = ell.ell2cart(beta1, lamb1)
+
+    alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000, all_points=True)
+    points_num = []
+    for beta, lamb in zip(betas, lambs):
+        points_num.append(ell.ell2cart(beta, lamb))
+    points_num = np.array(points_num)
+
+    alpha0, alpha1, s, points = gha2_ES(ell, P0, P1)
+    print(s_num)
+    print(s)
+    print(alpha0)
+    print(alpha1)
+    print(s - s_num)
+    show_points(points, points_num, P0, P1)

GHA/bessel.py

@@ -1,25 +0,0 @@
-from numpy import *
-import scipy as sp
-
-
-def gha1(re, phi_p1, lambda_p1, A_p1, s):
-    psi_p1 = re.phi2psi(phi_p1)
-    A_0 = arcsin(cos(psi_p1) * sin(A_p1))
-    temp = sin(psi_p1) / cos(A_0)
-    sigma_p1 = arcsin(sin(psi_p1) / cos(A_0))
-
-    sqrt_sigma = lambda sigma: sqrt(1 + re.e_ ** 2 * cos(A_0) ** 2 * sin(sigma) ** 2)
-    int_sqrt_sigma = lambda sigma: sp.integrate.quad(sqrt_sigma, sigma_p1, sigma)[0]
-
-    f_sigma_p2i = lambda sigma_p2i: (int_sqrt_sigma(sigma_p2i) - s / re.b)
-
-    sigma_p2_0 = sigma_p1 + s / re.a
-    sigma_p2 = sp.optimize.newton(f_sigma_p2i, sigma_p2_0)
-    psi_p2 = arcsin(cos(A_0) * sin(sigma_p2))
-    phi_p2 = re.psi2phi(psi_p2)
-    A_p2 = arcsin(sin(A_0) / cos(psi_p2))
-
-    f_d_lambda = lambda sigma: sin(A_0) * sqrt_sigma(sigma) / (1 - cos(A_0)**2 * sin(sigma)**2)
-    d_lambda = sqrt(1-re.e**2) * sp.integrate.quad(f_d_lambda, sigma_p1, sigma_p2)[0]
-    lambda_p2 = lambda_p1 + d_lambda
-    return phi_p2, lambda_p2, A_p2

GHA/gauss.py

@@ -1,107 +0,0 @@
-from numpy import sin, cos, pi, sqrt, tan, arcsin, arccos, arctan
-import ausgaben as aus
-
-
-def gha1(re, phi_p1, lambda_p1, A_p1, s, eps):
-    """
-    Berechnung der 1. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln
-    :param re: Klasse Ellipsoid
-    :param phi_p1: Breite Punkt 1
-    :param lambda_p1: Länge Punkt 1
-    :param A_p1: Azimut der geodätischen Linie in Punkt 1
-    :param s: Strecke zu Punkt 2
-    :param eps: Abbruchkriterium für Winkelgrößen
-    :return: Breite, Länge, Azimut von Punkt 2
-    """
-
-    t = lambda phi: tan(phi)
-    eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2)
-
-    F1 = lambda A, phi, s: 1 + (2 + 3 * t(phi)**2 + 2 * eta(phi)**2) / (24 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \
-                           + (t(phi)**2 - 1) * eta(phi) ** 2 / (8 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2
-
-    F2 = lambda A, phi, s: 1 + t(phi) ** 2 / (24 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \
-                           - (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / (24 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2
-
-    F3 = lambda A, phi, s: 1 + (1 + eta(phi)**2) / (12 * re.N(phi) ** 2) * sin(A) ** 2 * s ** 2 \
-                           + (3 + 8 * eta(phi)**2) / (24 * re.N(phi) ** 2) * cos(A) ** 2 * s ** 2
-
-    phi_p2_i = lambda A, phi: phi_p1 + cos(A) / re.M(phi) * s * F1(A, phi, s)
-    lambda_p2_i = lambda A, phi: lambda_p1 + sin(A) / (re.N(phi) * cos(phi)) * s * F2(A, phi, s)
-    A_p2_i = lambda A, phi: A_p1 + sin(A) * tan(phi) / re.N(phi) * s * F3(A, phi, s)
-
-    phi_p0_i = lambda phi2: (phi_p1 + phi2) / 2
-    A_p1_i = lambda A2: (A_p1 + A2) / 2
-
-    phi_p0 = []
-    A_p0 = []
-    phi_p2 = []
-    lambda_p2 = []
-    A_p2 = []
-
-    # 1. Näherung für P2
-    phi_p2.append(phi_p1 + cos(A_p1) / re.M(phi_p1) * s)
-    lambda_p2.append(lambda_p1 + sin(A_p1) / (re.N(phi_p1) * cos(phi_p1)) * s)
-    A_p2.append(A_p1 + sin(A_p1) * tan(phi_p1) / re.N(phi_p1) * s)
-
-    while True:
-        # Berechnug P0 durch Mittelbildung
-        phi_p0.append(phi_p0_i(phi_p2[-1]))
-        A_p0.append(A_p1_i(A_p2[-1]))
-        # Berechnung P2
-        phi_p2.append(phi_p2_i(A_p0[-1], phi_p0[-1]))
-        lambda_p2.append(lambda_p2_i(A_p0[-1], phi_p0[-1]))
-        A_p2.append(A_p2_i(A_p0[-1], phi_p0[-1]))
-
-        # Abbruchkriterium
-        if abs(phi_p2[-2] - phi_p2[-1]) < eps and \
-           abs(lambda_p2[-2] - lambda_p2[-1]) < eps and \
-           abs(A_p2[-2] - A_p2[-1]) < eps:
-            break
-
-    nks = 5
-    for i in range(len(phi_p2)):
-        print(f"P2[{i}]: {aus.gms('phi', phi_p2[i], nks)}\t{aus.gms('lambda', lambda_p2[i], nks)}\t{aus.gms('A', A_p2[i], nks)}")
-        if i != len(phi_p2)-1:
-            print(f"P0[{i}]: {aus.gms('phi', phi_p0[i], nks)}\t\t\t\t\t\t\t\t{aus.gms('A', A_p0[i], nks)}")
-    return phi_p2, lambda_p2, A_p2
-
-
-def gha2(re, phi_p1, lambda_p1, phi_p2, lambda_p2):
-    """
-    Berechnung der 2. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln
-    :param re: Klasse Ellipsoid
-    :param phi_p1: Breite Punkt 1
-    :param lambda_p1: Länge Punkt 1
-    :param phi_p2: Breite Punkt 2
-    :param lambda_p2: Länge Punkt 2
-    :return: Länge der geodätischen Linie, Azimut von P1 nach P2, Azimut von P2 nach P1
-    """
-
-    t = lambda phi: tan(phi)
-    eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2)
-
-    phi_0 = (phi_p1 + phi_p2) / 2
-    d_phi = phi_p2 - phi_p1
-    d_lambda = lambda_p2 - lambda_p1
-
-    f_A = lambda phi: (2 + 3*t(phi)**2 + 2*eta(phi)**2) / 24
-    f_B = lambda phi: ((t(phi)**2 - 1) * eta(phi)**2) / 8
-    # f_C = lambda phi: (t(phi)**2) / 24
-    f_D = lambda phi: (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / 24
-
-    F1 = lambda phi: d_phi * re.M(phi) * (1 - f_A(phi) * d_lambda ** 2 * cos(phi) ** 2 -
-                                          f_B(phi) * d_phi ** 2 / re.V(phi) ** 4)
-    F2 = lambda phi: d_lambda * re.N(phi) * cos(phi) * (1 - 1 / 24 * d_lambda ** 2 * sin(phi) ** 2 +
-                                                           f_D(phi) * d_phi ** 2 / re.V(phi) ** 4)
-
-    s = sqrt(F1(phi_0) ** 2 + F2(phi_0) ** 2)
-    A_0 = arctan(F2(phi_0) / F1(phi_0))
-    d_A = d_lambda * sin(phi_0) * (1 + (1 + eta(phi_0) ** 2) / 12 * s ** 2 * sin(A_0) ** 2 / re.N(phi_0) ** 2 +
-                                      (3 + 8 * eta(phi_0) ** 2) / 24 * s ** 2 * cos(A_0) ** 2 / re.N(phi_0) ** 2)
-
-    A_p1 = A_0 - d_A / 2
-    A_p2 = A_0 + d_A / 2
-    A_p2_p1 = A_p2 + pi
-
-    return s, A_p1, A_p2_p1

GHA/rk.py

@@ -1,17 +0,0 @@
-import Numerische_Integration.num_int_runge_kutta as rk
-from numpy import sin, cos, tan
-import winkelumrechnungen as wu
-from ellipsoide import EllipsoidBiaxial
-
-def gha1(re, x0, y0, z0, A0, s, num):
-    phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0)
-
-    f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c
-    f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c)
-    f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c
-
-    funktionswerte = rk.verfahren([f_phi, f_lam, f_A],
-                                  [0, phi0, lamb0, A0],
-                                  s, num)
-    coords = re.ell2cart(funktionswerte[-1][1], funktionswerte[-1][2], h0)
-    return coords

GHA_triaxial/gha1_ana.py

@@ -0,0 +1,145 @@
+import math
+from math import comb
+from typing import Tuple
+
+import numpy as np
+from numpy import arctan2, cos, sin
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+from GHA_triaxial.utils import pq_para
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_0_2pi
+
+
+def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]:
+    """
+    Panou, Korakitits 2020, 5ff.
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha0: Azimut im Startpunkt
+    :param s: Strecke
+    :param maxM: maximale Ordnung
+    :return: Zwischenpunkt, Azimut im Zwischenpunkt
+    """
+    x, y, z = point
+
+    # S. 6
+    x_m = [x]
+    y_m = [y]
+    z_m = [z]
+
+    p, q = pq_para(ell, point)
+
+    # 48-50
+    x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0))
+    y_m.append(p[1] * sin(alpha0) + q[1] * cos(alpha0))
+    z_m.append(p[2] * sin(alpha0) + q[2] * cos(alpha0))
+
+    # 34
+    H_ = lambda p: np.sum([comb(p, p - i) * (x_m[p - i] * x_m[i] +
+                                             1 / (1 - ell.ee ** 2) ** 2 * y_m[p - i] * y_m[i] +
+                                             1 / (1 - ell.ex ** 2) ** 2 * z_m[p - i] * z_m[i])
+                           for i in range(0, p + 1)])
+
+    # 35
+    h_ = lambda q: np.sum([comb(q, q-j) * (x_m[q-j+1] * x_m[j+1] +
+                                           1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] +
+                                           1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1])
+                           for j in range(0, q+1)])
+
+    # 31
+    hH_ = lambda t: 1/H_(0) * (h_(t) - np.sum([comb(t, l-1) * H_(t+1-l) * hH_t[l-1] for l in range(1, t+1)]))
+
+    # 28-30
+    x_ = lambda m: - np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)])
+    y_ = lambda m: -1 / (1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)])
+    z_ = lambda m: -1 / (1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)])
+
+    hH_t = []
+    a_m = []
+    b_m = []
+    c_m = []
+    for m in range(0, maxM+1):
+        if m >= 2:
+            hH_t.append(hH_(m-2))
+            x_m.append(x_(m))
+            y_m.append(y_(m))
+            z_m.append(z_(m))
+        fact_m = math.factorial(m)
+
+        # 22-24
+        a_m.append(x_m[m] / fact_m)
+        b_m.append(y_m[m] / fact_m)
+        c_m.append(z_m[m] / fact_m)
+
+    # 19-21
+    x_s = 0
+    for a in reversed(a_m):
+        x_s = x_s * s + a
+    y_s = 0
+    for b in reversed(b_m):
+        y_s = y_s * s + b
+    z_s = 0
+    for c in reversed(c_m):
+        z_s = z_s * s + c
+
+    p1 = np.array([x_s, y_s, z_s])
+    p_s, q_s = pq_para(ell, p1)
+
+    # 57-59
+    dx_s = 0
+    for i, a in reversed(list(enumerate(a_m[1:], start=1))):
+        dx_s = dx_s * s + i * a
+
+    dy_s = 0
+    for i, b in reversed(list(enumerate(b_m[1:], start=1))):
+        dy_s = dy_s * s + i * b
+
+    dz_s = 0
+    for i, c in reversed(list(enumerate(c_m[1:], start=1))):
+        dz_s = dz_s * s + i * c
+
+    # 52-53
+    sigma = np.array([dx_s, dy_s, dz_s])
+    P = float(p_s @ sigma)
+    Q = float(q_s @ sigma)
+
+    # 51
+    alpha1 = arctan2(P, Q)
+
+    if alpha1 < 0:
+        alpha1 += 2 * np.pi
+
+    return p1, wrap_0_2pi(alpha1)
+
+
+def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]:
+    """
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha0: Azimut im Startpunkt
+    :param s: Strecke
+    :param maxM: maximale Ordnung
+    :param maxPartCircum: maximale Aufteilung (1/x halber Ellipsoidumfang)
+    :return: Zielpunkt, Azimut im Zielpunkt
+    """
+    if s > np.pi / maxPartCircum * ell.ax:
+        s /= 2
+        point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum)
+        point_end, alpha_end = gha1_ana(ell, point_step, alpha_step, s, maxM, maxPartCircum)
+    else:
+        point_end, alpha_end = gha1_ana_step(ell, point, alpha0, s, maxM)
+
+    _, _, h = ell.cart2geod(point_end, "ligas3")
+    if h > 1e-5:
+        raise Exception("GHA1_ana: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
+
+    return point_end, wrap_0_2pi(alpha_end)
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    p0 = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(20))
+    p1, alpha1 = gha1_ana(ell, p0, wu.deg2rad(36), 200000, 70)
+    print(p1, wu.rad2gms(alpha1))

GHA_triaxial/gha1_approx.py

@@ -0,0 +1,108 @@
+from typing import Tuple
+
+import numpy as np
+import plotly.graph_objects as go
+from numpy import cos, sin
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+from GHA_triaxial.utils import louville_constant, pq_ell
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_0_2pi
+
+
+def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) \
+        -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray, NDArray]:
+    """
+    Berechung einer Näherungslösung der ersten Hauptaufgabe
+    :param ell: Ellipsoid
+    :param p0: Anfangspunkt
+    :param alpha0: Azimut im Anfangspunkt
+    :param s: Strecke bis zum Endpunkt
+    :param ds: Länge einzelner Streckenelemente
+    :param all_points: Ausgabe aller Punkte als Array?
+    :return: Endpunkt, Azimut im Endpunkt, optional alle Punkte
+    """
+    l0 = louville_constant(ell, p0, alpha0)
+    points = [p0]
+    alphas = [alpha0]
+    s_curr = 0.0
+    last_sigma = None
+    last_p = None
+
+    while s_curr < s:
+        ds_step = min(ds, s - s_curr)
+        if ds_step < 1e-8:
+            break
+
+        p1 = points[-1]
+        alpha1 = alphas[-1]
+
+        p, q = pq_ell(ell, p1)
+        if last_p is not None and np.dot(p, last_p) < 0:
+            p = -p
+            q = -q
+        last_p = p
+        sigma = p * sin(alpha1) + q * cos(alpha1)
+        if last_sigma is not None and np.dot(sigma, last_sigma) < 0:
+            sigma = -sigma
+            alpha1 += np.pi
+            alpha1 = wrap_0_2pi(alpha1)
+        p2 = p1 + ds_step * sigma
+        p2 = ell.point_onto_ellipsoid(p2)
+
+        dalpha = 1e-9
+        l2 = louville_constant(ell, p2, alpha1)
+        dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
+        if abs(dl_dalpha) < 1e-20:
+            alpha2 = alpha1 + 0
+        else:
+            alpha2 = alpha1 + (l0 - l2) / dl_dalpha
+        points.append(p2)
+        alphas.append(wrap_0_2pi(alpha2))
+
+        ds_step = np.linalg.norm(p2 - p1)
+        s_curr += ds_step
+        last_sigma = sigma
+        pass
+
+    if all_points:
+        return points[-1], alphas[-1], np.array(points), np.array(alphas)
+    else:
+        return points[-1], alphas[-1]
+
+def show_points(points: NDArray, p0: NDArray, p1: NDArray):
+    """
+    Anzeigen der Punkte
+    :param points: Array aller approximierten Punkte
+    :param p0: Startpunkt
+    :param p1: wahrer Endpunkt
+    """
+    fig = go.Figure()
+
+    fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2],
+                      mode='lines', line=dict(color="red", width=3), name="Approx")
+    fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
+                      mode='markers', marker=dict(color="green"), name="P0")
+    fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
+                      mode='markers', marker=dict(color="green"), name="P1")
+
+    fig.update_layout(
+        scene=dict(xaxis_title='X [km]',
+                   yaxis_title='Y [km]',
+                   zaxis_title='Z [km]',
+                   aspectmode='data'),
+        title="CHAMP")
+
+    fig.show()
+
+
+if __name__ == '__main__':
+    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    P0 = ell.ell2cart(wu.deg2rad(15), wu.deg2rad(15))
+    alpha0 = wu.deg2rad(270)
+    s = 1
+    P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=0.1, all_points=True)
+    # P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=40, maxPartCircum=32)
+    # print(np.linalg.norm(P1_app - P1_ana))
+    # show_points(points, P0, P0)

GHA_triaxial/gha1_num.py

@@ -0,0 +1,133 @@
+from typing import Callable, List, Tuple
+
+import numpy as np
+from numpy import arctan2, cos, sin
+from numpy.typing import NDArray
+
+import GHA_triaxial.numeric_examples_karney as ne_karney
+import runge_kutta as rk
+import winkelumrechnungen as wu
+from GHA_triaxial.gha1_ana import gha1_ana
+from GHA_triaxial.utils import alpha_ell2para, pq_ell
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_0_2pi
+
+
+def buildODE(ell: EllipsoidTriaxial) -> Callable:
+    """
+    Aufbau des DGL-Systems
+    :param ell: Ellipsoid
+    :return: DGL-System
+    """
+
+    def ODE(s: float, v: NDArray) -> NDArray:
+        """
+        DGL-System
+        :param s: unabhängige Variable
+        :param v: abhängige Variablen
+        :return: Ableitungen der abhängigen Variablen
+        """
+        x, dxds, y, dyds, z, dzds = v
+
+        H = ell.func_H(np.array([x, y, z]))
+        h = dxds ** 2 + 1 / (1 - ell.ee ** 2) * dyds ** 2 + 1 / (1 - ell.ex ** 2) * dzds ** 2
+
+        ddx = -(h / H) * x
+        ddy = -(h / H) * y / (1 - ell.ee ** 2)
+        ddz = -(h / H) * z / (1 - ell.ex ** 2)
+
+        return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
+
+    return ODE
+
+def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, num: int, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, List]:
+    """
+    Panou, Korakitits 2019
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha0: Azimut im Startpunkt
+    :param s: Strecke
+    :param num: Anzahl Zwischenpunkte
+    :param all_points: Ausgabe aller Punkte?
+    :return: Zielpunkt, Azimut im Zielpunkt (, alle Punkte)
+    """
+    phi, lam, _ = ell.cart2geod(point, "ligas3")
+    p0 = ell.geod2cart(phi, lam, 0)
+    x0, y0, z0 = p0
+
+    p, q = pq_ell(ell, p0)
+    dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0)
+    dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0)
+    dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0)
+
+    v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0])
+
+    ode = buildODE(ell)
+
+    _, werte = rk.rk4(ode, 0, v_init, s, num)
+    x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1]
+
+    point1 = np.array([x1, y1, z1])
+
+    p1, q1 = pq_ell(ell, point1)
+    sigma = np.array([dx1ds, dy1ds, dz1ds])
+    P = float(p1 @ sigma)
+    Q = float(q1 @ sigma)
+
+    alpha1 = arctan2(P, Q)
+
+    alpha1 = wrap_0_2pi(alpha1)
+
+    _, _, h = ell.cart2geod(point1, "ligas3")
+    if h > 1e-5:
+        raise Exception("GHA1_num: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
+
+    if all_points:
+        return point1, alpha1, werte
+    else:
+        return point1, alpha1
+
+if __name__ == "__main__":
+
+    # ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
+    # diffs_panou = []
+    # examples_panou = ne_panou.get_random_examples(5)
+    # for example in examples_panou:
+    #     beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example
+    #     P0 = ell.ell2cart(beta0, lamb0)
+    #
+    #     P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 100)
+    #     beta1_num, lamb1_num = ell.cart2ell(P1_num)
+    #
+    #     _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
+    #     P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 60)
+    #     beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
+    #     diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana)))
+    # diffs_panou = np.array(diffs_panou)
+    # mask_360 = (diffs_panou > 359) & (diffs_panou < 361)
+    # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
+    # print(diffs_panou)
+
+    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name("KarneyTest2024")
+    diffs_karney = []
+    # examples_karney = ne_karney.get_examples((30499, 30500, 40500))
+    examples_karney = ne_karney.get_random_examples(20)
+    for example in examples_karney:
+        beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
+        P0 = ell.ell2cart(beta0, lamb0)
+
+        P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 10000)
+        beta1_num, lamb1_num = ell.cart2ell(P1_num)
+
+        try:
+            _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
+            P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 45, maxPartCircum=32)
+            beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
+        except:
+            beta1_ana, lamb1_ana = np.inf, np.inf
+
+        diffs_karney.append((wu.rad2deg(abs(beta1-beta1_num)), wu.rad2deg(abs(lamb1-lamb1_num)), wu.rad2deg(abs(beta1-beta1_ana)), wu.rad2deg(abs(lamb1-lamb1_ana))))
+    diffs_karney = np.array(diffs_karney)
+    mask_360 = (diffs_karney > 359) & (diffs_karney < 361)
+    diffs_karney[mask_360] = np.abs(diffs_karney[mask_360] - 360)
+    print(diffs_karney)

GHA_triaxial/gha2_approx.py

@@ -0,0 +1,105 @@
+from typing import Tuple
+
+import numpy as np
+import plotly.graph_objects as go
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+from GHA_triaxial.gha2_num import gha2_num
+from GHA_triaxial.utils import sigma2alpha
+from ellipsoid_triaxial import EllipsoidTriaxial
+
+
+def gha2_approx(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]:
+    """
+    Numerische Approximation für die zweite Hauptaufgabe
+    :param ell: Ellipsoid
+    :param p0: Startpunkt
+    :param p1: Endpunkt
+    :param ds: maximales Streckenelement
+    :param all_points: Alle Punkte ausgeben?
+    :return:
+    """
+    points = np.array([p0, p1])
+    while True:
+        new_points = []
+
+        for i in range(len(points)-1):
+            new_points.append(points[i])
+
+            pi = points[i] + 1/2 * (points[i+1] - points[i])
+            pi = ell.point_onto_ellipsoid(pi)
+
+            new_points.append(pi)
+
+        new_points.append(points[-1])
+        points = np.array(new_points)
+
+        elements = np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)])
+
+        if np.average(elements) < ds:
+            break
+
+    p0i = ell.point_onto_ellipsoid(p0 + ds / 100 * (points[1] - p0) / np.linalg.norm(points[1] - p0))
+    sigma0 = (p0i - p0) / np.linalg.norm(p0i - p0)
+    alpha0 = sigma2alpha(ell, sigma0, p0)
+
+    p1i = ell.point_onto_ellipsoid(p1 - ds / 100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2]))
+    sigma1 = (p1 - p1i) / np.linalg.norm(p1 - p1i)
+    alpha1 = sigma2alpha(ell, sigma1, p1)
+
+    s = np.sum(np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)]))
+
+    if all_points:
+        return alpha0, alpha1, s, np.array(points)
+    else:
+        return alpha0, alpha1, s
+
+def show_points(points: NDArray, points_app: NDArray, p0: NDArray, p1: NDArray):
+    """
+    Anzeigen der Punkte
+    :param points: wahre Punkte der Linie
+    :param points_app: approximierte Punkte der Linie
+    :param p0: wahrer Startpunkt
+    :param p1: wahrer Endpunkt
+    """
+    fig = go.Figure()
+
+    fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2],
+                      mode='lines', line=dict(color="green", width=3), name="Analytisch")
+    fig.add_scatter3d(x=points_app[:, 0], y=points_app[:, 1], z=points_app[:, 2],
+                      mode='lines', line=dict(color="red", width=3), name="Approximiert")
+    fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
+                      mode='markers', marker=dict(color="black"), name="P0")
+    fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
+                      mode='markers', marker=dict(color="black"), name="P1")
+
+    fig.update_layout(
+        scene=dict(xaxis_title='X [km]',
+                   yaxis_title='Y [km]',
+                   zaxis_title='Z [km]',
+                   aspectmode='data'))
+
+    fig.show()
+
+
+if __name__ == '__main__':
+    ell = EllipsoidTriaxial.init_name("Bursa1970")
+
+    beta0, lamb0 = (0.1, 0.1)
+    P0 = ell.ell2cart(beta0, lamb0)
+    beta1, lamb1 = (0.3, np.pi)
+    P1 = ell.ell2cart(beta1, lamb1)
+
+    alpha0_app, alpha1_app, s_app, points = gha2_approx(ell, P0, P1, ds=100, all_points=True)
+    print("done")
+    alpha0, alpha1, s, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=10000, all_points=True)
+    points_ana = []
+    for beta, lamb in zip(betas, lambs):
+        points_ana.append(ell.ell2cart(beta, lamb))
+    points_ana = np.array(points_ana)
+
+    show_points(points_ana, points, P0, P1)
+    print(f"Differenz s: {s_app - s} m")
+    print(f"Differenz alpha0: {wu.rad2deg(alpha0_app - alpha0)}°")
+    print(f"Differenz alpha1: {wu.rad2deg(alpha1_app - alpha1)}°")

GHA_triaxial/gha2_num.py

@@ -0,0 +1,643 @@
+from typing import Tuple
+
+import numpy as np
+from numpy.typing import NDArray
+
+import GHA_triaxial.numeric_examples_karney as ne_karney
+import GHA_triaxial.numeric_examples_panou as ne_panou
+import ausgaben as aus
+import winkelumrechnungen as wu
+from ellipsoid_triaxial import EllipsoidTriaxial
+from runge_kutta import rk4, rk4_end, rk4_integral
+from utils_angle import cot, wrap_0_2pi, wrap_mpi_pi
+
+
+def norm_a(a: float) -> float:
+    a = float(a) % (2 * np.pi)
+    return a
+
+
+def azimut(E: float, G: float, dbeta_du: float, dlamb_du: float) -> float:
+    north = np.sqrt(E) * dbeta_du
+    east = np.sqrt(G) * dlamb_du
+    return norm_a(np.arctan2(east, north))
+
+
+def sph_azimuth(beta1, lam1, beta2, lam2):
+    dlam = wrap_mpi_pi(lam2 - lam1)
+    y = np.sin(dlam) * np.cos(beta2)
+    x = np.cos(beta1) * np.sin(beta2) - np.sin(beta1) * np.cos(beta2) * np.cos(dlam)
+    a = np.arctan2(y, x)
+    if a < 0:
+        a += 2 * np.pi
+    return a
+
+
+# Panou 2013
+def gha2_num(
+    ell: EllipsoidTriaxial,
+    beta_0: float,
+    lamb_0: float,
+    beta_1: float,
+    lamb_1: float,
+    n: int = 16000,
+    epsilon: float = 10**-12,
+    iter_max: int = 30,
+    all_points: bool = False,
+) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]:
+    """
+    :param ell: Ellipsoid
+    :param beta_0: Beta Punkt 0
+    :param lamb_0: Lambda Punkt 0
+    :param beta_1: Beta Punkt 1
+    :param lamb_1: Lambda Punkt 1
+    :param n: Anzahl Schritte
+    :param epsilon: Genauigkeit
+    :param iter_max: Maximale Iterationen
+    :param all_points: Ausgabe aller Punkte
+    :return: Azimut Startpunkt, Azumit Zielpunkt, Strecke
+    """
+
+    ax2 = float(ell.ax) * float(ell.ax)
+    ay2 = float(ell.ay) * float(ell.ay)
+    b2 = float(ell.b) * float(ell.b)
+    Ex2 = float(ell.Ex) * float(ell.Ex)
+    Ey2 = float(ell.Ey) * float(ell.Ey)
+    Ee2 = float(ell.Ee) * float(ell.Ee)
+    Ey4 = Ey2 * Ey2
+    Ee4 = Ee2 * Ee2
+    two_pi = 2.0 * np.pi
+
+    # Berechnung Koeffizienten, Gaußschen Fundamentalgrößen 1. Ordnung sowie deren Ableitungen
+    def BETA_LAMBDA(beta, lamb):
+        sb = np.sin(beta)
+        cb = np.cos(beta)
+        sl = np.sin(lamb)
+        cl = np.cos(lamb)
+
+        sb2 = sb * sb
+        cb2 = cb * cb
+        sl2 = sl * sl
+        cl2 = cl * cl
+
+        s2b = 2.0 * sb * cb
+        c2b = cb2 - sb2
+        s2l = 2.0 * sl * cl
+        c2l = cl2 - sl2
+
+        denB = Ex2 - Ey2 * sb2
+        denL = Ex2 - Ee2 * cl2
+
+        BETA = (ay2 * sb2 + b2 * cb2) / denB
+        LAMBDA = (ax2 * sl2 + ay2 * cl2) / denL
+
+        BETA_ = (ax2 * Ey2 * s2b) / (denB * denB)
+        LAMBDA_ = -(b2 * Ee2 * s2l) / (denL * denL)
+
+        BETA__ = (2.0 * ax2 * Ey4 * (s2b * s2b)) / (denB * denB * denB) + (2.0 * ax2 * Ey2 * c2b) / (
+            denB * denB
+        )
+        LAMBDA__ = (2.0 * b2 * Ee4 * (s2l * s2l)) / (denL * denL * denL) - (2.0 * b2 * Ee2 * s2l) / (
+            denL * denL
+        )
+
+        Q = Ey2 * cb2 + Ee2 * sl2
+
+        E = BETA * Q
+        G = LAMBDA * Q
+
+        E_beta = BETA_ * Q - BETA * Ey2 * s2b
+        E_lamb = BETA * Ee2 * s2l
+
+        G_beta = -LAMBDA * Ey2 * s2b
+        G_lamb = LAMBDA_ * Q + LAMBDA * Ee2 * s2l
+
+        E_beta_beta = BETA__ * Q - 2.0 * BETA_ * Ey2 * s2b - 2.0 * BETA * Ey2 * c2b
+        E_beta_lamb = BETA_ * Ee2 * s2l
+        E_lamb_lamb = 2.0 * BETA * Ee2 * c2l
+
+        G_beta_beta = -2.0 * LAMBDA * Ey2 * c2b
+        G_beta_lamb = -LAMBDA_ * Ey2 * s2b
+        G_lamb_lamb = LAMBDA__ * Q + 2.0 * LAMBDA_ * Ee2 * s2l + 2.0 * LAMBDA * Ee2 * c2l
+
+        return (
+            BETA,
+            LAMBDA,
+            E,
+            G,
+            BETA_,
+            LAMBDA_,
+            BETA__,
+            LAMBDA__,
+            E_beta,
+            E_lamb,
+            G_beta,
+            G_lamb,
+            E_beta_beta,
+            E_beta_lamb,
+            E_lamb_lamb,
+            G_beta_beta,
+            G_beta_lamb,
+            G_lamb_lamb,
+        )
+
+    # Berechnung der ODE Koeffizienten für Fall 1 (lambda_0 != lambda_1)
+    def p_coef(beta, lamb):
+        (
+            BETA,
+            LAMBDA,
+            E,
+            G,
+            BETA_,
+            LAMBDA_,
+            BETA__,
+            LAMBDA__,
+            E_beta,
+            E_lamb,
+            G_beta,
+            G_lamb,
+            E_beta_beta,
+            E_beta_lamb,
+            E_lamb_lamb,
+            G_beta_beta,
+            G_beta_lamb,
+            G_lamb_lamb,
+        ) = BETA_LAMBDA(beta, lamb)
+
+        p_3 = -0.5 * (E_lamb / G)
+        p_2 = (G_beta / G) - 0.5 * (E_beta / E)
+        p_1 = 0.5 * (G_lamb / G) - (E_lamb / E)
+        p_0 = 0.5 * (G_beta / E)
+
+        p_33 = -0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G**2))
+        p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G**2)) - 0.5 * (
+            (E * E_beta_beta - E_beta * E_beta) / (E**2)
+        )
+        p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G**2)) - (
+            (E * E_beta_lamb - E_beta * E_lamb) / (E**2)
+        )
+        p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E**2))
+
+        return BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00
+
+    # Berechnung der ODE Koeffizienten für Fall 2 (lambda_0 == lambda_1)
+    def q_coef(beta, lamb):
+        (
+            BETA,
+            LAMBDA,
+            E,
+            G,
+            BETA_,
+            LAMBDA_,
+            BETA__,
+            LAMBDA__,
+            E_beta,
+            E_lamb,
+            G_beta,
+            G_lamb,
+            E_beta_beta,
+            E_beta_lamb,
+            E_lamb_lamb,
+            G_beta_beta,
+            G_beta_lamb,
+            G_lamb_lamb,
+        ) = BETA_LAMBDA(beta, lamb)
+
+        q_3 = -0.5 * (G_beta / E)
+        q_2 = (E_lamb / E) - 0.5 * (G_lamb / G)
+        q_1 = 0.5 * (E_beta / E) - (G_beta / G)
+        q_0 = 0.5 * (E_lamb / G)
+
+        q_33 = -0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E**2))
+        q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E**2)) - 0.5 * (
+            (G * G_lamb_lamb - G_lamb * G_lamb) / (G**2)
+        )
+        q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E**2)) - (
+            (G * G_beta_lamb - G_beta * G_lamb) / (G**2)
+        )
+        q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G**2))
+
+        return BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00
+
+    def integrand_lambda(lamb, y):
+        beta = y[0]
+        beta_p = y[1]
+        (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
+        return np.sqrt(E * beta_p**2 + G)
+
+    def integrand_beta(beta, y):
+        lamb = y[0]
+        lamb_p = y[1]
+        (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
+        return np.sqrt(E + G * lamb_p**2)
+
+    def solve_lambda_branch(beta0, lamb0, beta1, lamb1_target, N_run, N_newt, it_max):
+        dlamb = float(lamb1_target - lamb0)
+        if abs(dlamb) < 1e-15:
+            return None
+
+        sgn = 1.0 if dlamb >= 0.0 else -1.0
+
+        def ode_lamb(lamb, v):
+            beta, beta_p, X3, X4 = v
+            (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
+
+            dbeta = beta_p
+            dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0
+            dX3 = X4
+            dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (
+                3 * p_3 * beta_p**2 + 2 * p_2 * beta_p + p_1
+            ) * X4
+            return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
+
+        alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1_target)
+        (_, _, E0, G0, *_) = BETA_LAMBDA(beta0, lamb0)
+        beta_p0_sph = np.sqrt(G0 / E0) * cot(alpha0_sph)
+
+        def solve_newton(beta_p0_init: float):
+            beta_p0 = float(beta_p0_init)
+            for _ in range(it_max):
+                v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
+                _, y_end = rk4_end(ode_lamb, lamb0, v0, dlamb, N_newt)
+
+                beta_end, _, X3_end, _ = y_end
+                delta = beta_end - beta1
+
+                if abs(delta) < epsilon:
+                    return True, beta_p0
+
+                if abs(X3_end) < 1e-20:
+                    return False, None
+
+                step = delta / X3_end
+                step = float(np.clip(step, -0.5, 0.5))
+                beta_p0 -= step
+
+            return False, None
+
+        seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
+        best = None
+
+        for seed in seeds:
+            ok, sol = solve_newton(seed)
+            if not ok:
+                continue
+            v0_sol = np.array([beta0, sol, 0.0, 1.0], dtype=float)
+            _, _, s_val = rk4_integral(ode_lamb, lamb0, v0_sol, dlamb, N_run, integrand_lambda)
+            if (best is None) or (s_val < best[0]):
+                best = (float(s_val), float(sol))
+
+        if best is None:
+            return None
+
+        return best[0], best[1], sgn, dlamb, ode_lamb
+
+    def solve_beta_branch(beta0, lamb0, beta1, lamb1, N_run, N_newt, it_max):
+        dbeta = float(beta1 - beta0)
+        if abs(dbeta) < 1e-15:
+            return None
+
+        sgn = 1.0 if dbeta >= 0.0 else -1.0
+
+        def ode_beta(beta, v):
+            lamb, lamb_p, Y3, Y4 = v
+            (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
+
+            dlamb = lamb_p
+            dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0
+            dY3 = Y4
+            dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + (
+                3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1
+            ) * Y4
+            return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
+
+        def solve_newton(lamb_p0_init: float):
+            lamb_p0 = float(lamb_p0_init)
+            for _ in range(it_max):
+                v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
+                _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N_newt)
+
+                lamb_end, _, Y3_end, _ = y_end
+                delta = lamb_end - lamb1
+
+                if abs(delta) < epsilon:
+                    return True, lamb_p0
+
+                if abs(Y3_end) < 1e-20:
+                    return False, None
+
+                step = delta / Y3_end
+                step = float(np.clip(step, -1.0, 1.0))
+                lamb_p0 -= step
+
+            return False, None
+
+        seeds = [0.0, 0.25, -0.25, 1.0, -1.0]
+        best = None
+
+        for seed in seeds:
+            ok, sol = solve_newton(seed)
+            if not ok:
+                continue
+            v0_sol = np.array([lamb0, sol, 0.0, 1.0], dtype=float)
+            _, _, s_val = rk4_integral(ode_beta, beta0, v0_sol, dbeta, N_run, integrand_beta)
+            if (best is None) or (s_val < best[0]):
+                best = (float(s_val), float(sol))
+
+        if best is None:
+            return None
+
+        return best[0], best[1], sgn, dbeta, ode_beta
+
+    lamb0 = float(wrap_mpi_pi(lamb_0))
+    lamb1 = float(wrap_mpi_pi(lamb_1))
+    beta0 = float(beta_0)
+    beta1 = float(beta_1)
+
+    N_full = int(n)
+    if N_full < 2:
+        N_full = 2
+
+    if all_points:
+        N_fast = min(2000, max(400, N_full // 10))
+    else:
+        N_fast = min(1500, max(300, N_full // 12))
+
+    k0 = int(np.round((lamb0 - lamb1) / two_pi))
+    lamb_targets = []
+    for dk in (-1, 0, 1):
+        lt = lamb1 + two_pi * float(k0 + dk)
+        dl = lt - lamb0
+        if abs(dl) <= np.pi + 1e-12:
+            lamb_targets.append(float(lt))
+    if not lamb_targets:
+        lamb_targets = [float(lamb1 + two_pi * float(k0))]
+
+    best_fast = None
+
+    for lt in lamb_targets:
+        if abs(lt - lamb0) >= 1e-15:
+            res = solve_lambda_branch(beta0, lamb0, beta1, lt, N_fast, min(N_fast, 800), min(iter_max, 12))
+            if res is None:
+                continue
+            s_fast, beta_p0_fast, sgn_fast, dlamb_fast, _ = res
+            cand = ("lambda", s_fast, lt, beta_p0_fast, sgn_fast, dlamb_fast)
+        else:
+            res = solve_beta_branch(beta0, lamb0, beta1, lamb1, N_fast, min(N_fast, 800), min(iter_max, 12))
+            if res is None:
+                continue
+            s_fast, lamb_p0_fast, sgn_fast, dbeta_fast, _ = res
+            cand = ("beta", s_fast, lt, lamb_p0_fast, sgn_fast, dbeta_fast)
+
+        if (best_fast is None) or (cand[1] < best_fast[1]):
+            best_fast = cand
+
+    if best_fast is None:
+        if abs(lamb1 - lamb0) >= 1e-15:
+            best_fast = ("lambda", 0.0, lamb1, None, 1.0, float(lamb1 - lamb0))
+        else:
+            best_fast = ("beta", 0.0, lamb1, None, 1.0, float(beta1 - beta0))
+
+    if best_fast[0] == "lambda":
+        lt = float(best_fast[2])
+        dlamb = float(lt - lamb0)
+        sgn = 1.0 if dlamb >= 0.0 else -1.0
+
+        def ode_lamb(lamb, v):
+            beta, beta_p, X3, X4 = v
+            (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
+
+            dbeta = beta_p
+            dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0
+            dX3 = X4
+            dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (
+                3 * p_3 * beta_p**2 + 2 * p_2 * beta_p + p_1
+            ) * X4
+            return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
+
+        alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lt)
+        (_, _, E0, G0, *_) = BETA_LAMBDA(beta0, lamb0)
+        beta_p0_sph = np.sqrt(G0 / E0) * cot(alpha0_sph)
+
+        beta_p0_init = best_fast[3]
+        if beta_p0_init is None:
+            beta_p0_init = beta_p0_sph
+        beta_p0_init = float(beta_p0_init)
+
+        N_newton = min(N_full, 4000)
+
+        def solve_newton_refine(beta_p0_init: float):
+            beta_p0 = float(beta_p0_init)
+            for _ in range(iter_max):
+                v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
+                _, y_end = rk4_end(ode_lamb, lamb0, v0, dlamb, N_newton)
+
+                beta_end, _, X3_end, _ = y_end
+                delta = beta_end - beta1
+
+                if abs(delta) < epsilon:
+                    return True, beta_p0
+
+                if abs(X3_end) < 1e-20:
+                    return False, None
+
+                step = delta / X3_end
+                step = float(np.clip(step, -0.5, 0.5))
+                beta_p0 -= step
+
+            return False, None
+
+        ok, beta_p0_sol = solve_newton_refine(beta_p0_init)
+
+        if not ok:
+            seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
+            best = None
+            for seed in seeds:
+                ok_s, sol_s = solve_newton_refine(seed)
+                if not ok_s:
+                    continue
+                v0_s = np.array([beta0, sol_s, 0.0, 1.0], dtype=float)
+                _, _, s_s = rk4_integral(ode_lamb, lamb0, v0_s, dlamb, min(N_full, 2000), integrand_lambda)
+                if (best is None) or (s_s < best[0]):
+                    best = (float(s_s), float(sol_s))
+            if best is None:
+                raise RuntimeError("GHA2_num: Keine Startwert-Variante konvergiert (lambda-Fall)")
+            beta_p0_sol = best[1]
+
+        beta_p0 = float(beta_p0_sol)
+        v0_final = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
+
+        if all_points:
+            lamb_list, states = rk4(ode_lamb, lamb0, v0_final, dlamb, N_full, False)
+            lamb_arr = np.array(lamb_list, dtype=float)
+            beta_arr = np.array([st[0] for st in states], dtype=float)
+            beta_p_arr = np.array([st[1] for st in states], dtype=float)
+
+            (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+            (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+
+            alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p_arr[0] * sgn, dlamb_du=1.0 * sgn)
+            alpha_1 = azimut(E_end, G_end, dbeta_du=beta_p_arr[-1] * sgn, dlamb_du=1.0 * sgn)
+
+            integrand = np.zeros(N_full + 1, dtype=float)
+            for i in range(N_full + 1):
+                (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
+                integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi)
+
+            h = abs(dlamb) / N_full
+            if N_full % 2 == 0:
+                S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(
+                    integrand[2:-1:2]
+                )
+                s = h / 3.0 * S
+            else:
+                s = np.trapz(integrand, dx=h)
+
+            return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr
+
+        _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N_full, integrand_lambda)
+        beta_end, beta_p_end, _, _ = y_end
+
+        (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
+        alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p0 * sgn, dlamb_du=1.0 * sgn)
+
+        (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), float(lamb0 + dlamb))
+        alpha_1 = azimut(E_end, G_end, dbeta_du=float(beta_p_end) * sgn, dlamb_du=1.0 * sgn)
+
+        return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s)
+
+    # Fall 2 (lambda_0 == lambda_1)
+    N = int(n)
+    dbeta = float(beta_1 - beta_0)
+
+    if abs(dbeta) < 1e-15:
+        if all_points:
+            return 0.0, 0.0, 0.0, np.array([]), np.array([])
+        return 0.0, 0.0, 0.0
+
+    sgn = 1.0 if dbeta >= 0.0 else -1.0
+
+    def ode_beta(beta, v):
+        lamb, lamb_p, Y3, Y4 = v
+        (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
+
+        dlamb = lamb_p
+        dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0
+        dY3 = Y4
+        dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + (
+            3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1
+        ) * Y4
+        return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
+
+    lamb_p0 = float(best_fast[3]) if (best_fast[0] == "beta" and best_fast[3] is not None) else 0.0
+
+    for _ in range(iter_max):
+        v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
+        _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N)
+
+        lamb_end, _, Y3_end, _ = y_end
+        delta = lamb_end - lamb1
+
+        if abs(delta) < epsilon:
+            break
+
+        if abs(Y3_end) < 1e-20:
+            raise RuntimeError("GHA2_num: Ableitung ~ 0 im beta-Fall")
+
+        step = delta / Y3_end
+        step = float(np.clip(step, -1.0, 1.0))
+        lamb_p0 -= step
+
+    v0_final = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
+
+    if all_points:
+        beta_list, states = rk4(ode_beta, beta0, v0_final, dbeta, N, False)
+        beta_arr = np.array(beta_list, dtype=float)
+        lamb_arr = np.array([st[0] for st in states], dtype=float)
+        lamb_p_arr = np.array([st[1] for st in states], dtype=float)
+
+        (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+        (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+
+        alpha_0 = azimut(E_start, G_start, dbeta_du=1.0 * sgn, dlamb_du=lamb_p_arr[0] * sgn)
+        alpha_1 = azimut(E_end, G_end, dbeta_du=1.0 * sgn, dlamb_du=lamb_p_arr[-1] * sgn)
+
+        integrand = np.zeros(N + 1, dtype=float)
+        for i in range(N + 1):
+            (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
+            integrand[i] = np.sqrt(Ei + Gi * lamb_p_arr[i] ** 2)
+
+        h = abs(dbeta) / N
+        if N % 2 == 0:
+            S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(
+                integrand[2:-1:2]
+            )
+            s = h / 3.0 * S
+        else:
+            s = np.trapz(integrand, dx=h)
+
+        return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr
+
+    _, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta)
+    lamb_end, lamb_p_end, _, _ = y_end
+
+    (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
+    alpha_0 = azimut(E_start, G_start, dbeta_du=1.0 * sgn, dlamb_du=lamb_p0 * sgn)
+
+    (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, float(lamb_end))
+    alpha_1 = azimut(E_end, G_end, dbeta_du=1.0 * sgn, dlamb_du=float(lamb_p_end) * sgn)
+
+    return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s)
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    beta1 = np.deg2rad(75)
+    lamb1 = np.deg2rad(-90)
+    beta2 = np.deg2rad(75)
+    lamb2 = np.deg2rad(66)
+    a0, a1, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=100)
+    print(aus.gms("a0", a0, 4))
+    print(aus.gms("a1", a1, 4))
+    print("s: ", s)
+    # print(aus.gms("a2", a2, 4))
+    # print(s)
+    cart1 = ell.para2cart(0, 0)
+    cart2 = ell.para2cart(0.4, 1.4)
+    beta1, lamb1 = ell.cart2ell(cart1)
+    beta2, lamb2 = ell.cart2ell(cart2)
+
+    a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
+    print(s)
+
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    diffs_panou = []
+    examples_panou = ne_panou.get_random_examples(4)
+    for example in examples_panou:
+        beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example
+        P0 = ell.ell2cart(beta0, lamb0)
+        try:
+            alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
+            diffs_panou.append(
+                (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num)))
+        except:
+            print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
+    diffs_panou = np.array(diffs_panou)
+    print(diffs_panou)
+
+    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    diffs_karney = []
+    # examples_karney = ne_karney.get_examples((30500, 40500))
+    examples_karney = ne_karney.get_random_examples(2)
+    for example in examples_karney:
+        beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example
+
+        try:
+            alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
+            diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num)))
+        except:
+            print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
+    diffs_karney = np.array(diffs_karney)
+    print(diffs_karney)
+
+    pass

GHA_triaxial/numeric_examples_karney.py

@@ -0,0 +1,118 @@
+import random
+from typing import List
+
+import winkelumrechnungen as wu
+from GHA_triaxial.utils import jacobi_konstante
+from ellipsoid_triaxial import EllipsoidTriaxial
+
+ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+file_path = r"Karney_2024_Testset.txt"
+
+def line2example(line: str) -> List:
+    """
+    Line-String in Liste umwandeln
+    :param line: Line-String
+    :return: Liste mit Zahlenwerten
+    """
+    split = line.split()
+    example = [float(value) for value in split[:7]]
+    for i, value in enumerate(example):
+        if i < 6:
+            example[i] = wu.deg2rad(value)
+            # example[i] = value
+    return example
+
+def get_random_examples(num: int, seed: int = None) -> List:
+    """
+    Rückgabe zufälliger Beispiele
+    beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12
+    :param num: Anzahl zufälliger Beispiele
+    :param seed: Random-Seed
+    :return: Liste mit Beispielen
+    """
+    if seed is not None:
+        random.seed(seed)
+    with open(file_path) as datei:
+        lines = datei.readlines()
+    examples = []
+    for i in range(num):
+        example = line2example(lines[random.randint(0, len(lines) - 1)])
+        examples.append(example)
+    return examples
+
+def get_examples(l_i: List) -> List:
+    """
+    Rückgabe ausgewählter Beispiele
+    beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12
+    :param l_i: Liste von Indizes
+    :return: Liste mit Beispielen
+    """
+    with open(file_path) as datei:
+        lines = datei.readlines()
+    examples = []
+    for i in l_i:
+        example = line2example(lines[i])
+        examples.append(example)
+    return examples
+
+
+def get_random_examples_gamma(group: str, num: int, seed: int = None, length: str = None) -> List:
+    """
+    Zufällige Beispiele aus Karney in Gruppen nach Einteilung anhand der Jacobi-Konstanten
+    :param group: Gruppe
+    :param num: Anzahl
+    :param seed: Random-Seed
+    :param length: long oder short, sond egal
+    :return: Liste mit Beispielen
+    """
+    eps = 1e-20
+    long_short = 2
+    if seed is not None:
+        random.seed(seed)
+    with open(file_path) as datei:
+        lines = datei.readlines()
+    examples = []
+    i = 0
+    while len(examples) < num and i < len(lines):
+        example = line2example(lines[random.randint(0, len(lines) - 1)])
+        if example in examples:
+            continue
+        i += 1
+
+        beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
+        gamma = jacobi_konstante(beta0, lamb0, alpha0_ell, ell)
+
+        if group not in ["a", "b", "c", "d", "e", "de"]:
+            break
+        elif group == "a" and not 1 >= gamma >= 0.01:
+            continue
+        elif group == "b" and not 0.01 > gamma > eps:
+            continue
+        elif group == "c" and not abs(gamma) <= eps:
+            continue
+        elif group == "d" and not -eps > gamma > -1e-17:
+            continue
+        elif group == "e" and not -1e-17 >= gamma >= -1:
+            continue
+        elif group == "de" and not -eps > gamma > -1:
+            continue
+
+        if length == "short":
+            if example[6] < long_short:
+                examples.append(example)
+        elif length == "long":
+            if example[6] >= long_short:
+                examples.append(example)
+        else:
+            examples.append(example)
+
+    return examples
+
+
+if __name__ == "__main__":
+    examples_a = get_random_examples_gamma("a", 10, 42)
+    examples_b = get_random_examples_gamma("b", 10, 42)
+    examples_c = get_random_examples_gamma("c", 10, 42)
+    examples_d = get_random_examples_gamma("d", 10, 42)
+    examples_e = get_random_examples_gamma("e", 10, 42)
+    pass

GHA_triaxial/numeric_examples_panou.py

@@ -0,0 +1,154 @@
+from typing import Tuple, List
+
+import winkelumrechnungen as wu
+import random
+
+table1 = [
+    (wu.deg2rad(0),  wu.deg2rad(0),    wu.deg2rad(0),   wu.deg2rad(90),   1.00000000000,
+     wu.gms2rad([90, 0, 0.0000]),      wu.gms2rad([90, 0, 0.0000]),       10018754.9569),
+
+    (wu.deg2rad(1),  wu.deg2rad(0),    wu.deg2rad(-80), wu.deg2rad(5),    0.05883743460,
+     wu.gms2rad([179, 7, 12.2719]),    wu.gms2rad([174, 40, 13.8487]),    8947130.7221),
+
+    (wu.deg2rad(5), wu.deg2rad(0),     wu.deg2rad(-60), wu.deg2rad(40),    0.34128138370,
+     wu.gms2rad([160, 13, 24.5001]),   wu.gms2rad([137, 26, 47.0036]),    8004762.4330),
+
+    (wu.deg2rad(30), wu.deg2rad(0),    wu.deg2rad(-30), wu.deg2rad(175),  0.86632464962,
+     wu.gms2rad([91, 7, 30.9337]),     wu.gms2rad([91, 7, 30.8672]),      19547128.7971),
+
+    (wu.deg2rad(60), wu.deg2rad(0),    wu.deg2rad(60), wu.deg2rad(175),   0.06207487624,
+     wu.gms2rad([2, 52, 26.2393]),     wu.gms2rad([177, 4, 13.6373]),     6705715.1610),
+
+    (wu.deg2rad(75), wu.deg2rad(0),    wu.deg2rad(80), wu.deg2rad(120),   0.11708984898,
+     wu.gms2rad([23, 20, 34.7823]),    wu.gms2rad([140, 55, 32.6385]),    2482501.2608),
+
+    (wu.deg2rad(80), wu.deg2rad(0),    wu.deg2rad(60), wu.deg2rad(90),    0.17478427424,
+     wu.gms2rad([72, 26, 50.4024]),    wu.gms2rad([159, 38, 30.3547]),    3519745.1283)
+]
+
+table2 = [
+    (wu.deg2rad(0), wu.deg2rad(-90),   wu.deg2rad(0), wu.deg2rad(89.5),   1.00000000000,
+     wu.gms2rad([90, 0, 0.0000]),      wu.gms2rad([90, 0, 0.0000]),       19981849.8629),
+
+    (wu.deg2rad(1), wu.deg2rad(-90),   wu.deg2rad(1), wu.deg2rad(89.5),   0.18979826428,
+     wu.gms2rad([10, 56, 33.6952]),    wu.gms2rad([169, 3, 26.4359]),     19776667.0342),
+
+    (wu.deg2rad(5), wu.deg2rad(-90),   wu.deg2rad(5), wu.deg2rad(89),     0.09398403161,
+     wu.gms2rad([5, 24, 48.3899]),     wu.gms2rad([174, 35, 12.6880]),    18889165.0873),
+
+    (wu.deg2rad(30), wu.deg2rad(-90),  wu.deg2rad(30), wu.deg2rad(86),    0.06004022935,
+     wu.gms2rad([3, 58, 23.8038]),     wu.gms2rad([176, 2, 7.2825]),      13331814.6078),
+
+    (wu.deg2rad(60), wu.deg2rad(-90),  wu.deg2rad(60), wu.deg2rad(78),    0.06076096484,
+     wu.gms2rad([6, 56, 46.4585]),     wu.gms2rad([173, 11, 5.9592]),     6637321.6350),
+
+    (wu.deg2rad(75), wu.deg2rad(-90),  wu.deg2rad(75), wu.deg2rad(66),    0.05805851008,
+     wu.gms2rad([12, 40, 34.9009]),    wu.gms2rad([168, 20, 26.7339]),    3267941.2812),
+
+    (wu.deg2rad(80), wu.deg2rad(-90),  wu.deg2rad(80), wu.deg2rad(55),    0.05817384452,
+     wu.gms2rad([18, 35, 40.7848]),    wu.gms2rad([164, 25, 34.0017]),    2132316.9048)
+]
+
+table3 = [
+    (wu.deg2rad(0), wu.deg2rad(0.5),   wu.deg2rad(80), wu.deg2rad(0.5),   0.05680316848,
+     wu.gms2rad([0, 0, 16.0757]),      wu.gms2rad([0, 1, 32.5762]),       8831874.3717),
+
+    (wu.deg2rad(-1), wu.deg2rad(5),    wu.deg2rad(75), wu.deg2rad(5),     0.05659149555,
+     wu.gms2rad([0, -1, 47.2105]),     wu.gms2rad([0, 6, 54.0958]),       8405370.4947),
+
+    (wu.deg2rad(-5), wu.deg2rad(30),   wu.deg2rad(60), wu.deg2rad(30),    0.04921108945,
+     wu.gms2rad([0, -4, 22.3516]),     wu.gms2rad([0, 8, 42.0756]),       7204083.8568),
+
+    (wu.deg2rad(-30), wu.deg2rad(45),  wu.deg2rad(30), wu.deg2rad(45),    0.04017812574,
+     wu.gms2rad([0, -3, 41.2461]),     wu.gms2rad([0, 3, 41.2461]),       6652788.1287),
+
+    (wu.deg2rad(-60), wu.deg2rad(60),  wu.deg2rad(5), wu.deg2rad(60),     0.02843082609,
+     wu.gms2rad([0, -8, 40.4575]),     wu.gms2rad([0, 4, 22.1675]),       7213412.4477),
+
+    (wu.deg2rad(-75), wu.deg2rad(85),  wu.deg2rad(1), wu.deg2rad(85),     0.00497802414,
+     wu.gms2rad([0, -6, 44.6115]),     wu.gms2rad([0, 1, 47.0474]),       8442938.5899),
+
+    (wu.deg2rad(-80),wu.deg2rad(89.5), wu.deg2rad(0), wu.deg2rad(89.5),   0.00050178253,
+     wu.gms2rad([0, -1, 27.9705]),     wu.gms2rad([0, 0, 16.0490]),       8888783.7815)
+]
+
+# table4 = [
+#     (wu.deg2rad(0), wu.deg2rad(0),     wu.deg2rad(0), wu.deg2rad(90),     1.00000000000,
+#      wu.gms2rad([90, 0, 0.0000]),      wu.gms2rad([90, 0, 0.0000]),       10018754.1714),
+#
+#     (wu.deg2rad(1), wu.deg2rad(0),     wu.deg2rad(0), wu.deg2rad(179.5),  0.30320665822,
+#      wu.gms2rad([17, 39, 11.0942]),    wu.gms2rad([162, 20, 58.9032]),    19884417.8083),
+#
+#     (wu.deg2rad(5), wu.deg2rad(0),     wu.deg2rad(-80), wu.deg2rad(170),  0.03104258442,
+#      wu.gms2rad([178, 12, 51.5083]),   wu.gms2rad([10, 17, 52.6423]),     11652530.7514),
+#
+#     (wu.deg2rad(30), wu.deg2rad(0),    wu.deg2rad(-75), wu.deg2rad(120),  0.24135347134,
+#      wu.gms2rad([163, 49, 4.4615]),    wu.gms2rad([68, 49, 50.9617]),     14057886.8752),
+#
+#     (wu.deg2rad(60), wu.deg2rad(0),    wu.deg2rad(-60), wu.deg2rad(40),   0.19408499032,
+#      wu.gms2rad([157, 9, 33.5589]),    wu.gms2rad([157, 9, 33.5589]),     13767414.8267),
+#
+#     (wu.deg2rad(75), wu.deg2rad(0),    wu.deg2rad(-30), wu.deg2rad(0.5),  0.00202789418,
+#      wu.gms2rad([179, 33, 3.8613]),    wu.gms2rad([179, 51, 57.0077]),    11661713.4496),
+#
+#     (wu.deg2rad(80), wu.deg2rad(0),    wu.deg2rad(-5), wu.deg2rad(120),   0.15201222384,
+#      wu.gms2rad([61, 5, 33.9600]),     wu.gms2rad([171, 13, 22.0148]),    11105138.2902),
+#
+#     (wu.deg2rad(0), wu.deg2rad(0),     wu.deg2rad(60), wu.deg2rad(0),     0.00000000000,
+#      wu.gms2rad([0, 0, 0.0000]),       wu.gms2rad([0, 0, 0.0000]),     6663348.2060)
+# ]
+
+tables = [table1, table2, table3]
+
+def get_example(table: int, example: int) -> Tuple:
+    """
+    Rückgabe eines Beispiels
+    :param table: Tabellen-Nummer
+    :param example: Beispiel-Nummer
+    :return: Bespiel
+    """
+    table -= 1
+    example -= 1
+    tables = get_tables()
+    beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = tables[table][example]
+    return beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s
+
+def get_tables() -> List:
+    """
+    Rückgabe aller Tabellen
+    :return: Alle Tabellen
+    """
+    sorted_tables = []
+    for table in tables:
+        sorted_tables.append([])
+        for example in table:
+            beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example
+            sorted_tables[-1].append((beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s))
+    return sorted_tables
+
+def get_random_examples(num: int, seed: int = None) -> List:
+    """
+    Rückgabe zufälliger Beispiele
+    :param num: Anzahl Beispiele
+    :param seed: Random-Seed
+    :return:
+    """
+    if seed is not None:
+        random.seed(seed)
+
+    examples = []
+    for i in range(num):
+        table = random.randint(1, 3)
+        if table == 4:
+            example = random.randint(1, 8)
+        else:
+            example = random.randint(1, 7)
+        example = get_example(table, example)
+        examples.append(example)
+    return examples
+
+
+if __name__ == "__main__":
+    # test = get_example(1, 4)
+    examples = get_random_examples(5)
+    pass

GHA_triaxial/panou.py

@@ -1,192 +0,0 @@
-import numpy as np
-import ellipsoide
-import Numerische_Integration.num_int_runge_kutta as rk
-import winkelumrechnungen as wu
-import ausgaben as aus
-import GHA.rk as ghark
-from scipy.special import factorial as fact
-from math import comb
-
-# Panou, Korakitits 2019
-
-def p_q(ell: ellipsoide.EllipsoidTriaxial, x, y, z):
-    H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2
-
-    n = np.array([x / np.sqrt(H), y / ((1 - ell.ee ** 2) * np.sqrt(H)), z / ((1 - ell.ex ** 2) * np.sqrt(H))])
-
-    beta, lamb, u = ell.cart2ell(x, y, z)
-    carts = ell.ell2cart(beta, lamb, u)
-    B = ell.Ex ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(beta) ** 2
-    L = ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2
-
-    c1 = x ** 2 + y ** 2 + z ** 2 - (ell.ax ** 2 + ell.ay ** 2 + ell.b ** 2)
-    c0 = (ell.ax ** 2 * ell.ay ** 2 + ell.ax ** 2 * ell.b ** 2 + ell.ay ** 2 * ell.b ** 2 -
-          (ell.ay ** 2 + ell.b ** 2) * x ** 2 - (ell.ax ** 2 + ell.b ** 2) * y ** 2 - (
-                      ell.ax ** 2 + ell.ay ** 2) * z ** 2)
-    t2 = (-c1 + np.sqrt(c1**2 - 4*c0)) / 2
-    t1 = c0 / t2
-    t2e = ell.ax ** 2 * np.sin(lamb) ** 2 + ell.ay ** 2 * np.cos(lamb) ** 2
-    t1e = ell.ay ** 2 * np.sin(beta) ** 2 + ell.b ** 2 * np.cos(beta) ** 2
-
-    F = ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2
-    p1 = -np.sqrt(L / (F * t2)) * ell.ax / ell.Ex * np.sqrt(B) * np.sin(lamb)
-    p2 = np.sqrt(L / (F * t2)) * ell.ay * np.cos(beta) * np.cos(lamb)
-    p3 = 1 / np.sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * np.sin(beta) * np.sin(2 * lamb)
-    # p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * ell.ax / (ell.Ex * ell.Ee) * np.sqrt(B) * np.sqrt(t2 - ell.ay ** 2)
-    # p2 = np.sign(x) * np.sqrt(L / (F * t2)) * ell.ay / (ell.Ey * ell.Ee) * np.sqrt((ell.ay ** 2 - t1) * (ell.ax ** 2 - t2))
-    # p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * ell.b / (ell.Ex * ell.Ey) * np.sqrt(
-    #     (t1 - ell.b ** 2) * (t2 - ell.ay ** 2) * (ell.ax ** 2 - t2))
-    p = np.array([p1, p2, p3])
-    q = np.array([n[1]*p[2]-n[2]*p[1],
-                  n[2]*p[0]-n[0]*p[2],
-                  n[1]*p[1]-n[1]*p[0]])
-
-    return {"H": H, "n": n, "beta": beta, "lamb": lamb, "u": u, "B": B, "L": L, "c1": c1, "c0": c0, "t1": t1, "t2": t2,
-            "F": F, "p": p, "q": q}
-def gha1_num(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num):
-    values = p_q(ell, x, y, z)
-    H = values["H"]
-    p = values["p"]
-    q = values["q"]
-
-    dxds0 = p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0)
-    dyds0 = p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0)
-    dzds0 = p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0)
-
-    h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2
-
-    f1 = lambda s, x, dxds, y, dyds, z, dzds: dxds
-    f2 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * x
-    f3 = lambda s, x, dxds, y, dyds, z, dzds: dyds
-    f4 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * y/(1-ell.ee**2)
-    f5 = lambda s, x, dxds, y, dyds, z, dzds: dzds
-    f6 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * z/(1-ell.ex**2)
-
-    funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [0, x, dxds0, y, dyds0, z, dzds0], s, num)
-    return funktionswerte
-
-def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points):
-    constantValues = []
-    for point in points:
-        x = point[1]
-        dxds = point[2]
-        y = point[3]
-        dyds = point[4]
-        z = point[5]
-        dzds = point[6]
-
-        values = p_q(ell, x, y, z)
-        p = values["p"]
-        q = values["q"]
-        t1 = values["t1"]
-        t2 = values["t2"]
-
-        P = p[0]*dxds + p[1]*dyds + p[2]*dzds
-        Q = q[0]*dxds + q[1]*dyds + q[2]*dzds
-        alpha = np.arctan(P/Q)
-
-        c = ell.ay**2 - (t1 * np.sin(alpha)**2 + t2 * np.cos(alpha)**2)
-        constantValues.append(c)
-    pass
-
-
-def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM):
-    """
-    Panou, Korakitits 2020, 5ff.
-    :param ell:
-    :param x:
-    :param y:
-    :param z:
-    :param alpha0:
-    :param s:
-    :param maxM:
-    :return:
-    """
-    x_m = [x]
-    y_m = [y]
-    z_m = [z]
-
-    # erste Ableitungen (7-8)
-    sqrtH = np.sqrt(ell.H(x, y, z))
-    n = np.array([x / sqrtH,
-                  y / ((1-ell.ee**2) * sqrtH),
-                  z / ((1-ell.ex**2) * sqrtH)])
-    u, v = ell.cart2para(x, y, z)
-    G = np.sqrt(1 - ell.ex**2 * np.cos(u)**2 - ell.ee**2 * np.sin(u)**2 * np.sin(v)**2)
-    q = np.array([-1/G * np.sin(u) * np.cos(v),
-                  -1/G * np.sqrt(1-ell.ee**2) * np.sin(u) * np.sin(v),
-                  1/G * np.sqrt(1-ell.ex**2) * np.cos(u)])
-    p = np.array([q[1]*n[2] - q[2]*n[1],
-                  q[2]*n[0] - q[0]*n[2],
-                  q[1]*n[1] - q[1]*n[0]])
-    x_m.append(p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0))
-    y_m.append(p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0))
-    z_m.append(p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0))
-
-    # H Ableitungen (7)
-    H_ = lambda p: np.sum([comb(p, i) * (x_m[p - i] * x_m[i] +
-                                         1 / (1-ell.ee**2) ** 2 * y_m[p-i] * y_m[i] +
-                                         1 / (1-ell.ex**2) ** 2 * z_m[p-i] * z_m[i]) for i in range(0, p+1)])
-
-    # h Ableitungen (7)
-    h_ = lambda q: np.sum([comb(q, j) * (x_m[q-j+1] * x_m[j+1] +
-                                         1 / (1 - ell.ee ** 2) ** 2 * y_m[q-j+1] * y_m[j+1] +
-                                         1 / (1 - ell.ex ** 2) ** 2 * z_m[q-j+1] * z_m[j+1]) for j in range(0, q + 1)])
-
-    # h/H Ableitungen (6)
-    hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) *
-                                np.sum([H_(t+1-l)/(fact(t+1-l)*fact(l-1))*hH_t[l-1] for l in range(1, t+1)]))
-
-    # xm, ym, zm Ableitungen (6)
-    x_ = lambda m: -np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)])
-    y_ = lambda m: -1/(1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)])
-    z_ = lambda m: -1/(1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)])
-
-    hH_t = []
-    a_m = []
-    b_m = []
-    c_m = []
-    for m in range(2, maxM+1):
-        hH_t.append(hH_(m-2))
-        x_m.append(x_(m))
-        a_m.append(x_m[m] / fact(m))
-        y_m.append(y_(m))
-        b_m.append(y_m[m] / fact(m))
-        z_m.append(z_(m))
-        c_m.append(z_m[m] / fact(m))
-
-    # am, bm, cm (6)
-    x_s = 0
-    for a in a_m:
-        x_s = x_s * s + a
-    y_s = 0
-    for b in b_m:
-        y_s = y_s * s + b
-    z_s = 0
-    for c in c_m:
-        z_s = z_s * s + c
-
-    return x_s, y_s, z_s
-    pass
-
-
-if __name__ == "__main__":
-    ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978")
-    ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial")
-    re = ellipsoide.EllipsoidBiaxial.init_name("Bessel")
-    x0 = 5672455.1954766
-    y0 = 2698193.7242382686
-    z0 = 1103177.6450055107
-    alpha0 = wu.gms2rad([20, 0, 0])
-    s = 100
-    num = 100
-    werteTri = gha1_num(ellbi, x0, y0, z0, alpha0, s, num)
-    print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8))
-    print(np.sqrt((x0-werteTri[-1][1])**2+(y0-werteTri[-1][3])**2+(z0-werteTri[-1][5])**2))
-    checkLiouville(ell, werteTri)
-    werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num)
-    print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8))
-    print(np.sqrt((x0-werteBi[0])**2+(y0-werteBi[1])**2+(z0-werteBi[2])**2))
-
-    werteAna = gha1_ana(ell, x0, y0, z0, alpha0, s, 7)
-    print(aus.xyz(werteAna[0], werteAna[1], werteAna[2], 8))

GHA_triaxial/utils.py

@@ -0,0 +1,201 @@
+from __future__ import annotations
+
+from typing import Tuple
+
+import numpy as np
+from numpy import arctan2, cos, sin, sqrt
+from numpy.typing import NDArray
+
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_0_2pi
+
+
+def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float:
+    """
+    Berechnung des Azimuts an einem Punkt anhand der Ableitung zu den kartesischen Koordinaten
+    :param ell: Ellipsoid
+    :param sigma: Ableitungsvektor ver kartesischen Koordinaten
+    :param point: Punkt
+    :return: Azimuts
+    """
+    p, q = pq_ell(ell, point)
+    P = float(p @ sigma)
+    Q = float(q @ sigma)
+
+    alpha = arctan2(P, Q)
+    return wrap_0_2pi(alpha)
+
+
+def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]:
+    """
+    Umrechnung des Azimuts bezogen auf parametrische Koordinaten zu ellipsoidischen
+    :param ell: Ellipsoid
+    :param u: parametrische Breite
+    :param v: parametrische Länge
+    :param alpha_para: Azimut bezogen auf parametrische Koordinaten
+    :return: Azimut bezogen auf ellipsoidische Koordinaten
+    """
+    point = ell.para2cart(u, v)
+    beta, lamb = ell.para2ell(u, v)
+
+    p_para, q_para = pq_para(ell, point)
+    sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
+
+    p_ell, q_ell = pq_ell(ell, point)
+    alpha_ell = arctan2(p_ell @ sigma_para, q_ell @ sigma_para)
+    sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell)
+
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
+        raise Exception("alpha_para2ell: Differenz in den Richtungsableitungen")
+
+    return beta, lamb, wrap_0_2pi(alpha_ell)
+
+
+def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]:
+    """
+    Umrechnung des Azimuts bezogen auf ellipsoidische Koordinaten zu parametrischen
+    :param ell: Ellipsoid
+    :param beta: ellipsoidische Breite
+    :param lamb: ellipsoidische Länge
+    :param alpha_ell: Azimut bezogen auf ellipsoidische Koordinaten
+    :return: Azimut bezogen auf parametrische Koordinaten
+    """
+    point = ell.ell2cart(beta, lamb)
+    u, v = ell.ell2para(beta, lamb)
+
+    p_ell, q_ell = pq_ell(ell, point)
+    sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell)
+
+    p_para, q_para = pq_para(ell, point)
+    alpha_para = arctan2(p_para @ sigma_ell, q_para @ sigma_ell)
+    sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
+
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
+        raise Exception("alpha_ell2para: Differenz in den Richtungsableitungen")
+
+    return u, v, wrap_0_2pi(alpha_para)
+
+
+def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha_ell: float) -> NDArray:
+    """
+    Berechnung der Richtungsableitungen in kartesischen Koordinaten aus ellipsoidischem Azimut
+    Panou (2019) [6]
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha_ell: ellipsoidischer Azimut
+    :return: Richtungsableitungen in kartesischen Koordinaten
+    """
+    p, q = pq_ell(ell, point)
+    sigma = p * sin(alpha_ell) + q * cos(alpha_ell)
+    return sigma
+
+
+def func_sigma_para(ell: EllipsoidTriaxial, point: NDArray, alpha_para: float) -> NDArray:
+    """
+    Berechnung der Richtungsableitungen in kartesischen Koordinaten aus parametischem Azimut
+    Panou, Korakitis (2019) [6]
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha_para: parametrischer Azimut
+    :return: Richtungsableitungen in kartesischen Koordinaten
+    """
+    p, q = pq_para(ell, point)
+    sigma = p * sin(alpha_para) + q * cos(alpha_para)
+    return sigma
+
+
+def louville_constant(ell: EllipsoidTriaxial, p0: NDArray, alpha_ell: float) -> float:
+    """
+    Berechnung der Louville Konstanten
+    Panou, Korakitis (2019) [6]
+    :param ell: Ellipsoid
+    :param p0: Punkt in kartesischen Koordinaten
+    :param alpha_ell: ellipsoidischer Azimut
+    :return:
+    """
+    beta, lamb = ell.cart2ell(p0)
+    l = ell.Ey**2 * cos(beta)**2 * sin(alpha_ell)**2 - ell.Ee**2 * sin(lamb)**2 * cos(alpha_ell)**2
+    return l
+
+
+def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
+    """
+    Berechnung von p (Tangente entlang konstantem beta) und q (Tangente entlang konstantem lambda)
+    Panou, Korakitits (2019) [5f.]
+    :param ell: Ellipsoid
+    :param point: Punkt
+    :return: p und q
+    """
+    n = ell.func_n(point)
+
+    beta, lamb = ell.cart2ell(point)
+    if abs(cos(beta)) < 1e-15 and abs(np.sin(lamb)) < 1e-15:
+        if beta > 0:
+            p = np.array([0, -1, 0])
+        else:
+            p = np.array([0, 1, 0])
+    else:
+        B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2
+        L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2
+
+        _, t2 = ell.func_t12(point)
+
+        F = ell.Ey ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(lamb) ** 2
+        p1 = -sqrt(L / (F * t2)) * ell.ax / ell.Ex * sqrt(B) * sin(lamb)
+        p2 = sqrt(L / (F * t2)) * ell.ay * cos(beta) * cos(lamb)
+        p3 = 1 / sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * sin(beta) * sin(2 * lamb)
+        p = np.array([p1, p2, p3])
+        p = p / np.linalg.norm(p)
+    q = np.array([n[1] * p[2] - n[2] * p[1],
+                  n[2] * p[0] - n[0] * p[2],
+                  n[0] * p[1] - n[1] * p[0]])
+    q = q / np.linalg.norm(q)
+
+    return p, q
+
+
+def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
+    """
+    Berechnung von p (Tangente entlang konstantem u) und q (Tangente entlang konstantem v)
+    Panou, Korakitits (2020)
+    :param ell: Ellipsoid
+    :param point: Punkt
+    :return: p und q
+    """
+    n = ell.func_n(point)
+    u, v = ell.cart2para(point)
+
+    # 41-47
+    G = sqrt(1 - ell.ex ** 2 * cos(u) ** 2 - ell.ee ** 2 * sin(u) ** 2 * sin(v) ** 2)
+    q = np.array([-1 / G * sin(u) * cos(v),
+                  -1 / G * sqrt(1 - ell.ee ** 2) * sin(u) * sin(v),
+                  1 / G * sqrt(1 - ell.ex ** 2) * cos(u)])
+    p = np.array([q[1] * n[2] - q[2] * n[1],
+                  q[2] * n[0] - q[0] * n[2],
+                  q[0] * n[1] - q[1] * n[0]])
+
+    p = p / np.linalg.norm(p)
+    q = q / np.linalg.norm(q)
+
+    return p, q
+
+
+def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTriaxial) -> float:
+    """
+    Jacobi-Konstante nach Karney (2025), Gl. (14)
+    :param beta: Beta Koordinate
+    :param omega: Omega Koordinate
+    :param alpha: Azimut alpha
+    :param ell: Ellipsoid
+    :return: Jacobi-Konstante
+    """
+    gamma_jacobi = float((ell.k ** 2) * (np.cos(beta) ** 2) * (np.sin(alpha) ** 2) - (ell.k_ ** 2) * (np.sin(omega) ** 2) * (np.cos(alpha) ** 2))
+    return gamma_jacobi
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    alpha_para = 0
+    u, v = ell.ell2para(np.pi/2, 0)
+    alpha_ell = alpha_para2ell(ell, u, v, alpha_para)
+    pass

Numerische_Integration/num_int_runge_kutta.py

@@ -1,62 +0,0 @@
-def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) -> list:
-    """
-    Runge-Kutta-Verfahren für ein beliebiges DGLS
-    :param funktionen: Liste mit allen Funktionen
-    :param startwerte: Liste mit allen Startwerten der Variablen
-    :param weite: gesamte Weite über die integriert werden soll
-    :param schritte: Anzahl der Schritte über die gesamte Weite
-    :return: Liste mit Listen für alle Wertepaare
-    """
-    h = weite / schritte
-    werte = [startwerte]
-    for i in range(schritte):
-
-        zuschlaege_grob = zuschlaege(funktionen, werte[-1], h)
-        werte_grob = [werte[-1][j] if j == 0 else werte[-1][j] + zuschlaege_grob[j - 1]
-                      for j in range(len(startwerte))]
-
-        zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2)
-        werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1]
-                        for j in range(len(startwerte))]
-
-        zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2)
-        werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1]
-                        for j in range(len(startwerte))]
-
-        werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j])
-                      for j in range(len(startwerte))]
-
-        werte.append(werte_korr)
-    return werte
-
-
-def zuschlaege(funktionen: list, startwerte: list, h: float) -> list:
-    """
-    Berechnung der Zuschläge eines einzelnen Schritts
-    :param funktionen: Liste mit allen Funktionen
-    :param startwerte: Liste mit allen Startwerten der Variablen
-    :param h: Schrittweite
-    :return: Liste mit Zuschlägen für die einzelnen Variablen
-    """
-    werte = [wert for wert in startwerte]
-
-    k1 = [h * funktion(*werte) for funktion in funktionen]
-
-    werte = [startwerte[i] + (h / 2 if i == 0 else k1[i - 1] / 2)
-             for i in range(len(startwerte))]
-
-    k2 = [h * funktion(*werte) for funktion in funktionen]
-
-    werte = [startwerte[i] + (h / 2 if i == 0 else k2[i - 1] / 2)
-             for i in range(len(startwerte))]
-
-    k3 = [h * funktion(*werte) for funktion in funktionen]
-
-    werte = [startwerte[i] + (h if i == 0 else k3[i - 1])
-             for i in range(len(startwerte))]
-
-    k4 = [h * funktion(*werte) for funktion in funktionen]
-
-    k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))]
-
-    return k_

Numerische_Integration/num_int_trapezformel.py

@@ -1,46 +0,0 @@
-def f(xi, yi):
-    ys = xi + yi
-    return ys
-
-
-# Gegeben:
-x0 = 0
-y0 = 0
-h = 0.2
-xmax = 0.4
-Ey = 0.0001
-
-x = [x0]
-y = [[y0]]
-
-n = 0
-while x[-1] < xmax:
-    x.append(x[-1]+h)
-    fn = f(x[n], y[n][-1])
-    if n == 0:
-        y_neu = y[n][-1] + h * fn
-    else:
-        y_neu = y[n-1][-1] + 2*h * fn
-    y.append([y_neu])
-    dy = 1
-    while dy > Ey:
-        y_neu = y[n][-1] + h/2 * (fn + f(x[n+1], y[n+1][-1]))
-        y[-1].append(y_neu)
-        dy = abs(y[-1][-2]-y[-1][-1])
-    n += 1
-
-print(x)
-print(y)
-
-werte = []
-for i in range(len(x)):
-    werte.append((x[i], y[i][-1]))
-
-for paar in werte:
-    print(f"({round(paar[0], 5)}, {round(paar[1], 5)})")
-
-integral = 0
-for i in range(len(werte)-1):
-    integral += (werte[i+1][1]+werte[i][1])/2 * (werte[i+1][0]-werte[i][0])
-
-print(f"Integral = {round(integral, 5)}")

Numerische_Integration/rk_DGLS_1Ordnung.py

@@ -1,7 +0,0 @@
-import num_int_runge_kutta as rk
-
-f = lambda ti, xi, yi: yi - ti
-g = lambda ti, xi, yi: xi + yi
-
-funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 3)
-print(funktionswerte)

Numerische_Integration/rk_DGL_1Ordnung.py

@@ -1,6 +0,0 @@
-import num_int_runge_kutta as rk
-
-f = lambda xi, yi: xi + yi
-
-funktionswerte = rk.verfahren([f], [0, 0], 1, 5)
-print(funktionswerte)

Numerische_Integration/rk_DGL_2Ordnung.py

@@ -1,10 +0,0 @@
-import numpy as np
-import num_int_runge_kutta as rk
-
-f = lambda ti, ui, phii: -4 * np.sin(phii)
-g = lambda ti, ui, phii: ui
-
-funktionswerte = rk.verfahren([f, g], [0, 1, 0], 0.6, 3)
-
-for wert in funktionswerte:
-    print(f"t = {round(wert[0],1)}s -> phi = {round(wert[2],5)}, phip = {round(wert[1],5)}, v = {round(2.45 * wert[1],5)}")

Numerische_Integration/rk_DGL_2Ordnung_2.py

@@ -1,7 +0,0 @@
-import num_int_runge_kutta as rk
-
-f = lambda xi, yi, ui: ui
-g = lambda xi, yi, ui: 4 * ui - 4 * yi
-
-funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 2)
-print(funktionswerte)

T.py

@@ -1 +0,0 @@
-print("T")

ausgaben.py

@@ -10,7 +10,7 @@ def xyz(x: float, y: float, z: float, stellen: int) -> str:
     :param stellen: Anzahl Nachkommastellen
     :return: String zur Ausgabe der Koordinaten
     """
-    return f"""x = {(round(x,stellen))} m  y = {(round(y,stellen))} m  z = {(round(z,stellen))} m"""
+    return f"""x = {(round(x, stellen))} m  y = {(round(y, stellen))} m  z = {(round(z, stellen))} m"""
 
 
 def gms(name: str, rad: float, stellen: int) -> str:
@@ -21,5 +21,5 @@ def gms(name: str, rad: float, stellen: int) -> str:
     :param stellen: Anzahl Nachkommastellen
     :return: String zur Ausgabe des Winkels
     """
-    gms = wu.rad2gms(rad)
-    return f"{name} = {int(gms[0])}° {int(gms[1])}' {round(gms[2],stellen):.{stellen}f}''"
+    values = wu.rad2gms(rad)
+    return f"{name} = {int(values[0])}° {int(values[1])}' {round(values[2], stellen):.{stellen}f}''"

ellipsoid_triaxial.py

@@ -0,0 +1,733 @@
+import math
+from typing import Tuple
+
+import numpy as np
+from numpy import arccos, arctan, arctan2, cos, pi, sin, sqrt
+from numpy.typing import NDArray
+
+import jacobian_Ligas
+from utils_angle import wrap_mhalfpi_halfpi, wrap_mpi_pi
+
+
+class EllipsoidTriaxial:
+    """
+    Klasse für dreiachsige Ellipsoide
+    Parameter: Formparameter
+    Funktionen: Koordinatenumrechnungen
+    """
+    def __init__(self, ax: float, ay: float, b: float):
+        self.ax = ax
+        self.ay = ay
+        self.b = b
+        self.ex = sqrt((self.ax**2 - self.b**2) / self.ax**2)
+        self.ey = sqrt((self.ay**2 - self.b**2) / self.ay**2)
+        self.ee = sqrt((self.ax**2 - self.ay**2) / self.ax**2)
+        self.ex_ = sqrt((self.ax**2 - self.b**2) / self.b**2)
+        self.ey_ = sqrt((self.ay**2 - self.b**2) / self.b**2)
+        self.ee_ = sqrt((self.ax**2 - self.ay**2) / self.ay**2)
+        self.Ex = sqrt(self.ax**2 - self.b**2)
+        self.Ey = sqrt(self.ay**2 - self.b**2)
+        self.Ee = sqrt(self.ax**2 - self.ay**2)
+        nenner = sqrt(max(self.ax * self.ax - self.b * self.b, 0.0))
+        self.k = sqrt(max(self.ay * self.ay - self.b * self.b, 0.0)) / nenner
+        self.k_ = sqrt(max(self.ax * self.ax - self.ay * self.ay, 0.0)) / nenner
+        self.e = sqrt(max(self.ax * self.ax - self.b * self.b, 0.0)) / self.ay
+
+    @classmethod
+    def init_name(cls, name: str) -> EllipsoidTriaxial:
+        """
+        Mögliche Ellipsoide: BursaSima1980round, KarneyTest2024, Fiction, BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972,
+        Bursa1970
+        Panou et al (2020)
+        :param name: Name des dreiachsigen Ellipsoids
+        :return: dreiachsiger Ellipsoid
+        """
+        if name == "BursaFialova1993":
+            ax = 6378171.36
+            ay = 6378101.61
+            b = 6356751.84
+            return cls(ax, ay, b)
+        elif name == "BursaSima1980":
+            ax = 6378172
+            ay = 6378102.7
+            b = 6356752.6
+            return cls(ax, ay, b)
+        elif name == "BursaSima1980round":
+            # Panou 2013
+            ax = 6378172
+            ay = 6378103
+            b = 6356753
+            return cls(ax, ay, b)
+        elif name == "Eitschberger1978":
+            ax = 6378173.435
+            ay = 6378103.9
+            b = 6356754.4
+            return cls(ax, ay, b)
+        elif name == "Bursa1972":
+            ax = 6378173
+            ay = 6378104
+            b = 6356754
+            return cls(ax, ay, b)
+        elif name == "Bursa1970":
+            ax = 6378173
+            ay = 6378105
+            b = 6356754
+            return cls(ax, ay, b)
+        elif name == "Fiction":
+            ax = 6000000
+            ay = 4000000
+            b = 2000000
+            return cls(ax, ay, b)
+        elif name == "KarneyTest2024":
+            ax = sqrt(2)
+            ay = 1
+            b = 1 / sqrt(2)
+            return cls(ax, ay, b)
+        else:
+            raise Exception(f"EllipsoidTriaxial.init_name: Name {name} unbekannt")
+
+    def func_H(self, point: NDArray) -> float:
+        """
+        Berechnung H
+        Panou, Korakitis 2019 [43]
+        :param point: Punkt
+        :return: H
+        """
+        x, y, z = point
+        return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2
+
+    def func_n(self, point: NDArray, H: float = None) -> NDArray:
+        """
+        Berechnung normalen Vektor
+        Panou, Korakitis 2019 [9-12]
+        :param point: Punkt
+        :param H:
+        :return:
+        """
+        if H is None:
+            H = self.func_H(point)
+        sqrtH = sqrt(H)
+        x, y, z = point
+        return np.array([x / sqrtH,
+                         y / ((1 - self.ee ** 2) * sqrtH),
+                         z / ((1 - self.ex ** 2) * sqrtH)])
+
+    def func_t12(self, point: NDArray) -> Tuple[float, float]:
+        """
+        Berechnung Wurzeln
+        Panou, Korakitis 2019 [9-12]
+        :param point: Punkt
+        :return: Wurzeln t1, t2
+        """
+        x, y, z = point
+        c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
+        c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
+              (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
+                      self.ax ** 2 + self.ay ** 2) * z ** 2)
+        if c1 ** 2 - 4 * c0 < -1e-9:
+            raise Exception("t1, t2: Negativer Wurzelterm")
+        elif c1 ** 2 - 4 * c0 < 0:
+            t2 = 0
+        else:
+            t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
+        if t2 == 0:
+            t2 = 1e-18
+        t1 = c0 / t2
+        return t1, t2
+
+    def ellu2cart(self, beta: float, lamb: float, u: float) -> NDArray:
+        """
+        Panou 2014 12ff.
+        :param beta: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :param u: radiale Koordinate entlang der kleinen Halbachse
+        :return: Punkt in kartesischen Koordinaten
+        """
+        x = sqrt(u**2 + self.Ex**2) * sqrt(cos(beta)**2 + self.Ee**2/self.Ex**2 * sin(beta)**2) * cos(lamb)
+        y = sqrt(u**2 + self.Ey**2) * cos(beta) * sin(lamb)
+        z = u * sin(beta) * sqrt(1 - self.Ee**2/self.Ex**2 * cos(lamb)**2)
+
+        return np.array([x, y, z])
+
+    def cart2ellu(self, point: NDArray) -> Tuple[float, float, float]:
+        """
+        Panou 2014 15ff.
+        :param point: Punkt in kartesischen Koordinaten
+        :return: ellipsoidische Breite, ellipsoidische Länge, radiale Koordinate entlang der kleinen Halbachse
+        """
+        x, y, z = point
+        c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2
+        c1 = (self.ax**2 * self.ay**2 + self.ax**2 * self.b**2 + self.ay**2 * self.b**2 -
+              (self.ay**2+self.b**2) * x**2 - (self.ax**2 + self.b**2) * y**2 - (self.ax**2 + self.ay**2) * z**2)
+        c0 = (self.ax**2 * self.ay**2 * self.b**2 - self.ay**2 * self.b**2 * x**2 -
+              self.ax**2 * self.b**2 * y**2 - self.ax**2 * self.ay**2 * z**2)
+
+        p = (c2**2 - 3*c1) / 9
+        q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54
+        omega = arccos(q / sqrt(p**3))
+
+        s1 = 2 * sqrt(p) * cos(omega/3) - c2/3
+        s2 = 2 * sqrt(p) * cos(omega/3 - 2*pi/3) - c2/3
+        s3 = 2 * sqrt(p) * cos(omega/3 - 4*pi/3) - c2/3
+
+        beta = arctan(sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
+        if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7:
+            lamb = arctan(sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
+        else:
+            lamb = 0
+        u = sqrt(self.b**2 + s1)
+
+        return beta, lamb, u
+
+    def ell2cart(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray:
+        """
+        Panou, Korakitis 2019 2
+        :param beta: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :return: Punkt in kartesischen Koordinaten
+        """
+        beta = np.asarray(beta, dtype=float)
+        lamb = np.asarray(lamb, dtype=float)
+
+        beta, lamb = np.broadcast_arrays(beta, lamb)
+
+        beta = np.where(
+            np.isclose(np.abs(beta), pi / 2, atol=1e-15),
+            beta * 8999999999999999 / 9000000000000000,
+            beta
+        )
+        B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2
+        L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2
+
+        x = self.ax / self.Ex * sqrt(B) * cos(lamb)
+        y = self.ay * cos(beta) * sin(lamb)
+        z = self.b / self.Ex * sin(beta) * sqrt(L)
+
+        xyz = np.stack((x, y, z), axis=-1)
+
+        # Pole
+        mask_south = beta == -pi / 2
+        mask_north = beta == pi / 2
+        xyz[mask_south] = np.array([0, 0, -self.b])
+        xyz[mask_north] = np.array([0, 0, self.b])
+
+        # Äquator
+        mask_eq = beta == 0
+        xyz[mask_eq & (lamb == -pi / 2)] = np.array([0, -self.ay, 0])
+        xyz[mask_eq & (lamb == pi / 2)] = np.array([0, self.ay, 0])
+        xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0])
+        xyz[mask_eq & (lamb == pi)] = np.array([-self.ax, 0, 0])
+
+        return xyz
+
+    def ell2cart_bektas(self, beta: float | NDArray, omega: float | NDArray) -> NDArray:
+        """
+        Bektas 2015
+        :param beta: ellipsoidische Breite
+        :param omega: ellipsoidische Länge
+        :return: Punkt in kartesischen Koordinaten
+        """
+        x = self.ax * cos(omega) * sqrt((self.ax**2 - self.ay**2 * sin(beta)**2 - self.b**2 * cos(beta)**2) / (self.ax**2 - self.b**2))
+        y = self.ay * cos(beta) * sin(omega)
+        z = self.b * sin(beta) * sqrt((self.ax**2 * sin(omega)**2 + self.ay**2 * cos(omega)**2 - self.b**2) / (self.ax**2 - self.b**2))
+
+        return np.array([x, y, z])
+
+    def ell2cart_karney(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray:
+        """
+        Karney 2025 Geographic Lib
+        :param beta: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :return: Punkt in kartesischen Koordinaten
+        """
+        k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2)
+        k_ = sqrt(self.ax**2 - self.ay**2) / sqrt(self.ax**2 - self.b**2)
+        X = self.ax * cos(lamb) * sqrt(k**2*cos(beta)**2+k_**2)
+        Y = self.ay * cos(beta) * sin(lamb)
+        Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2)
+        return np.array([X, Y, Z])
+
+    def cart2ell_yFake(self, point: NDArray, delta_y: float = 1e-4) -> Tuple[float, float]:
+        """
+        Bei Fehlschlagen von cart2ell
+        :param point: Punkt in kartesischen Koordinaten
+        :param delta_y: Startwert für Suche nach kleinstmöglichem delta_y
+        :return: ellipsoidische Breite und Länge
+        """
+        x, y, z = point
+        best_delta = np.inf
+        while True:
+            try:
+                y1 = y - delta_y
+                beta1, lamb1 = self.cart2ell(np.array([x, y1, z]), noFake=True)
+                point1 = self.ell2cart(beta1, lamb1)
+
+                y2 = y + delta_y
+                beta2, lamb2 = self.cart2ell(np.array([x, y2, z]), noFake=True)
+                point2 = self.ell2cart(beta2, lamb2)
+
+                pointM = (point1 + point2) / 2
+
+                actual_delta = np.linalg.norm(point - pointM)
+            except:
+                actual_delta = np.inf
+
+            if actual_delta < best_delta:
+                best_delta = actual_delta
+                delta_y /= 10
+            else:
+                delta_y *= 10
+
+                y1 = y - delta_y
+                beta1, lamb1 = self.cart2ell(np.array([x, y1, z]), noFake=True)
+
+                return beta1, lamb1
+
+    def cart2ell(self, point: NDArray, eps: float = 1e-12, maxI: int = 100, noFake: bool = False) -> Tuple[float, float]:
+        """
+        Panou, Korakitis 2019 3f. (num)
+        :param point:  Punkt in kartesischen Koordinaten
+        :param eps: zu erreichende Genauigkeit
+        :param maxI: maximale Anzahl Iterationen
+        :param noFake: y numerisch anpassen?
+        :return: ellipsoidische Breite und Länge
+        """
+        x, y, z = point
+        beta, lamb = self.cart2ell_panou(point)
+        delta_ell = np.array([np.inf, np.inf]).T
+        tiny = 1e-30
+
+        try:
+            i = 0
+            while np.linalg.norm(delta_ell) > eps and i < maxI:
+                x0, y0, z0 = self.ell2cart(beta, lamb)
+                delta_l = np.array([x-x0, y-y0, z-z0]).T
+
+                B = max(self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2, tiny)
+                L = max(self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2, tiny)
+
+                J = np.array([[(-self.ax * self.Ey ** 2) / (2 * self.Ex) * sin(2 * beta) / sqrt(B) * cos(lamb),
+                               -self.ax / self.Ex * sqrt(B) * sin(lamb)],
+                              [-self.ay * sin(beta) * sin(lamb),
+                               self.ay * cos(beta) * cos(lamb)],
+                              [self.b / self.Ex * cos(beta) * sqrt(L),
+                               (self.b * self.Ee ** 2) / (2 * self.Ex) * sin(beta) * sin(2 * lamb) / sqrt(L)]])
+
+                N = J.T @ J
+                det = N[0, 0] * N[1, 1] - N[0, 1] * N[1, 0]
+                N_inv = 1 / det * np.array([[N[1, 1], -N[0, 1]], [-N[1, 0], N[0, 0]]])
+                delta_ell = N_inv @ J.T @ delta_l
+                # delta_ell, *_ = np.linalg.lstsq(J, delta_l, rcond=None)
+
+                beta += delta_ell[0]
+                lamb += delta_ell[1]
+                i += 1
+
+            if i == maxI:
+                raise Exception("Umrechnung cart2ell: nicht konvergiert")
+
+            point_n = self.ell2cart(beta, lamb)
+            delta_r = np.linalg.norm(point - point_n, axis=-1)
+            if delta_r > 1e-6:
+                raise Exception("Umrechnung cart2ell: Punktdifferenz")
+
+            return wrap_mhalfpi_halfpi(beta), wrap_mpi_pi(lamb)
+
+        except Exception as e:
+            # Wenn die Berechnung fehlschlägt auf Grund von sehr kleinem y, solange anpassen, bis Umrechnung ohne Fehler
+            delta_y = 10 ** math.floor(math.log10(abs(self.ay/1000)))
+            if abs(y) < delta_y and not noFake:
+                return self.cart2ell_yFake(point, delta_y)
+            else:
+                raise e
+
+    def cart2ell_panou(self, point: NDArray) -> Tuple[float, float]:
+        """
+        Panou, Korakitis 2019 2f. (analytisch -> Näherung)
+        :param point: Punkt in kartesischen Koordinaten
+        :return: ellipsoidische Breite und Länge
+        """
+        x, y, z = point
+
+        eps = 1e-9
+
+        if abs(x) < eps and abs(y) < eps:  # Punkt in der z-Achse
+            beta = pi / 2 if z > 0 else -pi / 2
+            lamb = 0.0
+            return beta, lamb
+
+        elif abs(x) < eps and abs(z) < eps:  # Punkt in der y-Achse
+            beta = 0.0
+            lamb = pi / 2 if y > 0 else -pi / 2
+            return beta, lamb
+
+        elif abs(y) < eps and abs(z) < eps:  # Punkt in der x-Achse
+            beta = 0.0
+            lamb = 0.0 if x > 0 else pi
+            return beta, lamb
+
+        # ---- Allgemeiner Fall -----
+
+        t1, t2 = self.func_t12(point)
+
+        num_beta = max(t1 - self.b ** 2, 0)
+        den_beta = max(self.ay ** 2 - t1, 1e-30)
+        num_lamb = max(t2 - self.ay ** 2, 0)
+        den_lamb = max(self.ax ** 2 - t2, 1e-30)
+
+        beta = arctan(sqrt(num_beta / den_beta))
+        lamb = arctan(sqrt(num_lamb / den_lamb))
+
+        if z < 0:
+            beta = -beta
+
+        if y < 0:
+            lamb = -lamb
+
+        if x < 0:
+            lamb = pi - lamb
+
+        if abs(x) < eps:
+            lamb = -pi/2 if y < 0 else pi/2
+
+        elif abs(y) < eps:
+            lamb = 0 if x > 0 else pi
+
+        elif abs(z) < eps:
+            beta = 0
+
+        return beta, lamb
+
+    def cart2ell_bektas(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]:
+        """
+        Bektas 2015
+        :param point:  Punkt in kartesischen Koordinaten
+        :param eps: zu erreichende Genauigkeit
+        :param maxI: maximale Anzahl Iterationen
+        :return: ellipsoidische Breite und Länge
+        """
+        x, y, z = point
+        phi, lamb = self.cart2para(point)
+        p = sqrt((self.ax**2 - self.ay**2) / (self.ax**2 - self.b**2))
+        d_phi = np.inf
+        d_lamb = np.inf
+
+        i = 0
+        while d_phi > eps and d_lamb > eps and i < maxI:
+            lamb_new = arctan2(self.ax * y * sqrt((p**2-1) * sin(phi)**2 + 1), self.ay * x * cos(phi))
+            phi_new = arctan2(self.ay * z * sin(lamb), self.b * y * sqrt(1 - p**2 * cos(lamb)**2))
+            d_phi = abs(phi_new - phi)
+            phi = phi_new
+            d_lamb = abs(lamb_new - lamb)
+            lamb = lamb_new
+            i += 1
+
+        if i == maxI:
+            raise Exception("Umrechnung cart2ell: nicht konvergiert")
+
+        return phi, lamb
+
+    def geod2cart(self, phi: float | NDArray, lamb: float | NDArray, h: float) -> NDArray:
+        """
+        Ligas 2012, 250
+        :param phi: geodätische Breite
+        :param lamb: geodätische Länge
+        :param h: geodätische Höhe
+        :return: kartesische Koordinaten
+        """
+        v = self.ax / sqrt(1 - self.ex**2*sin(phi)**2-self.ee**2*cos(phi)**2*sin(lamb)**2)
+        xG = (v + h) * cos(phi) * cos(lamb)
+        yG = (v * (1-self.ee**2) + h) * cos(phi) * sin(lamb)
+        zG = (v * (1-self.ex**2) + h) * sin(phi)
+        return np.array([xG, yG, zG])
+
+    def cart2geod(self, point: NDArray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]:
+        """
+        Ligas 2012
+        :param mode: ligas1, ligas2, oder ligas3
+        :param point: Punkt in kartesischen Koordinaten
+        :param maxIter: maximale Anzahl Iterationen
+        :param maxLoa: Level of Accuracy, das erreicht werden soll
+        :return: geodätische Breite, Länge, Höhe
+        """
+        xG, yG, zG = point
+
+        eps = 1e-9
+
+        if abs(xG) < eps and abs(yG) < eps:  # Punkt in der z-Achse
+            phi = pi / 2 if zG > 0 else -pi / 2
+            lamb = 0.0
+            h = abs(zG) - self.b
+            return phi, lamb, h
+
+        elif abs(xG) < eps and abs(zG) < eps:  # Punkt in der y-Achse
+            phi = 0.0
+            lamb = pi / 2 if yG > 0 else -pi / 2
+            h = abs(yG) - self.ay
+            return phi, lamb, h
+
+        elif abs(yG) < eps and abs(zG) < eps:  # Punkt in der x-Achse
+            phi = 0.0
+            lamb = 0.0 if xG > 0 else pi
+            h = abs(xG) - self.ax
+            return phi, lamb, h
+
+        rG = sqrt(xG ** 2 + yG ** 2 + zG ** 2)
+        pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
+
+        E = 1 / self.ax**2
+        F = 1 / self.ay**2
+        G = 1 / self.b**2
+
+        i = 0
+        loa = np.inf
+
+        while i < maxIter and loa > maxLoa:
+            if mode == "ligas1":
+                invJ, fxE = jacobian_Ligas.case1(E, F, G, np.array([xG, yG, zG]), pE)
+            elif mode == "ligas2":
+                invJ, fxE = jacobian_Ligas.case2(E, F, G, np.array([xG, yG, zG]), pE)
+            elif mode == "ligas3":
+                invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE)
+            else:
+                raise Exception(f"cart2geod: Modus {mode} nicht bekannt")
+            pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1)
+            pEi = pEi.reshape(1, -1).flatten()
+            loa = sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
+            pE = pEi
+            i += 1
+
+        if i == maxIter and loa > maxLoa:
+            act_mode = int(mode[-1])
+            new_mode = 3 if act_mode == 1 else act_mode - 1
+            phi, lamb, h = self.cart2geod(point, f"ligas{new_mode}", maxIter, maxLoa)
+
+        else:
+            phi = arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2))
+            lamb = arctan(1/(1-self.ee**2) * pE[1]/pE[0])
+            h = np.sign(zG - pE[2]) * np.sign(pE[2]) * sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2)
+            if h < -self.ax:
+                act_mode = int(mode[-1])
+                new_mode = 3 if act_mode == 1 else act_mode - 1
+                phi, lamb, h = self.cart2geod(point, f"ligas{new_mode}", maxIter, maxLoa)
+            else:
+                if xG < 0 and yG < 0:
+                    lamb += -pi
+
+                elif xG < 0:
+                    lamb += pi
+
+                if abs(zG) < eps:
+                    phi = 0
+        wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb)
+        return wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb), h
+
+    def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray:
+        """
+        Panou, Korakitits 2020, 4
+        :param u: parametrische Breite
+        :param v: parametrische Länge
+        :return: Punkt in kartesischen Koordinaten
+        """
+        x = self.ax * cos(u) * cos(v)
+        y = self.ay * cos(u) * sin(v)
+        z = self.b * sin(u)
+        z = np.broadcast_to(z, np.shape(x))
+        return np.array([x, y, z])
+
+    def cart2para(self, point: NDArray) -> Tuple[float, float]:
+        """
+        Panou, Korakitits 2020, 4
+        :param point: Punkt in kartesischen Koordinaten
+        :return: parametrische Breite, Länge
+        """
+        x, y, z = point
+
+        u_check1 = z*sqrt(1 - self.ee**2)
+        u_check2 = sqrt(x**2 * (1-self.ee**2) + y**2) * sqrt(1-self.ex**2)
+        if u_check1 <= u_check2:
+            u = arctan2(u_check1, u_check2)
+        else:
+            u = pi/2 - arctan2(u_check2, u_check1)
+
+        v_check1 = y
+        v_check2 = x*sqrt(1-self.ee**2)
+        v_factor = sqrt(x**2*(1-self.ee**2)+y**2)
+        if v_check1 <= v_check2:
+            v = 2 * arctan2(v_check1, v_check2 + v_factor)
+        else:
+            v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor)
+        wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v)
+        return wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v)
+
+    def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]:
+        """
+        Umrechung von ellipsoidischen in parametrische Koordinaten (über kartesische Koordinaten)
+        :param beta: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :return: parametrische Breite, Länge
+        """
+        cart = self.ell2cart(beta, lamb)
+        return self.cart2para(cart)
+
+    def para2ell(self, u: float, v: float) -> Tuple[float, float]:
+        """
+        Umrechung von parametrischen in ellipsoidische Koordinaten (über kartesische Koordinaten)
+        :param u: parametrische Breite
+        :param v: parametrische Länge
+        :return: ellipsoidische Breite, Länge
+        """
+        cart = self.para2cart(u, v)
+        return self.cart2ell(cart)
+
+    def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]:
+        """
+        Umrechung von parametrischen in geodätische Koordinaten (über kartesische Koordinaten)
+        :param u: parametrische Breite
+        :param v: parametrische Länge
+        :param mode: ligas1, ligas2, oder ligas3
+        :param maxIter: maximale Anzahl Iterationen
+        :param maxLoa: Level of Accuracy, das erreicht werden soll
+        :return: geodätische Breite, Länge, Höhe
+        """
+        cart = self.para2cart(u, v)
+        return self.cart2geod(cart, mode, maxIter, maxLoa)
+
+    def geod2para(self, phi: float, lamb: float, h: float) -> Tuple[float, float]:
+        """
+        Umrechung von geodätischen in parametrische Koordinaten (über kartesische Koordinaten)
+        :param phi: geodätische Breite
+        :param lamb: geodätische Länge
+        :param h: geodätische Höhe
+        :return: parametrische Breite, Länge
+        """
+        cart = self.geod2cart(phi, lamb, h)
+        return self.cart2para(cart)
+
+    def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]:
+        """
+        Umrechung von ellipsoidischen in geodätische Koordinaten (über kartesische Koordinaten)
+        :param beta: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :param mode: ligas1, ligas2, oder ligas3
+        :param maxIter: maximale Anzahl Iterationen
+        :param maxLoa: Level of Accuracy, das erreicht werden soll
+        :return: geodätische Breite, Länge, Höhe
+        """
+        cart = self.ell2cart(beta, lamb)
+        return self.cart2geod(cart, mode, maxIter, maxLoa)
+
+    def geod2ell(self, phi: float, lamb: float, h: float) -> Tuple[float, float]:
+        """
+        Umrechung von geodätischen in ellipsoidische Koordinaten (über kartesische Koordinaten)
+        :param phi: geodätische Breite
+        :param lamb: geodätische Länge
+        :param h: geodätische Höhe
+        :return: ellipsoidische Breite, Länge
+        """
+        cart = self.geod2cart(phi, lamb, h)
+        return self.cart2ell(cart)
+
+    def point_on(self, point: NDArray) -> bool:
+        """
+        Test, ob ein Punkt auf dem Ellipsoid liegt
+        :param point: kartesische 3D-Koordinaten
+        :return: Punkt auf dem Ellispoid?
+        """
+        value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2
+        if abs(1-value) < 1e-6:
+            return True
+        else:
+            return False
+
+    def point_onto_ellipsoid(self, point: NDArray) -> NDArray:
+        """
+        Berechnung des Lotpunktes entlang der Normalkrümmung auf einem Ellipsoiden
+        :param point: Punkt in kartesischen Koordinaten, der gelotet werden soll
+        :return: Lotpunkt in kartesischen Koordinaten
+        """
+        phi, lamb, h = self.cart2geod(point, "ligas3")
+        p = self. geod2cart(phi, lamb, 0)
+        return p
+
+    def cart_ellh(self, point: NDArray, h: float) -> NDArray:
+        """
+        Punkt auf Ellipsoid hoch loten
+        :param point: Punkt auf dem Ellipsoid
+        :param h: Höhe über dem Ellipsoid
+        :return: hochgeloteter Punkt
+        """
+        phi, lamb, _ = self.cart2geod(point, "ligas3")
+        pointH = self. geod2cart(phi, lamb, h)
+        return pointH
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    # cart = ell.ell2cart(pi/2, 0)
+    # print(cart)
+    # cart = ell.ell2cart(pi/2*8999999999999999/9000000000000000, 0)
+    # print(cart)
+    elli = ell.cart2ell(np.array([0, 0.0, 1/sqrt(2)]))
+    print(elli)
+
+    # ell = EllipsoidTriaxial.init_name("BursaSima1980")
+    # diff_list = []
+    # diffs_para = []
+    # diffs_ell = []
+    # diffs_geod = []
+    # points = []
+    # for v_deg in range(-180, 181, 5):
+    #     for u_deg in range(-90, 91, 5):
+    #         v = wu.deg2rad(v_deg)
+    #         u = wu.deg2rad(u_deg)
+    #         point = ell.para2cart(u, v)
+    #         points.append(point)
+    #
+    #         elli = ell.cart2ell(point)
+    #         cart_elli = ell.ell2cart(elli[0], elli[1])
+    #         diff_ell = np.linalg.norm(point - cart_elli, axis=-1)
+    #
+    #         para = ell.cart2para(point)
+    #         cart_para = ell.para2cart(para[0], para[1])
+    #         diff_para = np.linalg.norm(point - cart_para, axis=-1)
+    #
+    #         geod = ell.cart2geod(point, "ligas3")
+    #         cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
+    #         diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1)
+    #
+    #         diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3])
+    #         diffs_ell.append([diff_ell])
+    #         diffs_para.append([diff_para])
+    #         diffs_geod.append([diff_geod3])
+    #
+    # diff_list = np.array(diff_list)
+    # diffs_ell = np.array(diffs_ell)
+    # diffs_para = np.array(diffs_para)
+    # diffs_geod = np.array(diffs_geod)
+    #
+    # pass
+    #
+    # points = np.array(points)
+    # fig = plt.figure()
+    # ax = fig.add_subplot(projection='3d')
+    #
+    # sc = ax.scatter(
+    #     points[:, 0],
+    #     points[:, 1],
+    #     points[:, 2],
+    #     c=diffs_ell,  # Farbcode = diff
+    #     cmap='viridis',  # Colormap
+    #     s=10 + 20 * diffs_ell,  # optional: Größe abhängig vom diff
+    #     alpha=0.8
+    # )
+    #
+    # # Farbskala
+    # cbar = plt.colorbar(sc)
+    # cbar.set_label("diff")
+    #
+    # ax.set_xlabel("X")
+    # ax.set_ylabel("Y")
+    # ax.set_zlabel("Z")
+    #
+    # plt.show()

ellipsoide.py

@@ -1,315 +0,0 @@
-import numpy as np
-import winkelumrechnungen as wu
-import ausgaben as aus
-import jacobian_Ligas
-
-
-class EllipsoidBiaxial:
-    def __init__(self, a: float, b: float):
-        self.a = a
-        self.b = b
-        self.c = a ** 2 / b
-        self.e = np.sqrt(a ** 2 - b ** 2) / a
-        self.e_ = np.sqrt(a ** 2 - b ** 2) / b
-
-    @classmethod
-    def init_name(cls, name: str):
-        if name == "Bessel":
-            a = 6377397.15508
-            b = 6356078.96290
-            return cls(a, b)
-        elif name == "Hayford":
-            a = 6378388
-            f = 1/297
-            b = a - a * f
-            return cls(a, b)
-        elif name == "Krassowski":
-            a = 6378245
-            f = 298.3
-            b = a - a * f
-            return cls(a, b)
-        elif name == "WGS84":
-            a = 6378137
-            f = 298.257223563
-            b = a - a * f
-            return cls(a, b)
-
-    @classmethod
-    def init_af(cls, a: float, f: float):
-        b = a - a * f
-        return cls(a, b)
-
-    V = lambda self, phi: np.sqrt(1 + self.e_ ** 2 * np.cos(phi) ** 2)
-    M = lambda self, phi: self.c / self.V(phi) ** 3
-    N = lambda self, phi: self.c / self.V(phi)
-
-    beta2psi = lambda self, beta: np.arctan(self.a / self.b * np.tan(beta))
-    beta2phi = lambda self, beta: np.arctan(self.a ** 2 / self.b ** 2 * np.tan(beta))
-
-    psi2beta = lambda self, psi: np.arctan(self.b / self.a * np.tan(psi))
-    psi2phi = lambda self, psi: np.arctan(self.a / self.b * np.tan(psi))
-
-    phi2beta = lambda self, phi: np.arctan(self.b ** 2 / self.a ** 2 * np.tan(phi))
-    phi2psi = lambda self, phi: np.arctan(self.b / self.a * np.tan(phi))
-
-    phi2p = lambda self, phi: self.N(phi) * np.cos(phi)
-
-    def cart2ell(self, Eh, Ephi, x, y, z):
-        p = np.sqrt(x**2+y**2)
-        # print(f"p = {round(p, 5)} m")
-
-        lamb = np.arctan(y/x)
-
-        phi_null = np.arctan(z/p*(1-self.e**2)**-1)
-
-        hi = [0]
-        phii = [phi_null]
-
-        i = 0
-
-        while True:
-            N = self.a*(1-self.e**2*np.sin(phii[i])**2)**(-1/2)
-            h = p/np.cos(phii[i])-N
-            phi = np.arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1))
-            hi.append(h)
-            phii.append(phi)
-            dh = abs(hi[i]-h)
-            dphi = abs(phii[i]-phi)
-            i = i+1
-            if dh < Eh:
-                if dphi < Ephi:
-                    break
-        for i in range(len(phii)):
-            # print(f"P3[{i}]: {aus.gms('phi', phii[i], 5)}\th = {round(hi[i], 5)} m")
-            pass
-        return phi, lamb, h
-
-    def ell2cart(self, phi, lamb, h):
-        W = np.sqrt(1 - self.e**2 * np.sin(phi)**2)
-        N = self.a / W
-        x = (N+h) * np.cos(phi) * np.cos(lamb)
-        y = (N+h) * np.cos(phi) * np.sin(lamb)
-        z = (N * (1-self.e**2) + h) * np.sin(lamb)
-        return x, y, z
-
-class EllipsoidTriaxial:
-    def __init__(self, ax: float, ay: float, b: float):
-        self.ax = ax
-        self.ay = ay
-        self.b = b
-        self.ex = np.sqrt((self.ax**2 - self.b**2) / self.ax**2)
-        self.ey = np.sqrt((self.ay**2 - self.b**2) / self.ay**2)
-        self.ee = np.sqrt((self.ax**2 - self.ay**2) / self.ax**2)
-        self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2)
-        self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2)
-        self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2)
-        self.Ex = np.sqrt(self.ax**2 - self.b**2)
-        self.Ey = np.sqrt(self.ay**2 - self.b**2)
-        self.Ee = np.sqrt(self.ax**2 - self.ay**2)
-
-    @classmethod
-    def init_name(cls, name: str):
-        if name == "BursaFialova1993":
-            ax = 6378171.36
-            ay = 6378101.61
-            b = 6356751.84
-            return cls(ax, ay, b)
-        elif name == "BursaSima1980":
-            ax = 6378172
-            ay = 6378102.7
-            b = 6356752.6
-            return cls(ax, ay, b)
-        elif name == "Eitschberger1978":
-            ax = 6378173.435
-            ay = 6378103.9
-            b = 6356754.4
-            return cls(ax, ay, b)
-        elif name == "Bursa1972":
-            ax = 6378173
-            ay = 6378104
-            b = 6356754
-            return cls(ax, ay, b)
-        elif name == "Bursa1970":
-            ax = 6378173
-            ay = 6378105
-            b = 6356754
-            return cls(ax, ay, b)
-        elif name == "Bessel-biaxial":
-            ax = 6377397.15509
-            ay = 6377397.15508
-            b = 6356078.96290
-            return cls(ax, ay, b)
-
-    def ell2cart(self, beta, lamb, u):
-        """
-        Panou 2014 12ff.
-        :param beta: ellipsoidische Breite
-        :param lamb: ellipsoidische Länge
-        :param u: Höhe
-        :return: kartesische Koordinaten
-        """
-        s1 = u**2 - self.b**2
-        s2 = -self.ay**2 * np.sin(beta)**2 - self.b**2 * np.cos(beta)**2
-        s3 = -self.ax**2 * np.sin(lamb)**2 - self.ay**2 * np.cos(lamb)**2
-        # print(s1, s2, s3)
-        xe = np.sqrt(((self.ax**2+s1) * (self.ax**2+s2) * (self.ax**2+s3)) /
-                     ((self.ax**2-self.ay**2) * (self.ax**2-self.b**2)))
-        ye = np.sqrt(((self.ay**2+s1) * (self.ay**2+s2) * (self.ay**2+s3)) /
-                     ((self.ay**2-self.ax**2) * (self.ay**2-self.b**2)))
-        ze = np.sqrt(((self.b**2+s1) * (self.b**2+s2) * (self.b**2+s3)) /
-                     ((self.b**2-self.ax**2) * (self.b**2-self.ax**2)))
-
-        x = np.sqrt(u**2 + self.Ex**2) * np.sqrt(np.cos(beta)**2 + self.Ee**2/self.Ex**2 * np.sin(beta)**2) * np.cos(lamb)
-        y = np.sqrt(u**2 + self.Ey**2) * np.cos(beta) * np.sin(lamb)
-        z = u * np.sin(beta) * np.sqrt(1 - self.Ee**2/self.Ex**2 * np.cos(lamb)**2)
-
-        return x, y, z
-
-    def cart2ell(self, x, y, z):
-        """
-        Panou 2014 15ff.
-        :param x:
-        :param y:
-        :param z:
-        :return:
-        """
-        c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2
-        c1 = (self.ax**2 * self.ay**2 + self.ax**2 * self.b**2 + self.ay**2 * self.b**2 -
-              (self.ay**2+self.b**2) * x**2 - (self.ax**2 + self.b**2) * y**2 - (self.ax**2 + self.ay**2) * z**2)
-        c0 = (self.ax**2 * self.ay**2 * self.b**2 - self.ay**2 * self.b**2 * x**2 -
-              self.ax**2 * self.b**2 * y**2 - self.ax**2 * self.ay**2 * z**2)
-
-        p = (c2**2 - 3*c1) / 9
-        q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54
-        omega = np.arccos(q / np.sqrt(p**3))
-
-        s1 = 2 * np.sqrt(p) * np.cos(omega/3) - c2/3
-        s2 = 2 * np.sqrt(p) * np.cos(omega/3 - 2*np.pi/3) - c2/3
-        s3 = 2 * np.sqrt(p) * np.cos(omega/3 - 4*np.pi/3) - c2/3
-        # print(s1, s2, s3)
-
-        beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
-        lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
-        u = np.sqrt(self.b**2 + s1)
-
-        return beta, lamb, u
-
-    def cart2geod(self, mode: str, xG, yG, zG, maxIter=30, maxLoa=0.005):
-        """
-        Ligas 2012
-        :param mode:
-        :param xG:
-        :param yG:
-        :param zG:
-        :param maxIter:
-        :param maxLoa:
-        :return:
-        """
-        rG = np.sqrt(xG**2 + yG**2 + zG**2)
-        pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
-
-        E = 1 / self.ax**2
-        F = 1 / self.ay**2
-        G = 1 / self.b**2
-
-        i = 0
-        loa = np.inf
-
-        while i < maxIter and loa > maxLoa:
-            if mode == "ligas1":
-                invJ, fxE = jacobian_Ligas.case1(E, F, G, np.array([xG, yG, zG]), pE)
-            elif mode == "ligas2":
-                invJ, fxE = jacobian_Ligas.case2(E, F, G, np.array([xG, yG, zG]), pE)
-            elif mode == "ligas3":
-                invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE)
-            pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1)
-            pEi = pEi.reshape(1, -1).flatten()
-            loa = np.sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
-            pE = pEi
-            i += 1
-
-        phi = np.arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / np.sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2))
-        lamb = np.arctan(1/(1-self.ee**2) * pE[1]/pE[0])
-        h = np.sign(zG-pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0]-xG)**2 + (pE[1]-yG)**2 + (pE[2]-zG)**2)
-
-        return phi, lamb, h
-
-    def geod2cart(self, phi, lamb, h):
-        """
-        Ligas 2012, 250
-        :param phi:
-        :param lamb:
-        :param h:
-        :return:
-        """
-        v = self.ax / np.sqrt(1 - self.ex**2*np.sin(phi)**2-self.ee**2*np.cos(phi)**2*np.sin(lamb)**2)
-        xG = (v + h) * np.cos(phi) * np.cos(lamb)
-        yG = (v * (1-self.ee**2) + h) * np.cos(phi) * np.sin(lamb)
-        zG = (v * (1-self.ex**2) + h) * np.sin(phi)
-        return xG, yG, zG
-
-    def para2cart(self, u, v):
-        """
-        Panou, Korakitits 2020, 4
-        :param u:
-        :param v:
-        :return:
-        """
-        x = self.ax * np.cos(u) * np.cos(v)
-        y = self.ay * np.cos(u) * np.cos(v)
-        z = self.b * np.sin(u)
-
-    def cart2para(self, x, y, z):
-        """
-        Panou, Korakitits 2020, 4
-        :param x:
-        :param y:
-        :param z:
-        :return:
-        """
-        u_check1 = z*np.sqrt(1 - self.ee**2)
-        u_check2 = np.sqrt(x**2 * (1-self.ee**2) + y**2) * np.sqrt(1-self.ex**2)
-        if u_check1 <= u_check2:
-            u = np.arctan(u_check1 / u_check2)
-        else:
-            u = np.pi/2 * np.arctan(u_check2 / u_check1)
-
-        v_check1 = y
-        v_check2 = x*np.sqrt(1-self.ee**2)
-        if v_check1 <= v_check2:
-            v = 2 * np.arctan(v_check1 / (v_check2 + np.sqrt(x**2*(1-self.ee**2)+y**2)))
-        else:
-            v = np.pi/2 - 2 * np.arctan(v_check2 / (v_check1 + np.sqrt(x**2*(1-self.ee**2)+y**2)))
-        return u, v
-
-    def H(self, x, y, z):
-        return x**2 + y**2/(1-self.ee**2)**2 + z**2/(1-self.ex**2)**2
-
-
-if __name__ == "__main__":
-    ellips = EllipsoidTriaxial.init_name("Eitschberger1978")
-
-    carts = ellips.ell2cart(wu.deg2rad(10), wu.deg2rad(30), 6378172)
-    ells = ellips.cart2ell(carts[0], carts[1], carts[2])
-    print(aus.gms("beta", ells[0], 3), aus.gms("lambda", ells[1], 3), "u =", ells[2])
-    ells2 = ellips.cart2ell(5712200, 2663400, 1106000)
-    carts2 = ellips.ell2cart(ells2[0], ells2[1], ells2[2])
-    print(aus.xyz(carts2[0], carts2[1], carts2[2], 10))
-
-    # stellen = 20
-    # geod1 = ellips.cart2geod("ligas1", 5712200, 2663400, 1106000)
-    # print(aus.gms("phi", geod1[0], stellen), aus.gms("lambda", geod1[1], stellen), "h =", geod1[2])
-    # geod2 = ellips.cart2geod("ligas2", 5712200, 2663400, 1106000)
-    # print(aus.gms("phi", geod2[0], stellen), aus.gms("lambda", geod2[1], stellen), "h =", geod2[2])
-    # geod3 = ellips.cart2geod("ligas3", 5712200, 2663400, 1106000)
-    # print(aus.gms("phi", geod3[0], stellen), aus.gms("lambda", geod3[1], stellen), "h =", geod3[2])
-    # cart1 = ellips.geod2cart(geod1[0], geod1[1], geod1[2])
-    # print(aus.xyz(cart1[0], cart1[1], cart1[2], 10))
-    # cart2 = ellips.geod2cart(geod2[0], geod2[1], geod2[2])
-    # print(aus.xyz(cart2[0], cart2[1], cart2[2], 10))
-    # cart3 = ellips.geod2cart(geod3[0], geod3[1], geod3[2])
-    # print(aus.xyz(cart3[0], cart3[1], cart3[2], 10))
-
-    # test_cart = ellips.geod2cart(0.175, 0.444, 100)
-    # print(aus.xyz(test_cart[0], test_cart[1], test_cart[2], 10))
-    pass

jacobian_Ligas.py

@@ -1,6 +1,19 @@
-import numpy as np
+from typing import Tuple
 
-def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+import numpy as np
+from numpy.typing import NDArray
+
+
+def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
+    """
+    Aufstellen des Gleichungssystem für den ersten Fall
+    :param E: Konstante E
+    :param F: Konstante F
+    :param G: Konstante G
+    :param pG: Punkt über dem Ellipsoid
+    :param pE: Punkt auf dem Ellipsoid
+    :return: inverse Jacobi-Matrix, Gleichungssystem
+    """
     j11 = 2 * E * pE[0]
     j12 = 2 * F * pE[1]
     j13 = 2 * G * pE[2]
@@ -23,7 +36,16 @@ def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
 
     return invJ, fxE
 
-def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+def case2(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
+    """
+    Aufstellen des Gleichungssystem für den zweiten Fall
+    :param E: Konstante E
+    :param F: Konstante F
+    :param G: Konstante G
+    :param pG: Punkt über dem Ellipsoid
+    :param pE: Punkt auf dem Ellipsoid
+    :return: inverse Jacobi-Matrix, Gleichungssystem
+    """
     j11 = 2 * E * pE[0]
     j12 = 2 * F * pE[1]
     j13 = 2 * G * pE[2]
@@ -35,17 +57,29 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
     j33 = (pE[1] - pG[1]) * G - F * pE[1]
 
     detJ = j11 * j22 * j33 - j21 * j12 * j33 + j21 * j13 * j32
-    invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22],
-                              [-j21*j33, j11*j33, j13*j21],
-                              [j21*j32, -j11*j32, j11*j22-j12*j21]])
+    if detJ == 0:
+        invJ, fxE = case3(E, F, G, pG, pE)
+    else:
+        invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22],
+                                  [-j21*j33, j11*j33, j13*j21],
+                                  [j21*j32, -j11*j32, j11*j22-j12*j21]])
 
-    fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
-                    (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
-                    (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
+        fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
+                        (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
+                        (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
 
     return invJ, fxE
 
-def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+def case3(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
+    """
+    Aufstellen des Gleichungssystem für den dritten Fall
+    :param E: Konstante E
+    :param F: Konstante F
+    :param G: Konstante G
+    :param pG: Punkt über dem Ellipsoid
+    :param pE: Punkt auf dem Ellipsoid
+    :return: inverse Jacobi-Matrix, Gleichungssystem
+    """
     j11 = 2 * E * pE[0]
     j12 = 2 * F * pE[1]
     j13 = 2 * G * pE[2]
@@ -57,13 +91,15 @@ def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
     j33 = (pE[1] - pG[1]) * G - F * pE[1]
 
     detJ = -j11 * j23 * j32 - j21 * j12 * j33 + j21 * j13 * j32
+    if detJ == 0:
+        invJ, fxE = case2(E, F, G, pG, pE)
+    else:
+        invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23],
+                                  [-j21*j33, j11*j33, -(j11*j23-j13*j21)],
+                                  [j21*j32, -j11*j32, -j12*j21]])
 
-    invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23],
-                              [-j21*j33, j11*j33, -(j11*j23-j13*j21)],
-                              [j21*j32, -j11*j32, -j12*j21]])
-
-    fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
-                    (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0],
-                    (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
+        fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
+                        (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0],
+                        (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
 
     return invJ, fxE

nicht abgeben/GHA_biaxial/bessel.py

@@ -0,0 +1,36 @@
+from typing import Tuple
+
+import scipy as sp
+from ellipsoid_biaxial import EllipsoidBiaxial
+from numpy import *
+
+
+def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0:float, s: float) -> Tuple[float, float, float]:
+    """
+    Berechnung der 1.GHA auf einem Rotationsellipsoid nach Bessel
+    :param re:
+    :param phi0:
+    :param lamb0:
+    :param alpha0:
+    :param s:
+    :return:
+    """
+    psi0 = re.phi2psi(phi0)
+    clairant = arcsin(cos(psi0) * sin(alpha0))
+    sigma0 = arcsin(sin(psi0) / cos(clairant))
+
+    sqrt_sigma = lambda sigma: sqrt(1 + re.e_ ** 2 * cos(clairant) ** 2 * sin(sigma) ** 2)
+    int_sqrt_sigma = lambda sigma: sp.integrate.quad(sqrt_sigma, sigma0, sigma)[0]
+
+    f_sigma1_i = lambda sigma1_i: (int_sqrt_sigma(sigma1_i) - s / re.b)
+
+    sigma1_0 = sigma0 + s / re.a
+    sigma1 = sp.optimize.newton(f_sigma1_i, sigma1_0)
+    psi1 = arcsin(cos(clairant) * sin(sigma1))
+    phi1 = re.psi2phi(psi1)
+    alpha1 = arcsin(sin(clairant) / cos(psi1))
+
+    f_d_lambda = lambda sigma: sin(clairant) * sqrt_sigma(sigma) / (1 - cos(clairant) ** 2 * sin(sigma) ** 2)
+    d_lambda = sqrt(1-re.e**2) * sp.integrate.quad(f_d_lambda, sigma0, sigma1)[0]
+    lamb1 = lamb0 + d_lambda
+    return phi1, lamb1, alpha1

nicht abgeben/GHA_biaxial/gauss.py

@@ -0,0 +1,101 @@
+from typing import Tuple
+
+from ellipsoid_biaxial import EllipsoidBiaxial
+from numpy import arctan, cos, sin, sqrt, tan
+
+
+def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, eps: float = 1e-12) -> Tuple[float, float, float]:
+    """
+    Berechnung der 1. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln
+    :param re: Klasse Ellipsoid
+    :param phi0: Breite Punkt 0
+    :param lamb0: Länge Punkt 0
+    :param alpha0: Azimut der geodätischen Linie in Punkt 1
+    :param s: Strecke zu Punkt 1
+    :param eps: Abbruchkriterium für Winkelgrößen
+    :return: Breite, Länge, Azimut von Punkt 21
+    """
+    t = lambda phi: tan(phi)
+    eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2)
+
+    F1 = lambda alpha, phi, s: 1 + (2 + 3 * t(phi) ** 2 + 2 * eta(phi) ** 2) / (24 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \
+                               + (t(phi)**2 - 1) * eta(phi) ** 2 / (8 * re.N(phi) ** 2) * cos(alpha) ** 2 * s ** 2
+
+    F2 = lambda alpha, phi, s: 1 + t(phi) ** 2 / (24 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \
+                               - (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / (24 * re.N(phi) ** 2) * cos(alpha) ** 2 * s ** 2
+
+    F3 = lambda alpha, phi, s: 1 + (1 + eta(phi) ** 2) / (12 * re.N(phi) ** 2) * sin(alpha) ** 2 * s ** 2 \
+                               + (3 + 8 * eta(phi)**2) / (24 * re.N(phi) ** 2) * cos(alpha) ** 2 * s ** 2
+
+    phi1_i = lambda alpha, phi: phi0 + cos(alpha) / re.M(phi) * s * F1(alpha, phi, s)
+    lamb1_i = lambda alpha, phi: lamb0 + sin(alpha) / (re.N(phi) * cos(phi)) * s * F2(alpha, phi, s)
+    alpha1_i = lambda alpha, phi: alpha0 + sin(alpha) * tan(phi) / re.N(phi) * s * F3(alpha, phi, s)
+
+    phi_m_i = lambda phi1: (phi0 + phi1) / 2
+    alpha_m_i = lambda alpha1: (alpha0 + alpha1) / 2
+
+    phi_m = []
+    alpha_m = []
+    phi1 = []
+    lamb1 = []
+    alpha1 = []
+
+    # 1. Näherung für P1
+    phi1.append(phi0 + cos(alpha0) / re.M(phi0) * s)
+    lamb1.append(lamb0 + sin(alpha0) / (re.N(phi0) * cos(phi0)) * s)
+    alpha1.append(alpha0 + sin(alpha0) * tan(phi0) / re.N(phi0) * s)
+
+    while True:
+        # Berechnug P_m durch Mittelbildung
+        phi_m.append(phi_m_i(phi1[-1]))
+        alpha_m.append(alpha_m_i(alpha1[-1]))
+        # Berechnung P1
+        phi1.append(phi1_i(alpha_m[-1], phi_m[-1]))
+        lamb1.append(lamb1_i(alpha_m[-1], phi_m[-1]))
+        alpha1.append(alpha1_i(alpha_m[-1], phi_m[-1]))
+
+        # Abbruchkriterium
+        if abs(phi1[-2] - phi1[-1]) < eps and \
+           abs(lamb1[-2] - lamb1[-1]) < eps and \
+           abs(alpha1[-2] - alpha1[-1]) < eps:
+            break
+    return phi1[-1], lamb1[-1], alpha1[-1]
+
+
+def gha2(re: EllipsoidBiaxial, phi0: float, lamb0: float, phi1: float, lamb1: float):
+    """
+    Berechnung der 2. Geodätische Hauptaufgabe nach Gauß´schen Mittelbreitenformeln
+    :param re: Klasse Ellipsoid
+    :param phi0: Breite Punkt 1
+    :param lamb0: Länge Punkt 1
+    :param phi1: Breite Punkt 2
+    :param lamb1: Länge Punkt 2
+    :return: Länge der geodätischen Linie, Azimut von P1 nach P2, Azimut von P2 nach P1
+    """
+
+    t = lambda phi: tan(phi)
+    eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2)
+
+    phi_0 = (phi0 + phi1) / 2
+    d_phi = phi1 - phi0
+    d_lambda = lamb1 - lamb0
+
+    f_A = lambda phi: (2 + 3*t(phi)**2 + 2*eta(phi)**2) / 24
+    f_B = lambda phi: ((t(phi)**2 - 1) * eta(phi)**2) / 8
+    # f_C = lambda phi: (t(phi)**2) / 24
+    f_D = lambda phi: (1 + eta(phi)**2 - 9 * eta(phi)**2 * t(phi)**2) / 24
+
+    F1 = lambda phi: d_phi * re.M(phi) * (1 - f_A(phi) * d_lambda ** 2 * cos(phi) ** 2 -
+                                          f_B(phi) * d_phi ** 2 / re.V(phi) ** 4)
+    F2 = lambda phi: d_lambda * re.N(phi) * cos(phi) * (1 - 1 / 24 * d_lambda ** 2 * sin(phi) ** 2 +
+                                                           f_D(phi) * d_phi ** 2 / re.V(phi) ** 4)
+
+    s = sqrt(F1(phi_0) ** 2 + F2(phi_0) ** 2)
+    A_0 = arctan(F2(phi_0) / F1(phi_0))
+    d_A = d_lambda * sin(phi_0) * (1 + (1 + eta(phi_0) ** 2) / 12 * s ** 2 * sin(A_0) ** 2 / re.N(phi_0) ** 2 +
+                                      (3 + 8 * eta(phi_0) ** 2) / 24 * s ** 2 * cos(A_0) ** 2 / re.N(phi_0) ** 2)
+
+    alpha0 = A_0 - d_A / 2
+    alpha1 = A_0 + d_A / 2
+
+    return alpha0, alpha1, s

nicht abgeben/GHA_biaxial/rk.py

@@ -0,0 +1,33 @@
+from typing import Tuple
+
+import numpy as np
+from ellipsoid_biaxial import EllipsoidBiaxial
+from numpy import cos, sin, tan
+from numpy.typing import NDArray
+
+import runge_kutta as rk
+
+
+def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, num: int) -> Tuple[float, float, float]:
+    """
+    Berechnung der 1. GHA auf einem Rotationsellipsoid mittels RK4
+    :param re:
+    :param phi0:
+    :param lamb0:
+    :param alpha0:
+    :param s:
+    :param num:
+    :return:
+    """
+    def buildODE():
+        def ODE(s: float, v: NDArray):
+            phi, lam, A = v
+            V = re.V(phi)
+            dphi = cos(A) * V ** 3 / re.c
+            dlam = sin(A) * V / (cos(phi) * re.c)
+            dA = tan(phi) * sin(A) * V / re.c
+            return np.array([dphi, dlam, dA])
+        return ODE
+
+    _, funktionswerte = rk.rk4(buildODE(), 0, np.array([phi0, lamb0, alpha0]), s, num)
+    return funktionswerte[-1][0], funktionswerte[-1][1], funktionswerte[-1][2]

nicht abgeben/Tests/alpha_conversion_test.ipynb

@@ -0,0 +1,101 @@
+{
+ "cells": [
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ],
+   "id": "a78faf7f4883772f",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "%reload_ext autoreload\n",
+    "%autoreload 2\n",
+    "import numpy as np\n",
+    "\n",
+    "import winkelumrechnungen as wu\n",
+    "from GHA_triaxial.utils import alpha_ell2para, alpha_para2ell\n",
+    "from ellipsoid_triaxial import EllipsoidTriaxial"
+   ],
+   "id": "46aa84a937fea491",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "ell = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n",
+    "diffs = []\n",
+    "for beta_deg in range(-90, 91, 15):\n",
+    "    for lamb_deg in range(-180, 180, 15):\n",
+    "        for alpha_deg in range(0, 360, 15):\n",
+    "            beta = wu.deg2rad(beta_deg)\n",
+    "            lamb = wu.deg2rad(lamb_deg)\n",
+    "            u, v = ell.ell2para(beta, lamb)\n",
+    "            alpha = wu.deg2rad(alpha_deg)\n",
+    "\n",
+    "            alpha_para_1, *_ = alpha_ell2para(ell, beta, lamb, alpha)\n",
+    "            alpha_ell_1, *_ = alpha_para2ell(ell, u, v, alpha_para_1)\n",
+    "            diff_1 = wu.deg2rad(abs(alpha_ell_1 - alpha))/3600\n",
+    "\n",
+    "            alpha_ell_2, *_ = alpha_para2ell(ell, u, v, alpha)\n",
+    "            alpha_para_2, *_ = alpha_ell2para(ell, beta, lamb, alpha_ell_2)\n",
+    "            diff_2 = wu.deg2rad(abs(alpha_para_2 - alpha))/3600\n",
+    "\n",
+    "            diffs.append((beta_deg, lamb_deg, alpha_deg, diff_1, diff_2))\n",
+    "diffs = np.array(diffs)"
+   ],
+   "id": "82fc6cbbe7d5abcb",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "i_max_ell = np.argmax(diffs[:, 3])\n",
+    "max_ell = diffs[i_max_ell, 3]\n",
+    "point_max_ell = diffs[i_max_ell, :3]\n",
+    "\n",
+    "i_max_para = np.argmax(diffs[:, 4])\n",
+    "max_para = diffs[i_max_para, 4]\n",
+    "point_max_para = diffs[i_max_para, :4]\n",
+    "\n",
+    "print(f'Für elliptisches Alpha = {point_max_ell[2]}° und beta = {point_max_ell[0]}°, lamb = {point_max_ell[1]}°: diff = {max_ell}\"')\n",
+    "print(f'Für parametrisches Alpha = {point_max_para[2]}° und beta = {point_max_para[0]}°, lamb = {point_max_para[1]}°: diff = {max_ell}\"')\n",
+    "pass"
+   ],
+   "id": "97b5b8c9ca5377ab",
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

nicht abgeben/Tests/conversions_test.ipynb

@@ -0,0 +1,246 @@
+{
+ "cells": [
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ],
+   "id": "746c5b9e4c0226e7",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "cell_type": "code",
+   "id": "initial_id",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "%reload_ext autoreload\n",
+    "%autoreload 2\n",
+    "from itertools import product\n",
+    "\n",
+    "import numpy as np\n",
+    "import pandas as pd\n",
+    "import plotly.graph_objects as go\n",
+    "\n",
+    "import winkelumrechnungen as wu\n",
+    "from ellipsoid_triaxial import EllipsoidTriaxial"
+   ],
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# ellips = \"KarneyTest2024\"\n",
+    "ellips = \"BursaSima1980\"\n",
+    "# ellips = \"Fiction\"\n",
+    "ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(ellips)"
+   ],
+   "id": "7b05ca89fcd7b331",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "def deg_range(start, stop, step):\n",
+    "    return [float(x) for x in range(start, stop + step, step)]\n",
+    "\n",
+    "def asymptotic_range(start, direction=\"up\", max_decimals=4):\n",
+    "    values = []\n",
+    "    for d in range(0, max_decimals + 1):\n",
+    "        step = 10 ** -d\n",
+    "        if direction == \"up\":\n",
+    "            values.append(start + (1 - step))\n",
+    "        else:\n",
+    "            values.append(start - (1 - step))\n",
+    "    return values"
+   ],
+   "id": "61a6b14fef0180ad",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "beta_5_85 = deg_range(5, 85, 5)\n",
+    "lambda_5_85 = deg_range(5, 85, 5)\n",
+    "beta_5_90 = deg_range(5, 90, 5)\n",
+    "lambda_5_90 = deg_range(5, 90, 5)\n",
+    "beta_0_90 = deg_range(0, 90, 5)\n",
+    "lambda_0_90 = deg_range(0, 90, 5)\n",
+    "beta_90 = [90.0]\n",
+    "lambda_90 = [90.0]\n",
+    "beta_0 = [0.0]\n",
+    "lambda_0 = [0.0]\n",
+    "beta_asym_89 = asymptotic_range(89.0, direction=\"up\")\n",
+    "lambda_asym_0 = asymptotic_range(1.0, direction=\"down\")"
+   ],
+   "id": "f7184980a4b930b7",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "groups = {\n",
+    "    1: list(product(beta_5_85, lambda_5_85)),\n",
+    "    2: list(product(beta_0, lambda_0_90)),\n",
+    "    3: list(product(beta_5_85, lambda_0)),\n",
+    "    4: list(product(beta_90, lambda_5_90)),\n",
+    "    5: list(product(beta_asym_89, lambda_asym_0)),\n",
+    "    6: list(product(beta_5_85, lambda_90)),\n",
+    "    7: list(product(lambda_asym_0, lambda_0_90)),\n",
+    "    8: list(product(beta_0_90, lambda_asym_0)),\n",
+    "    9: list(product(beta_asym_89, lambda_0_90)),\n",
+    "    10: list(product(beta_0_90, beta_asym_89)),\n",
+    "}"
+   ],
+   "id": "cea9fd9cce6a4fd1",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "for nr, points in groups.items():\n",
+    "    points_cart = []\n",
+    "    for point in points:\n",
+    "        beta, lamb = point\n",
+    "        cart = ell.ell2cart(wu.deg2rad(beta), wu.deg2rad(lamb))\n",
+    "        points_cart.append(cart)\n",
+    "    groups[nr] = points_cart"
+   ],
+   "id": "17a6a130782a89ce",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "results = {}\n",
+    "\n",
+    "for nr, points in groups.items():\n",
+    "    group_results = {\"ell\": [],\n",
+    "                     \"para\": [],\n",
+    "                     \"geod\": []}\n",
+    "    for point in points:\n",
+    "        elli = ell.cart2ell(point)\n",
+    "        cart_elli = ell.ell2cart(elli[0], elli[1])\n",
+    "        group_results[\"ell\"].append(np.linalg.norm(point - cart_elli, axis=-1))\n",
+    "\n",
+    "        para = ell.cart2para(point)\n",
+    "        cart_para = ell.para2cart(para[0], para[1])\n",
+    "        group_results[\"para\"].append(np.linalg.norm(point - cart_para, axis=-1))\n",
+    "\n",
+    "        geod = ell.cart2geod(point, \"ligas3\")\n",
+    "        cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])\n",
+    "        group_results[\"geod\"].append(np.linalg.norm(point - cart_geod, axis=-1))\n",
+    "\n",
+    "    group_results[\"ell\"] = np.array(group_results[\"ell\"])\n",
+    "    group_results[\"para\"] = np.array(group_results[\"para\"])\n",
+    "    group_results[\"geod\"] = np.array(group_results[\"geod\"])\n",
+    "    results[nr] = group_results"
+   ],
+   "id": "c3298ea233bca274",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# with open(f\"conversion_results_{ellips}.pkl\", \"wb\") as f:\n",
+    "#     pickle.dump(results, f)"
+   ],
+   "id": "e1285860be416ad3",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "# with open(f\"conversion_results_{ellips}.pkl\", \"rb\") as f:\n",
+    "#     results = pickle.load(f)"
+   ],
+   "id": "d26720e34595ccbc",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "df = pd.DataFrame({\n",
+    "    \"Gruppe\": [nr for nr in results.keys()],\n",
+    "    \"max_Δr_ell\": [f\"{max(result[\"ell\"]):.3g}\" for result in results.values()],\n",
+    "    \"max_Δr_para\": [f\"{max(result[\"para\"]):.3g}\" for result in results.values()],\n",
+    "    \"max_Δr_geod\": [f\"{max(result[\"geod\"]):.3g}\" for result in results.values()]\n",
+    "})"
+   ],
+   "id": "4e2e55e4699ec81e",
+   "outputs": [],
+   "execution_count": null
+  },
+  {
+   "metadata": {},
+   "cell_type": "code",
+   "source": [
+    "fig = go.Figure(data=[go.Table(\n",
+    "    header=dict(\n",
+    "        values=list(df.columns),\n",
+    "        fill_color=\"lightgrey\",\n",
+    "        align=\"left\"\n",
+    "    ),\n",
+    "    cells=dict(\n",
+    "        values=[df[col] for col in df.columns],\n",
+    "        align=\"left\"\n",
+    "    )\n",
+    ")])\n",
+    "fig.update_layout(\n",
+    "    template=\"simple_white\",\n",
+    "    width=650,\n",
+    "    height=len(groups)*20+80,\n",
+    "    margin=dict(l=20, r=20, t=20, b=20))\n",
+    "\n",
+    "fig.show()\n",
+    "# fig.write_image(f\"conversion_results_{ellips}.png\", width=650, height=len(groups)*20+80, scale=2)"
+   ],
+   "id": "c2fa82afef2d6e0e",
+   "outputs": [],
+   "execution_count": null
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 2
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython2",
+   "version": "2.7.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

nicht abgeben/ellipsoid_biaxial.py

@@ -0,0 +1,126 @@
+from typing import Tuple
+
+import numpy as np
+from numpy import arctan2, cos, sin, sqrt
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+
+
+class EllipsoidBiaxial:
+    """
+    Klasse für Rotationsellipdoide
+    """
+    def __init__(self, a: float, b: float):
+        self.a = a
+        self.b = b
+        self.c = a ** 2 / b
+        self.e = sqrt(a ** 2 - b ** 2) / a
+        self.e_ = sqrt(a ** 2 - b ** 2) / b
+
+    @classmethod
+    def init_name(cls, name: str) -> EllipsoidBiaxial:
+        """
+        Erstellen eines Rotationsellipdoids nach Namen
+        :param name: Name des Rotationsellipsoids
+        :return: Rotationsellipsoid
+        """
+        if name == "Bessel":
+            a = 6377397.15508
+            b = 6356078.96290
+            return cls(a, b)
+        elif name == "Hayford":
+            a = 6378388
+            f = 1/297
+            b = a - a * f
+            return cls(a, b)
+        elif name == "Krassowski":
+            a = 6378245
+            f = 298.3
+            b = a - a * f
+            return cls(a, b)
+        elif name == "WGS84":
+            a = 6378137
+            f = 298.257223563
+            b = a - a * f
+            return cls(a, b)
+        else:
+            raise Exception(f"EllipsoidBiaxial.init_name: Name {name} unbekannt")
+
+    @classmethod
+    def init_af(cls, a: float, f: float) -> EllipsoidBiaxial:
+        """
+        Erstellen eines Rotationsellipdoids aus der großen Halbachse und der Abplattung
+        :param a: große Halbachse
+        :param f: großen Halbachse
+        :return: Rotationsellipsoid
+        """
+        b = a - a * f
+        return cls(a, b)
+
+    V = lambda self, phi: sqrt(1 + self.e_ ** 2 * cos(phi) ** 2)
+    M = lambda self, phi: self.c / self.V(phi) ** 3
+    N = lambda self, phi: self.c / self.V(phi)
+
+    beta2psi = lambda self, beta: arctan2(self.a * sin(beta), self.b * cos(beta))
+    beta2phi = lambda self, beta: arctan2(self.a ** 2 * sin(beta), self.b ** 2 * cos(beta))
+
+    psi2beta = lambda self, psi: arctan2(self.b * sin(psi), self.a * cos(psi))
+    psi2phi = lambda self, psi: arctan2(self.a * sin(psi), self.b * cos(psi))
+
+    phi2beta = lambda self, phi: arctan2(self.b**2 * sin(phi), self.a**2 * cos(phi))
+    phi2psi = lambda self, phi: arctan2(self.b * sin(phi), self.a * cos(phi))
+
+    phi2p = lambda self, phi: self.N(phi) * cos(phi)
+
+    def bi_cart2ell(self, point: NDArray, Eh: float = 0.001, Ephi: float = wu.gms2rad([0, 0, 0.001])) -> Tuple[float, float, float]:
+        """
+        Umrechnung von kartesischen in ellipsoidische Koordinaten auf einem Rotationsellipsoid
+        # TODO: Quelle
+        :param point: Punkt in kartesischen Koordinaten
+        :param Eh: Grenzwert für die Höhe
+        :param Ephi: Grenzwert für die Breite
+        :return: ellipsoidische Breite, Länge, geodätische Höhe
+        """
+        x, y, z = point
+
+        lamb = arctan2(y, x)
+
+        p = sqrt(x**2+y**2)
+
+        phi_null = arctan2(z, p*(1 - self.e**2))
+
+        hi = [0]
+        phii = [phi_null]
+
+        i = 0
+
+        while True:
+            N = self.a / sqrt(1 - self.e**2 * sin(phii[i])**2)
+            h = p / cos(phii[i]) - N
+            phi = arctan2(z,  p * (1-(self.e**2*N) / (N+h)))
+            hi.append(h)
+            phii.append(phi)
+            dh = abs(hi[i]-h)
+            dphi = abs(phii[i]-phi)
+            i += 1
+            if dh < Eh:
+                if dphi < Ephi:
+                    break
+        return phi, lamb, h
+
+    def bi_ell2cart(self, phi: float, lamb: float, h: float) -> NDArray:
+        """
+        Umrechnung von ellipsoidischen in kartesische Koordinaten auf einem Rotationsellipsoid
+        # TODO: Quelle
+        :param phi: ellipsoidische Breite
+        :param lamb: ellipsoidische Länge
+        :param h: geodätische Höhe
+        :return: Punkt in kartesischen Koordinaten
+        """
+        W = sqrt(1 - self.e**2 * sin(phi)**2)
+        N = self.a / W
+        x = (N+h) * cos(phi) * cos(lamb)
+        y = (N+h) * cos(phi) * sin(lamb)
+        z = (N * (1-self.e**2) + h) * sin(phi)
+        return np.array([x, y, z])

nicht abgeben/kugel.py

@@ -0,0 +1,100 @@
+from typing import Tuple
+
+import numpy as np
+from numpy import arccos, arcsin, arctan2, cos, pi, sin, sqrt
+from numpy.typing import NDArray
+
+import winkelumrechnungen as wu
+
+
+def cart2sph(point: NDArray) -> Tuple[float, float, float]:
+    """
+    Umrechnung von kartesischen in sphärische Koordinaten
+    # TODO: Quelle
+    :param point: Punkt in kartesischen Koordinaten
+    :return: Radius, Breite, Länge
+    """
+    x, y, z = point
+    r = sqrt(x**2 + y**2 + z**2)
+    phi = arctan2(z, sqrt(x**2 + y**2))
+    lamb = arctan2(y, x)
+
+    return r, phi, lamb
+
+
+def sph2cart(r: float, phi: float, lamb: float) -> NDArray:
+    """
+    Umrechnung von sphärischen in kartesische Koordinaten
+    # TODO: Quelle
+    :param r: Radius
+    :param phi: Breite
+    :param lamb: Länge
+    :return: Punkt in kartesischen Koordinaten
+    """
+    x = r * cos(phi) * cos(lamb)
+    y = r * cos(phi) * sin(lamb)
+    z = r * sin(phi)
+
+    return np.array([x, y, z])
+
+
+def gha1(R: float, phi0: float, lamb0: float, s: float, alpha0: float) -> Tuple[float, float]:
+    """
+    Berechnung der 1. GHA auf der Kugel
+    # TODO: Quelle
+    :param R: Radius
+    :param phi0: Breite des Startpunktes
+    :param lamb0: Länge des Startpunktes
+    :param s: Strecke
+    :param alpha0: Azimut
+    :return: Breite, Länge des Zielpunktes
+    """
+    s_ = s / R
+
+    lamb1 = lamb0 + arctan2(sin(s_) * sin(alpha0),
+                            cos(phi0) * cos(s_) - sin(phi0) * sin(s_) * cos(alpha0))
+
+    phi1 = arcsin(sin(phi0) * cos(s_) + cos(phi0) * sin(s_) * cos(alpha0))
+
+    return phi1, lamb1
+
+
+def gha2(R: float, phi0: float, lamb0: float, phi1: float, lamb1: float) -> Tuple[float, float, float]:
+    """
+    Berechnung der 2. GHA auf der Kugel
+    # TODO: Quelle
+    :param R: Radius
+    :param phi0: Breite des Startpunktes
+    :param lamb0: Länge des Startpunktes
+    :param phi1: Breite des Zielpunktes
+    :param lamb1: Länge des Zielpunktes
+    :return: Azimut im Startpunkt, Azimut im Zielpunkt, Strecke
+    """
+    s_ = arccos(sin(phi0) * sin(phi1) + cos(phi0) * cos(phi1) * cos(lamb1 - lamb0))
+    s = R * s_
+
+    alpha0 = arctan2(cos(phi1) * sin(lamb1 - lamb0),
+                     cos(phi0) * sin(phi1) - sin(phi0) * cos(phi1) * cos(lamb1 - lamb0))
+
+    alpha1 = arctan2(-cos(phi0) * sin(lamb1 - lamb0),
+                     cos(phi1) * sin(phi0) - sin(phi1) * cos(phi0) * cos(lamb1 - lamb0))
+    if alpha1 < 0:
+        alpha1 += 2 * pi
+
+    return alpha0, alpha1, s
+
+
+if __name__ == "__main__":
+    R = 6378815.904                                             # Bern
+    phi0 = wu.deg2rad(10)
+    lamb0 = wu.deg2rad(40)
+    alpha0 = wu.deg2rad(100)
+    s = 10000
+    phi1, lamb1 = gha1(R, phi0, lamb0, s, alpha0)
+    alpha0_g, alpha1, s_g = gha2(R, phi0, lamb0, phi1, lamb1)
+    phi1 = wu.rad2deg(phi1)
+    lamb1 = wu.rad2deg(lamb1)
+    alpha0_g = wu.rad2deg(alpha0_g)
+    alpha1 = wu.rad2deg(alpha1)
+
+    pass

nicht abgeben/test.py

@@ -0,0 +1,8 @@
+import numpy as np
+
+import ellipsoid_triaxial
+
+ell = ellipsoid_triaxial.EllipsoidTriaxial.init_name("KarneyTest2024")
+
+cart = ell.para2cart(0, np.pi/2)
+print(cart)

phi_ellipse.py

@@ -1,55 +0,0 @@
-import winkelumrechnungen as wu
-
-
-def polyapp_tscheby_hayford(s: float) -> float:
-    """
-    Berechnung der ellipsoidisch geodätischen Breite.
-    Polynomapproximation mittels Tschebyscheff-Polynomen.
-    Auf dem Hayford-Ellipsoid.
-    :param s: Strecke auf einer Ellipse vom Äquator aus
-    :type s: float
-    :return: ellipsoidisch geodätische Breite
-    :rtype: float
-    """
-    c0 = 1
-    c1 = -0.00837809325
-    c2 = 0.00428127367
-    c3 = -0.00114523986
-    c4 = 0.00023219707
-    c5 = -0.00004421222
-    c6 = 0.00000570244
-    alpha = wu.gms2rad([0, 0, 325643.97199])
-    s90 = 10002288.2990
-
-    xi = s/s90
-
-    phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) +
-                        c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6))
-    return phi
-
-
-def polyapp_tscheby_bessel(s: float) -> float:
-    """
-    Berechnung der ellipsoidisch geodätischen Breite.
-    Polynomapproximation mittels Tschebyscheff-Polynomen.
-    Auf dem Bessel-Ellipsoid.
-    :param s: Strecke auf einer Ellipse vom Äquator aus
-    :type s: float
-    :return: ellipsoidisch geodätische Breite
-    :rtype: float
-    """
-    c0 = 1
-    c1 = -0.00831729565
-    c2 = 0.00424914906
-    c3 = -0.00113566119
-    c4 = 0.00022976983
-    c5 = -0.00004363980
-    c6 = 0.00000562025
-    alpha = wu.gms2rad([0, 0, 325632.08677])
-    s90 = 10000855.7644
-
-    xi = s/s90
-
-    phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) +
-                        c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6))
-    return phi

requirements.txt

@@ -0,0 +1,7 @@
+numpy~=2.3.4
+plotly~=6.4.0
+pandas~=2.3.3
+scipy~=1.16.3
+dash-bootstrap-components~=2.0.4
+dash~=4.0.0
+matplotlib~=3.10.7

runge_kutta.py

@@ -0,0 +1,155 @@
+from typing import Callable
+
+import numpy as np
+from numpy.typing import NDArray
+
+
+def rk4(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]:
+    """
+    Standard Runge-Kutta Verfahren 4. Ordnung
+    :param ode: ODE-System als Funktion
+    :param t0: Startwert der unabhängigen Variable
+    :param v0: Startwerte
+    :param weite: Integrationsweite
+    :param schritte: Schrittzahl
+    :param fein: Fein-Rechnung?
+    :return: Variable und Funktionswerte an jedem Stützpunkt
+    """
+    h = weite/schritte
+
+    t_list = [t0]
+    werte = [v0]
+
+    for _ in range(schritte):
+        t = t_list[-1]
+        v = werte[-1]
+
+        if not fein:
+            v_next = rk4_step(ode, t, v, h)
+
+        else:
+            v_grob = rk4_step(ode, t, v, h)
+            v_half = rk4_step(ode, t, v, 0.5 * h)
+            v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+            v_next = v_fein + (v_fein - v_grob) / 15.0
+
+        t_list.append(t + h)
+        werte.append(v_next)
+
+    return t_list, werte
+
+def rk4_step(ode: Callable, t: float, v: NDArray, h: float) -> NDArray:
+    """
+    Ein Schritt des Runge-Kutta Verfahrens 4. Ordnung
+    :param ode: ODE-System als Funktion
+    :param t: unabhängige Variable
+    :param v: abhängige Variablen
+    :param h: Schrittweite
+    :return: abhängige Variablen nach einem Schritt
+    """
+    k1 = ode(t, v)
+    k2 = ode(t + 0.5 * h, v + 0.5 * h * k1)
+    k3 = ode(t + 0.5 * h, v + 0.5 * h * k2)
+    k4 = ode(t + h, v + h * k3)
+    return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+
+def rk4_end(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, fein: bool = False):
+    """
+    Standard Runge-Kutta Verfahren 4. Ordnung, nur Ausgabe der letzten Variablenwerte
+    :param ode: ODE-System als Funktion
+    :param t0: Startwert der unabhängigen Variable
+    :param v0: Startwerte
+    :param weite: Integrationsweite
+    :param schritte: Schrittzahl
+    :param fein: Fein-Rechnung?
+    :return: Variable und Funktionswerte am letzten Stützpunkt
+    """
+    h = weite / schritte
+    t = float(t0)
+    v = np.array(v0, dtype=float, copy=True)
+
+    for _ in range(schritte):
+        if not fein:
+            v_next = rk4_step(ode, t, v, h)
+        else:
+            v_grob = rk4_step(ode, t, v, h)
+            v_half = rk4_step(ode, t, v, 0.5 * h)
+            v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+            v_next = v_fein + (v_fein - v_grob) / 15.0
+
+        t += h
+        v = v_next
+
+    return t, v
+
+# RK4 mit Simpson bzw. Trapez
+def rk4_integral(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, integrand_at: Callable, fein: bool = False, simpson: bool = True):
+    """
+    Runge-Kutta Verfahren 4. Ordnung mit Simpson bzw. Trapez
+    :param ode: ODE-System als Funktion
+    :param t0: Startwert der unabhängigen Variable
+    :param v0: Startwerte
+    :param weite: Integrationsweite
+    :param integrand_at: Funktion
+    :param schritte: Schrittzahl
+    :param fein: Fein-Rechnung?
+    :param simpson: Simpson? Wenn nein, dann Trapez
+    :return: Variable und Funktionswerte am letzten Stützpunkt
+    """
+    h = weite / schritte
+    habs = abs(h)
+
+    t = float(t0)
+    v = np.array(v0, dtype=float, copy=True)
+
+    if simpson and (schritte % 2 == 0):
+        f0 = float(integrand_at(t, v))
+        odd_sum = 0.0
+        even_sum = 0.0
+        fN = None
+
+        for i in range(1, schritte + 1):
+            if not fein:
+                v_next = rk4_step(ode, t, v, h)
+            else:
+                v_grob = rk4_step(ode, t, v, h)
+                v_half = rk4_step(ode, t, v, 0.5 * h)
+                v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+                v_next = v_fein + (v_fein - v_grob) / 15.0
+
+            t += h
+            v = v_next
+
+            fi = float(integrand_at(t, v))
+            if i == schritte:
+                fN = fi
+            elif i % 2 == 1:
+                odd_sum += fi
+            else:
+                even_sum += fi
+
+        S = f0 + fN + 4.0 * odd_sum + 2.0 * even_sum
+        s = (habs / 3.0) * S
+        return t, v, s
+
+    f_prev = float(integrand_at(t, v))
+    acc = 0.0
+
+    for _ in range(schritte):
+        if not fein:
+            v_next = rk4_step(ode, t, v, h)
+        else:
+            v_grob = rk4_step(ode, t, v, h)
+            v_half = rk4_step(ode, t, v, 0.5 * h)
+            v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+            v_next = v_fein + (v_fein - v_grob) / 15.0
+
+        t += h
+        v = v_next
+
+        f_cur = float(integrand_at(t, v))
+        acc += 0.5 * (f_prev + f_cur)
+        f_prev = f_cur
+
+    s = habs * acc
+    return t, v, s

s_ellipse.py

@@ -1,78 +0,0 @@
-import numpy as np
-
-
-def reihenentwicklung(es: float, c: float, phi: float) -> float:
-    """
-    Berechnung der Strecke auf einer Ellipse.
-    Reihenentwicklung.
-    :param es: zweite numerische Exzentrizität
-    :type es: float
-    :param c: Polkrümmungshalbmesser
-    :type c: float
-    :param phi: ellipsoidisch geodästische Breite in Radiant
-    :type phi: float
-    :return: Strecke auf einer Ellipse vom Äquator aus
-    :rtype: float
-    """
-    Ass = 1 - 3/4*es**2 + 45/64*es**4 - 175/256*es**6 + 11025/16384*es**8
-    Bss = - 3/4*es**2 + 15/16*es**4 - 525/512*es**6 + 2205/2048*es**8
-    Css = 15/64*es**4 - 105/256*es**6 + 2205/4096*es**8
-    Dss = - 35/512*es**6 + 315/2048*es**8
-    print(f"A'' = {round(Ass, 10):.10f}\nB'' = {round(Bss, 10):.10f}\nC'' = {round(Css, 10):.10f}\nD'' = {round(Dss, 10):.10f}")
-
-    s = c * (Ass*phi + 1/2*Bss * np.sin(2*phi) + 1/4*Css * np.sin(4*phi) + 1/6*Dss * np.sin(6*phi))
-    return s
-
-
-def polyapp_tscheby_hayford(phi: float) -> float:
-    """
-    Berechnung der Strecke auf einer Ellipse.
-    Polynomapproximation mittels Tschebyscheff-Polynomen.
-    Auf dem Hayford-Ellipsoid.
-    :param phi: ellipsoidisch geodästische Breite in Radiant
-    :type phi: float
-    :return: Strecke auf einer Ellipse vom Äquator aus
-    :rtype: float
-    """
-    c1s = 0.00829376218
-    c2s = -0.00398963425
-    c3s = 0.00084200710
-    c4s = -0.0000648906
-    c5s = -0.00001075680
-    c6s = 0.00000396474
-    c7s = -0.00000046347
-    alpha = 9951793.0123
-
-    xi = 2 / np.pi * phi
-    xis = xi ^ 2
-
-    s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3
-                      + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7)
-    return s
-
-
-def polyapp_tscheby_bessel(phi: float) -> float:
-    """
-    Berechnung der Strecke auf einer Ellipse.
-    Polynomapproximation mittels Tschebyscheff-Polynomen.
-    Auf dem Bessel-Ellipsoid.
-    :param phi: ellipsoidisch geodästische Breite in Radiant
-    :type phi: float
-    :return: Strecke auf einer Ellipse vom Äquator aus
-    :rtype: float
-    """
-    c1s = 0.00823417717
-    c2s = -0.00396170744
-    c3s = 0.00083680249
-    c4s = -0.00006488462
-    c5s = -0.00001053242
-    c6s = 0.00000390854
-    c7s = -0.00000045768
-    alpha = 9950730.8876
-
-    xi = 2 / np.pi * phi
-    xis = xi ** 2
-
-    s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3
-                      + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7)
-    return s

test.py

@@ -1,29 +0,0 @@
-import numpy as np
-from scipy.special import factorial as fact
-from math import comb
-
-J = np.array([
-    [2, 3, 0],
-    [0, 3, 0],
-    [6, 0, 4]
-])
-
-xi = np.array([1, 2, 3])
-xi_col = xi.reshape(-1, 1)
-print(xi_col)
-xi_row = xi_col.reshape(1, -1).flatten()
-print(xi_row)
-
-# Spaltenvektor-Variante
-res_col = xi[:, None] - J @ xi[:, None]
-
-# Zeilenvektor-Variante
-res_row = xi[None, :] - xi[None, :] @ J
-
-print("Spaltenvektor:")
-print(res_col[0,0])
-print("Zeilenvektor:")
-print(res_row)
-t = 5
-l = 2
-print(fact(t+1-l) / (fact(t+1-l) * fact(l-1)), comb(l-1, t+1-l))

utils_angle.py

@@ -0,0 +1,54 @@
+import numpy as np
+
+import winkelumrechnungen as wu
+
+
+def arccot(x: float) -> float:
+    """
+    Berechnung von arccot eines Winkels
+    :param x: Winkel
+    :return: arccot(Winkel)
+    """
+    return np.arctan2(1.0, x)
+
+
+def cot(x: float) -> float:
+    """
+    Berechnung von cot eines Winkels
+    :param x: Winkel
+    :return: cot(Winkel)
+    """
+    return np.cos(x) / np.sin(x)
+
+
+def wrap_mpi_pi(x: float) -> float:
+    """
+    Wrap eines Winkels in den Wertebereich [-π, π)
+    :param x: Winkel
+    :return: Winkel in [-π, π)
+    """
+    return (x + np.pi) % (2 * np.pi) - np.pi
+
+
+def wrap_mhalfpi_halfpi(x: float) -> float:
+    """
+    Wrap eines Winkels in den Wertebereich [-π/2, π/2)
+    :param x: Winkel
+    :return: Winkel in [-π/2, π/2)
+    """
+    return (x + np.pi / 2) % np.pi - np.pi / 2
+
+
+def wrap_0_2pi(x: float) -> float:
+    """
+    Wrap eines Winkels in den Wertebereich [0, 2π)
+    :param x: Winkel
+    :return: Winkel in [0, 2π)
+    """
+    return x % (2 * np.pi)
+
+
+if __name__ == "__main__":
+    print(wu.rad2deg(wrap_mhalfpi_halfpi(wu.deg2rad(181))))
+    print(wu.rad2deg(wrap_0_2pi(wu.deg2rad(181))))
+    print(wu.rad2deg(wrap_mpi_pi(wu.deg2rad(181))))

winkelumrechnungen.py

@@ -1,4 +1,5 @@
 from numpy import *
+import numpy as np
 
 
 def deg2gms(deg: float) -> list:
@@ -9,14 +10,11 @@ def deg2gms(deg: float) -> list:
     :return: Winkel in Grad-Minuten-Sekunden
     :rtype: list
     """
-    gra = deg // 1
-    min = gra % 1
-    gra = gra // 1
-    min *= 60
-    sek = min % 1
-    min = min // 1
-    sek *= 60
-    return [gra, min, sek]
+    gra = int(deg)
+    minu_f = (deg - gra) * 60
+    minu = int(minu_f)
+    sek = (minu_f - minu) * 60
+    return [gra, minu, sek]
 
 
 def deg2gra(deg: float) -> float:
@@ -30,13 +28,13 @@ def deg2gra(deg: float) -> float:
     return deg * 10/9
 
 
-def deg2rad(deg: float) -> float:
+def deg2rad(deg: float | np.ndarray) -> float | np.ndarray:
     """
     Umrechnung von Grad in Radiant
     :param deg: Winkel in Grad
-    :type deg: float
+    :type deg: float or np.ndarray
     :return: Winkel in Radiant
-    :rtype: float
+    :rtype: float or np.ndarray
     """
     return deg * pi / 180
 
@@ -50,14 +48,11 @@ def gra2gms(gra: float) -> list:
     :rtype: list
     """
     deg = gra2deg(gra)
-    gra = deg // 1
-    min = gra % 1
-    gra = gra // 1
-    min *= 60
-    sek = min % 1
-    min = min // 1
-    sek *= 60
-    return [gra, min, sek]
+    gra = int(deg)
+    minu_f = (deg - gra) * 60
+    minu = int(minu_f)
+    sek = (minu_f - minu) * 60
+    return [gra, minu, sek]
 
 
 def gra2rad(gra: float) -> float:
@@ -113,13 +108,11 @@ def rad2gms(rad: float) -> list:
     :rtype: list
     """
     deg = rad2deg(rad)
-    min = deg % 1
-    gra = deg // 1
-    min *= 60
-    sek = min % 1
-    min = min // 1
-    sek *= 60
-    return [gra, min, sek]
+    gra = int(deg)
+    minu_f = (deg - gra) * 60
+    minu = int(minu_f)
+    sek = (minu_f - minu) * 60
+    return [gra, minu, sek]
 
 
 def gms2rad(gms: list) -> float: