diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 734dc22..c898cfc 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -1,382 +1,510 @@ 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 -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) - 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) - - # 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 - - # 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)) - - 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 = 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 * 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) - - # 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) - - 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, 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) - -def p_coef(beta, lamb): - - (BETA, LAMBDA, E, G, - BETA_, LAMBDA_, BETA__, LAMBDA__, - E_beta, E_lamb, G_beta, G_lamb, - E_beta_beta, E_beta_lamb, E_lamb_lamb, - G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb) - - p_3 = - 0.5 * (E_lamb / G) - p_2 = (G_beta / G) - 0.5 * (E_beta / E) - p_1 = 0.5 * (G_lamb / G) - (E_lamb / E) - p_0 = 0.5 * (G_beta / E) - - p_33 = - 0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G ** 2)) - p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G ** 2)) - 0.5 * ((E * E_beta_beta - E_beta * E_beta) / (E ** 2)) - p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E ** 2)) - - return (BETA, LAMBDA, E, G, - p_3, p_2, p_1, p_0, - p_33, p_22, p_11, p_00) - -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_, LAMBDA_, BETA__, LAMBDA__, - E_beta, E_lamb, G_beta, G_lamb, - E_beta_beta, E_beta_lamb, E_lamb_lamb, - G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb) - - q_3 = - 0.5 * (G_beta / E) - q_2 = (E_lamb / E) - 0.5 * (G_lamb / G) - q_1 = 0.5 * (E_beta / E) - (G_beta / G) - q_0 = 0.5 * (E_lamb / G) - - q_33 = - 0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E ** 2)) - q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E ** 2)) - 0.5 * ((G * G_lamb_lamb - G_lamb * G_lamb) / (G ** 2)) - q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G ** 2)) - - return (BETA, LAMBDA, E, G, - q_3, q_2, q_1, q_0, - q_33, q_22, q_11, q_00) - -def buildODEbeta(): - def ODE(beta, v): - lamb, lamb_p, Y3, Y4 = v - - (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 - # 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_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 + ) + 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 + ) + + 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 + + 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 * (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_ * (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 * 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__ * (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) + ) + + 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, + ) + + def p_coef(beta, lamb): + ( + BETA, + LAMBDA, + E, + G, + BETA_, + LAMBDA_, + BETA__, + LAMBDA__, + E_beta, + E_lamb, + G_beta, + G_lamb, + E_beta_beta, + E_beta_lamb, + E_lamb_lamb, + G_beta_beta, + G_beta_lamb, + G_lamb_lamb, + ) = BETA_LAMBDA(beta, lamb) + + p_3 = -0.5 * (E_lamb / G) + p_2 = (G_beta / G) - 0.5 * (E_beta / E) + p_1 = 0.5 * (G_lamb / G) - (E_lamb / E) + p_0 = 0.5 * (G_beta / E) + + p_33 = -0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G**2)) + p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G**2)) - 0.5 * ( + (E * E_beta_beta - E_beta * E_beta) / (E**2) + ) + p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G**2)) - ( + (E * E_beta_lamb - E_beta * E_lamb) / (E**2) + ) + p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E**2)) + + return (BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) + + def q_coef(beta, lamb): + ( + BETA, + LAMBDA, + E, + G, + BETA_, + LAMBDA_, + BETA__, + LAMBDA__, + E_beta, + E_lamb, + G_beta, + G_lamb, + E_beta_beta, + E_beta_lamb, + E_lamb_lamb, + G_beta_beta, + G_beta_lamb, + G_lamb_lamb, + ) = BETA_LAMBDA(beta, lamb) + + q_3 = -0.5 * (G_beta / E) + q_2 = (E_lamb / E) - 0.5 * (G_lamb / G) + q_1 = 0.5 * (E_beta / E) - (G_beta / G) + q_0 = 0.5 * (E_lamb / G) + + q_33 = -0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E**2)) + q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E**2)) - 0.5 * ( + (G * G_lamb_lamb - G_lamb * G_lamb) / (G**2) + ) + q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E**2)) - ( + (G * G_beta_lamb - G_beta * G_lamb) / (G**2) + ) + q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G**2)) + + return (BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) + + def rk4_last(f, t0, y0, dt, N): + h = dt / N + 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 + + def rk4_last_with_integral(f, t0, y0, dt, N, integrand_at): + + h = dt / N + habs = abs(h) + t = t0 + y = np.array(y0, dtype=float, copy=True) + + if N % 2 == 0: + # Simpson streaming + f0 = integrand_at(t, y) + odd_sum = 0.0 + even_sum = 0.0 + + 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: - beta_0 = 0.0 - else: - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) - beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph) + def buildODElamb(): + def ODE(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) + + 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(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 + 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) - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) - beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) + ok, beta_p0_sol = solve_newton(beta_p0_sph) - guesses = [ - beta_p0_sph, - 0.5 * beta_p0_sph, - 2.0 * beta_p0_sph, - -beta_p0_sph, - -0.5 * beta_p0_sph, - ] + if not ok: - best = None + 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: - ok, beta_p0_sol, lamb_list_cand, states_cand = solve_newton(g) - if not ok: - continue + for g in candidates: + ok_g, beta_p0_sol_g = solve_newton(g) + if not ok_g: + 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) + startwerte_g = np.array([beta_1, beta_p0_sol_g, 0.0, 1.0], dtype=float) + _, _, s_quick = rk4_last_with_integral( + 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]): + best = (s_quick, beta_p0_sol_g) + + if best is None: + raise RuntimeError("Keine Startwert-Variante konvergiert.") + + beta_p0_sol = best[1] + + beta_0 = beta_p0_sol + startwerte_final = np.array([beta_1, beta_0, 0.0, 1.0], dtype=float) + + if all_points: + lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte_final, dlamb, N, 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) + + (_, _, 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]) + + 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, *_) = BETA_LAMBDA(beta_arr_c[i], lamb_arr_c[i]) - integrand[i] = np.sqrt(Ei * beta_p_arr_c[i] ** 2 + Gi) + (_, _, 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: - 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 + 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_cand = np.trapz(integrand, dx=h) + s = 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.") - - s_best, beta_0, lamb_list, werte = best - - 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) - - 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] - - (_, _, 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) - 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 - if N % 2 == 0: - S = integrand[0] + integrand[-1] \ - + 4.0 * np.sum(integrand[1:-1:2]) \ - + 2.0 * np.sum(integrand[2:-1:2]) - s = h / 3.0 * S - else: - s = np.trapz(integrand, dx=h) - - 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: - return alpha_1, alpha_2, s - else: # lamb_1 == lamb_2 - N = n - dbeta = beta_2 - beta_1 + # all_points == False (schnell) + _, y_end, s = rk4_last_with_integral(ode_lamb, lamb_1, startwerte_final, dlamb, N, integrand_lambda) + beta_end, beta_p_end, _, _ = y_end - if abs(dbeta) < 1e-15: - if all_points: - return 0, 0, 0, np.array([]), np.array([]) - else: - return 0, 0, 0 + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) + (_, _, E2, G2, *_) = BETA_LAMBDA(beta_2, lamb_2) - lamb_0 = 0 + alpha_1 = arccot(np.sqrt(E1 / G1) * beta_0) + alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_end) - ode_beta = buildODEbeta() + alpha_1 = normalize_alpha_0_pi(float(alpha_1)) + alpha_2 = normalize_alpha_0_pi(float(alpha_2)) - for i in range(iter_max): - startwerte = [lamb_1, lamb_0, 0.0, 1.0] + return alpha_1, alpha_2, s - beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False) + N = n + dbeta = beta_2 - beta_1 - beta_end = beta_list[-1] - lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1] + 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 - d_lamb_end_d_lambda0 = Y3_end - delta = lamb_end - lamb_2 + # 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) - if abs(delta) < epsilon: - break + 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 - if abs(d_lamb_end_d_lambda0) < 1e-20: - raise RuntimeError("Abbruch (Ableitung ~ 0).") + return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float) - max_step = 1.0 - step = delta / d_lamb_end_d_lambda0 - if abs(step) > max_step: - step = np.sign(step) * max_step + return ODE - lamb_0 = lamb_0 - step + ode_beta = buildODEbeta() - beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False) + # Newton auf lambda'_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_arr = np.zeros(N + 1) - beta_arr = np.array(beta_list) - lamb_arr = np.zeros(N + 1) - lambda_p_arr = np.zeros(N + 1) + lamb_end, lamb_p_end, Y3_end, _ = states[-1] + delta = lamb_end - lamb_2 - 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] + if abs(delta) < epsilon: + break + + if abs(Y3_end) < 1e-20: + raise RuntimeError("Abbruch (Ableitung ~ 0).") + + step = delta / Y3_end + max_step = 1.0 + if abs(step) > max_step: + step = np.sign(step) * max_step + lamb_0 = lamb_0 - step + + startwerte_final = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float) + + if all_points: + beta_list, states = rk.rk4(ode_beta, beta_1, startwerte_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) # 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]) + (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) - 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]) + 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) * 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, - *_) = 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 alpha_1, alpha_2, s, beta_arr, lamb_arr + + # 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") - # # 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)) + ell = EllipsoidTriaxial.init_name("BursaSima1980round") + 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, n=5000) + print(aus.gms("a1", a1, 4)) # # print(aus.gms("a2", a2, 4)) # # print(s) # cart1 = ell.para2cart(0, 0)