Ellipsoid-Auswahl angepasst

revert Delete Tests/gha_resultsKarney.pkl

ES/Hansen_ES_CMA.py

@@ -1,7 +1,8 @@
 import numpy as np
+from numpy.typing import NDArray
 
 
-def felli(x):
+def felli(x: NDArray) -> float:
     N = x.shape[0]
     if N < 2:
         raise ValueError("dimension must be greater than one")
@@ -11,8 +12,7 @@ def felli(x):
 
 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):
-
+          bestEver=np.inf, noImproveGen=0, absTolImprove=1e-12, maxNoImproveGen=100, sigmaImprove=1e-12):
     if func_kwargs is None:
         func_kwargs = {}
 
@@ -27,7 +27,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
         N = xmean.shape[0]
 
     if stopeval is None:
-        stopeval = int(1e3 * N**2)
+        stopeval = int(1e3 * N ** 2)
 
     # Strategy parameter setting: Selection
     lambda_ = 4 + int(np.floor(3 * np.log(N)))
@@ -37,14 +37,14 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
     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)
+    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)
+    c1 = 2 / ((N + 1.3) ** 2 + mueff)
     cmu = min(1 - c1,
-              2 * (mueff - 2 + 1 / mueff) / ((N + 2)**2 + 2 * mueff))
+              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
@@ -54,7 +54,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
     D = np.eye(N)
     C = B @ D @ (B @ D).T
     eigeneval = 0
-    chiN = np.sqrt(N) * (1 - 1/(4*N) + 1/(21 * N**2))
+    chiN = np.sqrt(N) * (1 - 1 / (4 * N) + 1 / (21 * N ** 2))
 
     # Generation Loop
     counteval = 0
@@ -64,7 +64,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
 
     gen = 0
 
-    print(f'    [CMA-ES] Start: lambda = {lambda_}, sigma ={round(sigma, 6)}, stopeval = {stopeval}')
+    # print(f'    [CMA-ES] Start: lambda = {lambda_}, sigma ={round(sigma, 6)}, stopeval = {stopeval}')
 
     while counteval < stopeval:
         gen += 1
@@ -91,27 +91,24 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
             bestEver = fbest
             noImproveGen = 0
         else:
-            noImproveGen = noImproveGen + 1
+            noImproveGen += 1
 
-
-        if gen == 1 or gen%50==0:
-            print(f'    [CMA-ES] Gen {gen}, best = {round(fbest, 6)}, sigma = {sigma:.3g}')
+        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.')
+            # 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}')
+            # 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 = 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)
@@ -140,17 +137,17 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
         # 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")
+            # print('    [CMA-ES] stopfitness erreicht.')
+            # print("warning: flat fitness, consider reformulating the objective")
             break
 
-        #print(f"{counteval}: {arfitness[0]}")
+        # print(f"{counteval}: {arfitness[0]}")
 
-    #Final Message
-    #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)}')
+    # print(f'    [CMA-ES] Ende: Gen = {gen}, best = {round(bestValue, 6)}')
     return xmin
 
 

ES/gha1_ES.py

@@ -1,26 +1,17 @@
 from __future__ import annotations
-from typing import List, Optional, Tuple
+
+from typing import List, Tuple
+
 import numpy as np
-from ellipsoide import EllipsoidTriaxial
-from GHA_triaxial.gha1_ana import gha1_ana
-from GHA_triaxial.gha1_approx import gha1_approx
-from Hansen_ES_CMA import escma
-from utils_angle import wrap_to_pi
 from numpy.typing import NDArray
 
-
-def ellipsoid_formparameter(ell: EllipsoidTriaxial):
-    """
-    Berechnet die Formparameter des dreiachsigen Ellipsoiden nach Karney (2025), Gl. (2)
-    :param ell: Ellipsoid
-    :return: e, k und k'
-    """
-    nenner = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0))
-    k = np.sqrt(max(ell.ay * ell.ay - ell.b * ell.b, 0.0)) / nenner
-    k_ = np.sqrt(max(ell.ax * ell.ax - ell.ay * ell.ay, 0.0)) / nenner
-    e = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0)) / ell.ay
-
-    return e, k, k_
+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) \
@@ -37,7 +28,7 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
                 R (XYZ) = Punkt in XYZ
     """
     # Berechnungshilfen
-    omega = wrap_to_pi(omega)
+    omega = wrap_mpi_pi(omega)
     cb = np.cos(beta)
     sb = np.sin(beta)
     co = np.cos(omega)
@@ -58,45 +49,32 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
     # --- Ableitungen - Karney Gl. (5a,b,c)---
     # E = ∂R/∂ω
     dX_dw = -ell.ax * so * Sx / D
-    dY_dw =  ell.ay * cb * co
+    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
+    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 = np.linalg.norm(E)
-    Nn = np.linalg.norm(N)
-    Un = np.linalg.norm(U)
+    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 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
-    """
-    e, k, k_ = ellipsoid_formparameter(ell)
-    gamma_jacobi = float((k ** 2) * (np.cos(beta) ** 2) * (np.sin(alpha) ** 2) - (k_ ** 2) * (np.sin(omega) ** 2) * (np.cos(alpha) ** 2))
-
-    return gamma_jacobi
-
-
 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
@@ -110,11 +88,12 @@ def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_
     """
     v = (P_curr - P_prev).astype(float)
     vT = v - float(np.dot(v, U_hat_curr)) * U_hat_curr
-    vT_hat = vT / np.linalg.norm(vT)
+    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_to_pi(float(np.arctan2(sE, sN)))
+    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,
@@ -136,10 +115,21 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
     # 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|
-    d_beta = ds * float(np.cos(alpha_i)) / Nn_i
-    d_omega = ds * float(np.sin(alpha_i)) / En_i
+
+    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_to_pi(omega_i + d_omega)
+    omega_pred = wrap_mpi_pi(omega_i + d_omega)
 
     xmean = np.array([beta_pred, omega_pred], dtype=float)
 
@@ -156,10 +146,10 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
         :return: Fitnesswert (f)
         """
         beta = x[0]
-        omega = wrap_to_pi(x[1])
+        omega = wrap_mpi_pi(x[1])
 
-        P = ell.ell2cart(beta, omega) # in kartesischer Koordinaten
-        d = float(np.linalg.norm(P - P_i)) # Distanz zwischen
+        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
@@ -178,19 +168,19 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
 
         return f
 
-
-    xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES
+    xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0)  # Aufruf CMA-ES
 
     beta_best = xb[0]
-    omega_best = wrap_to_pi(xb[1])
-    P_best = ell.ell2cart(beta_best, omega_best)
+    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):
+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
@@ -199,15 +189,16 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
     :param alpha0: Azimut Startkoordinate
     :param s_total: Gesamtstrecke
     :param maxSegLen: maximale Segmentlänge
-    :return: Zielpunkt Pk und Azimut am Zielpunkt
+    :param all_points: Alle Punkte ausgeben?
+    :return: Zielpunkt Pk, Azimut am Zielpunkt und Punktliste
     """
     beta = float(beta0)
-    omega = wrap_to_pi(float(omega0))
-    alpha = wrap_to_pi(float(alpha0))
+    omega = wrap_mpi_pi(float(omega0))
+    alpha = wrap_mpi_pi(float(alpha0))
 
-    gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0
+    gamma0 = jacobi_konstante(beta, omega, alpha, ell)  # Referenz-γ0
 
-    points: List[NDArray] = [ell.ell2cart(beta, omega)]
+    P_all: List[NDArray] = [ell.ell2cart_karney(beta, omega)]
     alpha_end: List[float] = [alpha]
 
     s_acc = 0.0
@@ -216,35 +207,41 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
     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}")
+        # 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
-        points.append(P)
-        alpha_end.append(alpha)
+        P_all.append(P)
+        alpha_end.append(wrap_mpi_pi(alpha))
         if step > nsteps_est + 50:
-            raise RuntimeError("Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
-    Pk = points[-1]
-    alpha1 = alpha_end[-1]
+            raise RuntimeError("GHA1_ES: Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
+    Pk = P_all[-1]
+    alpha1 = float(alpha_end[-1])
 
-    return Pk, alpha1
+    if all_points:
+        return Pk, alpha1, np.array(P_all)
+    else:
+        return Pk, alpha1
 
 
 if __name__ == "__main__":
     ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    s = 188891.650873
-    alpha0 = 70/(180/np.pi)
-
-    P0 = ell.ell2cart(5/(180/np.pi), -90/(180/np.pi))
+    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 = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=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

@@ -1,14 +1,13 @@
+from typing import Tuple
+
 import numpy as np
-from Hansen_ES_CMA import escma
-from ellipsoide import EllipsoidTriaxial
-from numpy.typing import NDArray
 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
-
-ell_ES: EllipsoidTriaxial = None
-P_left: NDArray = None
-P_right: NDArray = None
+from ellipsoid_triaxial import EllipsoidTriaxial
 
 
 def Sehne(P1: NDArray, P2: NDArray) -> float:
@@ -19,57 +18,57 @@ def Sehne(P1: NDArray, P2: NDArray) -> float:
     :return: Bogenlänge s
     """
     R12 = P2-P1
-    s = np.linalg.norm(R12)
+    s = float(np.linalg.norm(R12))
+
     return s
 
 
