From 39641b52939657c2818b58ddef8a61ecf11df750 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 15:52:22 +0100 Subject: [PATCH 1/9] ell --- GHA_triaxial/gha2_num.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 734dc22..8a9f313 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -17,7 +17,7 @@ def sph_azimuth(beta1, lam1, beta2, lam2): a += 2 * np.pi return a -def BETA_LAMBDA(beta, lamb): +def BETA_LAMBDA(ell, 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) @@ -63,7 +63,7 @@ def p_coef(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) + G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(ell, beta, lamb) p_3 = - 0.5 * (E_lamb / G) p_2 = (G_beta / G) - 0.5 * (E_beta / E) @@ -102,7 +102,7 @@ def q_coef(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) + G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(ell, beta, lamb) q_3 = - 0.5 * (G_beta / E) q_2 = (E_lamb / E) - 0.5 * (G_lamb / G) @@ -164,7 +164,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float if abs(dlamb) < 1e-15: beta_0 = 0.0 else: - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph) ode_lamb = buildODElamb() @@ -196,7 +196,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float return False, None, None, None alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) - (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) guesses = [ @@ -220,7 +220,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): - (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr_c[i], lamb_arr_c[i]) + (_, _, Ei, Gi, *_) = BETA_LAMBDA(ell, beta_arr_c[i], lamb_arr_c[i]) integrand[i] = np.sqrt(Ei * beta_p_arr_c[i] ** 2 + Gi) h = abs(dlamb) / N @@ -253,9 +253,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float beta_p_arr[i] = state[1] (_, _, E1, G1, - *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) (_, _, E2, G2, - *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) + *_) = BETA_LAMBDA(ell, 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]) @@ -263,7 +263,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) + *_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i]) integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi) h = abs(dlamb) / N @@ -341,9 +341,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float # Azimute (BETA1, LAMBDA1, E1, G1, - *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) (BETA2, LAMBDA2, E2, G2, - *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) + *_) = BETA_LAMBDA(ell, 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]) @@ -351,7 +351,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) + *_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i]) integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2) h = abs(dbeta) / N From 09ae06e9b2ebfca90ca2d6b149769635b8a0a49f Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 15:56:09 +0100 Subject: [PATCH 2/9] ell2 --- GHA_triaxial/gha2_num.py | 267 ++++++++++++++++++++------------------- 1 file changed, 140 insertions(+), 127 deletions(-) diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 8a9f313..c243f4a 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -17,124 +17,6 @@ def sph_azimuth(beta1, lam1, beta2, lam2): a += 2 * np.pi return a -def BETA_LAMBDA(ell, 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(ell, 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(ell, 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, @@ -155,6 +37,137 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float """ # 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) + + # 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 if lamb_1 != lamb_2: N = n @@ -164,7 +177,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float if abs(dlamb) < 1e-15: beta_0 = 0.0 else: - (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph) ode_lamb = buildODElamb() @@ -196,7 +209,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float return False, None, None, None alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2) - (_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1) + (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1) beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) guesses = [ @@ -220,7 +233,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): - (_, _, Ei, Gi, *_) = BETA_LAMBDA(ell, beta_arr_c[i], lamb_arr_c[i]) + (_, _, 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 @@ -253,9 +266,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float beta_p_arr[i] = state[1] (_, _, E1, G1, - *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (_, _, E2, G2, - *_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1]) + *_) = 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]) @@ -263,7 +276,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(ell, 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 @@ -341,9 +354,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float # Azimute (BETA1, LAMBDA1, E1, G1, - *_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0]) + *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (BETA2, LAMBDA2, E2, G2, - *_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1]) + *_) = 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]) @@ -351,7 +364,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float integrand = np.zeros(N + 1) for i in range(N + 1): (_, _, Ei, Gi, - *_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i]) + *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2) h = abs(dbeta) / N From ac1436a7f72022c81b32aa8ee5e34f6cf857bef9 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 15:58:11 +0100 Subject: [PATCH 3/9] stopfitness=1e-14 --- Hansen_ES_CMA.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Hansen_ES_CMA.py b/Hansen_ES_CMA.py index d90b602..6a7311f 100644 --- a/Hansen_ES_CMA.py +++ b/Hansen_ES_CMA.py @@ -9,7 +9,7 @@ def felli(x): return float(np.sum((1e6 ** exponents) * (x ** 2))) -def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-10, stopeval=None, +def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=None, func_args=(), func_kwargs=None, seed=None, bestEver = np.inf, noImproveGen = 0, absTolImprove = 1e-12, maxNoImproveGen = 100, sigmaImprove = 1e-12): From 36e243d6d431000899d89b1afae34f101f23e7d3 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 16:54:28 +0100 Subject: [PATCH 4/9] Anpassungen --- GHA_triaxial/gha2_ES.py | 10 ++++------ Hansen_ES_CMA.py | 5 +++-- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/GHA_triaxial/gha2_ES.py b/GHA_triaxial/gha2_ES.py index d403b02..3ed29aa 100644 --- a/GHA_triaxial/gha2_ES.py +++ b/GHA_triaxial/gha2_ES.py @@ -49,7 +49,7 @@ def midpoint_fitness(x: tuple) -> float: return f -def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, stopeval: int = 2000, maxIter: int = 10000, all_points: bool = False): +def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, maxIter: int = 10000, all_points: bool = False): """ Berechnen der 2. GHA mithilfe der CMA-ES. Die CMA-ES optimiert sukzessive den Mittelpunkt zwischen Start- und Zielpunkt. Der Abbruch der Berechnung erfolgt, wenn alle Segmentlängen <= maxSegLen sind. @@ -58,7 +58,6 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = :param P0: Startpunkt :param Pk: Zielpunkt :param maxSegLen: maximale Segmentlänge - :param stopeval: maximale Durchläufe der CMA-ES :param maxIter: maximale Durchläufe der Mittelpunktsgenerierung :param all_points: Ergebnisliste mit allen Punkte, die wahlweise mit ausgegeben werden kann :return: Richtungswinkel des Start- und Zielpunktes und Gesamtlänge @@ -67,7 +66,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = ell_ES = ell R0 = (ell.ax + ell.ay + ell.b) / 3 if maxSegLen is None: - maxSegLen = R0 * 1 / (637.4) # 10km Segment bei mittleren Erdradius + maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius sigma_uv_nom = 1e-3 * (maxSegLen / R0) # ~1e-5 @@ -101,8 +100,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = sigmaStep = sigma_uv_nom * (Sehne(A, B) / maxSegLen) - u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep, stopfitness=-np.inf, - stopeval=stopeval) + u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) P_next = ell.para2cart(u, v) new_points.append(P_next) @@ -171,7 +169,7 @@ if __name__ == '__main__': beta1, lamb1 = (0.7, 0.3) P1 = ell.ell2cart(beta1, lamb1) - alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=10000, all_points=True) + alpha0, alpha1, s_num, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000, all_points=True) points_num = [] for beta, lamb in zip(betas, lambs): points_num.append(ell.ell2cart(beta, lamb)) diff --git a/Hansen_ES_CMA.py b/Hansen_ES_CMA.py index 6a7311f..055e3ee 100644 --- a/Hansen_ES_CMA.py +++ b/Hansen_ES_CMA.py @@ -9,8 +9,8 @@ def felli(x): return float(np.sum((1e6 ** exponents) * (x ** 2))) -def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=None, - func_args=(), func_kwargs=None, seed=None, +def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000, + func_args=(), func_kwargs=None, seed=0, bestEver = np.inf, noImproveGen = 0, absTolImprove = 1e-12, maxNoImproveGen = 100, sigmaImprove = 1e-12): if func_kwargs is None: @@ -142,6 +142,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=None sigma = sigma * np.exp(0.2 + cs / damps) print(' [CMA-ES] stopfitness erreicht.') #print("warning: flat fitness, consider reformulating the objective") + break #print(f"{counteval}: {arfitness[0]}") From 82444208d71c4ca2bf841081ecf0dc9d0d2852df Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 19:20:09 +0100 Subject: [PATCH 5/9] gha1_ES --- GHA_triaxial/gha1_ES.py | 228 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 GHA_triaxial/gha1_ES.py diff --git a/GHA_triaxial/gha1_ES.py b/GHA_triaxial/gha1_ES.py new file mode 100644 index 0000000..c3ffe11 --- /dev/null +++ b/GHA_triaxial/gha1_ES.py @@ -0,0 +1,228 @@ +from __future__ import annotations +from typing import List, Optional, Tuple +import numpy as np + +from GHA_triaxial.utils import pq_ell +from ellipsoide import EllipsoidTriaxial +from GHA_triaxial.gha1_ana import gha1_ana +from GHA_triaxial.gha1_approx import gha1_approx +from Hansen_ES_CMA import escma +from utils_angle import wrap_to_pi +from numpy.typing import NDArray + +ell: EllipsoidTriaxial = None + + +def ellipsoid_formparameter(ell: EllipsoidTriaxial): + """ + Berechnet die Formparameter des dreiachsigen Ellipsoiden nach Karney (2025), Gl. (2) + :param ell: Ellipsoid + :return: e, k und k' + """ + nenner = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0)) + k = np.sqrt(max(ell.ay * ell.ay - ell.b * ell.b, 0.0)) / nenner + k_ = np.sqrt(max(ell.ax * ell.ax - ell.ay * ell.ay, 0.0)) / nenner + e = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0)) / ell.ay + return e, k, k_ + +def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \ + -> Tuple[NDArray, NDArray, NDArray, float, float, NDArray]: + """ + Analytische ENU-Basis in ellipsoidalen Koordinaten (β, ω) nach Karney. + R(β, ω) gemäß Karney (2025), Eq. (4). :contentReference[oaicite:2]{index=2} + Die Ableitungen E=∂R/∂ω und N=∂R/∂β ergeben sich durch direktes Ableiten. + + Rückgabe: + E_hat, N_hat, U_hat, En=||E||, Nn=||N||, R (XYZ) + """ + + omega = wrap_to_pi(omega) + + cb = np.cos(beta) + sb = np.sin(beta) + co = np.cos(omega) + so = np.sin(omega) + + # D = sqrt(a^2 - c^2) + D2 = ell.ax*ell.ax - ell.b*ell.b + if D2 <= 0: + raise ValueError("Degenerierter Fall a^2 - c^2 <= 0 (nahe Kugel).") + D = np.sqrt(D2) + + # Sx = sqrt(a^2 - b^2 sin^2β - c^2 cos^2β) + Sx2 = ell.ax*ell.ax - ell.ay*ell.ay*(sb*sb) - ell.b*ell.b*(cb*cb) + if Sx2 < 0: # numerische Schutzklemme + Sx2 = 0.0 + Sx = np.sqrt(Sx2) + + # Sz = sqrt(a^2 sin^2ω + b^2 cos^2ω - c^2) + Sz2 = ell.ax*ell.ax*(so*so) + ell.ay*ell.ay*(co*co) - ell.b*ell.b + if Sz2 < 0: + Sz2 = 0.0 + Sz = np.sqrt(Sz2) + + # Karney Eq. (4) + X = ell.ax * co * Sx / D + Y = ell.ay * cb * so + Z = ell.b * sb * Sz / D + R = np.array([X, Y, Z], dtype=float) + + # --- Ableitungen --- + # E = ∂R/∂ω + dX_dw = -ell.ax * so * Sx / D + dY_dw = ell.ay * cb * co + dZ_dw = ell.b * sb * (so * co * (ell.ax*ell.ax - ell.ay*ell.ay) / Sz) / D + E = np.array([dX_dw, dY_dw, dZ_dw], dtype=float) + + # N = ∂R/∂β + dX_db = ell.ax * co * (sb * cb * (ell.b*ell.b - ell.ay*ell.ay) / Sx) / D + dY_db = -ell.ay * sb * so + dZ_db = ell.b * cb * Sz / D + N = np.array([dX_db, dY_db, dZ_db], dtype=float) + + # U ~ Grad(x^2/a^2 + y^2/b^2 + z^2/c^2 - 1) + U = np.array([X/(ell.ax*ell.ax), Y/(ell.ay*ell.ay), Z/(ell.b*ell.b)], dtype=float) + + En = np.linalg.norm(E) + Nn = np.linalg.norm(N) + Un = np.linalg.norm(U) + + N_hat = N / Nn + E_hat = E / En + U_hat = U / Un + + return E_hat, N_hat, U_hat, En, Nn, R + + + +def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTriaxial) -> float: + """ + Jacobi-Konstante nach Karney (2025), Gl. (14): γ = k^2 cos^2β sin^2α − k'^2 sin^2ω cos^2α + :param beta: Beta Koordinate + :param omega: Omega Koordinate + :param alpha: Azimut alpha + :param ell: Ellipsoid + :return: Jacobi-Konstante + """ + e, k, k_ = ellipsoid_formparameter(ell) + cb = np.cos(beta) + so = np.sin(omega) + sa = np.sin(alpha) + ca = np.cos(alpha) + return float((k ** 2) * (cb ** 2) * (sa ** 2) - (k_ ** 2) * (so ** 2) * (ca ** 2)) + + +def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_hat_curr: NDArray, U_hat_curr: NDArray): + v = (P_curr - P_prev).astype(float) + vT = v - float(np.dot(v, U_hat_curr)) * U_hat_curr + vT_hat = vT / np.linalg.norm(vT) + sE = float(np.dot(vT_hat, E_hat_curr)) + sN = float(np.dot(vT_hat, N_hat_curr)) + return wrap_to_pi(float(np.arctan2(sE, sN))) + + +def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: float, gamma0: float, + ell: EllipsoidTriaxial, maxSegLen: float = 10000.0, sigma0: float = None): + + # Startbasis (für Predictor + optionales alpha_start) + E_i, N_i, U_i, En_i, Nn_i, P_i = ENU_beta_omega(beta_i, omega_i, ell) + + # Predictor: dβ ≈ ds cosα / |N|, dω ≈ ds sinα / |E| + d_beta = ds * float(np.cos(alpha_target)) / Nn_i + d_omega = ds * float(np.sin(alpha_target)) / En_i + + beta_pred = beta_i + d_beta + omega_pred = wrap_to_pi(omega_i + d_omega) + + xmean = np.array([beta_pred, omega_pred], dtype=float) + + if sigma0 is None: + R0 = (ell.ax + ell.ay + ell.b) / 3 + sigma0 = 1e-3 * (ds / R0) + + def fitness(x: NDArray) -> float: + beta = x[0] + omega = wrap_to_pi(x[1]) + + P = ell.ell2cart(beta, omega) + d = float(np.linalg.norm(P - P_i)) + + # length penalty + J_len = ((d - ds) / ds) ** 2 + if d > maxSegLen * 1.02: + J_len += 1e3 * ((d / maxSegLen) - 1.02) ** 2 + w_len = 1.0 + + # alpha at end, computed using previous point (for Jacobi gamma) + E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta, omega, ell) + alpha_end = azimuth_at_ESpoint(P_i, P, E_j, N_j, U_j) + + # Jacobi gamma at candidate/end + g_end = jacobi_konstante(beta, omega, alpha_end, ell) + J_gamma = (g_end - gamma0) ** 2 + w_gamma = 10 + + return float(w_len * J_len + w_gamma * J_gamma) + + # Aufruf CMA-ES + xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) + + beta_best = float(np.clip(float(xb[0]), -0.499999 * np.pi, 0.499999 * np.pi)) + omega_best = wrap_to_pi(float(xb[1])) + P_best = ell.ell2cart(beta_best, omega_best) + E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta_best, omega_best, ell) + alpha_end = azimuth_at_ESpoint(P_i, P_best, E_j, N_j, U_j) + + return beta_best, omega_best, P_best, alpha_end + + +def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float): + + beta = float(beta0) + omega = wrap_to_pi(float(omega0)) + alpha = wrap_to_pi(float(alpha0)) + + gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0 + + points: List[NDArray] = [ell.ell2cart(beta, omega)] + + s_acc = 0.0 + step = 0 + nsteps_est = int(np.ceil(s_total / maxSegLen)) + + while s_acc < s_total - 1e-9: + step += 1 + ds = min(maxSegLen, s_total - s_acc) + + print(f"[GHA1-ES] Step {step}/{nsteps_est} ds={ds:.3f} m s_acc={s_acc:.3f} m beta={beta:.6f} omega={omega:.6f} alpha={alpha:.6f}") + + beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_target=alpha, ds=ds, gamma0=gamma0, + ell=ell, maxSegLen=maxSegLen) + + s_acc += ds + points.append(P) + + if step > nsteps_est + 50: + raise RuntimeError("Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.") + + Pk = points[-1] + + return Pk + + +if __name__ == "__main__": + ell = EllipsoidTriaxial.init_name("BursaSima1980round") + s = 1888916.50873 + alpha0 = 70/(180/np.pi) + + P0 = ell.ell2cart(5/(180/np.pi), -90/(180/np.pi)) + point1, alpha1 = gha1_ana(ell, P0, alpha0=alpha0, s=s, maxM=100, maxPartCircum=32) + point1app, alpha1app = gha1_approx(ell, P0, alpha0=alpha0, s=s, ds=1000) + + res = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=1000.0) + + print(point1) + print(res) + # print("alpha1 (am Endpunkt):", res.alpha1) + print(res - point1) + print(point1app - point1, "approx") \ No newline at end of file From e11edd2a43704cf4bc6395c68ffe7e76385dd2e1 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 21:01:06 +0100 Subject: [PATCH 6/9] henreick --- GHA_triaxial/gha1_ES.py | 107 ++++++++++++++++++++++------------------ 1 file changed, 58 insertions(+), 49 deletions(-) diff --git a/GHA_triaxial/gha1_ES.py b/GHA_triaxial/gha1_ES.py index c3ffe11..ac00280 100644 --- a/GHA_triaxial/gha1_ES.py +++ b/GHA_triaxial/gha1_ES.py @@ -1,8 +1,6 @@ from __future__ import annotations from typing import List, Optional, Tuple import numpy as np - -from GHA_triaxial.utils import pq_ell from ellipsoide import EllipsoidTriaxial from GHA_triaxial.gha1_ana import gha1_ana from GHA_triaxial.gha1_approx import gha1_approx @@ -10,8 +8,6 @@ from Hansen_ES_CMA import escma from utils_angle import wrap_to_pi from numpy.typing import NDArray -ell: EllipsoidTriaxial = None - def ellipsoid_formparameter(ell: EllipsoidTriaxial): """ @@ -23,51 +19,43 @@ def ellipsoid_formparameter(ell: EllipsoidTriaxial): k = np.sqrt(max(ell.ay * ell.ay - ell.b * ell.b, 0.0)) / nenner k_ = np.sqrt(max(ell.ax * ell.ax - ell.ay * ell.ay, 0.0)) / nenner e = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0)) / ell.ay + return e, k, k_ + def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \ -> Tuple[NDArray, NDArray, NDArray, float, float, NDArray]: """ - Analytische ENU-Basis in ellipsoidalen Koordinaten (β, ω) nach Karney. - R(β, ω) gemäß Karney (2025), Eq. (4). :contentReference[oaicite:2]{index=2} - Die Ableitungen E=∂R/∂ω und N=∂R/∂β ergeben sich durch direktes Ableiten. - - Rückgabe: - E_hat, N_hat, U_hat, En=||E||, Nn=||N||, R (XYZ) + Analytische ENU-Basis in ellipsoidische Koordinaten (β, ω) nach Karney (2025), S. 2 + :param beta: Beta Koordinate + :param omega: Omega Koordinate + :param ell: Ellipsoid + :return: E_hat = Einheitsrichtung entlang wachsendem ω (East) + N_hat = Einheitsrichtung entlang wachsendem β (North) + U_hat = Einheitsnormale (Up) + En & Nn = Längen der unnormierten Ableitungen + R (XYZ) = Punkt in XYZ """ - + # Berechnungshilfen omega = wrap_to_pi(omega) - cb = np.cos(beta) sb = np.sin(beta) co = np.cos(omega) so = np.sin(omega) - # D = sqrt(a^2 - c^2) - D2 = ell.ax*ell.ax - ell.b*ell.b - if D2 <= 0: - raise ValueError("Degenerierter Fall a^2 - c^2 <= 0 (nahe Kugel).") - D = np.sqrt(D2) - + D = np.sqrt(ell.ax*ell.ax - ell.b*ell.b) # Sx = sqrt(a^2 - b^2 sin^2β - c^2 cos^2β) - Sx2 = ell.ax*ell.ax - ell.ay*ell.ay*(sb*sb) - ell.b*ell.b*(cb*cb) - if Sx2 < 0: # numerische Schutzklemme - Sx2 = 0.0 - Sx = np.sqrt(Sx2) - + Sx = np.sqrt(ell.ax*ell.ax - ell.ay*ell.ay*(sb*sb) - ell.b*ell.b*(cb*cb)) # Sz = sqrt(a^2 sin^2ω + b^2 cos^2ω - c^2) - Sz2 = ell.ax*ell.ax*(so*so) + ell.ay*ell.ay*(co*co) - ell.b*ell.b - if Sz2 < 0: - Sz2 = 0.0 - Sz = np.sqrt(Sz2) + Sz = np.sqrt(ell.ax*ell.ax*(so*so) + ell.ay*ell.ay*(co*co) - ell.b*ell.b) - # Karney Eq. (4) + # Karney Gl. (4) X = ell.ax * co * Sx / D Y = ell.ay * cb * so Z = ell.b * sb * Sz / D R = np.array([X, Y, Z], dtype=float) - # --- Ableitungen --- + # --- Ableitungen - Karney Gl. (5a,b,c)--- # E = ∂R/∂ω dX_dw = -ell.ax * so * Sx / D dY_dw = ell.ay * cb * co @@ -80,7 +68,7 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \ dZ_db = ell.b * cb * Sz / D N = np.array([dX_db, dY_db, dZ_db], dtype=float) - # U ~ Grad(x^2/a^2 + y^2/b^2 + z^2/c^2 - 1) + # U = Grad(x^2/a^2 + y^2/b^2 + z^2/c^2 - 1) U = np.array([X/(ell.ax*ell.ax), Y/(ell.ay*ell.ay), Z/(ell.b*ell.b)], dtype=float) En = np.linalg.norm(E) @@ -94,10 +82,9 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \ return E_hat, N_hat, U_hat, En, Nn, R - def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTriaxial) -> float: """ - Jacobi-Konstante nach Karney (2025), Gl. (14): γ = k^2 cos^2β sin^2α − k'^2 sin^2ω cos^2α + Jacobi-Konstante nach Karney (2025), Gl. (14) :param beta: Beta Koordinate :param omega: Omega Koordinate :param alpha: Azimut alpha @@ -105,28 +92,48 @@ def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTria :return: Jacobi-Konstante """ e, k, k_ = ellipsoid_formparameter(ell) - cb = np.cos(beta) - so = np.sin(omega) - sa = np.sin(alpha) - ca = np.cos(alpha) - return float((k ** 2) * (cb ** 2) * (sa ** 2) - (k_ ** 2) * (so ** 2) * (ca ** 2)) + gamma_jacobi = float((k ** 2) * (np.cos(beta) ** 2) * (np.sin(alpha) ** 2) - (k_ ** 2) * (np.sin(omega) ** 2) * (np.cos(alpha) ** 2)) + + return gamma_jacobi -def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_hat_curr: NDArray, U_hat_curr: NDArray): +def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_hat_curr: NDArray, U_hat_curr: NDArray) -> float: + """ + Berechnet das Azimut in der lokalen Tangentialebene am aktuellen Punkt P_curr, gemessen + an der Bewegungsrichtung vom vorherigen Punkt P_prev nach P_curr. + :param P_prev: vorheriger Punkt + :param P_curr: aktueller Punkt + :param E_hat_curr: Einheitsvektor der lokalen Tangentialrichtung am Punkt P_curr + :param N_hat_curr: Einheitsvektor der lokalen Tangentialrichtung am Punkt P_curr + :param U_hat_curr: Einheitsnormalenvektor am Punkt P_curr + :return: Azimut in Radiant + """ v = (P_curr - P_prev).astype(float) vT = v - float(np.dot(v, U_hat_curr)) * U_hat_curr vT_hat = vT / np.linalg.norm(vT) sE = float(np.dot(vT_hat, E_hat_curr)) sN = float(np.dot(vT_hat, N_hat_curr)) + return wrap_to_pi(float(np.arctan2(sE, sN))) -def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: float, gamma0: float, - ell: EllipsoidTriaxial, maxSegLen: float = 10000.0, sigma0: float = None): +def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: float, gamma0: float, + ell: EllipsoidTriaxial, maxSegLen: float = 10000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]: + """ + + :param beta_i: + :param omega_i: + :param alpha_target: + :param ds: + :param gamma0: + :param ell: + :param maxSegLen: + :param sigma0: + :return: + """ # Startbasis (für Predictor + optionales alpha_start) E_i, N_i, U_i, En_i, Nn_i, P_i = ENU_beta_omega(beta_i, omega_i, ell) - # Predictor: dβ ≈ ds cosα / |N|, dω ≈ ds sinα / |E| d_beta = ds * float(np.cos(alpha_target)) / Nn_i d_omega = ds * float(np.sin(alpha_target)) / En_i @@ -164,8 +171,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: return float(w_len * J_len + w_gamma * J_gamma) - # Aufruf CMA-ES - xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) + xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES beta_best = float(np.clip(float(xb[0]), -0.499999 * np.pi, 0.499999 * np.pi)) omega_best = wrap_to_pi(float(xb[1])) @@ -176,6 +182,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: return beta_best, omega_best, P_best, alpha_end + def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float): beta = float(beta0) @@ -185,6 +192,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0 points: List[NDArray] = [ell.ell2cart(beta, omega)] + alpha_end: List[float] = [alpha] s_acc = 0.0 step = 0 @@ -198,31 +206,32 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_target=alpha, ds=ds, gamma0=gamma0, ell=ell, maxSegLen=maxSegLen) - s_acc += ds points.append(P) - + alpha_end.append(alpha) if step > nsteps_est + 50: raise RuntimeError("Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.") - Pk = points[-1] + alpha1 = alpha_end[-1] + + return Pk, alpha1 - return Pk if __name__ == "__main__": ell = EllipsoidTriaxial.init_name("BursaSima1980round") - s = 1888916.50873 + s = 188891.650873 alpha0 = 70/(180/np.pi) P0 = ell.ell2cart(5/(180/np.pi), -90/(180/np.pi)) point1, alpha1 = gha1_ana(ell, P0, alpha0=alpha0, s=s, maxM=100, maxPartCircum=32) point1app, alpha1app = gha1_approx(ell, P0, alpha0=alpha0, s=s, ds=1000) - res = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=1000.0) + res, alpha = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=1000.0) print(point1) print(res) + print(alpha) # print("alpha1 (am Endpunkt):", res.alpha1) print(res - point1) print(point1app - point1, "approx") \ No newline at end of file From f75672ec36a66e54274607d4391dc66edf897211 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 21:03:59 +0100 Subject: [PATCH 7/9] henreick --- GHA_triaxial/gha1_ES.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/GHA_triaxial/gha1_ES.py b/GHA_triaxial/gha1_ES.py index ac00280..06b11af 100644 --- a/GHA_triaxial/gha1_ES.py +++ b/GHA_triaxial/gha1_ES.py @@ -117,9 +117,8 @@ def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_ return wrap_to_pi(float(np.arctan2(sE, sN))) - def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: float, gamma0: float, - ell: EllipsoidTriaxial, maxSegLen: float = 10000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]: + ell: EllipsoidTriaxial, maxSegLen: float = 1000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]: """ :param beta_i: @@ -183,8 +182,17 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: -def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float): +def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000): + """ + :param ell: + :param beta0: + :param omega0: + :param alpha0: + :param s_total: + :param maxSegLen: + :return: + """ beta = float(beta0) omega = wrap_to_pi(float(omega0)) alpha = wrap_to_pi(float(alpha0)) @@ -227,7 +235,7 @@ if __name__ == "__main__": point1, alpha1 = gha1_ana(ell, P0, alpha0=alpha0, s=s, maxM=100, maxPartCircum=32) point1app, alpha1app = gha1_approx(ell, P0, alpha0=alpha0, s=s, ds=1000) - res, alpha = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=1000.0) + res, alpha = gha1_ES(ell, beta0=5/(180/np.pi), omega0=-90/(180/np.pi), alpha0=alpha0, s_total=s, maxSegLen=1000) print(point1) print(res) From e894c7089a1fd992e740fb25a5f2772837daa413 Mon Sep 17 00:00:00 2001 From: Nico Date: Thu, 5 Feb 2026 21:29:20 +0100 Subject: [PATCH 8/9] final --- GHA_triaxial/gha1_ES.py | 85 ++++++++++++++++++++++------------------- 1 file changed, 45 insertions(+), 40 deletions(-) diff --git a/GHA_triaxial/gha1_ES.py b/GHA_triaxial/gha1_ES.py index 06b11af..fbdb159 100644 --- a/GHA_triaxial/gha1_ES.py +++ b/GHA_triaxial/gha1_ES.py @@ -117,26 +117,27 @@ def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_ return wrap_to_pi(float(np.arctan2(sE, sN))) -def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: float, gamma0: float, - ell: EllipsoidTriaxial, maxSegLen: float = 1000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]: +def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float, gamma0: float, + ell: EllipsoidTriaxial, maxSegLen: float = 1000.0, sigma0: float = None) -> Tuple[float, float, NDArray, float]: """ - - :param beta_i: - :param omega_i: - :param alpha_target: - :param ds: - :param gamma0: - :param ell: - :param maxSegLen: + Berechnung der 1. GHA mithilfe der CMA-ES. + Die CMA-ES optimiert sukzessive einen Punkt, der maxSegLen vom vorherigen Punkt entfernt und zusätzlich auf der + geodätischen Linien liegt. Somit entsteht ein Geodäten ähnlicher Polygonzug auf der Oberfläche des dreiachsigen Ellipsoids. + :param beta_i: Beta Koordinate am Punkt i + :param omega_i: Omega Koordinate am Punkt i + :param alpha_i: Azimut am Punkt i + :param ds: Gesamtlänge + :param gamma0: Jacobi-Konstante am Startpunkt + :param ell: Ellipsoid + :param maxSegLen: maximale Segmentlänge :param sigma0: :return: """ - # Startbasis (für Predictor + optionales alpha_start) + # Startbasis E_i, N_i, U_i, En_i, Nn_i, P_i = ENU_beta_omega(beta_i, omega_i, ell) - # Predictor: dβ ≈ ds cosα / |N|, dω ≈ ds sinα / |E| - d_beta = ds * float(np.cos(alpha_target)) / Nn_i - d_omega = ds * float(np.sin(alpha_target)) / En_i - + # Prediktor: dβ ≈ ds cosα / |N|, dω ≈ ds sinα / |E| + d_beta = ds * float(np.cos(alpha_i)) / Nn_i + d_omega = ds * float(np.sin(alpha_i)) / En_i beta_pred = beta_i + d_beta omega_pred = wrap_to_pi(omega_i + d_omega) @@ -144,36 +145,44 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: if sigma0 is None: R0 = (ell.ax + ell.ay + ell.b) / 3 - sigma0 = 1e-3 * (ds / R0) + sigma0 = 1e-5 * (ds / R0) def fitness(x: NDArray) -> float: + """ + Fitnessfunktion: Fitnesscheck erfolgt anhand der Segmentlänge und der Jacobi-Konstante. + Die Segmentlänge muss möglichst gut zum Sollwert passen. Die Jacobi-Konstante am Punkt x muss zur + Jacobi-Konstanten am Startpunkt passen, damit der Polygonzug auf derselben geodätischen Linie bleibt. + :param x: Koordinate in beta, lambda aus der CMA-ES + :return: Fitnesswert (f) + """ beta = x[0] omega = wrap_to_pi(x[1]) - P = ell.ell2cart(beta, omega) - d = float(np.linalg.norm(P - P_i)) + P = ell.ell2cart(beta, omega) # in kartesischer Koordinaten + d = float(np.linalg.norm(P - P_i)) # Distanz zwischen - # length penalty + # maxSegLen einhalten J_len = ((d - ds) / ds) ** 2 - if d > maxSegLen * 1.02: - J_len += 1e3 * ((d / maxSegLen) - 1.02) ** 2 w_len = 1.0 - # alpha at end, computed using previous point (for Jacobi gamma) + # Azimut für Jacobi-Konstante E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta, omega, ell) alpha_end = azimuth_at_ESpoint(P_i, P, E_j, N_j, U_j) - # Jacobi gamma at candidate/end + # Jacobi-Konstante g_end = jacobi_konstante(beta, omega, alpha_end, ell) J_gamma = (g_end - gamma0) ** 2 w_gamma = 10 - return float(w_len * J_len + w_gamma * J_gamma) + f = float(w_len * J_len + w_gamma * J_gamma) + + return f + xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES - beta_best = float(np.clip(float(xb[0]), -0.499999 * np.pi, 0.499999 * np.pi)) - omega_best = wrap_to_pi(float(xb[1])) + beta_best = xb[0] + omega_best = wrap_to_pi(xb[1]) P_best = ell.ell2cart(beta_best, omega_best) E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta_best, omega_best, ell) alpha_end = azimuth_at_ESpoint(P_i, P_best, E_j, N_j, U_j) @@ -181,17 +190,16 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_target: float, ds: return beta_best, omega_best, P_best, alpha_end - def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000): """ - - :param ell: - :param beta0: - :param omega0: - :param alpha0: - :param s_total: - :param maxSegLen: - :return: + Aufruf der 1. GHA mittels CMA-ES + :param ell: Ellipsoid + :param beta0: Beta Startkoordinate + :param omega0: Omega Startkoordinate + :param alpha0: Azimut Startkoordinate + :param s_total: Gesamtstrecke + :param maxSegLen: maximale Segmentlänge + :return: Zielpunkt Pk und Azimut am Zielpunkt """ beta = float(beta0) omega = wrap_to_pi(float(omega0)) @@ -205,15 +213,13 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_acc = 0.0 step = 0 nsteps_est = int(np.ceil(s_total / maxSegLen)) - while s_acc < s_total - 1e-9: step += 1 ds = min(maxSegLen, s_total - s_acc) - print(f"[GHA1-ES] Step {step}/{nsteps_est} ds={ds:.3f} m s_acc={s_acc:.3f} m beta={beta:.6f} omega={omega:.6f} alpha={alpha:.6f}") - beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_target=alpha, ds=ds, gamma0=gamma0, - ell=ell, maxSegLen=maxSegLen) + beta, omega, P, alpha = optimize_next_point(beta_i=beta, omega_i=omega, alpha_i=alpha, ds=ds, gamma0=gamma0, + ell=ell, maxSegLen=maxSegLen) s_acc += ds points.append(P) alpha_end.append(alpha) @@ -225,7 +231,6 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, return Pk, alpha1 - if __name__ == "__main__": ell = EllipsoidTriaxial.init_name("BursaSima1980round") s = 188891.650873 From a864cd9279baa4a6aa5bbd172232379d43b31830 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Thu, 5 Feb 2026 21:35:48 +0100 Subject: [PATCH 9/9] ES1 in Tests --- Tests/algorithms_test.ipynb | 2401 +++++------------------------------ 1 file changed, 296 insertions(+), 2105 deletions(-) diff --git a/Tests/algorithms_test.ipynb b/Tests/algorithms_test.ipynb index bc8e2b5..88bd2d5 100644 --- a/Tests/algorithms_test.ipynb +++ b/Tests/algorithms_test.ipynb @@ -30,17 +30,22 @@ "%autoreload 2\n", "import time\n", "from numpy import nan\n", + "import numpy as np\n", + "import math\n", "import winkelumrechnungen as wu\n", "import os\n", "from contextlib import contextmanager, redirect_stdout, redirect_stderr\n", "import plotly.graph_objects as go\n", "import warnings\n", + "import pickle\n", + "import random\n", "\n", "from ellipsoide import EllipsoidTriaxial\n", "from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n", "\n", "from GHA_triaxial.gha1_num import gha1_num\n", "from GHA_triaxial.gha1_ana import gha1_ana\n", + "from GHA_triaxial.gha1_ES import gha1_ES\n", "from GHA_triaxial.gha1_approx import gha1_approx\n", "\n", "from GHA_triaxial.gha2_num import gha2_num\n", @@ -82,23 +87,17 @@ }, "cell_type": "code", "source": [ - "# steps_gha1_num = [2000, 5000, 10000, 20000]\n", - "# maxM_gha1_ana = [20, 50]\n", - "# parts_gha1_ana = [4, 8]\n", - "# dsPart_gha1_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", - "#\n", - "# steps_gha2_num = [2000, 5000, 10000, 20000]\n", - "# dsPart_gha2_ES = [600, 1250, 6000] # entspricht bei der Erde ca. 10000, 5000, 1000\n", - "# dsPart_gha2_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", + "# dsPart = [60, 125, 600, 1250, 6000, 60000] entspricht bei der Erde ca. 100km, 50km, 10km, 5km, 1km, 100m\n", "\n", - "steps_gha1_num = [2000, 5000, 10000]\n", + "steps_gha1_num = [20, 50, 100, 200, 500, 1000, 5000]\n", "maxM_gha1_ana = [20, 50]\n", - "parts_gha1_ana = [2, 8, 32]\n", - "dsPart_gha1_approx = [600, 1250, 6000, 60000]\n", + "parts_gha1_ana = [2, 8]\n", + "dsPart_gha1_ES = [60, 600, 1250]\n", + "dsPart_gha1_approx = [600, 1250, 6000]\n", "\n", - "steps_gha2_num = [100, 1000]\n", - "dsPart_gha2_ES = [20, 100]\n", - "dsPart_gha2_approx = [600, 1250, 6000, 60000]" + "steps_gha2_num = [200, 500, 1000]\n", + "dsPart_gha2_ES = [60, 600, 1250]\n", + "dsPart_gha2_approx = [600, 1250, 6000]" ], "id": "96093cdde03f8d57", "outputs": [], @@ -113,11 +112,12 @@ }, "cell_type": "code", "source": [ - "test = \"Karney\"\n", - "# test = \"Panou\"\n", + "# test = \"Karney\"\n", + "test = \"Panou\"\n", + "# test = \"Random\"\n", "if test == \"Karney\":\n", " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n", - " examples = get_examples_karney(10, 42)\n", + " examples = get_examples_karney(4, seed=42)\n", "elif test == \"Panou\":\n", " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n", " tables = get_tables_panou()\n", @@ -126,153 +126,267 @@ " for i, table in enumerate(tables):\n", " for example in table:\n", " table_indices.append(i+1)\n", - " examples.append(example)" + " examples.append(example)\n", + "elif test == \"Random\":\n", + " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n", + " examples = [] # beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s\n", + " random.seed(42)\n", + " for _ in range(10):\n", + " beta0 = wu.deg2rad(random.randint(-90, 90))\n", + " lamb0 = wu.deg2rad(random.randint(-179, 180))\n", + " alpha0_ell = wu.deg2rad(random.randint(0, 359))\n", + " s = random.randint(10000, int(np.pi*ell.b))\n", + " examples.append([beta0, lamb0, alpha0_ell, s])\n", + " pass" ], "id": "6e384cc01c2dbe", "outputs": [], "execution_count": 17 }, { - "metadata": { - "jupyter": { - "is_executing": true - }, - "ExecuteTime": { - "start_time": "2026-02-04T18:11:18.785718Z" - } - }, + "metadata": {}, "cell_type": "code", + "outputs": [], + "execution_count": null, "source": [ - "results = {}\n", - "for i, example in enumerate(examples):\n", - " print(f\"----- Beispiel {i+1}/{len(examples)}\")\n", - " example_results = {}\n", + "if test != \"Random\":\n", + " results = {}\n", + " for i, example in enumerate(examples):\n", + " print(f\"----- Beispiel {i+1}/{len(examples)}\")\n", + " example_results = {}\n", "\n", - " beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n", - " P0 = ell.ell2cart(beta0, lamb0)\n", - " P1 = ell.ell2cart(beta1, lamb1)\n", - " _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n", + " beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n", + " P0 = ell.ell2cart(beta0, lamb0)\n", + " P1 = ell.ell2cart(beta1, lamb1)\n", + " _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n", "\n", - " for steps in steps_gha1_num:\n", - " start = time.perf_counter()\n", - " try:\n", - " P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n", - " end = time.perf_counter()\n", - " beta1_num, lamb1_num = ell.cart2ell(P1_num)\n", - " d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600\n", - " d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600\n", - " d_time = end - start\n", - " example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n", - "\n", - " for maxM in maxM_gha1_ana:\n", - " for parts in parts_gha1_ana:\n", + " for steps in steps_gha1_num:\n", " start = time.perf_counter()\n", " try:\n", - " P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n", + " P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n", + " print(f\"GHA1_num_{steps}\", P1_num)\n", " end = time.perf_counter()\n", - " beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", - " _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", - " d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n", - " d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n", + " beta1_num, lamb1_num = ell.cart2ell(P1_num)\n", + " d_beta1 = abs(beta1_num - beta1)\n", + " d_lamb1 = abs(lamb1_num - lamb1)\n", + " d_alpha1 = abs(alpha1_num_1 - alpha1_ell)\n", " d_time = end - start\n", - " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", " except Exception as e:\n", " print(e)\n", - " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n", + " example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n", "\n", - " for dsPart in dsPart_gha1_approx:\n", - " ds = ell.ax/dsPart\n", - " start = time.perf_counter()\n", - " try:\n", - " P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n", - " end = time.perf_counter()\n", - " beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n", - " d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n", - " d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", - " d_time = end - start\n", - " example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + " for maxM in maxM_gha1_ana:\n", + " for parts in parts_gha1_ana:\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n", + " print(f\"GHA1_ana_{maxM}_{parts}\", P1_ana)\n", + " end = time.perf_counter()\n", + " beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", + " _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", + " d_beta1 = abs(beta1_ana - beta1)\n", + " d_lamb1 = abs(lamb1_ana - lamb1)\n", + " d_alpha1 = abs(alpha1_ana_ell - alpha1_ell)\n", + " d_time = end - start\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n", "\n", - " for steps in steps_gha2_num:\n", - " start = time.perf_counter()\n", - " try:\n", - " with warnings.catch_warnings():\n", - " warnings.simplefilter(\"ignore\", RuntimeWarning)\n", - " alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n", - " end = time.perf_counter()\n", - " d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n", - " d_s = abs(s_num - s) / 1000\n", - " d_time = end - start\n", - " example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n", + " for dsPart in dsPart_gha1_ES:\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_ES, alpha1_ES = gha1_ES(ell, beta0, lamb0, alpha0_ell, s, maxSegLen=dsPart)\n", + " end = time.perf_counter()\n", + " beta1_ES, lamb1_ES = ell.cart2ell(P1_ES)\n", + " d_beta1 = abs(beta1_ES - beta1)\n", + " d_lamb1 = abs(lamb1_ES - lamb1)\n", + " d_alpha1 = abs(alpha1_ES - alpha1_ell)\n", + " d_time = end - start\n", + " example_results[f\"GHA1_ES_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_ES_{dsPart}\"] = (nan, nan, nan, nan)\n", "\n", - " for dsPart in dsPart_gha2_ES:\n", - " ds = ell.ax/dsPart\n", - " start = time.perf_counter()\n", - " try:\n", - " with suppress_print():\n", - " alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n", - " end = time.perf_counter()\n", - " d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n", - " d_s = abs(s_ES - s) / 1000\n", - " d_time = end - start\n", - " example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n", + " for dsPart in dsPart_gha1_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n", + " end = time.perf_counter()\n", + " beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n", + " d_beta1 = abs(beta1_approx - beta1)\n", + " d_lamb1 = abs(lamb1_approx - lamb1)\n", + " d_alpha1 = abs(alpha1_approx - alpha1_ell)\n", + " d_time = end - start\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", "\n", - " for dsPart in dsPart_gha2_approx:\n", - " ds = ell.ax/dsPart\n", - " start = time.perf_counter()\n", - " try:\n", - " alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n", - " end = time.perf_counter()\n", - " d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", - " d_s = abs(s_approx - s) / 1000\n", - " d_time = end - start\n", - " example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + " # ----------------------------------------------\n", "\n", - " results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results" + " for steps in steps_gha2_num:\n", + " start = time.perf_counter()\n", + " try:\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\", RuntimeWarning)\n", + " alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n", + " print(alpha0_num, alpha1_num_2, s_num)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(alpha0_num - alpha0_ell)\n", + " d_alpha1 = abs(alpha1_num_2 - alpha1_ell)\n", + " d_s = abs(s_num - s)\n", + " d_time = end - start\n", + " example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha2_ES:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " with suppress_print():\n", + " alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(alpha0_ES - alpha0_ell)\n", + " d_alpha1 = abs(alpha1_ES - alpha1_ell)\n", + " d_s = abs(s_ES - s)\n", + " d_time = end - start\n", + " example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha2_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(alpha0_approx - alpha0_ell)\n", + " d_alpha1 = abs(alpha1_approx - alpha1_ell)\n", + " d_s = abs(s_approx - s)\n", + " d_time = end - start\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results" ], - "id": "ef8849908ed3231e", - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "----- Beispiel 1/10\n", - "Keine Multi-Start-Variante konvergiert.\n", - "Keine Multi-Start-Variante konvergiert.\n", - "----- Beispiel 2/10\n", - "Keine Multi-Start-Variante konvergiert.\n", - "Keine Multi-Start-Variante konvergiert.\n", - "----- Beispiel 3/10\n", - "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", - "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", - "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", - "Fehler in der Umrechnung cart2ell\n", - "Keine Multi-Start-Variante konvergiert.\n", - "Keine Multi-Start-Variante konvergiert.\n" - ] - } + "id": "2a87e028089a215" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "if test == \"Random\":\n", + " results = {}\n", + " for i, example in enumerate(examples):\n", + " print(f\"----- Beispiel {i+1}/{len(examples)}\")\n", + " example_results = {}\n", + "\n", + " beta0, lamb0, alpha0_ell, s = example\n", + " P0 = ell.ell2cart(beta0, lamb0)\n", + " _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n", + "\n", + " # P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=60, maxPartCircum=4)\n", + " # beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", + " # _, _, alpha1_ana = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", + "\n", + " P1, alpha1_ell = gha1_num(ell, P0, alpha0_ell, s, num=10000)\n", + " beta1, lamb1 = ell.cart2ell(P1)\n", + "\n", + " for maxM in maxM_gha1_ana:\n", + " for parts in parts_gha1_ana:\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n", + " end = time.perf_counter()\n", + " beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", + " _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", + " d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n", + " d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n", + " d_time = end - start\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha1_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n", + " end = time.perf_counter()\n", + " beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n", + " d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n", + " d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", + " d_time = end - start\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " # for steps in steps_gha2_num:\n", + " # start = time.perf_counter()\n", + " # try:\n", + " # with warnings.catch_warnings():\n", + " # warnings.simplefilter(\"ignore\", RuntimeWarning)\n", + " # alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n", + " # end = time.perf_counter()\n", + " # d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n", + " # d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n", + " # d_s = abs(s_num - s) / 1000\n", + " # d_time = end - start\n", + " # example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " # except Exception as e:\n", + " # print(e)\n", + " # example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n", + " #\n", + " # for dsPart in dsPart_gha2_ES:\n", + " # ds = ell.ax/dsPart\n", + " # start = time.perf_counter()\n", + " # try:\n", + " # with suppress_print():\n", + " # alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n", + " # end = time.perf_counter()\n", + " # d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n", + " # d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n", + " # d_s = abs(s_ES - s) / 1000\n", + " # d_time = end - start\n", + " # example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " # except Exception as e:\n", + " # print(e)\n", + " # example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha2_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", + " d_s = abs(s_approx - s) / 1000\n", + " d_time = end - start\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results\n", + "\n", + "\n" ], - "execution_count": null + "id": "d8868b5f0e6f9b42" }, { "metadata": { @@ -307,13 +421,41 @@ "execution_count": 8 }, { - "metadata": { - "ExecuteTime": { - "end_time": "2026-02-04T18:01:30.114488Z", - "start_time": "2026-02-04T18:01:25.723109Z" - } - }, + "metadata": {}, "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": [ + "metrics_gha1 = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n", + "metrics_gha2 = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n", + "listed_results = {}\n", + "for example, example_metrics in results.items():\n", + " for method, method_metrics in example_metrics.items():\n", + " if \"GHA1\" in method:\n", + " if method not in listed_results.keys():\n", + " listed_results[method] = {metric: [] for metric in metrics_gha1}\n", + " for i, metric in enumerate(method_metrics):\n", + " if '[\"]' in metrics_gha1[i]:\n", + " listed_results[method][metrics_gha1[i]].append(wu.rad2deg(metric)*3600)\n", + " else:\n", + " listed_results[method][metrics_gha1[i]].append(metric)\n", + " if \"GHA2\" in method:\n", + " if method not in listed_results.keys():\n", + " listed_results[method] = {metric: [] for metric in metrics_gha2}\n", + " for i, metric in enumerate(method_metrics):\n", + " if '[\"]' in metrics_gha2[i]:\n", + " listed_results[method][metrics_gha2[i]].append(wu.rad2deg(metric)*3600)\n", + " else:\n", + " listed_results[method][metrics_gha2[i]].append(metric)\n", + " pass" + ], + "id": "8fd6f220440f4494" + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, "source": [ "def format_max(values, is_angle=False):\n", " arr = np.array(values, dtype=float)\n", @@ -322,10 +464,6 @@ " maxi = np.nanmax(np.abs(arr))\n", " if maxi is None or (isinstance(maxi,float) and (math.isnan(maxi))):\n", " return \"nan\"\n", - " # i = 2\n", - " # while maxi > np.nanmean(np.abs(arr)):\n", - " # maxi = np.sort(np.abs(arr[~np.isnan(arr)]))[-i]\n", - " # i += 1\n", " if is_angle:\n", " maxi = wu.rad2deg(maxi)*3600\n", " if f\"{maxi:.3g}\" == 0:\n", @@ -334,7 +472,6 @@ "\n", "\n", "def build_max_table(gha_prefix, title, group_value = None):\n", - " # metrics\n", " if gha_prefix==\"GHA1\":\n", " metrics = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n", " angle_mask = [True, True, True, False]\n", @@ -342,12 +479,11 @@ " metrics = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n", " angle_mask = [True, True, False, False]\n", "\n", - " # collect keys in this group\n", " if group_value is None:\n", - " example_keys = [example_key for example_key in results.keys()]\n", + " example_keys = [example_key for example_key in list(results.keys())]\n", " else:\n", " example_keys = [example_key for example_key, group_index in zip(results.keys(), table_indices) if group_index==group_value]\n", - " # all variant keys (inner dict keys) matching prefix\n", + "\n", " algorithms = sorted({algorithm for example_key in example_keys for algorithm in results[example_key].keys() if algorithm.startswith(gha_prefix)})\n", "\n", " header = [\"Algorithmus\", \"Parameter\", \"NaN\"] + list(metrics)\n", @@ -402,1950 +538,7 @@ "for fig in figs:\n", " fig.show()" ], - "id": "b46d57fc0d794e28", - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "C:\\Users\\moell\\AppData\\Local\\Temp\\ipykernel_25588\\2480185505.py:5: RuntimeWarning:\n", - "\n", - "All-NaN slice encountered\n", - "\n" - ] - }, - { - "data": { - "application/vnd.plotly.v1+json": { - "data": [ - { - "cells": { - "align": "center", - "values": [ - [ - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "ana", - "approx", - "approx", - "num", - "num" - ], - [ - "20_16", - "20_2", - "20_24", - "20_4", - "20_8", - "40_16", - "40_2", - "40_24", - "40_4", - "40_8", - "60_16", - "60_2", - "60_24", - "60_4", - "60_8", - "1250", - "600", - "2000", - "5000" - ], - [ - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 - ], - [ - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "5.1e+03", - "7.16e+03", - "7.15e+03", - "3.67e+03", - "3.67e+03" - ], - [ - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.06e+04", - "1.67e+04", - "1.67e+04", - "3.99e+03", - "3.99e+03" - ], - [ - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "1.55e+04", - "2.47e+03", - "2.46e+03", - "1.82e+03", - "1.82e+03" - ], - [ - "0.02", - "0.00746", - "0.0353", - "0.0123", - "0.0103", - "0.0868", - "0.0558", - "0.201", - "0.0561", - "0.0435", - "0.251", - "0.0799", - "0.495", - "0.0816", - "0.126", - "0.731", - "0.294", - "0.0847", - "0.231" - ] - ] - }, - "header": { - "align": "center", - "fill": { - "color": "lightgrey" - }, - "font": { - "size": 13 - }, - "values": [ - "Algorithmus", - "Parameter", - "NaN", - "dBeta [\"]", - "dLambda [\"]", - "dAlpha1 [\"]", - "time [s]" - ] - }, - "type": "table" - } - ], - "layout": { - "template": { - "data": { - "barpolar": [ - { - "marker": { - "line": { - "color": "white", - "width": 0.5 - }, - "pattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - } - }, - "type": "barpolar" - } - ], - "bar": [ - { - "error_x": { - "color": "rgb(36,36,36)" - }, - "error_y": { - "color": "rgb(36,36,36)" - }, - "marker": { - "line": { - "color": "white", - "width": 0.5 - }, - "pattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - } - }, - "type": "bar" - } - ], - "carpet": [ - { - "aaxis": { - "endlinecolor": "rgb(36,36,36)", - "gridcolor": "white", - "linecolor": "white", - "minorgridcolor": "white", - "startlinecolor": "rgb(36,36,36)" - }, - "baxis": { - "endlinecolor": "rgb(36,36,36)", - "gridcolor": "white", - "linecolor": "white", - "minorgridcolor": "white", - "startlinecolor": "rgb(36,36,36)" - }, - "type": "carpet" - } - ], - "choropleth": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "choropleth" - } - ], - "contourcarpet": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "contourcarpet" - } - ], - "contour": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "contour" - } - ], - "heatmap": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "heatmap" - } - ], - "histogram2dcontour": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "histogram2dcontour" - } - ], - "histogram2d": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "histogram2d" - } - ], - "histogram": [ - { - "marker": { - "line": { - "color": "white", - "width": 0.6 - } - }, - "type": "histogram" - } - ], - "mesh3d": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "mesh3d" - } - ], - "parcoords": [ - { - "line": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "parcoords" - } - ], - "pie": [ - { - "automargin": true, - "type": "pie" - } - ], - "scatter3d": [ - { - "line": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatter3d" - } - ], - "scattercarpet": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattercarpet" - } - ], - "scattergeo": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattergeo" - } - ], - "scattergl": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattergl" - } - ], - "scattermapbox": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattermapbox" - } - ], - "scattermap": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattermap" - } - ], - "scatterpolargl": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterpolargl" - } - ], - "scatterpolar": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterpolar" - } - ], - "scatter": [ - { - "fillpattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - }, - "type": "scatter" - } - ], - "scatterternary": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterternary" - } - ], - "surface": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "surface" - } - ], - "table": [ - { - "cells": { - "fill": { - "color": "rgb(237,237,237)" - }, - "line": { - "color": "white" - } - }, - "header": { - "fill": { - "color": "rgb(217,217,217)" - }, - "line": { - "color": "white" - } - }, - "type": "table" - } - ] - }, - "layout": { - "annotationdefaults": { - "arrowhead": 0, - "arrowwidth": 1 - }, - "autotypenumbers": "strict", - "coloraxis": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "colorscale": { - "diverging": [ - [ - 0.0, - "rgb(103,0,31)" - ], - [ - 0.1, - "rgb(178,24,43)" - ], - [ - 0.2, - "rgb(214,96,77)" - ], - [ - 0.3, - "rgb(244,165,130)" - ], - [ - 0.4, - "rgb(253,219,199)" - ], - [ - 0.5, - "rgb(247,247,247)" - ], - [ - 0.6, - "rgb(209,229,240)" - ], - [ - 0.7, - "rgb(146,197,222)" - ], - [ - 0.8, - "rgb(67,147,195)" - ], - [ - 0.9, - "rgb(33,102,172)" - ], - [ - 1.0, - "rgb(5,48,97)" - ] - ], - "sequential": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "sequentialminus": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ] - }, - "colorway": [ - "#1F77B4", - "#FF7F0E", - "#2CA02C", - "#D62728", - "#9467BD", - "#8C564B", - "#E377C2", - "#7F7F7F", - "#BCBD22", - "#17BECF" - ], - "font": { - "color": "rgb(36,36,36)" - }, - "geo": { - "bgcolor": "white", - "lakecolor": "white", - "landcolor": "white", - "showlakes": true, - "showland": true, - "subunitcolor": "white" - }, - "hoverlabel": { - "align": "left" - }, - "hovermode": "closest", - "mapbox": { - "style": "light" - }, - "paper_bgcolor": "white", - "plot_bgcolor": "white", - "polar": { - "angularaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "bgcolor": "white", - "radialaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - } - }, - "scene": { - "xaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "yaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "zaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - } - }, - "shapedefaults": { - "fillcolor": "black", - "line": { - "width": 0 - }, - "opacity": 0.3 - }, - "ternary": { - "aaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "baxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "bgcolor": "white", - "caxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - } - }, - "title": { - "x": 0.05 - }, - "xaxis": { - "automargin": true, - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside", - "title": { - "standoff": 15 - }, - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "yaxis": { - "automargin": true, - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside", - "title": { - "standoff": 15 - }, - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - } - } - }, - "margin": { - "l": 20, - "r": 20, - "t": 60, - "b": 20 - }, - "title": { - "text": "Karney - GHA1" - }, - "width": 800, - "height": 280 - }, - "config": { - "plotlyServerURL": "https://plot.ly" - } - } - }, - "metadata": {}, - "output_type": "display_data", - "jetTransient": { - "display_id": null - } - }, - { - "data": { - "application/vnd.plotly.v1+json": { - "data": [ - { - "cells": { - "align": "center", - "values": [ - [ - "ES", - "approx", - "num" - ], - [ - "20", - "600", - "2000" - ], - [ - 0, - 0, - 2 - ], - [ - "1.07e+04", - "9.95e+03", - "nan" - ], - [ - "4.35e+03", - "2.77e+03", - "nan" - ], - [ - "0.000886", - "0.000888", - "nan" - ], - [ - "1.47", - "0.108", - "nan" - ] - ] - }, - "header": { - "align": "center", - "fill": { - "color": "lightgrey" - }, - "font": { - "size": 13 - }, - "values": [ - "Algorithmus", - "Parameter", - "NaN", - "dAlpha0 [\"]", - "dAlpha1 [\"]", - "dStrecke [m]", - "time [s]" - ] - }, - "type": "table" - } - ], - "layout": { - "template": { - "data": { - "barpolar": [ - { - "marker": { - "line": { - "color": "white", - "width": 0.5 - }, - "pattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - } - }, - "type": "barpolar" - } - ], - "bar": [ - { - "error_x": { - "color": "rgb(36,36,36)" - }, - "error_y": { - "color": "rgb(36,36,36)" - }, - "marker": { - "line": { - "color": "white", - "width": 0.5 - }, - "pattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - } - }, - "type": "bar" - } - ], - "carpet": [ - { - "aaxis": { - "endlinecolor": "rgb(36,36,36)", - "gridcolor": "white", - "linecolor": "white", - "minorgridcolor": "white", - "startlinecolor": "rgb(36,36,36)" - }, - "baxis": { - "endlinecolor": "rgb(36,36,36)", - "gridcolor": "white", - "linecolor": "white", - "minorgridcolor": "white", - "startlinecolor": "rgb(36,36,36)" - }, - "type": "carpet" - } - ], - "choropleth": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "choropleth" - } - ], - "contourcarpet": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "contourcarpet" - } - ], - "contour": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "contour" - } - ], - "heatmap": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "heatmap" - } - ], - "histogram2dcontour": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "histogram2dcontour" - } - ], - "histogram2d": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "histogram2d" - } - ], - "histogram": [ - { - "marker": { - "line": { - "color": "white", - "width": 0.6 - } - }, - "type": "histogram" - } - ], - "mesh3d": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "type": "mesh3d" - } - ], - "parcoords": [ - { - "line": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "parcoords" - } - ], - "pie": [ - { - "automargin": true, - "type": "pie" - } - ], - "scatter3d": [ - { - "line": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatter3d" - } - ], - "scattercarpet": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattercarpet" - } - ], - "scattergeo": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattergeo" - } - ], - "scattergl": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattergl" - } - ], - "scattermapbox": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattermapbox" - } - ], - "scattermap": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scattermap" - } - ], - "scatterpolargl": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterpolargl" - } - ], - "scatterpolar": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterpolar" - } - ], - "scatter": [ - { - "fillpattern": { - "fillmode": "overlay", - "size": 10, - "solidity": 0.2 - }, - "type": "scatter" - } - ], - "scatterternary": [ - { - "marker": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "type": "scatterternary" - } - ], - "surface": [ - { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - }, - "colorscale": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "type": "surface" - } - ], - "table": [ - { - "cells": { - "fill": { - "color": "rgb(237,237,237)" - }, - "line": { - "color": "white" - } - }, - "header": { - "fill": { - "color": "rgb(217,217,217)" - }, - "line": { - "color": "white" - } - }, - "type": "table" - } - ] - }, - "layout": { - "annotationdefaults": { - "arrowhead": 0, - "arrowwidth": 1 - }, - "autotypenumbers": "strict", - "coloraxis": { - "colorbar": { - "outlinewidth": 1, - "tickcolor": "rgb(36,36,36)", - "ticks": "outside" - } - }, - "colorscale": { - "diverging": [ - [ - 0.0, - "rgb(103,0,31)" - ], - [ - 0.1, - "rgb(178,24,43)" - ], - [ - 0.2, - "rgb(214,96,77)" - ], - [ - 0.3, - "rgb(244,165,130)" - ], - [ - 0.4, - "rgb(253,219,199)" - ], - [ - 0.5, - "rgb(247,247,247)" - ], - [ - 0.6, - "rgb(209,229,240)" - ], - [ - 0.7, - "rgb(146,197,222)" - ], - [ - 0.8, - "rgb(67,147,195)" - ], - [ - 0.9, - "rgb(33,102,172)" - ], - [ - 1.0, - "rgb(5,48,97)" - ] - ], - "sequential": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ], - "sequentialminus": [ - [ - 0.0, - "#440154" - ], - [ - 0.1111111111111111, - "#482878" - ], - [ - 0.2222222222222222, - "#3e4989" - ], - [ - 0.3333333333333333, - "#31688e" - ], - [ - 0.4444444444444444, - "#26828e" - ], - [ - 0.5555555555555556, - "#1f9e89" - ], - [ - 0.6666666666666666, - "#35b779" - ], - [ - 0.7777777777777778, - "#6ece58" - ], - [ - 0.8888888888888888, - "#b5de2b" - ], - [ - 1.0, - "#fde725" - ] - ] - }, - "colorway": [ - "#1F77B4", - "#FF7F0E", - "#2CA02C", - "#D62728", - "#9467BD", - "#8C564B", - "#E377C2", - "#7F7F7F", - "#BCBD22", - "#17BECF" - ], - "font": { - "color": "rgb(36,36,36)" - }, - "geo": { - "bgcolor": "white", - "lakecolor": "white", - "landcolor": "white", - "showlakes": true, - "showland": true, - "subunitcolor": "white" - }, - "hoverlabel": { - "align": "left" - }, - "hovermode": "closest", - "mapbox": { - "style": "light" - }, - "paper_bgcolor": "white", - "plot_bgcolor": "white", - "polar": { - "angularaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "bgcolor": "white", - "radialaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - } - }, - "scene": { - "xaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "yaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "zaxis": { - "backgroundcolor": "white", - "gridcolor": "rgb(232,232,232)", - "gridwidth": 2, - "linecolor": "rgb(36,36,36)", - "showbackground": true, - "showgrid": false, - "showline": true, - "ticks": "outside", - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - } - }, - "shapedefaults": { - "fillcolor": "black", - "line": { - "width": 0 - }, - "opacity": 0.3 - }, - "ternary": { - "aaxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "baxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - }, - "bgcolor": "white", - "caxis": { - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside" - } - }, - "title": { - "x": 0.05 - }, - "xaxis": { - "automargin": true, - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside", - "title": { - "standoff": 15 - }, - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - }, - "yaxis": { - "automargin": true, - "gridcolor": "rgb(232,232,232)", - "linecolor": "rgb(36,36,36)", - "showgrid": false, - "showline": true, - "ticks": "outside", - "title": { - "standoff": 15 - }, - "zeroline": false, - "zerolinecolor": "rgb(36,36,36)" - } - } - }, - "margin": { - "l": 20, - "r": 20, - "t": 60, - "b": 20 - }, - "title": { - "text": "Karney - GHA2" - }, - "width": 800, - "height": 280 - }, - "config": { - "plotlyServerURL": "https://plot.ly" - } - } - }, - "metadata": {}, - "output_type": "display_data", - "jetTransient": { - "display_id": null - } - } - ], - "execution_count": 10 + "id": "f1fb73407f5b1c8a" }, { "metadata": { @@ -2356,8 +549,6 @@ }, "cell_type": "code", "source": [ - "import numpy as np\n", - "import math\n", "from collections import defaultdict\n", "\n", "def to_latex_sci(x, sig=3):\n", @@ -2574,13 +765,13 @@ "\n", "\n", "# --- Beispielaufruf ---\n", - "example_keys = list(results.keys()) # oder gefiltert auf Panou Gruppe etc.\n", + "example_keys = list(key for i, key in enumerate(results.keys()) if table_indices[i] == 3) # oder gefiltert auf Panou Gruppe etc.\n", "latex = build_latex_table_from_results(\n", " results,\n", - " gha_prefix=\"GHA2\",\n", + " gha_prefix=\"GHA1\",\n", " example_keys=example_keys,\n", - " caption=r\"Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid\",\n", - " label=\"tab:results_algorithms\",\n", + " caption=r\"Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid (Gruppe 3)\",\n", + " label=\"tab:results_gha1_g1\",\n", " include_nan_col=False, # True, wenn du NaN-Spalte willst\n", " wu=wu # wichtig, falls Winkelwerte rad sind\n", ")\n",