diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 0548ec2..eb54c5d 100644 --- a/GHA_triaxial/gha2_num.py +++ b/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) / ( - ell.Ex**2 - ell.Ey**2 * np.sin(beta) ** 2 + sb = np.sin(beta) + cb = np.cos(beta) + sl = np.sin(lamb) + cl = np.cos(lamb) + + sb2 = sb * sb + cb2 = cb * cb + sl2 = sl * sl + cl2 = cl * cl + + s2b = 2.0 * sb * cb + c2b = cb2 - sb2 + s2l = 2.0 * sl * cl + c2l = cl2 - sl2 + + denB = Ex2 - Ey2 * sb2 + denL = Ex2 - Ee2 * cl2 + + BETA = (ay2 * sb2 + b2 * cb2) / denB + LAMBDA = (ax2 * sl2 + ay2 * cl2) / denL + + BETA_ = (ax2 * Ey2 * s2b) / (denB * denB) + LAMBDA_ = -(b2 * Ee2 * s2l) / (denL * denL) + + BETA__ = (2.0 * ax2 * Ey4 * (s2b * s2b)) / (denB * denB * denB) + (2.0 * ax2 * Ey2 * c2b) / ( + denB * denB ) - LAMBDA = (ell.ax**2 * np.sin(lamb) ** 2 + ell.ay**2 * np.cos(lamb) ** 2) / ( - ell.Ex**2 - ell.Ee**2 * np.cos(lamb) ** 2 + LAMBDA__ = (2.0 * b2 * Ee4 * (s2l * s2l)) / (denL * denL * denL) - (2.0 * b2 * Ee2 * s2l) / ( + denL * denL ) - 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 + Q = Ey2 * cb2 + Ee2 * sl2 - 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 - ) + E = BETA * Q + G = LAMBDA * Q - E = BETA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2) - G = LAMBDA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2) + E_beta = BETA_ * Q - BETA * Ey2 * s2b + E_lamb = BETA * Ee2 * s2l - 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) + G_beta = -LAMBDA * Ey2 * s2b + G_lamb = LAMBDA_ * Q + LAMBDA * Ee2 * s2l - 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) - ) + 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 - 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) - - 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) - ) + 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, @@ -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) - lamb_1 = wrap_to_pi(lamb_1) + def solve_lambda_branch(beta0, lamb0, beta1, lamb1_target, N_run, N_newt, it_max): + dlamb = float(lamb1_target - lamb0) + if abs(dlamb) < 1e-15: + return None - # 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] - N_quick = min(N, 2000) + seeds = [beta_p0_sph, -beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph] best = None - for g in candidates: - ok_g, sol = solve_newton(g) - if not ok_g: + for seed in seeds: + ok_s, sol_s = solve_newton_refine(seed) + if not ok_s: continue - v0_g = np.array([beta0, sol, 0.0, 1.0], dtype=float) - _, _, s_quick = rk4_integral(ode_lamb, lamb0, v0_g, dlamb, N_quick, integrand_lambda) - if (best is None) or (s_quick < best[0]): - best = (s_quick, sol) + 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] @@ -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 + 1, dtype=float) - for i in range(N + 1): + 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 - if N % 2 == 0: + 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] ) @@ -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)