-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)
-    """
-    global ell_ES, P_left, P_right
-
-    u, v = x
-    P_middle = ell_ES.para2cart(u, v)
-    d1 = Sehne(P_left, P_middle)
-    d2 = Sehne(P_middle, P_right)
-    base = d1 + d2
-
-    # midpoint penalty (dimensionslos)
-    # relative Differenz, skaliert stabil ü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
-
-
-def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, maxIter: int = 10000, all_points: bool = False):
+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: Parameter des triaxialen Ellipsoids
+    :param ell: Ellipsoid
     :param P0: Startpunkt
     :param Pk: Zielpunkt
     :param maxSegLen: maximale Segmentlänge
-    :param maxIter: maximale Durchläufe der Mittelpunktsgenerierung
-    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben werden kann
-    :return: Richtungswinkel des Start- und Zielpunktes und Gesamtlänge
+    :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben wird
+    :return: Richtungswinkel in RAD des Start- und Zielpunktes und Gesamtlänge
     """
-    global ell_ES
-    ell_ES = ell
+    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
+        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
@@ -82,48 +81,47 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
 
         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)
+            # print(dAB)
 
             if dAB > maxSegLen:
-                global P_left, P_right
+                # global P_left, P_right
                 P_left, P_right = A, B
-                Au, Av = ell_ES.cart2para(A)
-                Bu, Bv = ell_ES.cart2para(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)
+                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("Abbruch: maximale Iterationen überschritten.")
-
+                    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")
+        # 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_ES.point_onto_ellipsoid(P0 + 10.0 * (points[1] - P0) / np.linalg.norm(points[1] - P0))
+        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_ES, sigma0, P0)
+        alpha0 = sigma2alpha(ell, sigma0, P0)
 
-        p1i = ell_ES.point_onto_ellipsoid(Pk - 10.0 * (Pk - points[-2]) / np.linalg.norm(Pk - points[-2]))
+        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_ES, sigma1, Pk)
+        alpha1 = sigma2alpha(ell, sigma1, Pk)
     else:
         alpha0 = None
         alpha1 = None
@@ -166,7 +164,7 @@ if __name__ == '__main__':
 
     beta0, lamb0 = (0.2, 0.1)
     P0 = ell.ell2cart(beta0, lamb0)
-    beta1, lamb1 = (0.7, 0.3)
+    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)
@@ -175,8 +173,10 @@ if __name__ == '__main__':
         points_num.append(ell.ell2cart(beta, lamb))
     points_num = np.array(points_num)
 
-    alpha0, alpha1, s, points = gha2_ES(ell, P0, P1, all_points=True)
+    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_triaxial/gha1_ana.py

@@ -1,13 +1,15 @@
+import math
 from math import comb
 from typing import Tuple
 
 import numpy as np
-from numpy import sin, cos, arctan2
-from numpy._typing import NDArray
-from scipy.special import factorial as fact
+from numpy import arctan2, cos, sin
+from numpy.typing import NDArray
 
-from ellipsoide import EllipsoidTriaxial
+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]:
@@ -64,7 +66,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
             x_m.append(x_(m))
             y_m.append(y_(m))
             z_m.append(z_(m))
-        fact_m = fact(m)
+        fact_m = math.factorial(m)
 
         # 22-24
         a_m.append(x_m[m] / fact_m)
@@ -109,10 +111,10 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
     if alpha1 < 0:
         alpha1 += 2 * np.pi
 
-    return p1, alpha1
+    return p1, wrap_0_2pi(alpha1)
 
 
-def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 16) -> Tuple[NDArray, float]:
+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
@@ -131,6 +133,13 @@ def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, ma
 
     _, _, h = ell.cart2geod(point_end, "ligas3")
     if h > 1e-5:
-        raise Exception("Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
+        raise Exception("GHA1_ana: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
 
-    return point_end, alpha_end
+    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

@@ -1,11 +1,18 @@
-import numpy as np
-from ellipsoide import EllipsoidTriaxial
-from GHA_triaxial.gha1_ana import gha1_ana
-from GHA_triaxial.utils import func_sigma_ell, louville_constant
-import plotly.graph_objects as go
-import winkelumrechnungen as wu
+from typing import Tuple
 
-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]:
+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
@@ -20,6 +27,8 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
     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)
@@ -29,22 +38,33 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
         p1 = points[-1]
         alpha1 = alphas[-1]
 
-        sigma = func_sigma_ell(ell, p1, alpha1)
+        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-6
+        dalpha = 1e-9
         l2 = louville_constant(ell, p2, alpha1)
         dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
-        alpha2 = alpha1 + (l0 - l2) / dl_dalpha
-
+        if abs(dl_dalpha) < 1e-20:
+            alpha2 = alpha1 + 0
+        else:
+            alpha2 = alpha1 + (l0 - l2) / dl_dalpha
         points.append(p2)
-        alphas.append(alpha2)
+        alphas.append(wrap_0_2pi(alpha2))
 
         ds_step = np.linalg.norm(p2 - p1)
         s_curr += ds_step
-        if s_curr > 10000000:
-            pass
+        last_sigma = sigma
+        pass
 
     if all_points:
         return points[-1], alphas[-1], np.array(points), np.array(alphas)
@@ -78,11 +98,11 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray):
 
 
 if __name__ == '__main__':
-    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    P0 = ell.para2cart(0.2, 0.3)
-    alpha0 = wu.deg2rad(35)
-    s = 13000000
-    P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
-    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16)
-    show_points(points, P0, P1_ana)
-    print(np.linalg.norm(P1_app - P1_ana))
+    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

@@ -1,15 +1,16 @@
+from typing import Callable, List, Tuple
+
 import numpy as np
-from numpy import sin, cos, arctan2
-import ellipsoide
-import runge_kutta as rk
-import winkelumrechnungen as wu
-import GHA_triaxial.numeric_examples_karney as ne_karney
-from GHA_triaxial.gha1_ana import gha1_ana
-from ellipsoide import EllipsoidTriaxial
-from typing import Callable, Tuple, List
+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:
@@ -75,8 +76,11 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu
 
     alpha1 = arctan2(P, Q)
 
-    if alpha1 < 0:
-        alpha1 += 2 * np.pi
+    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
@@ -104,7 +108,7 @@ if __name__ == "__main__":
     # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
     # print(diffs_panou)
 
-    ell: EllipsoidTriaxial = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
+    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)

GHA_triaxial/gha2_approx.py

@@ -1,12 +1,13 @@
-import numpy as np
-from ellipsoide import EllipsoidTriaxial
-from GHA_triaxial.gha2_num import gha2_num
-import plotly.graph_objects as go
-import winkelumrechnungen as wu
-from numpy.typing import NDArray
 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]:

GHA_triaxial/gha2_num.py

@@ -1,433 +1,643 @@
-import numpy as np
-from ellipsoide import EllipsoidTriaxial
-import runge_kutta as rk
 from typing import Tuple
+
+import numpy as np
 from numpy.typing import NDArray
 
-from utils_angle import arccot, cot, wrap_to_pi
+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):
-    # sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß)
-    dlam = wrap_to_pi(lam2 - lam1)
+    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)  # (-pi, pi]
+    a = np.arctan2(y, x)
     if a < 0:
         a += 2 * np.pi
     return a
 
 
 # Panou 2013
