Optimierung

GHA_triaxial/gha2_num.py

@@ -56,65 +56,67 @@ def gha2_num(
     :return: Azimut Startpunkt, Azumit Zielpunkt, Strecke
     """
 
+    ax2 = float(ell.ax) * float(ell.ax)
+    ay2 = float(ell.ay) * float(ell.ay)
+    b2 = float(ell.b) * float(ell.b)
+    Ex2 = float(ell.Ex) * float(ell.Ex)
+    Ey2 = float(ell.Ey) * float(ell.Ey)
+    Ee2 = float(ell.Ee) * float(ell.Ee)
+    Ey4 = Ey2 * Ey2
+    Ee4 = Ee2 * Ee2
+    two_pi = 2.0 * np.pi
+
     # Berechnung Koeffizienten, Gaußschen Fundamentalgrößen 1. Ordnung sowie deren Ableitungen
     def BETA_LAMBDA(beta, lamb):
-        BETA = (ell.ay**2 * np.sin(beta) ** 2 + ell.b**2 * np.cos(beta) ** 2) / (
+        sb = np.sin(beta)
-            ell.Ex**2 - ell.Ey**2 * np.sin(beta) ** 2
+        cb = np.cos(beta)
+        sl = np.sin(lamb)
+        cl = np.cos(lamb)
+
+        sb2 = sb * sb
+        cb2 = cb * cb
+        sl2 = sl * sl
+        cl2 = cl * cl
+
+        s2b = 2.0 * sb * cb
+        c2b = cb2 - sb2
+        s2l = 2.0 * sl * cl
+        c2l = cl2 - sl2
+
+        denB = Ex2 - Ey2 * sb2
+        denL = Ex2 - Ee2 * cl2
+
+        BETA = (ay2 * sb2 + b2 * cb2) / denB
+        LAMBDA = (ax2 * sl2 + ay2 * cl2) / denL
+
+        BETA_ = (ax2 * Ey2 * s2b) / (denB * denB)
+        LAMBDA_ = -(b2 * Ee2 * s2l) / (denL * denL)
+
+        BETA__ = (2.0 * ax2 * Ey4 * (s2b * s2b)) / (denB * denB * denB) + (2.0 * ax2 * Ey2 * c2b) / (
+            denB * denB
         )
-        LAMBDA = (ell.ax**2 * np.sin(lamb) ** 2 + ell.ay**2 * np.cos(lamb) ** 2) / (
+        LAMBDA__ = (2.0 * b2 * Ee4 * (s2l * s2l)) / (denL * denL * denL) - (2.0 * b2 * Ee2 * s2l) / (
-            ell.Ex**2 - ell.Ee**2 * np.cos(lamb) ** 2
+            denL * denL
         )
 
-        BETA_ = (ell.ax**2 * ell.Ey**2 * np.sin(2 * beta)) / (
+        Q = Ey2 * cb2 + Ee2 * sl2
-            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
 
-        BETA__ = (
+        E = BETA * Q
-            (2 * ell.ax**2 * ell.Ey**4 * np.sin(2 * beta) ** 2)
+        G = LAMBDA * Q
-            / (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
-        )
 
-        E = BETA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2)
+        E_beta = BETA_ * Q - BETA * Ey2 * s2b
-        G = LAMBDA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2)
+        E_lamb = BETA * Ee2 * s2l
 
-        E_beta = (
+        G_beta = -LAMBDA * Ey2 * s2b
-            BETA_ * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2)
+        G_lamb = LAMBDA_ * Q + LAMBDA * Ee2 * s2l
-            - BETA * ell.Ey**2 * np.sin(2 * beta)
-        )
-        E_lamb = BETA * ell.Ee**2 * np.sin(2 * lamb)
 
-        G_beta = -LAMBDA * ell.Ey**2 * np.sin(2 * beta)
+        E_beta_beta = BETA__ * Q - 2.0 * BETA_ * Ey2 * s2b - 2.0 * BETA * Ey2 * c2b
-        G_lamb = (
+        E_beta_lamb = BETA_ * Ee2 * s2l
-            LAMBDA_ * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2)
+        E_lamb_lamb = 2.0 * BETA * Ee2 * c2l
-            + LAMBDA * ell.Ee**2 * np.sin(2 * lamb)
-        )
 
-        E_beta_beta = (
+        G_beta_beta = -2.0 * LAMBDA * Ey2 * c2b
-            BETA__ * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2)
+        G_beta_lamb = -LAMBDA_ * Ey2 * s2b
-            - 2 * BETA_ * ell.Ey**2 * np.sin(2 * beta)
+        G_lamb_lamb = LAMBDA__ * Q + 2.0 * LAMBDA_ * Ee2 * s2l + 2.0 * LAMBDA * Ee2 * c2l
-            - 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(lamb) ** 2)
-            + 2 * LAMBDA_ * ell.Ee**2 * np.sin(2 * lamb)
-            + 2 * LAMBDA * ell.Ee**2 * np.cos(2 * lamb)
-        )
 
         return (
             BETA,
@@ -227,20 +229,13 @@ def gha2_num(
         (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
         return np.sqrt(E + G * lamb_p**2)
 
-    lamb_0 = wrap_to_pi(lamb_0)
+    def solve_lambda_branch(beta0, lamb0, beta1, lamb1_target, N_run, N_newt, it_max):
-    lamb_1 = wrap_to_pi(lamb_1)
+        dlamb = float(lamb1_target - lamb0)
+        if abs(dlamb) < 1e-15:
+            return None
 
-    # Fall 1 (lambda_0 != lambda_1)
-    if abs(lamb_1 - lamb_0) >= 1e-15:
-        N = int(n)
-        dlamb = wrap_to_pi(lamb_1 - lamb_0)
         sgn = 1.0 if dlamb >= 0.0 else -1.0
 
-        beta0 = float(beta_0)
-        lamb0 = float(lamb_0)
-        beta1 = float(beta_1)
-        lamb1 = float(lamb_1)
-
         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)
@@ -253,13 +248,183 @@ def gha2_num(
             ) * X4
             return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
 
-        alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1)
+        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)
 
-        N_newton = min(N, 4000)
-
         def solve_newton(beta_p0_init: float):
+            beta_p0 = float(beta_p0_init)
+            for _ in range(it_max):
+                v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
+                _, y_end = rk4_end(ode_lamb, lamb0, v0, dlamb, N_newt)
+
+                beta_end, _, X3_end, _ = y_end
+                delta = beta_end - beta1
+
+                if abs(delta) < epsilon:
+                    return True, beta_p0
+
+                if abs(X3_end) < 1e-20:
+                    return False, None
+
+                step = delta / X3_end
+                step = float(np.clip(step, -0.5, 0.5))
+                beta_p0 -= step
+
+            return False, None
+
+        seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
+        best = None
+
+        for seed in seeds:
+            ok, sol = solve_newton(seed)
+            if not ok:
+                continue
+            v0_sol = np.array([beta0, sol, 0.0, 1.0], dtype=float)
+            _, _, s_val = rk4_integral(ode_lamb, lamb0, v0_sol, dlamb, N_run, integrand_lambda)
+            if (best is None) or (s_val < best[0]):
+                best = (float(s_val), float(sol))
+
+        if best is None:
+            return None
+
+        return best[0], best[1], sgn, dlamb, ode_lamb
+
+    def solve_beta_branch(beta0, lamb0, beta1, lamb1, N_run, N_newt, it_max):
+        dbeta = float(beta1 - beta0)
+        if abs(dbeta) < 1e-15:
+            return None
+
+        sgn = 1.0 if dbeta >= 0.0 else -1.0
+
+        def ode_beta(beta, v):
+            lamb, lamb_p, Y3, Y4 = v
+            (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
+
+            dlamb = lamb_p
+            dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0
+            dY3 = Y4
+            dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + (
+                3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1
+            ) * Y4
+            return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
+
+        def solve_newton(lamb_p0_init: float):
+            lamb_p0 = float(lamb_p0_init)
+            for _ in range(it_max):
+                v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
+                _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N_newt)
+
+                lamb_end, _, Y3_end, _ = y_end
+                delta = lamb_end - lamb1
+
+                if abs(delta) < epsilon:
+                    return True, lamb_p0
+
+                if abs(Y3_end) < 1e-20:
+                    return False, None
+
+                step = delta / Y3_end
+                step = float(np.clip(step, -1.0, 1.0))
+                lamb_p0 -= step
+
+            return False, None
+
+        seeds = [0.0, 0.25, -0.25, 1.0, -1.0]
+        best = None
+
+        for seed in seeds:
+            ok, sol = solve_newton(seed)
+            if not ok:
+                continue
+            v0_sol = np.array([lamb0, sol, 0.0, 1.0], dtype=float)
+            _, _, s_val = rk4_integral(ode_beta, beta0, v0_sol, dbeta, N_run, integrand_beta)
+            if (best is None) or (s_val < best[0]):
+                best = (float(s_val), float(sol))
+
+        if best is None:
+            return None
+
+        return best[0], best[1], sgn, dbeta, ode_beta
+
+    lamb0 = float(wrap_to_pi(lamb_0))
+    lamb1 = float(wrap_to_pi(lamb_1))
+    beta0 = float(beta_0)
+    beta1 = float(beta_1)
+
+    N_full = int(n)
+    if N_full < 2:
+        N_full = 2
+
+    if all_points:
+        N_fast = min(2000, max(400, N_full // 10))
+    else:
+        N_fast = min(1500, max(300, N_full // 12))
+
+    k0 = int(np.round((lamb0 - lamb1) / two_pi))
+    lamb_targets = []
+    for dk in (-1, 0, 1):
+        lt = lamb1 + two_pi * float(k0 + dk)
+        dl = lt - lamb0
+        if abs(dl) <= np.pi + 1e-12:
+            lamb_targets.append(float(lt))
+    if not lamb_targets:
+        lamb_targets = [float(lamb1 + two_pi * float(k0))]
+
+    best_fast = None
+
+    for lt in lamb_targets:
+        if abs(lt - lamb0) >= 1e-15:
+            res = solve_lambda_branch(beta0, lamb0, beta1, lt, N_fast, min(N_fast, 800), min(iter_max, 12))
+            if res is None:
+                continue
+            s_fast, beta_p0_fast, sgn_fast, dlamb_fast, _ = res
+            cand = ("lambda", s_fast, lt, beta_p0_fast, sgn_fast, dlamb_fast)
+        else:
+            res = solve_beta_branch(beta0, lamb0, beta1, lamb1, N_fast, min(N_fast, 800), min(iter_max, 12))
+            if res is None:
+                continue
+            s_fast, lamb_p0_fast, sgn_fast, dbeta_fast, _ = res
+            cand = ("beta", s_fast, lt, lamb_p0_fast, sgn_fast, dbeta_fast)
+
+        if (best_fast is None) or (cand[1] < best_fast[1]):
+            best_fast = cand
+
+    if best_fast is None:
+        if abs(lamb1 - lamb0) >= 1e-15:
+            best_fast = ("lambda", 0.0, lamb1, None, 1.0, float(lamb1 - lamb0))
+        else:
+            best_fast = ("beta", 0.0, lamb1, None, 1.0, float(beta1 - beta0))
+
+    if best_fast[0] == "lambda":
+        lt = float(best_fast[2])
+        dlamb = float(lt - lamb0)
+        sgn = 1.0 if dlamb >= 0.0 else -1.0
+
+        def ode_lamb(lamb, v):
+            beta, beta_p, X3, X4 = v
+            (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
+
+            dbeta = beta_p
+            dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0
+            dX3 = X4
+            dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (
+                3 * p_3 * beta_p**2 + 2 * p_2 * beta_p + p_1
+            ) * X4
+            return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
+
+        alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lt)
+        (_, _, E0, G0, *_) = BETA_LAMBDA(beta0, lamb0)
+        beta_p0_sph = np.sqrt(G0 / E0) * cot(alpha0_sph)
+
+        beta_p0_init = best_fast[3]
+        if beta_p0_init is None:
+            beta_p0_init = beta_p0_sph
+        beta_p0_init = float(beta_p0_init)
+
+        N_newton = min(N_full, 4000)
+
+        def solve_newton_refine(beta_p0_init: float):
             beta_p0 = float(beta_p0_init)
             for _ in range(iter_max):
                 v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
@@ -275,25 +440,24 @@ def gha2_num(
                     return False, None
 
                 step = delta / X3_end
-                step = np.clip(step, -0.5, 0.5)
+                step = float(np.clip(step, -0.5, 0.5))
                 beta_p0 -= step
 
             return False, None
 
-        ok, beta_p0_sol = solve_newton(beta_p0_sph)
+        ok, beta_p0_sol = solve_newton_refine(beta_p0_init)
 
         if not ok:
-            candidates = [-beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
+            seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
-            N_quick = min(N, 2000)
             best = None
-            for g in candidates:
+            for seed in seeds:
-                ok_g, sol = solve_newton(g)
+                ok_s, sol_s = solve_newton_refine(seed)
-                if not ok_g:
+                if not ok_s:
                     continue
-                v0_g = np.array([beta0, sol, 0.0, 1.0], dtype=float)
+                v0_s = np.array([beta0, sol_s, 0.0, 1.0], dtype=float)
-                _, _, s_quick = rk4_integral(ode_lamb, lamb0, v0_g, dlamb, N_quick, integrand_lambda)
+                _, _, s_s = rk4_integral(ode_lamb, lamb0, v0_s, dlamb, min(N_full, 2000), integrand_lambda)
-                if (best is None) or (s_quick < best[0]):
+                if (best is None) or (s_s < best[0]):
-                    best = (s_quick, sol)
+                    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]
@@ -302,7 +466,7 @@ def gha2_num(
         v0_final = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
 
         if all_points:
-            lamb_list, states = rk4(ode_lamb, lamb0, v0_final, dlamb, N, False)
+            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)
@@ -313,14 +477,13 @@ def gha2_num(
             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)
 
-            # Distanz aus Arrays
+            integrand = np.zeros(N_full + 1, dtype=float)
-            integrand = np.zeros(N + 1, dtype=float)
+            for i in range(N_full + 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)
 
-            h = abs(dlamb) / N
+            h = abs(dlamb) / N_full
-            if N % 2 == 0:
+            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]
                 )
@@ -330,13 +493,13 @@ def gha2_num(
 
             return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr
 
-        _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N, integrand_lambda)
+        _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N_full, integrand_lambda)
         beta_end, beta_p_end, _, _ = y_end
 
         (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
         alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p0 * sgn, dlamb_du=1.0 * sgn)
 
-        (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), lamb1)
+        (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), float(lamb0 + dlamb))
         alpha_1 = azimut(E_end, G_end, dbeta_du=float(beta_p_end) * sgn, dlamb_du=1.0 * sgn)
 
         return float(alpha_0), float(alpha_1), float(s)
@@ -350,10 +513,6 @@ def gha2_num(
             return 0.0, 0.0, 0.0, np.array([]), np.array([])
         return 0.0, 0.0, 0.0
 
-    beta0 = float(beta_0)
-    lamb0 = float(lamb_0)
-    beta1 = float(beta_1)
-    lamb1 = float(lamb_1)
     sgn = 1.0 if dbeta >= 0.0 else -1.0
 
     def ode_beta(beta, v):
@@ -368,7 +527,8 @@ def gha2_num(
         ) * Y4
         return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
 
-    lamb_p0 = 0.0
+    lamb_p0 = float(best_fast[3]) if (best_fast[0] == "beta" and best_fast[3] is not None) else 0.0
+
     for _ in range(iter_max):
         v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
         _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N)
@@ -383,7 +543,7 @@ def gha2_num(
             raise RuntimeError("GHA2_num: Ableitung ~ 0 im beta-Fall")
 
         step = delta / Y3_end
-        step = np.clip(step, -1.0, 1.0)
+        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)