Merge remote-tracking branch 'origin/main'

# Conflicts:
#	GHA_triaxial/gha2_num.py

GHA_triaxial/gha2_num.py

@@ -1,102 +1,149 @@
 import numpy as np
 from ellipsoide import EllipsoidTriaxial
 import runge_kutta as rk
+import GHA_triaxial.numeric_examples_karney as ne_karney
+import GHA_triaxial.numeric_examples_panou as ne_panou
+import winkelumrechnungen as wu
 from typing import Tuple
 from numpy.typing import NDArray
+import ausgaben as aus
+from utils_angle import cot, wrap_to_pi
 
-from utils_angle import arccot, cot, wrap_to_pi
 
+def arccot(x):
+    x = np.asarray(x)
+    a = np.arctan2(1.0, x)
+    return np.where(x < 0.0, a - np.pi, a)
+
+def normalize_alpha_0_pi(alpha):
+    if alpha < 0.0:
+        alpha += np.pi
+    return alpha
 
 def sph_azimuth(beta1, lam1, beta2, lam2):
-    # sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß)
     dlam = wrap_to_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,
+def gha2_num(
-             n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False
+    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]:
-    """
 
-    :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 n: Anzahl Schritte
-    :param epsilon:
-    :param iter_max: Maximale Anzhal Iterationen
-    :param all_points:
-    :return:
-    """
 
-    # h_x, h_y, h_e entsprechen E_x, E_y, E_e
     def BETA_LAMBDA(beta, 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)
+            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)
+            ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2
+        )
 
