Merge remote-tracking branch 'origin/main'

GHA_triaxial/gha1_num.py

@@ -12,29 +12,7 @@ from numpy.typing import NDArray
 from GHA_triaxial.utils import alpha_ell2para, pq_ell
 
 
-def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, num: int, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, List]:
+def buildODE(ell: EllipsoidTriaxial) -> Callable:
-    """
-    Panou, Korakitits 2019
-    :param ell: Ellipsoid
-    :param point: Punkt in kartesischen Koordinaten
-    :param alpha0: Azimut im Startpunkt
-    :param s: Strecke
-    :param num: Anzahl Zwischenpunkte
-    :param all_points: Ausgabe aller Punkte?
-    :return: Zielpunkt, Azimut im Zielpunkt (, alle Punkte)
-    """
-    phi, lam, _ = ell.cart2geod(point, "ligas3")
-    p0 = ell.geod2cart(phi, lam, 0)
-    x0, y0, z0 = p0
-
-    p, q = pq_ell(ell, p0)
-    dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0)
-    dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0)
-    dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0)
-
-    v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0])
-
-    def buildODE(ell: EllipsoidTriaxial) -> Callable:
     """
     Aufbau des DGL-Systems
     :param ell: Ellipsoid
@@ -61,6 +39,28 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu
 
     return ODE
 
+def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, num: int, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, List]:
+    """
+    Panou, Korakitits 2019
+    :param ell: Ellipsoid
+    :param point: Punkt in kartesischen Koordinaten
+    :param alpha0: Azimut im Startpunkt
+    :param s: Strecke
+    :param num: Anzahl Zwischenpunkte
+    :param all_points: Ausgabe aller Punkte?
+    :return: Zielpunkt, Azimut im Zielpunkt (, alle Punkte)
+    """
+    phi, lam, _ = ell.cart2geod(point, "ligas3")
+    p0 = ell.geod2cart(phi, lam, 0)
+    x0, y0, z0 = p0
+
+    p, q = pq_ell(ell, p0)
+    dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0)
+    dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0)
+    dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0)
+
+    v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0])
+
     ode = buildODE(ell)
 
     _, werte = rk.rk4(ode, 0, v_init, s, num)

GHA_triaxial/gha2_num.py

@@ -1,42 +1,13 @@
 import numpy as np
 from ellipsoide import EllipsoidTriaxial
 import runge_kutta as rk
-import GHA_triaxial.numeric_examples_karney as ne_karney
-import GHA_triaxial.numeric_examples_panou as ne_panou
-import winkelumrechnungen as wu
 from typing import Tuple
 from numpy.typing import NDArray
 
-# Panou 2013
+from utils_angle import arccot, cot, wrap_to_pi
-def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float,
-             n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False
-             ) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]:
-    """
 
-    :param ell: triaxiales Ellipsoid
-    :param beta_1: reduzierte ellipsoidische Breite Punkt 1
-    :param lamb_1: elllipsoidische Länge Punkt 1
-    :param beta_2: reduzierte ellipsoidische Breite Punkt 2
-    :param lamb_2: elllipsoidische Länge Punkt 2
-    :param n: Anzahl Schritte
-    :param epsilon:
-    :param iter_max: Maximale Anzhal Iterationen
-    :param all_points:
-    :return:
-    """
 
-    # h_x, h_y, h_e entsprechen E_x, E_y, E_e
+def sph_azimuth(beta1, lam1, beta2, lam2):
-
-    def arccot(x):
-        return np.arctan2(1.0, x)
-
-    def cot(a):
-        return np.cos(a) / np.sin(a)
-
-    def wrap_to_pi(x):
-        return (x + np.pi) % (2 * np.pi) - np.pi
-
-    def sph_azimuth(beta1, lam1, beta2, lam2):
     # sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß)
     dlam = wrap_to_pi(lam2 - lam1)
     y = np.sin(dlam) * np.cos(beta2)
@@ -46,7 +17,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
         a += 2 * np.pi
     return a
 
-    def BETA_LAMBDA(beta, lamb):
+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)
@@ -86,7 +57,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
             E_beta_beta, E_beta_lamb, E_lamb_lamb,
             G_beta_beta, G_beta_lamb, G_lamb_lamb)
 
-    def p_coef(beta, lamb):
+def p_coef(beta, lamb):
 
     (BETA, LAMBDA, E, G,
      BETA_, LAMBDA_, BETA__, LAMBDA__,
@@ -108,7 +79,24 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
             p_3, p_2, p_1, p_0,
             p_33, p_22, p_11, p_00)
 
-    def q_coef(beta, lamb):
+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__,
@@ -130,49 +118,46 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
             q_3, q_2, q_1, q_0,
             q_33, q_22, q_11, q_00)
 
-    if lamb_1 != lamb_2:
+def buildODEbeta():
-        # def functions():
+    def ODE(beta, v):
-        #     def f_beta(lamb, beta, beta_p, X3, X4):
+        lamb, lamb_p, Y3, Y4 = v
-        #         return beta_p
-        #
-        #     def f_beta_p(lamb, beta, beta_p, X3, X4):
-        #         (BETA, LAMBDA, E, G,
-        #          p_3, p_2, p_1, p_0,
-        #          p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
-        #         return p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
-        #
-        #     def f_X3(lamb, beta, beta_p, X3, X4):
-        #         return X4
-        #
-        #     def f_X4(lamb, beta, beta_p, X3, X4):
-        #         (BETA, LAMBDA, E, G,
-        #          p_3, p_2, p_1, p_0,
-        #          p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
-        #         return (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 [f_beta, f_beta_p, f_X3, f_X4]
-
-        def buildODElamb():
-            def ODE(lamb, v):
-                beta, beta_p, X3, X4 = v
 
         (BETA, LAMBDA, E, G,
-                 p_3, p_2, p_1, p_0,
+         q_3, q_2, q_1, q_0,
-                 p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
+         q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
 
-                dbeta = beta_p
+        dlamb = lamb_p
-                dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
+        dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
-                dX3 = X4
+        dY3 = Y4
-                dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
+        dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
-                    (3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
+              (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
 
-                return np.array([dbeta, dbeta_p, dX3, dX4])
+        return np.array([dlamb, dlamb_p, dY3, dY4])
 
     return ODE
 
-        N = n
+# Panou 2013
+def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float,
+             n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False
+             ) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]:
+    """
 