-def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: 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]:
+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: triaxiales Ellipsoid
-    :param beta_1: reduzierte ellipsoidische Breite Punkt 1
-    :param lamb_1: elllipsoidische Länge Punkt 1
-    :param beta_2: reduzierte ellipsoidische Breite Punkt 2
-    :param lamb_2: elllipsoidische Länge Punkt 2
+    :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:
-    :param iter_max: Maximale Anzhal Iterationen
-    :param all_points:
-    :return:
+    :param epsilon: Genauigkeit
+    :param iter_max: Maximale Iterationen
+    :param all_points: Ausgabe aller Punkte
+    :return: Azimut Startpunkt, Azumit Zielpunkt, Strecke
     """
 
-    # h_x, h_y, h_e entsprechen E_x, E_y, E_e
+    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)
 
-        BETA = (ell.ay ** 2 * np.sin(beta) ** 2 + ell.b ** 2 * np.cos(beta) ** 2) / (
-                    ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2)
-        LAMBDA = (ell.ax ** 2 * np.sin(lamb) ** 2 + ell.ay ** 2 * np.cos(lamb) ** 2) / (
-                    ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2)
+        sb2 = sb * sb
+        cb2 = cb * cb
+        sl2 = sl * sl
+        cl2 = cl * cl
 
-        # Erste Ableitungen von ΒETA und LAMBDA
-        BETA_ = (ell.ax ** 2 * ell.Ey ** 2 * np.sin(2 * beta)) / (ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2
-        LAMBDA_ = - (ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 2
+        s2b = 2.0 * sb * cb
+        c2b = cb2 - sb2
+        s2l = 2.0 * sl * cl
+        c2l = cl2 - sl2
 
-        # Zweite Ableitungen von ΒETA und LAMBDA
-        BETA__ = ((2 * ell.ax ** 2 * ell.Ey ** 4 * np.sin(2 * beta) ** 2) / (
-                    ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 3) + (
-                             (2 * ell.ax ** 2 * ell.Ey ** 2 * np.cos(2 * beta)) / (
-                                 ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2)
-        LAMBDA__ = (((2 * ell.b ** 2 * ell.Ee ** 4 * np.sin(2 * lamb) ** 2) / (
-                    ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 3) -
-                    ((2 * ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (
-                                ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 2))
+        denB = Ex2 - Ey2 * sb2
+        denL = Ex2 - Ee2 * cl2
 
-        E = BETA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
-        F = 0
-        G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
+        BETA = (ay2 * sb2 + b2 * cb2) / denB
+        LAMBDA = (ax2 * sl2 + ay2 * cl2) / denL
 
-        # Erste Ableitungen von E und G
-        E_beta = BETA_ * (
-                    ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) - BETA * ell.Ey ** 2 * np.sin(
-            2 * beta)
-        E_lamb = BETA * ell.Ee ** 2 * np.sin(2 * lamb)
+        BETA_ = (ax2 * Ey2 * s2b) / (denB * denB)
+        LAMBDA_ = -(b2 * Ee2 * s2l) / (denL * denL)
 
-        G_beta = - LAMBDA * ell.Ey ** 2 * np.sin(2 * beta)
-        G_lamb = LAMBDA_ * (
-                    ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + LAMBDA * ell.Ee ** 2 * np.sin(
-            2 * lamb)
+        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
+        )
 
-        # Zweite Ableitungen von E und G
-        E_beta_beta = BETA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(
-            lamb) ** 2) - 2 * BETA_ * ell.Ey ** 2 * np.sin(2 * beta) - 2 * BETA * ell.Ey ** 2 * np.cos(2 * beta)
-        E_beta_lamb = BETA_ * ell.Ee ** 2 * np.sin(2 * lamb)
-        E_lamb_lamb = 2 * BETA * ell.Ee ** 2 * np.cos(2 * lamb)
+        Q = Ey2 * cb2 + Ee2 * sl2
 
-        G_beta_beta = - 2 * LAMBDA * ell.Ey ** 2 * np.cos(2 * beta)
-        G_beta_lamb = - LAMBDA_ * ell.Ey ** 2 * np.sin(2 * beta)
-        G_lamb_lamb = LAMBDA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(
-            lamb) ** 2) + 2 * LAMBDA_ * ell.Ee ** 2 * np.sin(2 * lamb) + 2 * LAMBDA * ell.Ee ** 2 * np.cos(2 * lamb)
+        E = BETA * Q
+        G = LAMBDA * Q
 
-        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)
+        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)
 
-        (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_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))
+        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)
-
-    def buildODElamb():
-        def ODE(lamb, v):
-            beta, beta_p, X3, X4 = v
-
-            (BETA, LAMBDA, E, G,
-             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])
-
-        return ODE
+        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)
 
-        (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_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))
+        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)
+        return BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00
 
-    def buildODEbeta():
-        def ODE(beta, v):
-            lamb, lamb_p, Y3, Y4 = v
+    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)
 
-            (BETA, LAMBDA, E, G,
-             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])
-
-        return ODE
-
-    if lamb_1 != lamb_2:
-        N = n
-        dlamb = lamb_2 - lamb_1
-        alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
+    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:
-            beta_0 = 0.0
-        else:
-            (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
-            beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph)
+            return None
 
-        ode_lamb = buildODElamb()
+        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)
 
-            for _ in range(iter_max):
-                startwerte = np.array([beta_1, beta_p0, 0.0, 1.0], dtype=float)
-                lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte, dlamb, N, False)
-
-                beta_end, beta_p_end, X3_end, X4_end = states[-1]
-                delta = beta_end - beta_2
+                beta_end, _, X3_end, _ = y_end
+                delta = beta_end - beta1
 
                 if abs(delta) < epsilon:
-                    return True, beta_p0, lamb_list, states
+                    return True, beta_p0
 
-                d_beta_end_d_beta0 = X3_end
-                if abs(d_beta_end_d_beta0) < 1e-20:
-                    return False, None, None, None
+                if abs(X3_end) < 1e-20:
+                    return False, None
 
-                step = delta / d_beta_end_d_beta0
-                max_step = 0.5
-                if abs(step) > max_step:
-                    step = np.sign(step) * max_step
+                step = delta / X3_end
+                step = float(np.clip(step, -0.5, 0.5))
+                beta_p0 -= step
 
-                beta_p0 = beta_p0 - step
-
-            return False, None, None, None
-
-        alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
-        (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
-        beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph)
-
-        guesses = [
-            beta_p0_sph,
-            0.5 * beta_p0_sph,
-            2.0 * beta_p0_sph,
-            -beta_p0_sph,
-            -0.5 * beta_p0_sph,
-        ]
+            return False, None
 
+        seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
         best = None
 
-        for g in guesses:
-            ok, beta_p0_sol, lamb_list_cand, states_cand = solve_newton(g)
+        for seed in seeds:
+            ok, sol = solve_newton(seed)
             if not ok:
                 continue
-
-            beta_arr_c = np.array([st[0] for st in states_cand], dtype=float)
-            beta_p_arr_c = np.array([st[1] for st in states_cand], dtype=float)
-            lamb_arr_c = np.array(lamb_list_cand, dtype=float)
-
-            integrand = np.zeros(N + 1)
-            for i in range(N + 1):
-                (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr_c[i], lamb_arr_c[i])
-                integrand[i] = np.sqrt(Ei * beta_p_arr_c[i] ** 2 + Gi)
-
-            h = abs(dlamb) / 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_cand = h / 3.0 * S
-            else:
-                s_cand = np.trapz(integrand, dx=h)
-
-            if (best is None) or (s_cand < best[0]):
-                best = (s_cand, beta_p0_sol, lamb_list_cand, states_cand)
+            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:
-            raise RuntimeError("Keine Multi-Start-Variante konvergiert.")
+            return None
 
-        s_best, beta_0, lamb_list, werte = best
+        return best[0], best[1], sgn, dlamb, ode_lamb
 
-        beta_arr = np.zeros(N + 1)
-        # lamb_arr = np.zeros(N + 1)
-        lamb_arr = np.array(lamb_list)
-        beta_p_arr = np.zeros(N + 1)
+    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
 
-        for i, state in enumerate(werte):
-            # lamb_arr[i] = state[0]
-            # beta_arr[i] = state[1]
-            # beta_p_arr[i] = state[2]
-            beta_arr[i] = state[0]
-            beta_p_arr[i] = state[1]
+        sgn = 1.0 if dbeta >= 0.0 else -1.0
 
-        (_, _, E1, G1,
-         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
-        (_, _, E2, G2,
-         *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+        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)
 
-        alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0])
-        alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1])
+            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)
 
-        integrand = np.zeros(N + 1)
-        for i in range(N + 1):
-            (_, _, Ei, Gi,
-             *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
-            integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi)
+        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)
 
-        h = abs(dlamb) / 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
+                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:
-            s = np.trapz(integrand, dx=h)
+            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)
 
-        beta0 = beta_arr[0]
-        lamb0 = lamb_arr[0]
-        c = np.sqrt(
-            (np.cos(beta0) ** 2 + (ell.Ee**2 / ell.Ex**2) * np.sin(beta0) ** 2) * np.sin(alpha_1) ** 2
-            + (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2
-        )
+        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:
-            return alpha_1, alpha_2, s, beta_arr, lamb_arr
-        else:
-            return alpha_1, alpha_2, s
+            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)
 
-    else:  # lamb_1 == lamb_2
-        N = n
-        dbeta = beta_2 - beta_1
+            (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+            (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
 
-        if abs(dbeta) < 1e-15:
-            if all_points:
-                return 0, 0, 0, np.array([]), np.array([])
+            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:
-                return 0, 0, 0
+                s = np.trapz(integrand, dx=h)
 
-        lamb_0 = 0
+            return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr
 
-        ode_beta = buildODEbeta()
+        _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N_full, integrand_lambda)
+        beta_end, beta_p_end, _, _ = y_end
 
-        for i in range(iter_max):
-            startwerte = [lamb_1, lamb_0, 0.0, 1.0]
+        (_, _, 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)
 
-            beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
+        (_, _, 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)
 
-            beta_end = beta_list[-1]
-            lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
+        return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s)
 
-            d_lamb_end_d_lambda0 = Y3_end
-            delta = lamb_end - lamb_2
+    # Fall 2 (lambda_0 == lambda_1)
+    N = int(n)
+    dbeta = float(beta_1 - beta_0)
 
-            if abs(delta) < epsilon:
-                break
+    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
 
-            if abs(d_lamb_end_d_lambda0) < 1e-20:
-                raise RuntimeError("Abbruch (Ableitung ~ 0).")
+    sgn = 1.0 if dbeta >= 0.0 else -1.0
 
-            max_step = 1.0
-            step = delta / d_lamb_end_d_lambda0
-            if abs(step) > max_step:
-                step = np.sign(step) * max_step
+    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)
 
-            lamb_0 = lamb_0 - step
+        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)
 
-        beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False)
+    lamb_p0 = float(best_fast[3]) if (best_fast[0] == "beta" and best_fast[3] is not None) else 0.0
 
-        # beta_arr = np.zeros(N + 1)
-        beta_arr = np.array(beta_list)
-        lamb_arr = np.zeros(N + 1)
-        lambda_p_arr = np.zeros(N + 1)
+    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)
 
-        for i, state in enumerate(werte):
-            # beta_arr[i] = state[0]
-            # lamb_arr[i] = state[1]
-            # lambda_p_arr[i] = state[2]
-            lamb_arr[i] = state[0]
-            lambda_p_arr[i] = state[1]
+        lamb_end, _, Y3_end, _ = y_end
+        delta = lamb_end - lamb1
 
-        # Azimute
-        (BETA1, LAMBDA1, E1, G1,
-         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
-        (BETA2, LAMBDA2, E2, G2,
-         *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+        if abs(delta) < epsilon:
+            break
 
-        alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lambda_p_arr[0])
-        alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lambda_p_arr[-1])
+        if abs(Y3_end) < 1e-20:
+            raise RuntimeError("GHA2_num: Ableitung ~ 0 im beta-Fall")
 
-        integrand = np.zeros(N + 1)
+        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 * lambda_p_arr[i] ** 2)
+            (_, _, 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 = 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)
 
-        if all_points:
-            return alpha_1, alpha_2, s, beta_arr, lamb_arr
-        else:
-            return alpha_1, alpha_2, s
+        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("Fiction")
-    # # beta1 = np.deg2rad(75)
-    # # lamb1 = np.deg2rad(-90)
-    # # beta2 = np.deg2rad(75)
-    # # lamb2 = np.deg2rad(66)
-    # # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
-    # # print(aus.gms("a1", a1, 4))
-    # # 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)
+    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)
 
-    # 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)
+    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

@@ -1,6 +1,12 @@
 import random
+from typing import List
+
 import winkelumrechnungen as wu
-from typing import List, Tuple
+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:
     """