-        # Erste Ableitungen von ΒETA und LAMBDA
+        BETA_ = (ell.ax ** 2 * ell.Ey ** 2 * np.sin(2 * beta)) / (
-        BETA_ = (ell.ax ** 2 * ell.Ey ** 2 * np.sin(2 * beta)) / (ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2
+            ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2
-        LAMBDA_ = - (ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 2
+        ) ** 2
+        LAMBDA_ = -(ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (
+            ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2
+        ) ** 2
 
-        # Zweite Ableitungen von ΒETA und LAMBDA
+        BETA__ = (
-        BETA__ = ((2 * ell.ax ** 2 * ell.Ey ** 4 * np.sin(2 * beta) ** 2) / (
+            (2 * ell.ax ** 2 * ell.Ey ** 4 * np.sin(2 * beta) ** 2)
-                    ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 3) + (
+            / (ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 3
-                             (2 * ell.ax ** 2 * ell.Ey ** 2 * np.cos(2 * beta)) / (
+            + (2 * ell.ax ** 2 * ell.Ey ** 2 * np.cos(2 * beta))
-                                 ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2) ** 2)
+            / (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) -
+        LAMBDA__ = (
-                    ((2 * ell.b ** 2 * ell.Ee ** 2 * np.sin(2 * lamb)) / (
+            (2 * ell.b ** 2 * ell.Ee ** 4 * np.sin(2 * lamb) ** 2)
-                                ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 2) ** 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
+        )
 
         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)
 
-        # Erste Ableitungen von E und G
+        E_beta = (
-        E_beta = BETA_ * (
+            BETA_ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
-                    ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) - BETA * ell.Ey ** 2 * np.sin(
+            - BETA * ell.Ey ** 2 * np.sin(2 * beta)
-            2 * beta)
+        )
         E_lamb = BETA * ell.Ee ** 2 * np.sin(2 * lamb)
 
         G_beta = -LAMBDA * ell.Ey ** 2 * np.sin(2 * beta)
-        G_lamb = LAMBDA_ * (
+        G_lamb = (
-                    ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + LAMBDA * ell.Ee ** 2 * np.sin(
+            LAMBDA_ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
-            2 * lamb)
+            + LAMBDA * ell.Ee ** 2 * np.sin(2 * lamb)
+        )
 
-        # Zweite Ableitungen von E und G
+        E_beta_beta = (
-        E_beta_beta = BETA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(
+            BETA__ * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
-            lamb) ** 2) - 2 * BETA_ * ell.Ey ** 2 * np.sin(2 * beta) - 2 * BETA * ell.Ey ** 2 * np.cos(2 * beta)
+            - 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)
 
         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(
+        G_lamb_lamb = (
-            lamb) ** 2) + 2 * LAMBDA_ * ell.Ee ** 2 * np.sin(2 * lamb) + 2 * LAMBDA * ell.Ee ** 2 * np.cos(2 * 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)
+        )
 
-        return (BETA, LAMBDA, E, G,
+        return (
-                BETA_, LAMBDA_, BETA__, LAMBDA__,
+            BETA,
-                E_beta, E_lamb, G_beta, G_lamb,
+            LAMBDA,
-                E_beta_beta, E_beta_lamb, E_lamb_lamb,
+            E,
-                G_beta_beta, G_beta_lamb, G_lamb_lamb)
+            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,
+        )
 
     def p_coef(beta, lamb):
-
+        (
-        (BETA, LAMBDA, E, G,
+            BETA,
-         BETA_, LAMBDA_, BETA__, LAMBDA__,
+            LAMBDA,
-         E_beta, E_lamb, G_beta, G_lamb,
+            E,
-         E_beta_beta, E_beta_lamb, E_lamb_lamb,
+            G,
-         G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
+            BETA_,
+            LAMBDA_,
+            BETA__,
+            LAMBDA__,
+            E_beta,
+            E_lamb,
+            G_beta,
+            G_lamb,
+            E_beta_beta,
+            E_beta_lamb,
+            E_lamb_lamb,
+            G_beta_beta,
+            G_beta_lamb,
+            G_lamb_lamb,
+        ) = BETA_LAMBDA(beta, lamb)
 
         p_3 = -0.5 * (E_lamb / G)
         p_2 = (G_beta / G) - 0.5 * (E_beta / E)
@@ -104,38 +151,37 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
         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_22 = ((G * G_beta_beta - G_beta * G_beta) / (G**2)) - 0.5 * (
-        p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2))
+            (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,
+        return (BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00)
-                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
 
     def q_coef(beta, lamb):
-
+        (
-        (BETA, LAMBDA, E, G,
+            BETA,
-         BETA_, LAMBDA_, BETA__, LAMBDA__,
+            LAMBDA,
-         E_beta, E_lamb, G_beta, G_lamb,
+            E,
-         E_beta_beta, E_beta_lamb, E_lamb_lamb,
+            G,
-         G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
+            BETA_,
+            LAMBDA_,
+            BETA__,
+            LAMBDA__,
+            E_beta,
+            E_lamb,
+            G_beta,
+            G_lamb,
+            E_beta_beta,
+            E_beta_lamb,
+            E_lamb_lamb,
+            G_beta_beta,
+            G_beta_lamb,
+            G_lamb_lamb,
+        ) = BETA_LAMBDA(beta, lamb)
 
         q_3 = -0.5 * (G_beta / E)
         q_2 = (E_lamb / E) - 0.5 * (G_lamb / G)
@@ -143,253 +189,321 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
         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_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E**2)) - 0.5 * (
-        q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2))
+            (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,
+        return (BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00)
-                q_3, q_2, q_1, q_0,
-                q_33, q_22, q_11, q_00)
 
-    def buildODEbeta():
+    def rk4_last(f, t0, y0, dt, N):
-        def ODE(beta, v):
+        h = dt / N
-            lamb, lamb_p, Y3, Y4 = v
+        t = t0
+        y = np.array(y0, dtype=float, copy=True)
+        for _ in range(N):
+            k1 = f(t, y)
+            k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
+            k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
+            k4 = f(t + h, y + h * k3)
+            y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+            t = t + h
+        return t, y
 
-            (BETA, LAMBDA, E, G,
+    def rk4_last_with_integral(f, t0, y0, dt, N, integrand_at):
-             q_3, q_2, q_1, q_0,
-             q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
 
-            dlamb = lamb_p
+        h = dt / N
-            dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
+        habs = abs(h)
-            dY3 = Y4
+        t = t0
-            dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
+        y = np.array(y0, dtype=float, copy=True)
-                  (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
 
-            return np.array([dlamb, dlamb_p, dY3, dY4])
+        if N % 2 == 0:
+            # Simpson streaming
+            f0 = integrand_at(t, y)
+            odd_sum = 0.0
+            even_sum = 0.0
 
-        return ODE
+            for i in range(1, N + 1):
+                k1 = f(t, y)
+                k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
+                k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
+                k4 = f(t + h, y + h * k3)
+                y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+                t = t + h
+
+                fi = integrand_at(t, y)
+                if i == N:
+                    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, y, s
+
+        # Trapez streaming
+        f_prev = integrand_at(t, y)
+        acc = 0.0
+        for _ in range(N):
+            k1 = f(t, y)
+            k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
+            k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
+            k4 = f(t + h, y + h * k3)
+            y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+            t = t + h
+            f_cur = integrand_at(t, y)
+            acc += 0.5 * (f_prev + f_cur)
+            f_prev = f_cur
+        s = habs * acc
+        return t, y, s
+
+    def integrand_lambda(lamb, y):
+        beta = y[0]
+        beta_p = y[1]
+        (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
+        return np.sqrt(E * beta_p**2 + G)
+
+    def integrand_beta(beta, y):
+        lamb = y[0]
+        lamb_p = y[1]
+        (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
+        return np.sqrt(E + G * lamb_p**2)
 
     if lamb_1 != lamb_2:
         N = n
         dlamb = lamb_2 - lamb_1
-        alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
 
-        if abs(dlamb) < 1e-15:
+        def buildODElamb():
-            beta_0 = 0.0
+            def ODE(lamb, v):
-        else:
+                beta, beta_p, X3, X4 = v
-            (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
+                (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
-            beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph)
+
+                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)
+
+            return ODE
 
         ode_lamb = buildODElamb()
 
+        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) if abs(dlamb) >= 1e-15 else 0.0
+
+
+        N_newton = min(N, 4000)
+
         def solve_newton(beta_p0_init: float):
             beta_p0 = float(beta_p0_init)
-
             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)
+                _, y_end = rk4_last(ode_lamb, lamb_1, startwerte, dlamb, N_newton)
+                beta_end, _, X3_end, _ = y_end
 
-                beta_end, beta_p_end, X3_end, X4_end = states[-1]
                 delta = beta_end - beta_2
-
                 if abs(delta) < epsilon:
-                    return True, beta_p0, lamb_list, states
+                    return True, beta_p0
 
-                d_beta_end_d_beta0 = X3_end
+                if abs(X3_end) < 1e-20:
-                if abs(d_beta_end_d_beta0) < 1e-20:
+                    return False, None
-                    return False, None, None, None
 
-                step = delta / d_beta_end_d_beta0
+                step = delta / X3_end
                 max_step = 0.5
                 if abs(step) > max_step:
                     step = np.sign(step) * max_step
-
                 beta_p0 = beta_p0 - step
 
-            return False, None, None, None
+            return False, None
 
-        alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
+        ok, beta_p0_sol = solve_newton(beta_p0_sph)
-        (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
-        beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph)
 
-        guesses = [
+        if not ok:
-            beta_p0_sph,
-            0.5 * beta_p0_sph,
-            2.0 * beta_p0_sph,
-            -beta_p0_sph,
-            -0.5 * beta_p0_sph,
-        ]
 
+            candidates = [-beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
+            N_quick = min(N, 2000)
             best = None
 
-        for g in guesses:
+            for g in candidates:
-            ok, beta_p0_sol, lamb_list_cand, states_cand = solve_newton(g)
+                ok_g, beta_p0_sol_g = solve_newton(g)
-            if not ok:
+                if not ok_g:
                     continue
 
-            beta_arr_c = np.array([st[0] for st in states_cand], dtype=float)
+                startwerte_g = np.array([beta_1, beta_p0_sol_g, 0.0, 1.0], dtype=float)
-            beta_p_arr_c = np.array([st[1] for st in states_cand], dtype=float)
+                _, _, s_quick = rk4_last_with_integral(
-            lamb_arr_c = np.array(lamb_list_cand, dtype=float)
+                    ode_lamb, lamb_1, startwerte_g, dlamb, N_quick, integrand_lambda
+                )
 
-            integrand = np.zeros(N + 1)
+                if (best is None) or (s_quick < best[0]):
-            for i in range(N + 1):
+                    best = (s_quick, beta_p0_sol_g)
-                (_, _, 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)
 
             if best is None:
-            raise RuntimeError("Keine Multi-Start-Variante konvergiert.")
+                raise RuntimeError("Keine Startwert-Variante konvergiert.")
 
-        s_best, beta_0, lamb_list, werte = best
+            beta_p0_sol = best[1]
 
-        beta_arr = np.zeros(N + 1)
+        beta_0 = beta_p0_sol
-        # lamb_arr = np.zeros(N + 1)
+        startwerte_final = np.array([beta_1, beta_0, 0.0, 1.0], dtype=float)
-        lamb_arr = np.array(lamb_list)
-        beta_p_arr = np.zeros(N + 1)
 
-        for i, state in enumerate(werte):
+        if all_points:
-            # lamb_arr[i] = state[0]
+            lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte_final, dlamb, N, False)
-            # beta_arr[i] = state[1]
-            # beta_p_arr[i] = state[2]
-            beta_arr[i] = state[0]
-            beta_p_arr[i] = state[1]
 
-        (_, _, E1, G1,
+            lamb_arr = np.array(lamb_list, dtype=float)
-         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+            beta_arr = np.array([st[0] for st in states], dtype=float)
-        (_, _, E2, G2,
+            beta_p_arr = np.array([st[1] for st in states], dtype=float)
-         *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+
+            (_, _, E1, G1, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+            (_, _, E2, G2, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
 
             alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0])
             alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1])
 
-        integrand = np.zeros(N + 1)
+            alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+            alpha_2 = normalize_alpha_0_pi(float(alpha_2))
+
+            # Distanz s aus Arrays (Simpson/Trapz)
+            integrand = np.zeros(N + 1, dtype=float)
             for i in range(N + 1):
-            (_, _, Ei, Gi,
+                (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
-             *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
                 integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi)
 
             h = abs(dlamb) / N
             if N % 2 == 0:
-            S = integrand[0] + integrand[-1] \
+                S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(integrand[2:-1:2])
-                + 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)
 
-        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 all_points:
             return alpha_1, alpha_2, s, beta_arr, lamb_arr
-        else:
+
+        _, y_end, s = rk4_last_with_integral(ode_lamb, lamb_1, startwerte_final, dlamb, N, integrand_lambda)
+        beta_end, beta_p_end, _, _ = y_end
+
+        (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
+        (_, _, E2, G2, *_) = BETA_LAMBDA(beta_2, lamb_2)
+
+        alpha_1 = arccot(np.sqrt(E1 / G1) * beta_0)
+        alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_end)
+
+        alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+        alpha_2 = normalize_alpha_0_pi(float(alpha_2))
+
         return alpha_1, alpha_2, s
 
-    else:  # lamb_1 == lamb_2
     N = n
     dbeta = beta_2 - beta_1
 
     if abs(dbeta) < 1e-15:
         if all_points:
-                return 0, 0, 0, np.array([]), np.array([])
+            return 0.0, 0.0, 0.0, np.array([]), np.array([])
-            else:
+        return 0.0, 0.0, 0.0
-                return 0, 0, 0
 
-        lamb_0 = 0
+    # ODE-System (lambda, lambda', Y3, Y4) in Abhängigkeit von beta
+    def buildODEbeta():
+        def ODE(beta, v):
+            lamb, lamb_p, Y3, Y4 = v
+            (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
+
+            dlamb = lamb_p
+            dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0
+            dY3 = Y4
+            dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + (
+                3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1
+            ) * Y4
+
+            return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
+
+        return ODE
 
     ode_beta = buildODEbeta()
 
-        for i in range(iter_max):
+    # Newton auf lambda'_0
-            startwerte = [lamb_1, lamb_0, 0.0, 1.0]
+    lamb_0 = 0.0
+    for _ in range(iter_max):
+        startwerte = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float)
+        beta_list, states = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
 
-            beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
+        lamb_end, lamb_p_end, Y3_end, _ = states[-1]
-
-            beta_end = beta_list[-1]
-            lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
-
-            d_lamb_end_d_lambda0 = Y3_end
         delta = lamb_end - lamb_2
 
         if abs(delta) < epsilon:
             break
 
-            if abs(d_lamb_end_d_lambda0) < 1e-20:
+        if abs(Y3_end) < 1e-20:
             raise RuntimeError("Abbruch (Ableitung ~ 0).")
 
+        step = delta / Y3_end
         max_step = 1.0
-            step = delta / d_lamb_end_d_lambda0
         if abs(step) > max_step:
             step = np.sign(step) * max_step
-
         lamb_0 = lamb_0 - step
 
-        beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False)
+    startwerte_final = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float)
 
-        # beta_arr = np.zeros(N + 1)
+    if all_points:
-        beta_arr = np.array(beta_list)
+        beta_list, states = rk.rk4(ode_beta, beta_1, startwerte_final, dbeta, N, False)
-        lamb_arr = np.zeros(N + 1)
-        lambda_p_arr = np.zeros(N + 1)
 
-        for i, state in enumerate(werte):
+        beta_arr = np.array(beta_list, dtype=float)
-            # beta_arr[i] = state[0]
+        lamb_arr = np.array([st[0] for st in states], dtype=float)
-            # lamb_arr[i] = state[1]
+        lamb_p_arr = np.array([st[1] for st in states], dtype=float)
-            # lambda_p_arr[i] = state[2]
-            lamb_arr[i] = state[0]
-            lambda_p_arr[i] = state[1]
 
         # Azimute
-        (BETA1, LAMBDA1, E1, G1,
+        (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
-         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+        (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
-        (BETA2, LAMBDA2, E2, G2,
-         *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
 
-        alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lambda_p_arr[0])
+        alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_p_arr[0])
-        alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lambda_p_arr[-1])
+        alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_arr[-1])
 
-        integrand = np.zeros(N + 1)
+        # optionaler Quadrantenfix (robust)
+        alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+        alpha_2 = normalize_alpha_0_pi(float(alpha_2))
+
+        # Distanz
+        integrand = np.zeros(N + 1, dtype=float)
         for i in range(N + 1):
-            (_, _, Ei, Gi,
+            (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
-             *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
+            integrand[i] = np.sqrt(Ei + Gi * lamb_p_arr[i] ** 2)
-            integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2)
 
         h = abs(dbeta) / N
         if N % 2 == 0:
-            S = integrand[0] + integrand[-1] \
+            S = integrand[0] + integrand[-1] + 4.0 * np.sum(integrand[1:-1:2]) + 2.0 * np.sum(integrand[2:-1:2])
-                + 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:
+
+    # all_points == False: streaming Integral
+    _, y_end, s = rk4_last_with_integral(ode_beta, beta_1, startwerte_final, dbeta, N, integrand_beta)
+    lamb_end, lamb_p_end, _, _ = y_end
+
+    (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_1, lamb_1)
+    (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_2, lamb_2)
+
+    alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_0)
+    alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_end)
+
+    alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+    alpha_2 = normalize_alpha_0_pi(float(alpha_2))
+
     return alpha_1, alpha_2, s
 
 
 if __name__ == "__main__":
-    # ell = EllipsoidTriaxial.init_name("Fiction")
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    # # beta1 = np.deg2rad(75)
+    beta1 = np.deg2rad(75)
-    # # lamb1 = np.deg2rad(-90)
+    lamb1 = np.deg2rad(-90)
-    # # beta2 = np.deg2rad(75)
+    beta2 = np.deg2rad(75)
-    # # lamb2 = np.deg2rad(66)
+    lamb2 = np.deg2rad(66)
-    # # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
+    a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
-    # # print(aus.gms("a1", a1, 4))
+    print(aus.gms("a1", a1, 4))
     # # print(aus.gms("a2", a2, 4))
     # # print(s)
     # cart1 = ell.para2cart(0, 0)

dashboard.py

@@ -1,4 +1,4 @@
-from dash import Dash, html, dcc, Input, Output, State, no_update
+from dash import Dash, dash, html, dcc, Input, Output, State, no_update, ctx
 import plotly.graph_objects as go
 import numpy as np
 import dash_bootstrap_components as dbc
@@ -12,6 +12,7 @@ import ausgaben as aus
 
 from GHA_triaxial.gha1_ana import gha1_ana
 from GHA_triaxial.gha1_num import gha1_num
+from GHA_triaxial.gha1_ES import gha1_ES
 from GHA_triaxial.gha1_approx import gha1_approx
 
 from GHA_triaxial.gha2_num import gha2_num
@@ -23,6 +24,7 @@ app = Dash(__name__, suppress_callback_exceptions=True, external_stylesheets=[db
 app.title = "Geodätische Hauptaufgaben"
 
 
+# Erzeugen der Eingabefelder
 def inputfeld(left_text, input_id, right_text="", width=200, min=None, max=None):
     return html.Div(
         children=[
@@ -42,6 +44,7 @@ def inputfeld(left_text, input_id, right_text="", width=200, min=None, max=None)
         style={"display": "flex", "alignItems": "center", "marginBottom": "10px"},
     )
 
+# Erzeugen der Checklisten inkl. Eingabefelder
 def method_row(label, cb_id, input_id=None, value="", info=""):
     base_row_style = {
         "display": "flex",
@@ -102,8 +105,6 @@ def method_row(label, cb_id, input_id=None, value="", info=""):
     return html.Div(children, style=base_row_style)
 
 
-
-
 def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
     fig = go.Figure()
 
@@ -228,11 +229,12 @@ def figure_points(fig, points):
         ))
     return fig
 
-def figure_lines(fig, line, color):
+def figure_lines(fig, line, name, color):
     """
 
     :param fig: plotly.graph_objects.Figure
     :param line: Punktliste [[x1,y1,z1], [x2,y2,z2]]
+    :param name: Name
     :param color: Farbe
     :return: plotly.graph_objects.Figure
     """
@@ -242,10 +244,12 @@ def figure_lines(fig, line, color):
         x=points[:, 0], y=points[:, 1], z=points[:, 2],
         mode="lines",
         line=dict(width=4, color=color),
-        name="Strecke", showlegend=False
+        name=name, showlegend=False
     ))
     return fig
 
+# HTML der beiden Tabs
+# Tab 1
 pane_gha1 = html.Div(
     [
         inputfeld("β₀", "input-GHA1-beta0", "°", min=-90, max=90),
@@ -255,6 +259,7 @@ pane_gha1 = html.Div(
 
         method_row("Analytisch", "cb-ana-1", "input-ana-1", ""),
         method_row("Numerisch", "cb-num-1", "input-num-n-1", "2000", info="Anzahl Schritte"),
+        method_row("Stochastisch (ES)", "cb-stoch-1", "input-stoch-n-1", "1000", info="Info"),
         method_row("Approximiert", "cb-approx-1", "input-approx-ds-1", "1000", info="Länge Streckensegment [m]"),
 
 
@@ -287,6 +292,7 @@ pane_gha1 = html.Div(
     style={"display": "block"},
 )
 
+# Tab2
 pane_gha2 = html.Div(
     [
         inputfeld("β₀", "input-GHA2-beta0", "°", min=-90, max=90),
@@ -295,7 +301,7 @@ pane_gha2 = html.Div(
         inputfeld("λ₁", "input-GHA2-lamb1", "°", min=-180, max=180),
 
         method_row("Numerisch", "cb-num-2", "input-num-n-2", "2000", info="Anzahl Schritte"),
-        method_row("Stochastisch", "cb-stoch-2", "input-stoch-2", ""),
+        method_row("Stochastisch", "cb-stoch-2", "input-stoch-n-2", "1000", info="Info"),
         method_row("Approximiert", "cb-approx-2", "input-approx-ds-2", "1000", info="Länge Streckensegment [m]"),
 
         html.Div(
@@ -326,7 +332,7 @@ pane_gha2 = html.Div(
 )
 
 
-# HTML-Elemente
+# HTML-Gerüst
 app.layout = html.Div(
     style={"fontFamily": "Arial", "padding": "10px", "width": "95%", "margin": "0 auto"},
     children=[
@@ -342,6 +348,7 @@ app.layout = html.Div(
             },
             children=[
 
+                # Linker Bereich
                 html.Div(
                     style={"flex": "1 1 420px", "maxWidth": "640px"},
                     children=[
@@ -387,6 +394,7 @@ app.layout = html.Div(
                     ],
                 ),
 
+                # Rechter Bereich
                 html.Div(
                     style={
                         "flex": "1 1 750px",
@@ -452,6 +460,21 @@ def switch_tabs(tab):
     show2 = {"display": "block"} if tab == "tab-GHA2" else {"display": "none"}
     return show1, show2
 
+# Funktionen zum Aktivieren der Eingabefelder innerhalb der Checklisten
+@app.callback(
+    Output("input-num-n-1", "disabled"),
+    Input("cb-num-1", "value"),
+)
+def toggle_ds(v):
+    return "on" not in (v or [])
+
+@app.callback(
+    Output("input-stoch-n-1", "disabled"),
+    Input("cb-stoch-1", "value"),
+)
+def toggle_ds(v):
+    return "on" not in (v or [])
+
 @app.callback(
     Output("input-approx-ds-1", "disabled"),
     Input("cb-approx-1", "value"),
@@ -459,9 +482,17 @@ def switch_tabs(tab):
 def toggle_ds(v):
     return "on" not in (v or [])
 
+
 @app.callback(
-    Output("input-num-n-1", "disabled"),
+    Output("input-num-n-2", "disabled"),
-    Input("cb-num-1", "value"),
+    Input("cb-num-2", "value"),
+)
+def toggle_ds(v):
+    return "on" not in (v or [])
+
+@app.callback(
+    Output("input-stoch-n-2", "disabled"),
+    Input("cb-stoch-2", "value"),
 )
 def toggle_ds(v):
     return "on" not in (v or [])
@@ -473,21 +504,9 @@ def toggle_ds(v):
 def toggle_ds(v):
     return "on" not in (v or [])
 
-@app.callback(
-    Output("input-num-n-2", "disabled"),
-    Input("cb-num-2", "value"),
-)
-def toggle_ds(v):
-    return "on" not in (v or [])
-
-@app.callback(
-    Output("input-stoch-2", "disabled"),
-    Input("cb-stoch-2", "value"),
-)
-def toggle_ds(v):
-    return "on" not in (v or [])
 
 
+# Abfrage ob Berechnungsverfahren gewählt
 @app.callback(
     Output("tabs-GHA1-out", "children"),
     Input("button-calc-gha1", "n_clicks"),
@@ -528,8 +547,7 @@ def compute_gha1_ana(n1, cb_ana, beta0, lamb0, s, a0, ax, ay, b):
     #if not method1:
         #return html.Span("Bitte Berechnungsverfahren wählen.", style={"color": "red"}), None
     if "on" not in (cb_ana or []):
-        return out, no_update
+        return "", None
-
 
     ell = EllipsoidTriaxial(ax, ay, b)
     beta_rad = wu.deg2rad(float(beta0))
@@ -576,8 +594,8 @@ def compute_gha1_ana(n1, cb_ana, beta0, lamb0, s, a0, ax, ay, b):
 def compute_gha1_num(n1, cb_num, n_in, beta0, lamb0, s, a0, ax, ay, b):
     if not n1:
         return no_update, no_update
-    if not n1 or "on" not in (cb_num or []):
+    if  "on" not in (cb_num or []):
-        return no_update, no_update
+        return "", None
 
     n_in = int(n_in) if n_in else 2000
 
@@ -614,6 +632,8 @@ def compute_gha1_num(n1, cb_num, n_in, beta0, lamb0, s, a0, ax, ay, b):
     Output("output-gha1-stoch", "children"),
     Output("store-gha1-stoch", "data"),
     Input("button-calc-gha1", "n_clicks"),
+    State("cb-stoch-1", "value"),
+    State("input-stoch-n-1", "value"),
     State("input-GHA1-beta0", "value"),
     State("input-GHA1-lamb0", "value"),
     State("input-GHA1-s", "value"),
@@ -624,11 +644,13 @@ def compute_gha1_num(n1, cb_num, n_in, beta0, lamb0, s, a0, ax, ay, b):
     State("method-checklist-1", "value"),
     prevent_initial_call=True,
 )
-def compute_gha1_stoch(n1, beta0, lamb0, s, a0, ax, ay, b, method1):
+def compute_gha1_stoch(n1, cb_stoch, n_in, beta0, lamb0, s, a0, ax, ay, b, method1):
     if not n1:
         return no_update, no_update
-    if "stochastisch" not in (method1 or []):
+    if "on" not in (cb_stoch or []):
-        return no_update, no_update
+        return "", None
+
+    n_in = int(n_in) if n_in else 100
 
     ell = EllipsoidTriaxial(ax, ay, b)
     beta_rad = wu.deg2rad(float(beta0))
@@ -636,7 +658,7 @@ def compute_gha1_stoch(n1, beta0, lamb0, s, a0, ax, ay, b, method1):
     alpha_rad = wu.deg2rad(float(a0))
     s_val = float(s)
 
-    betas, lambs, alphas, S_real = gha1_es(
+    betas, lambs, alphas, S_real = gha1_ES(
         beta_rad, lamb_rad, alpha_rad,
         s_val,
         10000,
@@ -682,7 +704,7 @@ def compute_gha1_approx(n1, cb_approx, ds_in, beta0, lamb0, s, a0, ax, ay, b):
     if not n1:
         return no_update, no_update
     if not n1 or "on" not in (cb_approx or []):
-        return no_update, no_update
+        return "", None
 
     ds_in = int(ds_in) if ds_in else 1000
 
@@ -740,7 +762,7 @@ def compute_gha2_num(n2, cb_num, n_in, beta0, lamb0, beta1, lamb1, ax, ay, b):
     if None in (beta0, lamb0, beta1, lamb1):
         return html.Span("Bitte β₀, λ₀, β₁ und λ₁ eingeben.", style={"color": "red"}), None
     if "on" not in (cb_num or []):
-        return no_update, no_update
+        return "", None
 
     n_in = int(n_in) if n_in else 2000
 
@@ -770,6 +792,7 @@ def compute_gha2_num(n2, cb_num, n_in, beta0, lamb0, beta1, lamb1, ax, ay, b):
     store = {
         "points": [("P0", P0, "black"), ("P1", P1, "black")],
         "polyline": polyline,
+        "name": "numerisch",
         "color": "#ff8c00",
     }
     return out, store
@@ -779,7 +802,7 @@ def compute_gha2_num(n2, cb_num, n_in, beta0, lamb0, beta1, lamb1, ax, ay, b):
     Output("store-gha2-stoch", "data"),
     Input("button-calc-gha2", "n_clicks"),
     State("cb-stoch-2", "value"),
-    State("input-stoch-2", "value"),
+    State("input-stoch-n-2", "value"),
     State("input-GHA2-beta0", "value"),
     State("input-GHA2-lamb0", "value"),
     State("input-GHA2-beta1", "value"),
@@ -789,11 +812,13 @@ def compute_gha2_num(n2, cb_num, n_in, beta0, lamb0, beta1, lamb1, ax, ay, b):
     State("input-b", "value"),
     prevent_initial_call=True,
 )
-def compute_gha2_stoch(n2, cb_stoch, n_stoch, beta0, lamb0, beta1, lamb1, ax, ay, b):
+def compute_gha2_stoch(n2, cb_stoch, n_in, beta0, lamb0, beta1, lamb1, ax, ay, b):
     if not n2:
         return no_update, no_update
     if "on" not in (cb_stoch or []):
-        return no_update, no_update
+        return "", None
+
+    n_in = int(n_in) if n_in else 1000
 
     ell = EllipsoidTriaxial(ax, ay, b)
 
@@ -812,8 +837,14 @@ def compute_gha2_stoch(n2, cb_stoch, n_stoch, beta0, lamb0, beta1, lamb1, ax, ay
         html.Span(f"{aus.gms('α₀', a0_stoch, 4)}, α₁ = {a1_stoch}, s = {s_stoch:.4f} m"),
     ])
 
+    store = {
+        "points": [("P0", P0, "black"), ("P1", P1, "black")],
+        "polyline": points,
+        "name": "stochastisch (ES)",
+        "color": "#d62728",
+    }
 
-    return out, {"points": None, "polyline": None, "color": "#d62728"}
+    return out, store
 
 @app.callback(
     Output("output-gha2-approx", "children"),
@@ -834,7 +865,7 @@ def compute_gha2_approx(n2, cb_approx, ds_in, beta0, lamb0, beta1, lamb1, ax, ay
     if not n2:
         return no_update, no_update
     if "on" not in (cb_approx or []):
-        return no_update, no_update
+        return "", None
 
     ds_in = int(ds_in) if ds_in else 1000
 
@@ -858,6 +889,7 @@ def compute_gha2_approx(n2, cb_approx, ds_in, beta0, lamb0, beta1, lamb1, ax, ay
     store = {
         "points": [("P0", P0, "black"), ("P1", P1, "black")],
         "polyline": points,
+        "name": "aproximiert",
         "color": "#00c2fc",
     }
 
@@ -870,6 +902,8 @@ def compute_gha2_approx(n2, cb_approx, ds_in, beta0, lamb0, beta1, lamb1, ax, ay
     Input("input-ay", "value"),
     Input("input-b", "value"),
     Input("dropdown-coors-type", "value"),
+    Input("button-calc-gha1", "n_clicks"),
+    Input("button-calc-gha2", "n_clicks"),
     Input("store-gha1-ana", "data"),
     Input("store-gha1-num", "data"),
     Input("store-gha1-stoch", "data"),
@@ -878,34 +912,48 @@ def compute_gha2_approx(n2, cb_approx, ds_in, beta0, lamb0, beta1, lamb1, ax, ay
     Input("store-gha2-stoch", "data"),
     Input("store-gha2-approx", "data"),
 )
-def render_all(ax, ay, b, coords_type, store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx, store_gha2_num, store_gha2_stoch, store_gha2_approx):
+def render_all(ax, ay, b, coords_type,
+               n_clicks_gha1, n_clicks_gha2,
+               store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx,
+               store_gha2_num, store_gha2_stoch, store_gha2_approx):
+
     if None in (ax, ay, b):
         return go.Figure()
 
+    n1 = n_clicks_gha1 or 0
+    n2 = n_clicks_gha2 or 0
+
     ell = EllipsoidTriaxial(ax, ay, b)
     fig = ellipsoid_figure(ell, title="")
     fig = figure_constant_lines(fig, ell, coords_type)
 
-    def add_from_store(fig, store):
+    def add_from_store(fig, store, expected_calc_id):
         if not store:
             return fig
+        if store.get("calc_id") != expected_calc_id:
+            return fig
         pts = store.get("points") or []
         if pts:
             fig = figure_points(fig, pts)
         line = store.get("polyline")
+        name = store.get("name", "")
         if line:
-            fig = figure_lines(fig, line, store.get("color", "#ff8c00"))
+            fig = figure_lines(fig, line, name, store.get("color", "#ff8c00"))
         return fig
 
-    for st in (store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx, store_gha2_num, store_gha2_stoch, store_gha2_approx):
+    for st in (store_gha1_ana, store_gha1_num, store_gha1_stoch, store_gha1_approx):
-        fig = add_from_store(fig, st)
+        fig = add_from_store(fig, st, n1)
+
+    for st in (store_gha2_num, store_gha2_stoch, store_gha2_approx):
+        fig = add_from_store(fig, st, n2)
 
     return fig
 
+
 if __name__ == "__main__":
+    # Automatisiertes Öffnen der Seite im Browser
     HOST = "127.0.0.1"
     PORT = 8050
-
     #Timer(1.0, webbrowser.open_new_tab(f"http://{HOST}:{PORT}/")).start
 
     app.run(host=HOST, port=PORT, debug=False)