+    :param ell: triaxiales Ellipsoid
+    :param beta_1: reduzierte ellipsoidische Breite Punkt 1
+    :param lamb_1: elllipsoidische Länge Punkt 1
+    :param beta_2: reduzierte ellipsoidische Breite Punkt 2
+    :param lamb_2: elllipsoidische Länge Punkt 2
+    :param n: Anzahl Schritte
+    :param epsilon:
+    :param iter_max: Maximale Anzhal Iterationen
+    :param all_points:
+    :return:
+    """
+
+    # h_x, h_y, h_e entsprechen E_x, E_y, E_e
+
+    if lamb_1 != lamb_2:
+        N = n
         dlamb = lamb_2 - lamb_1
         alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
 
@@ -182,10 +167,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
             (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
             beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph)
 
-        converged = False
-        iterations = 0
-
-        # funcs = functions()
         ode_lamb = buildODElamb()
 
         def solve_newton(beta_p0_init: float):
@@ -307,11 +288,10 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
             return alpha_1, alpha_2, s
 
     else:  # lamb_1 == lamb_2
-
         N = n
         dbeta = beta_2 - beta_1
 
-        if abs(dbeta) < 10**-15:
+        if abs(dbeta) < 1e-15:
             if all_points:
                 return 0, 0, 0, np.array([]), np.array([])
             else:
@@ -319,68 +299,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
 
         lamb_0 = 0
 
-        converged = False
-        iterations = 0
-
-        # def functions_beta():
-        #     def g_lamb(beta, lamb, lamb_p, Y3, Y4):
-        #         return lamb_p
-        #
-        #     def g_lamb_p(beta, lamb, lamb_p, Y3, Y4):
-        #         (BETA, LAMBDA, E, G,
-        #          q_3, q_2, q_1, q_0,
-        #          q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
-        #         return q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
-        #
-        #     def g_Y3(beta, lamb, lamb_p, Y3, Y4):
-        #         return Y4
-        #
-        #     def g_Y4(beta, lamb, lamb_p, Y3, Y4):
-        #         (BETA, LAMBDA, E, G,
-        #          q_3, q_2, q_1, q_0,
-        #          q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
-        #         return (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 [g_lamb, g_lamb_p, g_Y3, g_Y4]
-
-        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
-
-        # funcs_beta = functions_beta()
         ode_beta = buildODEbeta()
 
         for i in range(iter_max):
-            iterations = i + 1
-
             startwerte = [lamb_1, lamb_0, 0.0, 1.0]
 
-            # werte = rk.verfahren(funcs_beta, startwerte, dbeta, N, False)
             beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
 
             beta_end = beta_list[-1]
-            # beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
             lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
 
             d_lamb_end_d_lambda0 = Y3_end
             delta = lamb_end - lamb_2
 
             if abs(delta) < epsilon:
-                converged = True
                 break
 
             if abs(d_lamb_end_d_lambda0) < 1e-20:
@@ -393,7 +325,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
 
             lamb_0 = lamb_0 - step
 
-        # werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N, False)
         beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False)
 
         # beta_arr = np.zeros(N + 1)

utils_angle.py

@@ -0,0 +1,13 @@
+import numpy as np
+
+
+def arccot(x):
+    return np.arctan2(1.0, x)
+
+
+def cot(a):
+    return np.cos(a) / np.sin(a)
+
+
+def wrap_to_pi(x):
+    return (x + np.pi) % (2 * np.pi) - np.pi