@@ -26,7 +32,7 @@ def get_random_examples(num: int, seed: int = None) -> List:
     """
     if seed is not None:
         random.seed(seed)
-    with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
+    with open(file_path) as datei:
         lines = datei.readlines()
     examples = []
     for i in range(num):
@@ -41,7 +47,7 @@ def get_examples(l_i: List) -> List:
     :param l_i: Liste von Indizes
     :return: Liste mit Beispielen
     """
-    with open("Karney_2024_Testset.txt") as datei:
+    with open(file_path) as datei:
         lines = datei.readlines()
     examples = []
     for i in l_i:
@@ -50,5 +56,63 @@ def get_examples(l_i: List) -> List:
     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__":
-    get_random_examples(10)
+    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/utils.py

@@ -1,11 +1,13 @@
+from __future__ import annotations
+
 from typing import Tuple
 
 import numpy as np
-from numpy import arctan2, sin, cos, sqrt
-from numpy._typing import NDArray
+from numpy import arctan2, cos, sin, sqrt
 from numpy.typing import NDArray
 
-from ellipsoide import EllipsoidTriaxial
+from ellipsoid_triaxial import EllipsoidTriaxial
+from utils_angle import wrap_0_2pi
 
 
 def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float:
@@ -21,7 +23,7 @@ def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float
     Q = float(q @ sigma)
 
     alpha = arctan2(P, Q)
-    return alpha
+    return wrap_0_2pi(alpha)
 
 
 def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]:
@@ -43,10 +45,10 @@ def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float
     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-12:
-        raise Exception("Alpha Umrechnung fehlgeschlagen")
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
+        raise Exception("alpha_para2ell: Differenz in den Richtungsableitungen")
 
-    return beta, lamb, alpha_ell
+    return beta, lamb, wrap_0_2pi(alpha_ell)
 
 
 def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]:
@@ -68,10 +70,10 @@ def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell:
     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-9:
-        raise Exception("Alpha Umrechnung fehlgeschlagen")
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
+        raise Exception("alpha_ell2para: Differenz in den Richtungsableitungen")
 
-    return u, v, alpha_para
+    return u, v, wrap_0_2pi(alpha_para)
 
 
 def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha_ell: float) -> NDArray:
@@ -124,25 +126,26 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
     :param point: Punkt
     :return: p und q
     """
-    x, y, z = point
     n = ell.func_n(point)
 
     beta, lamb = ell.cart2ell(point)
-    B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2
-    L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2
+    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
 
-    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 + sqrt(c1 ** 2 - 4 * c0)) / 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)
+        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]])
@@ -171,13 +174,28 @@ def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
                   q[2] * n[0] - q[0] * n[2],
                   q[0] * n[1] - q[1] * n[0]])
 
-    t1 = np.dot(n, q)
-    t2 = np.dot(n, p)
-    t3 = np.dot(p, q)
-    if not (t1 < 1e-10 or t1 > 1-1e-10) and not (t2 < 1e-10 or t2 > 1-1e-10) and not (t3 < 1e-10 or t3 > 1-1e-10):
-        raise Exception("Fehler in den normierten Vektoren")
-
     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

Tests/algorithms_test.ipynb

@@ -1,846 +0,0 @@
-{
- "cells": [
-  {
-   "cell_type": "code",
-   "id": "initial_id",
-   "metadata": {
-    "collapsed": true,
-    "ExecuteTime": {
-     "end_time": "2026-02-04T17:56:25.185202Z",
-     "start_time": "2026-02-04T17:56:22.991833Z"
-    }
-   },
-   "source": [
-    "%load_ext autoreload\n",
-    "%autoreload 2"
-   ],
-   "outputs": [],
-   "execution_count": 1
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T17:56:35.920624Z",
-     "start_time": "2026-02-04T17:56:25.201916Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "%reload_ext autoreload\n",
-    "%autoreload 2\n",
-    "import time\n",
-    "from numpy import nan\n",
-    "import numpy as np\n",
-    "import math\n",
-    "import winkelumrechnungen as wu\n",
-    "import os\n",
-    "from contextlib import contextmanager, redirect_stdout, redirect_stderr\n",
-    "import plotly.graph_objects as go\n",
-    "import warnings\n",
-    "import pickle\n",
-    "import random\n",
-    "\n",
-    "from ellipsoide import EllipsoidTriaxial\n",
-    "from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n",
-    "\n",
-    "from GHA_triaxial.gha1_num import gha1_num\n",
-    "from GHA_triaxial.gha1_ana import gha1_ana\n",
-    "from GHA_triaxial.gha1_ES import gha1_ES\n",
-    "from GHA_triaxial.gha1_approx import gha1_approx\n",
-    "\n",
-    "from GHA_triaxial.gha2_num import gha2_num\n",
-    "from GHA_triaxial.gha2_ES import gha2_ES\n",
-    "from GHA_triaxial.gha2_approx import gha2_approx\n",
-    "\n",
-    "from GHA_triaxial.numeric_examples_panou import get_tables as get_tables_panou\n",
-    "from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney"
-   ],
-   "id": "961cb22764c5bcb9",
-   "outputs": [],
-   "execution_count": 2
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T17:56:36.470974Z",
-     "start_time": "2026-02-04T17:56:35.937646Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "@contextmanager\n",
-    "def suppress_print():\n",
-    "    with open(os.devnull, 'w') as fnull:\n",
-    "        with redirect_stdout(fnull), redirect_stderr(fnull):\n",
-    "            yield"
-   ],
-   "id": "e4740625a366ac13",
-   "outputs": [],
-   "execution_count": 3
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T18:11:14.405810Z",
-     "start_time": "2026-02-04T18:11:14.169636Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "# dsPart = [60, 125, 600, 1250, 6000, 60000] entspricht bei der Erde ca. 100km, 50km, 10km, 5km, 1km, 100m\n",
-    "\n",
-    "steps_gha1_num = [20, 50, 100, 200, 500, 1000, 5000]\n",
-    "maxM_gha1_ana = [20, 50]\n",
-    "parts_gha1_ana = [2, 8]\n",
-    "dsPart_gha1_ES = [60, 600, 1250]\n",
-    "dsPart_gha1_approx = [600, 1250, 6000]\n",
-    "\n",
-    "steps_gha2_num = [200, 500, 1000]\n",
-    "dsPart_gha2_ES = [60, 600, 1250]\n",
-    "dsPart_gha2_approx = [600, 1250, 6000]"
-   ],
-   "id": "96093cdde03f8d57",
-   "outputs": [],
-   "execution_count": 16
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T18:11:15.988345Z",
-     "start_time": "2026-02-04T18:11:15.774321Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "# test = \"Karney\"\n",
-    "test = \"Panou\"\n",
-    "# test = \"Random\"\n",
-    "if test == \"Karney\":\n",
-    "    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n",
-    "    examples = get_examples_karney(4, seed=42)\n",
-    "elif test == \"Panou\":\n",
-    "    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n",
-    "    tables = get_tables_panou()\n",
-    "    table_indices = []\n",
-    "    examples = []\n",
-    "    for i, table in enumerate(tables):\n",
-    "        for example in table:\n",
-    "            table_indices.append(i+1)\n",
-    "            examples.append(example)\n",
-    "elif test == \"Random\":\n",
-    "    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n",
-    "    examples = [] # beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s\n",
-    "    random.seed(42)\n",
-    "    for _ in range(10):\n",
-    "        beta0 = wu.deg2rad(random.randint(-90, 90))\n",
-    "        lamb0 = wu.deg2rad(random.randint(-179, 180))\n",
-    "        alpha0_ell = wu.deg2rad(random.randint(0, 359))\n",
-    "        s = random.randint(10000, int(np.pi*ell.b))\n",
-    "        examples.append([beta0, lamb0, alpha0_ell, s])\n",
-    "    pass"
-   ],
-   "id": "6e384cc01c2dbe",
-   "outputs": [],
-   "execution_count": 17
-  },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "outputs": [],
-   "execution_count": null,
-   "source": [
-    "if test != \"Random\":\n",
-    "    results = {}\n",
-    "    for i, example in enumerate(examples):\n",
-    "        print(f\"----- Beispiel {i+1}/{len(examples)}\")\n",
-    "        example_results = {}\n",
-    "\n",
-    "        beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n",
-    "        P0 = ell.ell2cart(beta0, lamb0)\n",
-    "        P1 = ell.ell2cart(beta1, lamb1)\n",
-    "        _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
-    "\n",
-    "        for steps in steps_gha1_num:\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n",
-    "                print(f\"GHA1_num_{steps}\", P1_num)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_num, lamb1_num = ell.cart2ell(P1_num)\n",
-    "                d_beta1 = abs(beta1_num - beta1)\n",
-    "                d_lamb1 = abs(lamb1_num - lamb1)\n",
-    "                d_alpha1 = abs(alpha1_num_1 - alpha1_ell)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for maxM in maxM_gha1_ana:\n",
-    "            for parts in parts_gha1_ana:\n",
-    "                start = time.perf_counter()\n",
-    "                try:\n",
-    "                    P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n",
-    "                    print(f\"GHA1_ana_{maxM}_{parts}\", P1_ana)\n",
-    "                    end = time.perf_counter()\n",
-    "                    beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
-    "                    _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
-    "                    d_beta1 = abs(beta1_ana - beta1)\n",
-    "                    d_lamb1 = abs(lamb1_ana - lamb1)\n",
-    "                    d_alpha1 = abs(alpha1_ana_ell - alpha1_ell)\n",
-    "                    d_time = end - start\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "                except Exception as e:\n",
-    "                    print(e)\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha1_ES:\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_ES, alpha1_ES = gha1_ES(ell, beta0, lamb0, alpha0_ell, s, maxSegLen=dsPart)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_ES, lamb1_ES = ell.cart2ell(P1_ES)\n",
-    "                d_beta1 = abs(beta1_ES - beta1)\n",
-    "                d_lamb1 = abs(lamb1_ES - lamb1)\n",
-    "                d_alpha1 = abs(alpha1_ES - alpha1_ell)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_ES_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_ES_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha1_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n",
-    "                d_beta1 = abs(beta1_approx - beta1)\n",
-    "                d_lamb1 = abs(lamb1_approx - lamb1)\n",
-    "                d_alpha1 = abs(alpha1_approx - alpha1_ell)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        # ----------------------------------------------\n",
-    "\n",
-    "        for steps in steps_gha2_num:\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                with warnings.catch_warnings():\n",
-    "                    warnings.simplefilter(\"ignore\", RuntimeWarning)\n",
-    "                    alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n",
-    "                    print(alpha0_num, alpha1_num_2, s_num)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(alpha0_num - alpha0_ell)\n",
-    "                d_alpha1 = abs(alpha1_num_2 - alpha1_ell)\n",
-    "                d_s = abs(s_num - s)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha2_ES:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                with suppress_print():\n",
-    "                    alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(alpha0_ES - alpha0_ell)\n",
-    "                d_alpha1 = abs(alpha1_ES - alpha1_ell)\n",
-    "                d_s = abs(s_ES - s)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha2_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(alpha0_approx - alpha0_ell)\n",
-    "                d_alpha1 = abs(alpha1_approx - alpha1_ell)\n",
-    "                d_s = abs(s_approx - s)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results"
-   ],
-   "id": "2a87e028089a215"
-  },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "outputs": [],
-   "execution_count": null,
-   "source": [
-    "if test == \"Random\":\n",
-    "    results = {}\n",
-    "    for i, example in enumerate(examples):\n",
-    "        print(f\"----- Beispiel {i+1}/{len(examples)}\")\n",
-    "        example_results = {}\n",
-    "\n",
-    "        beta0, lamb0, alpha0_ell, s = example\n",
-    "        P0 = ell.ell2cart(beta0, lamb0)\n",
-    "        _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
-    "\n",
-    "        # P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=60, maxPartCircum=4)\n",
-    "        # beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
-    "        # _, _, alpha1_ana = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
-    "\n",
-    "        P1, alpha1_ell = gha1_num(ell, P0, alpha0_ell, s, num=10000)\n",
-    "        beta1, lamb1 = ell.cart2ell(P1)\n",
-    "\n",
-    "        for maxM in maxM_gha1_ana:\n",
-    "            for parts in parts_gha1_ana:\n",
-    "                start = time.perf_counter()\n",
-    "                try:\n",
-    "                    P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n",
-    "                    end = time.perf_counter()\n",
-    "                    beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
-    "                    _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
-    "                    d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n",
-    "                    d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n",
-    "                    d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n",
-    "                    d_time = end - start\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "                except Exception as e:\n",
-    "                    print(e)\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha1_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n",
-    "                d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n",
-    "                d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        # for steps in steps_gha2_num:\n",
-    "        #     start = time.perf_counter()\n",
-    "        #     try:\n",
-    "        #         with warnings.catch_warnings():\n",
-    "        #             warnings.simplefilter(\"ignore\", RuntimeWarning)\n",
-    "        #             alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n",
-    "        #         end = time.perf_counter()\n",
-    "        #         d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n",
-    "        #         d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n",
-    "        #         d_s = abs(s_num - s) / 1000\n",
-    "        #         d_time = end - start\n",
-    "        #         example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "        #     except Exception as e:\n",
-    "        #         print(e)\n",
-    "        #         example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n",
-    "        #\n",
-    "        # for dsPart in dsPart_gha2_ES:\n",
-    "        #     ds = ell.ax/dsPart\n",
-    "        #     start = time.perf_counter()\n",
-    "        #     try:\n",
-    "        #         with suppress_print():\n",
-    "        #             alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n",
-    "        #         end = time.perf_counter()\n",
-    "        #         d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n",
-    "        #         d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n",
-    "        #         d_s = abs(s_ES - s) / 1000\n",
-    "        #         d_time = end - start\n",
-    "        #         example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "        #     except Exception as e:\n",
-    "        #         print(e)\n",
-    "        #         example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha2_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
-    "                d_s = abs(s_approx - s) / 1000\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results\n",
-    "\n",
-    "\n"
-   ],
-   "id": "d8868b5f0e6f9b42"
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T14:35:09.900827Z",
-     "start_time": "2026-02-04T14:35:09.690955Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "# with open(f\"gha_results{test}.pkl\", \"wb\") as f:\n",
-    "#     pickle.dump(results, f)"
-   ],
-   "id": "664105ea35d50a7b",
-   "outputs": [],
-   "execution_count": 7
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T14:35:10.113762Z",
-     "start_time": "2026-02-04T14:35:09.910687Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "# with open(f\"gha_results{test}.pkl\", \"rb\") as f:\n",
-    "#     results = pickle.load(f)"
-   ],
-   "id": "ad8aa9dcc8af4a05",
-   "outputs": [],
-   "execution_count": 8
-  },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "outputs": [],
-   "execution_count": null,
-   "source": [
-    "metrics_gha1 = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n",
-    "metrics_gha2 = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n",
-    "listed_results = {}\n",
-    "for example, example_metrics in results.items():\n",
-    "    for method, method_metrics in example_metrics.items():\n",
-    "        if \"GHA1\" in method:\n",
-    "            if method not in listed_results.keys():\n",
-    "                listed_results[method] = {metric: [] for metric in metrics_gha1}\n",
-    "            for i, metric in enumerate(method_metrics):\n",
-    "                if '[\"]' in metrics_gha1[i]:\n",
-    "                    listed_results[method][metrics_gha1[i]].append(wu.rad2deg(metric)*3600)\n",
-    "                else:\n",
-    "                    listed_results[method][metrics_gha1[i]].append(metric)\n",
-    "        if \"GHA2\" in method:\n",
-    "            if method not in listed_results.keys():\n",
-    "                listed_results[method] = {metric: [] for metric in metrics_gha2}\n",
-    "            for i, metric in enumerate(method_metrics):\n",
-    "                if '[\"]' in metrics_gha2[i]:\n",
-    "                    listed_results[method][metrics_gha2[i]].append(wu.rad2deg(metric)*3600)\n",
-    "                else:\n",
-    "                    listed_results[method][metrics_gha2[i]].append(metric)\n",
-    "        pass"
-   ],
-   "id": "8fd6f220440f4494"
-  },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "outputs": [],
-   "execution_count": null,
-   "source": [
-    "def format_max(values, is_angle=False):\n",
-    "    arr = np.array(values, dtype=float)\n",
-    "    if arr.size==0:\n",
-    "        return np.nan\n",
-    "    maxi = np.nanmax(np.abs(arr))\n",
-    "    if maxi is None or (isinstance(maxi,float) and (math.isnan(maxi))):\n",
-    "        return \"nan\"\n",
-    "    if is_angle:\n",
-    "        maxi = wu.rad2deg(maxi)*3600\n",
-    "    if f\"{maxi:.3g}\" == 0:\n",
-    "        pass\n",
-    "    return f\"{maxi:.3g}\"\n",
-    "\n",
-    "\n",
-    "def build_max_table(gha_prefix, title, group_value = None):\n",
-    "    if gha_prefix==\"GHA1\":\n",
-    "        metrics = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n",
-    "        angle_mask = [True, True, True, False]\n",
-    "    else:\n",
-    "        metrics = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n",
-    "        angle_mask = [True, True, False, False]\n",
-    "\n",
-    "    if group_value is None:\n",
-    "        example_keys = [example_key for example_key in list(results.keys())]\n",
-    "    else:\n",
-    "        example_keys = [example_key for example_key, group_index in zip(results.keys(), table_indices) if group_index==group_value]\n",
-    "\n",
-    "    algorithms = sorted({algorithm for example_key in example_keys for algorithm in results[example_key].keys() if algorithm.startswith(gha_prefix)})\n",
-    "\n",
-    "    header = [\"Algorithmus\", \"Parameter\", \"NaN\"] + list(metrics)\n",
-    "    cells = [[] for i in range(len(metrics) + 3)]\n",
-    "    for algorithm in algorithms:\n",
-    "\n",
-    "        ghaNr, variant, params = algorithm.split(\"_\", 2)\n",
-    "        cells[0].append(variant)\n",
-    "        cells[1].append(params)\n",
-    "        nan_values = []\n",
-    "        for example_key in example_keys:\n",
-    "            nan_values.append(results[example_key][algorithm][0])\n",
-    "        cells[2].append(np.sum(np.isnan(nan_values)))\n",
-    "        for i, metric in enumerate(metrics):\n",
-    "            values = []\n",
-    "            for example_key in example_keys:\n",
-    "                values.append(results[example_key][algorithm][i])\n",
-    "            cells[i+3].append(format_max(values, is_angle=angle_mask[i]))\n",
-    "\n",
-    "    header = dict(\n",
-    "        values=header,\n",
-    "        align=\"center\",\n",
-    "        fill_color=\"lightgrey\",\n",
-    "        font=dict(size=13)\n",
-    "    )\n",
-    "    cells = dict(\n",
-    "        values=cells,\n",
-    "        align=\"center\"\n",
-    "    )\n",
-    "\n",
-    "    fig = go.Figure(data=[go.Table(header=header, cells=cells)])\n",
-    "    fig.update_layout(title=title,\n",
-    "    template=\"simple_white\",\n",
-    "    width=800,\n",
-    "    height=280,\n",
-    "    margin=dict(l=20, r=20, t=60, b=20))\n",
-    "    return fig\n",
-    "\n",
-    "figs = []\n",
-    "if test == \"Panou\":\n",
-    "    for table_index in sorted(set(table_indices)):\n",
-    "        fig1 = build_max_table(\"GHA1\", f\"{test} - Gruppe {table_index} - GHA1\", table_index)\n",
-    "        fig2 = build_max_table(\"GHA2\", f\"{test} - Gruppe {table_index} - GHA2\", table_index)\n",
-    "        figs.append(fig1)\n",
-    "        figs.append(fig2)\n",
-    "elif test == \"Karney\":\n",
-    "    fig1 = build_max_table(\"GHA1\", f\"{test} - GHA1\")\n",
-    "    fig2 = build_max_table(\"GHA2\", f\"{test} - GHA2\")\n",
-    "    figs.append(fig1)\n",
-    "    figs.append(fig2)\n",
-    "\n",
-    "for fig in figs:\n",
-    "    fig.show()"
-   ],
-   "id": "f1fb73407f5b1c8a"
-  },
-  {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-02-04T18:08:09.603687Z",
-     "start_time": "2026-02-04T18:08:09.269698Z"
-    }
-   },
-   "cell_type": "code",
-   "source": [
-    "from collections import defaultdict\n",
-    "\n",
-    "def to_latex_sci(x, sig=3):\n",
-    "    \"\"\"\n",
-    "    Formatiert eine Zahl als LaTeX, z.B. 8.02e-08 -> $8.02\\\\cdot10^{-8}$.\n",
-    "    Für normale Zahlen -> $0.216$.\n",
-    "    Für nan -> nan (ohne $...$).\n",
-    "    \"\"\"\n",
-    "    if x is None:\n",
-    "        return \"nan\"\n",
-    "    try:\n",
-    "        xf = float(x)\n",
-    "    except Exception:\n",
-    "        return str(x)\n",
-    "\n",
-    "    if math.isnan(xf):\n",
-    "        return \"nan\"\n",
-    "\n",
-    "    # nahe 0 -> 0\n",
-    "    if xf == 0.0:\n",
-    "        return \"$0$\"\n",
-    "\n",
-    "    s = f\"{xf:.{sig}g}\"  # z.B. '8.02e-08' oder '0.216'\n",
-    "    if \"e\" in s or \"E\" in s:\n",
-    "        mant, exp = s.lower().split(\"e\")\n",
-    "        exp = int(exp)\n",
-    "        return f\"${mant}\\\\cdot10^{{{exp}}}$\"\n",
-    "    else:\n",
-    "        return f\"${s}$\"\n",
-    "\n",
-    "\n",
-    "def nan_count_for_algorithm(results, example_keys, algorithm):\n",
-    "    vals = []\n",
-    "    for k in example_keys:\n",
-    "        vals.extend(results[k][algorithm])\n",
-    "    return int(np.sum(np.isnan(np.array(vals, dtype=float))))\n",
-    "\n",
-    "\n",
-    "def max_abs_metric(results, example_keys, algorithm, metric_index, is_angle=False, wu=None):\n",
-    "    \"\"\"\n",
-    "    Echter Max(|...|) über alle example_keys.\n",
-    "    Wenn is_angle=True: Werte werden als rad angenommen und in Bogensekunden umgerechnet.\n",
-    "    \"\"\"\n",
-    "    arr = []\n",
-    "    for k in example_keys:\n",
-    "        arr.append(results[k][algorithm][metric_index])\n",
-    "    arr = np.array(arr, dtype=float)\n",
-    "\n",
-    "    if arr.size == 0:\n",
-    "        return np.nan\n",
-    "    m = np.nanmax(np.abs(arr))\n",
-    "    if is_angle:\n",
-    "        if wu is None:\n",
-    "            raise ValueError(\"wu wird benötigt für rad2deg, wenn is_angle=True\")\n",
-    "        m = wu.rad2deg(m) * 3600.0\n",
-    "    return m\n",
-    "\n",
-    "\n",
-    "def parse_variant_params(algorithm_key):\n",
-    "    \"\"\"\n",
-    "    Erwartet: GHA1_ana_20_4  -> variant='ana', params='20_4'\n",
-    "    Gibt zusätzlich eine hübsche Parameterdarstellung zurück.\n",
-    "    \"\"\"\n",
-    "    ghaNr, variant, params = algorithm_key.split(\"_\", 2)\n",
-    "\n",
-    "    # Standard: params 그대로\n",
-    "    pretty = params\n",
-    "\n",
-    "    # Für deine Tabellen: ana hat z.B. '20_4' -> '4, 20' (Segmentgröße, Ordnung)\n",
-    "    if variant == \"ana\":\n",
-    "        # bei dir: params = '20_4' oder '50_8' -> Ordnung_Segment\n",
-    "        # du willst aber: Segment, Ordnung -> '4, 20'\n",
-    "        a, b = params.split(\"_\")\n",
-    "        order = a\n",
-    "        seg = b\n",
-    "        pretty = f\"{seg}, {order}\"\n",
-    "    else:\n",
-    "        # num: '2000' bleibt '2000'\n",
-    "        # approx/ES: oft eine Zahl mit Punkt/Unterstrich? -> 그대로\n",
-    "        pretty = params.replace(\"_\", \", \")\n",
-    "\n",
-    "    return variant, params, pretty\n",
-    "\n",
-    "\n",
-    "def build_latex_table_from_results(\n",
-    "    results,\n",
-    "    gha_prefix=\"GHA1\",\n",
-    "    example_keys=None,\n",
-    "    caption=\"\",\n",
-    "    label=\"tab:results_algorithms\",\n",
-    "    include_nan_col=False,\n",
-    "    wu=None\n",
-    "):\n",
-    "    \"\"\"\n",
-    "    Erzeugt LaTeX Tabular (inkl. table-Umgebung) im gewünschten Stil.\n",
-    "    \"\"\"\n",
-    "\n",
-    "    if example_keys is None:\n",
-    "        example_keys = list(results.keys())\n",
-    "\n",
-    "    # Metriken & Winkel-Maske wie bei dir\n",
-    "    if gha_prefix == \"GHA1\":\n",
-    "        metric_headers = [\n",
-    "            r\"$\\max(|\\Delta \\beta|)$ [$''$]\",\n",
-    "            r\"$\\max(|\\Delta \\lambda|)$ [$''$]\",\n",
-    "            r\"$\\max(|\\Delta \\alpha_1|)$ [$''$]\",\n",
-    "            r\"time [s]\"\n",
-    "        ]\n",
-    "        angle_mask = [True, True, True, False]\n",
-    "    else:\n",
-    "        metric_headers = [\n",
-    "            r\"$\\max(|\\Delta \\alpha_0|)$ [$''$]\",\n",
-    "            r\"$\\max(|\\Delta \\alpha_1|)$ [$''$]\",\n",
-    "            r\"$\\max(|\\Delta s|)$ [m]\",\n",
-    "            r\"time [s]\"\n",
-    "        ]\n",
-    "        angle_mask = [True, True, False, False]\n",
-    "\n",
-    "    # Alle Algorithmen sammeln\n",
-    "    algorithms = sorted({\n",
-    "        alg for k in example_keys\n",
-    "        for alg in results[k].keys()\n",
-    "        if alg.startswith(gha_prefix)\n",
-    "    })\n",
-    "\n",
-    "    # Gruppieren nach variant (ana, num, approx, ES, ...)\n",
-    "    grouped = defaultdict(list)\n",
-    "    for alg in algorithms:\n",
-    "        variant, raw_params, pretty = parse_variant_params(alg)\n",
-    "        grouped[variant].append((alg, pretty))\n",
-    "\n",
-    "    # gewünschte Reihenfolge (falls vorhanden)\n",
-    "    variant_order = [\"ana\", \"num\", \"approx\", \"ES\"]\n",
-    "    ordered_variants = [v for v in variant_order if v in grouped] + [v for v in grouped.keys() if v not in variant_order]\n",
-    "\n",
-    "    # LaTeX Header\n",
-    "    cols = 2 + (1 if include_nan_col else 0) + len(metric_headers)\n",
-    "    colspec = \"|\" + \"|\".join([\"c\"] * cols) + \"|\"\n",
-    "\n",
-    "    header_cells = [\"Methode\", \"Parameterwerte\"]\n",
-    "    if include_nan_col:\n",
-    "        header_cells.append(\"NaN\")\n",
-    "    header_cells += metric_headers\n",
-    "\n",
-    "    lines = []\n",
-    "    lines.append(r\"\\begin{table}[H]\")\n",
-    "    lines.append(r\"\\centering\")\n",
-    "    if caption:\n",
-    "        lines.append(rf\"\\caption{{{caption}}}\")\n",
-    "    lines.append(rf\"\\label{{{label}}}\")\n",
-    "    lines.append(\"\")\n",
-    "    lines.append(rf\"\\begin{{tabular}}{{{colspec}}}\")\n",
-    "    lines.append(r\"\\hline\")\n",
-    "    lines.append(\" & \".join(header_cells) + r\" \\\\\")\n",
-    "    lines.append(r\"\\Xhline{1.5pt}\")\n",
-    "\n",
-    "    # Zeilen bauen\n",
-    "    for variant in ordered_variants:\n",
-    "        rows = grouped[variant]\n",
-    "        # für stable output: sort by pretty param (numerisch)\n",
-    "        # (du kannst das ändern, wenn du eine andere Reihenfolge willst)\n",
-    "        def sort_key(t):\n",
-    "            _, pretty = t\n",
-    "            # versuche numerisch zu sortieren\n",
-    "            try:\n",
-    "                parts = [float(p.strip()) for p in pretty.split(\",\")]\n",
-    "                return parts\n",
-    "            except Exception:\n",
-    "                return [pretty]\n",
-    "        rows = sorted(rows, key=sort_key)\n",
-    "\n",
-    "        multi_n = len(rows)\n",
-    "        method_name = {\n",
-    "            \"ana\": \"analytisch\",\n",
-    "            \"num\": \"numerisch\",\n",
-    "            \"approx\": \"approximiert\"\n",
-    "        }.get(variant, variant)\n",
-    "\n",
-    "        for idx, (alg, pretty_param) in enumerate(rows):\n",
-    "            row_cells = []\n",
-    "\n",
-    "            if multi_n > 1:\n",
-    "                if idx == 0:\n",
-    "                    row_cells.append(rf\"\\multirow{{{multi_n}}}{{*}}{{{method_name}}}\")\n",
-    "                else:\n",
-    "                    row_cells.append(\"\")  # Multirow fortsetzen\n",
-    "            else:\n",
-    "                row_cells.append(method_name)\n",
-    "\n",
-    "            row_cells.append(pretty_param)\n",
-    "\n",
-    "            if include_nan_col:\n",
-    "                nans = nan_count_for_algorithm(results, example_keys, alg)\n",
-    "                row_cells.append(str(nans))\n",
-    "\n",
-    "            # Metriken\n",
-    "            for mi in range(len(metric_headers)):\n",
-    "                m = max_abs_metric(results, example_keys, alg, mi, is_angle=angle_mask[mi], wu=wu)\n",
-    "                row_cells.append(to_latex_sci(m, sig=3))\n",
-    "\n",
-    "            # Zeile + Linienlogik wie in deinem Beispiel\n",
-    "            latex_row = \" & \".join(row_cells) + r\" \\\\\"\n",
-    "\n",
-    "            # nach letzter Zeile einer Methode eine \\hline (wie bei dir)\n",
-    "            if idx == multi_n - 1:\n",
-    "                latex_row += r\"\\hline\"\n",
-    "            lines.append(latex_row)\n",
-    "\n",
-    "    lines.append(r\"\\end{tabular}\")\n",
-    "    lines.append(\"\")\n",
-    "    lines.append(r\"\\end{table}\")\n",
-    "\n",
-    "    return \"\\n\".join(lines)\n",
-    "\n",
-    "\n",
-    "# --- Beispielaufruf ---\n",
-    "example_keys = list(key for i, key in enumerate(results.keys()) if table_indices[i] == 3)  # oder gefiltert auf Panou Gruppe etc.\n",
-    "latex = build_latex_table_from_results(\n",
-    "    results,\n",
-    "    gha_prefix=\"GHA1\",\n",
-    "    example_keys=example_keys,\n",
-    "    caption=r\"Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid (Gruppe 3)\",\n",
-    "    label=\"tab:results_gha1_g1\",\n",
-    "    include_nan_col=False,  # True, wenn du NaN-Spalte willst\n",
-    "    wu=wu  # wichtig, falls Winkelwerte rad sind\n",
-    ")\n",
-    "print(latex)"
-   ],
-   "id": "98de5e2a4f419483",
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\\begin{table}[H]\n",
-      "\\centering\n",
-      "\\caption{Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid}\n",
-      "\\label{tab:results_algorithms}\n",
-      "\n",
-      "\\begin{tabular}{|c|c|c|c|c|c|}\n",
-      "\\hline\n",
-      "Methode & Parameterwerte & $\\max(|\\Delta \\alpha_0|)$ [$''$] & $\\max(|\\Delta \\alpha_1|)$ [$''$] & $\\max(|\\Delta s|)$ [m] & time [s] \\\\\n",
-      "\\Xhline{1.5pt}\n",
-      "numerisch & 2000 & nan & nan & nan & nan \\\\\\hline\n",
-      "approximiert & 600 & $9.95\\cdot10^{3}$ & $2.77\\cdot10^{3}$ & $0.000888$ & $0.108$ \\\\\\hline\n",
-      "ES & 20 & $1.07\\cdot10^{4}$ & $4.35\\cdot10^{3}$ & $0.000886$ & $1.47$ \\\\\\hline\n",
-      "\\end{tabular}\n",
-      "\n",
-      "\\end{table}\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "C:\\Users\\moell\\AppData\\Local\\Temp\\ipykernel_25588\\1918125892.py:53: RuntimeWarning:\n",
-      "\n",
-      "All-NaN slice encountered\n",
-      "\n"
-     ]
-    }
-   ],
-   "execution_count": 15
-  },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "outputs": [],
-   "execution_count": null,
-   "source": "",
-   "id": "a431f99070cdee3c"
-  }
- ],
- "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
-}

Tests/test_biaxial.py

@@ -1,34 +0,0 @@
-import numpy as np
-from ellipsoide import EllipsoidBiaxial
-from GHA_biaxial.bessel import gha1 as gha1_bessel
-from GHA_biaxial.gauss import gha1 as gha1_gauss
-from GHA_biaxial.rk import gha1 as gha1_rk
-from GHA_biaxial.gauss import gha2 as gha2_gauss
-
-re = EllipsoidBiaxial.init_name("Bessel")
-
-# phi0 = 0.6
-# lamb0 = 1.2
-# alpha0 = 0.45
-# s = 123456
-#
-# values_bessel = gha1_bessel(re, phi0, lamb0, alpha0, s)
-# alpha1_bessel = values_bessel[-1]
-# p1_bessel = re.bi_ell2cart(values_bessel[0], values_bessel[1], 0)
-#
-# values_gauss1 = gha1_gauss(re, phi0, lamb0, alpha0, s)
-# alpha1_gauss1 = values_gauss1[-1]
-# p1_gauss = re.bi_ell2cart(values_gauss1[0], values_gauss1[1], 0)
-#
-# values_rk = gha1_rk(re, phi0, lamb0 , alpha0, s, 10000)
-# alpha1_rk = values_rk[-1]
-# p1_rk = re.bi_ell2cart(values_rk[0], values_rk[1], 0)
-#
-# alpha0_gauss, alpha1_gauss2, s_gauss = gha2_gauss(re, phi0, lamb0, values_gauss1[0], values_gauss1[1])
-
-phi0 = 0.6
-lamb0 = 1.2
-
-cart = re.bi_ell2cart(phi0, lamb0, 0)
-ell = re.bi_cart2ell(cart)
-pass

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

@@ -1,116 +1,20 @@
-import numpy as np
-from numpy import sin, cos, arctan, arctan2, sqrt, pi, arccos
-import winkelumrechnungen as wu
-import jacobian_Ligas
-import matplotlib.pyplot as plt
-from typing import Tuple
-from numpy.typing import NDArray
 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
 
-class EllipsoidBiaxial:
-    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
+import jacobian_Ligas
+from utils_angle import wrap_mhalfpi_halfpi, wrap_mpi_pi
 
-    @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: 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: np.arctan2(self.a * np.sin(beta), self.b * np.cos(beta))
-    beta2phi = lambda self, beta: np.arctan2(self.a ** 2 * np.sin(beta), self.b ** 2 * np.cos(beta))
-
-    psi2beta = lambda self, psi: np.arctan2(self.b * np.sin(psi), self.a * np.cos(psi))
-    psi2phi = lambda self, psi: np.arctan2(self.a * np.sin(psi), self.b * np.cos(psi))
-
-    phi2beta = lambda self, phi: np.arctan2(self.b**2 * np.sin(phi), self.a**2 * np.cos(phi))
-    phi2psi = lambda self, phi: np.arctan2(self.b * np.sin(phi), self.a * np.cos(phi))
-
-    phi2p = lambda self, phi: self.N(phi) * cos(phi)
-
-    def bi_cart2ell(self, point: NDArrayself, 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 = 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])
 
 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
@@ -124,14 +28,19 @@ class EllipsoidTriaxial:
         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):
+    def init_name(cls, name: str) -> EllipsoidTriaxial:
         """
-        Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, BursaSima1980round, Eitschberger1978, Bursa1972,
-        Bursa1970, BesselBiaxial, Fiction, KarneyTest2024
+        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
@@ -164,11 +73,6 @@ class EllipsoidTriaxial:
             ay = 6378105
             b = 6356754
             return cls(ax, ay, b)
-        elif name == "BesselBiaxial":
-            ax = 6377397.15509
-            ay = 6377397.15508
-            b = 6356078.96290
-            return cls(ax, ay, b)
         elif name == "Fiction":
             ax = 6000000
             ay = 4000000
@@ -179,6 +83,8 @@ class EllipsoidTriaxial:
             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:
         """
@@ -218,8 +124,10 @@ class EllipsoidTriaxial:
         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 < 0:
-            t2 = np.nan
+        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:
@@ -261,7 +169,6 @@ class EllipsoidTriaxial:
         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
-        # print(s1, s2, s3)
 
         beta = arctan(sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
         if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7:
@@ -284,6 +191,11 @@ class EllipsoidTriaxial:
 
         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
 
@@ -412,14 +324,14 @@ class EllipsoidTriaxial:
                 i += 1
 
             if i == maxI:
-                raise Exception("Umrechnung ist nicht konvergiert")
+                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("Fehler in der Umrechnung cart2ell")
+                raise Exception("Umrechnung cart2ell: Punktdifferenz")
 
-            return beta, lamb
+            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
@@ -511,7 +423,7 @@ class EllipsoidTriaxial:
             i += 1
 
         if i == maxI:
-            raise Exception("Umrechung ist nicht konvergiert")
+            raise Exception("Umrechnung cart2ell: nicht konvergiert")
 
         return phi, lamb
 
@@ -577,6 +489,8 @@ class EllipsoidTriaxial:
                 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)
@@ -598,15 +512,15 @@ class EllipsoidTriaxial:
                 phi, lamb, h = self.cart2geod(point, f"ligas{new_mode}", maxIter, maxLoa)
             else:
                 if xG < 0 and yG < 0:
-                    lamb = -pi + lamb
+                    lamb += -pi
 
                 elif xG < 0:
-                    lamb = pi + lamb
+                    lamb += pi
 
                 if abs(zG) < eps:
                     phi = 0
-
-        return phi, lamb, h
+        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:
         """
@@ -643,8 +557,8 @@ class EllipsoidTriaxial:
             v = 2 * arctan2(v_check1, v_check2 + v_factor)
         else:
             v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor)
-
-        return u, v
+        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]:
         """
@@ -749,63 +663,71 @@ class EllipsoidTriaxial:
 
 
 if __name__ == "__main__":
-    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)
+    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)
 
-            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()
+    # 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()

jacobian_Ligas.py

@@ -1,6 +1,8 @@
+from typing import Tuple
+
 import numpy as np
 from numpy.typing import NDArray
-from typing import Tuple
+
 
 def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
     """
@@ -34,7 +36,7 @@ def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArr
 
     return invJ, fxE
 
-def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, 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
@@ -68,7 +70,7 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple
 
     return invJ, fxE
 
-def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, 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

nicht abgeben/GHA_biaxial/bessel.py

@@ -1,10 +1,20 @@
-from numpy import *
-import scipy as sp
-from ellipsoide import EllipsoidBiaxial
 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))

nicht abgeben/GHA_biaxial/gauss.py

@@ -1,8 +1,8 @@
-from numpy import sin, cos, pi, sqrt, tan, arcsin, arccos, arctan
-import ausgaben as aus
-from ellipsoide import EllipsoidBiaxial
 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]:
     """

nicht abgeben/GHA_biaxial/rk.py

@@ -1,13 +1,26 @@
-import runge_kutta as rk
-from numpy import sin, cos, tan
-import winkelumrechnungen as wu
-from ellipsoide import EllipsoidBiaxial
+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, v):
+        def ODE(s: float, v: NDArray):
             phi, lam, A = v
             V = re.V(phi)
             dphi = cos(A) * V ** 3 / re.c

nicht abgeben/Tests/alpha_conversion_test.ipynb

@@ -1,56 +1,41 @@
 {
  "cells": [
   {
+   "metadata": {},
    "cell_type": "code",
-   "id": "initial_id",
-   "metadata": {
-    "collapsed": true,
-    "ExecuteTime": {
-     "end_time": "2026-01-20T15:30:31.978159Z",
-     "start_time": "2026-01-20T15:30:31.835157Z"
-    }
-   },
    "source": [
     "%load_ext autoreload\n",
     "%autoreload 2"
    ],
+   "id": "a78faf7f4883772f",
    "outputs": [],
-   "execution_count": 1
+   "execution_count": null
   },
   {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-01-20T15:30:33.910807Z",
-     "start_time": "2026-01-20T15:30:32.803089Z"
-    }
-   },
+   "metadata": {},
    "cell_type": "code",
    "source": [
     "%reload_ext autoreload\n",
     "%autoreload 2\n",
+    "import numpy as np\n",
+    "\n",
     "import winkelumrechnungen as wu\n",
-    "from ellipsoide import EllipsoidTriaxial\n",
-    "from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n",
-    "import numpy as np"
+    "from GHA_triaxial.utils import alpha_ell2para, alpha_para2ell\n",
+    "from ellipsoid_triaxial import EllipsoidTriaxial"
    ],
-   "id": "9ad815aea55574e3",
+   "id": "46aa84a937fea491",
    "outputs": [],
-   "execution_count": 2
+   "execution_count": null
   },
   {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-01-20T15:33:40.785362Z",
-     "start_time": "2026-01-20T15:33:34.296487Z"
-    }
-   },
+   "metadata": {},
    "cell_type": "code",
    "source": [
     "ell = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n",
     "diffs = []\n",
-    "for beta_deg in range(-180, 181, 45):\n",
-    "    for lamb_deg in range(-90, 91, 45):\n",
-    "        for alpha_deg in range(0, 360, 45):\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",
@@ -67,17 +52,12 @@
     "            diffs.append((beta_deg, lamb_deg, alpha_deg, diff_1, diff_2))\n",
     "diffs = np.array(diffs)"
    ],
-   "id": "98b9b220118deb3f",
+   "id": "82fc6cbbe7d5abcb",
    "outputs": [],
-   "execution_count": 6
+   "execution_count": null
   },
   {
-   "metadata": {
-    "ExecuteTime": {
-     "end_time": "2026-01-20T15:33:50.497990Z",
-     "start_time": "2026-01-20T15:33:50.261115Z"
-    }
-   },
+   "metadata": {},
    "cell_type": "code",
    "source": [
     "i_max_ell = np.argmax(diffs[:, 3])\n",
@@ -92,18 +72,9 @@
     "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": "3c74b65b0e85e3c2",
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Für elliptisches Alpha = 315.0° und beta = -90.0°, lamb = -90.0°: diff = 3.426945967752335e-05\"\n",
-      "Für parametrisches Alpha = 315.0° und beta = -90.0°, lamb = -90.0°: diff = 3.426945967752335e-05\"\n"
-     ]
-    }
-   ],
-   "execution_count": 7
+   "id": "97b5b8c9ca5377ab",
+   "outputs": [],
+   "execution_count": null
   }
  ],
  "metadata": {

nicht abgeben/Tests/conversions_test.ipynb

@@ -20,13 +20,14 @@
    "source": [
     "%reload_ext autoreload\n",
     "%autoreload 2\n",
-    "import pickle\n",
-    "import numpy as np\n",
-    "import winkelumrechnungen as wu\n",
     "from itertools import product\n",
+    "\n",
+    "import numpy as np\n",
     "import pandas as pd\n",
-    "from ellipsoide import EllipsoidTriaxial\n",
-    "import plotly.graph_objects as go"
+    "import plotly.graph_objects as go\n",
+    "\n",
+    "import winkelumrechnungen as wu\n",
+    "from ellipsoid_triaxial import EllipsoidTriaxial"
    ],
    "outputs": [],
    "execution_count": null

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

@@ -1,7 +1,9 @@
-import numpy as np
-from numpy import sqrt, arctan2, sin, cos, arcsin, arccos
-from numpy.typing import NDArray
 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
 
 
@@ -77,7 +79,7 @@ def gha2(R: float, phi0: float, lamb0: float, phi1: float, lamb1: float) -> Tupl
     alpha1 = arctan2(-cos(phi0) * sin(lamb1 - lamb0),
                      cos(phi1) * sin(phi0) - sin(phi1) * cos(phi0) * cos(lamb1 - lamb0))
     if alpha1 < 0:
-        alpha1 += 2 * np.pi
+        alpha1 += 2 * pi
 
     return alpha0, alpha1, s
 

nicht abgeben/plots.ipynb

@@ -9,10 +9,11 @@
    },
    "cell_type": "code",
    "source": [
-    "import plotly.graph_objects as go\n",
     "import numpy as np\n",
-    "from ellipsoide import EllipsoidTriaxial\n",
-    "import winkelumrechnungen as wu"
+    "import plotly.graph_objects as go\n",
+    "\n",
+    "import winkelumrechnungen as wu\n",
+    "from ellipsoid_triaxial import EllipsoidTriaxial"
    ],
    "id": "731173e4745cfe7c",
    "outputs": [],

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)

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

@@ -1,7 +1,10 @@
+from typing import Callable
+
 import numpy as np
+from numpy.typing import NDArray
 
 
-def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]:
+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
@@ -9,7 +12,7 @@ def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool
     :param v0: Startwerte
     :param weite: Integrationsweite
     :param schritte: Schrittzahl
-    :param fein:
+    :param fein: Fein-Rechnung?
     :return: Variable und Funktionswerte an jedem Stützpunkt
     """
     h = weite/schritte
@@ -35,9 +38,118 @@ def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool
 
     return t_list, werte
 
-def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray:
+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)
+    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

test.py

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

utils_angle.py

@@ -1,13 +1,54 @@
 import numpy as np
 
+import winkelumrechnungen as wu
 
-def arccot(x):
+
+def arccot(x: float) -> float:
+    """
+    Berechnung von arccot eines Winkels
+    :param x: Winkel
+    :return: arccot(Winkel)
+    """
     return np.arctan2(1.0, x)
 
 
-def cot(a):
-    return np.cos(a) / np.sin(a)
+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_to_pi(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))))