Algorithmen Test

GHA_triaxial/ES_gha2.py

@@ -128,14 +128,14 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, stepLenTarget: flo
         P_prev = P_next
 
     print('Maximale Schrittanzahl erreicht.')
-    P_all.append(P_end)
+    # P_all.append(P_end)
     totalLen += Bogenlaenge(P_prev, P_end)
 
-    p0i = ell.point_onto_ellipsoid(P0 + 10 * (P_all[1] - P0) / np.linalg.norm(P_all[1] - P0))
+    p0i = ell.point_onto_ellipsoid(P0 + stepLenTarget/1000 * (P_all[1] - P0) / np.linalg.norm(P_all[1] - P0))
     sigma0 = (p0i - P0) / np.linalg.norm(p0i - P0)
     alpha0 = sigma2alpha(ell_ES, sigma0, P0)
 
-    p1i = ell.point_onto_ellipsoid(Pk - 10 * (Pk - P_all[-2]) / np.linalg.norm(Pk - P_all[-2]))
+    p1i = ell.point_onto_ellipsoid(Pk - stepLenTarget/1000 * (Pk - P_all[-2]) / np.linalg.norm(Pk - P_all[-2]))
     sigma1 = (Pk - p1i) / np.linalg.norm(Pk - p1i)
     alpha1 = sigma2alpha(ell_ES, sigma1, Pk)
 

GHA_triaxial/approx_gha1.py

@@ -1,6 +1,6 @@
 import numpy as np
 from ellipsoide import EllipsoidTriaxial
-from .panou import louville_constant, func_sigma_ell, gha1_ana
+from GHA_triaxial.panou import louville_constant, func_sigma_ell, gha1_ana
 import plotly.graph_objects as go
 import winkelumrechnungen as wu
 
@@ -19,24 +19,31 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
     points = [p0]
     alphas = [alpha0]
     s_curr = 0.0
+
     while s_curr < s:
         ds_step = min(ds, s - s_curr)
         if ds_step < 1e-8:
             break
+
         p1 = points[-1]
         alpha1 = alphas[-1]
+
         sigma = func_sigma_ell(ell, p1, alpha1)
         p2 = p1 + ds_step * sigma
         p2 = ell.point_onto_ellipsoid(p2)
-        ds_step = np.linalg.norm(p2 - p1)
 
-        points.append(p2)
         dalpha = 1e-6
         l2 = louville_constant(ell, p2, alpha1)
         dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
         alpha2 = alpha1 + (l0 - l2) / dl_dalpha
+
+        points.append(p2)
         alphas.append(alpha2)
+
+        ds_step = np.linalg.norm(p2 - p1)
         s_curr += ds_step
+        if s_curr > 10000000:
+            pass
 
     if all_points:
         return points[-1], alphas[-1], np.array(points)
@@ -71,10 +78,10 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray):
 
 if __name__ == '__main__':
     ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    P0 = ell.para2cart(0, 0)
+    P0 = ell.para2cart(0.2, 0.3)
-    alpha0 = wu.deg2rad(90)
+    alpha0 = wu.deg2rad(35)
-    s = 1000000
+    s = 13000000
-    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=32)
+    P1_app, alpha1_app, points = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
-    P1_app, alpha1_app, points = gha1_approx(ell, P0, alpha0, s, ds=5000, all_points=True)
+    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16)
     show_points(points, P0, P1_ana)
     print(np.linalg.norm(P1_app - P1_ana))

GHA_triaxial/approx_gha1_2.py

@@ -0,0 +1,120 @@
+import numpy as np
+from ellipsoide import EllipsoidTriaxial
+from panou import louville_constant, func_sigma_ell, gha1_ana
+import plotly.graph_objects as go
+import winkelumrechnungen as wu
+from numpy import sin, cos, arccos
+
+def Bogenlaenge(P1: NDArray, P2: NDArray) -> float:
+    """
+    Berechnung der mittleren Bogenlänge zwischen zwei kartesischen Punkten
+    :param P1: kartesische Koordinate Punkt 1
+    :param P2: kartesische Koordinate Punkt 2
+    :return: Bogenlänge s
+    """
+    R1 = np.linalg.norm(P1)
+    R2 = np.linalg.norm(P2)
+    R = 0.5 * (R1 + R2)
+    if P1 @ P2 / (R1 * R2) > 1:
+        s = np.linalg.norm(P1 - P2)
+    else:
+        theta = arccos(P1 @ P2 / (R1 * R2))
+        s = float(R * theta)
+    return s
+
+def gha1_approx2(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]:
+    """
+    Berechung einer Näherungslösung der ersten Hauptaufgabe
+    :param ell: Ellipsoid
+    :param p0: Anfangspunkt
+    :param alpha0: Azimut im Anfangspunkt
+    :param s: Strecke bis zum Endpunkt
+    :param ds: Länge einzelner Streckenelemente
+    :param all_points: Ausgabe aller Punkte als Array?
+    :return: Endpunkt, Azimut im Endpunkt, optional alle Punkte
+    """
+    l0 = louville_constant(ell, p0, alpha0)
+    points = [p0]
+    alphas = [alpha0]
+    s_curr = 0.0
+
+    while s_curr < s:
+        ds_target = min(ds, s - s_curr)
+        if ds_target < 1e-8:
+            break
+
+        p1 = points[-1]
+        alpha1 = alphas[-1]
+        alpha1_mid = alphas[-1]
+        p2 = points[-1]
+        alpha2 = alphas[-1]
+
+        i = 0
+        while i < 2:
+            i += 1
+
+            sigma = func_sigma_ell(ell, p1, alpha1_mid)
+            p2_new = p1 + ds_target * sigma
+            p2_new = ell.point_onto_ellipsoid(p2_new)
+            p2 = p2_new
+
+            j = 0
+            while j < 2:
+                j += 1
+
+                dalpha = 1e-6
+                l2 = louville_constant(ell, p2, alpha2)
+                dl_dalpha = (louville_constant(ell, p2, alpha2 + dalpha) - l2) / dalpha
+                alpha2_new = alpha2 + (l0 - l2) / dl_dalpha
+                alpha2 = alpha2_new
+
+            alpha1_mid = (alpha1 + alpha2) / 2
+
+        points.append(p2)
+        alphas.append(alpha2)
+
+        ds_actual = np.linalg.norm(p2 - p1)
+        s_curr += ds_actual
+        if s_curr > 10000000:
+            pass
+
+    if all_points:
+        return points[-1], alphas[-1], np.array(points)
+    else:
+        return points[-1], alphas[-1]
+
+def show_points(points: NDArray, p0: NDArray, p1: NDArray):
+    """
+    Anzeigen der Punkte
+    :param points: Array aller approximierten Punkte
+    :param p0: Startpunkt
+    :param p1: wahrer Endpunkt
+    """
+    fig = go.Figure()
+
+    fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2],
+                      mode='lines', line=dict(color="red", width=3), name="Approx")
+    fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
+                      mode='markers', marker=dict(color="green"), name="P0")
+    fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
+                      mode='markers', marker=dict(color="green"), name="P1")
+
+    fig.update_layout(
+        scene=dict(xaxis_title='X [km]',
+                   yaxis_title='Y [km]',
+                   zaxis_title='Z [km]',
+                   aspectmode='data'),
+        title="CHAMP")
+
+    fig.show()
+
+
+if __name__ == '__main__':
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    P0 = ell.para2cart(0.2, 0.3)
+    alpha0 = wu.deg2rad(35)
+    s = 13000000
+    P1_app, alpha1_app, points = gha1_approx2(ell, P0, alpha0, s, ds=10000, all_points=True)
+    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16)
+    show_points(points, P0, P1_ana)
+    print(np.linalg.norm(P1_app - P1_ana))

GHA_triaxial/approx_gha2.py

@@ -9,7 +9,7 @@ from typing import Tuple
 from utils import sigma2alpha
 
 
-def gha2(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]:
+def gha2_approx(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]:
     """
     Numerische Approximation für die zweite Hauptaufgabe
     :param ell: Ellipsoid
@@ -83,15 +83,15 @@ def show_points(points: NDArray, points_app: NDArray, p0: NDArray, p1: NDArray):
 
 
 if __name__ == '__main__':
-    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
 
     beta0, lamb0 = (0.2, 0.1)
     P0 = ell.ell2cart(beta0, lamb0)
     beta1, lamb1 = (0.7, 0.3)
     P1 = ell.ell2cart(beta1, lamb1)
 
-    alpha0_app, alpha1_app, s_app, points = gha2(ell, P0, P1, ds=1e-4, all_points=True)
+    alpha0_app, alpha1_app, s_app, points = gha2_approx(ell, P0, P1, ds=1000, all_points=True)
-
+    print("done")
     alpha0, alpha1, s, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=5000, all_points=True)
     points_ana = []
     for beta, lamb in zip(betas, lambs):

GHA_triaxial/numeric_examples_karney.py

@@ -26,7 +26,7 @@ def get_random_examples(num: int, seed: int = None) -> List:
     """
     if seed is not None:
         random.seed(seed)
-    with open("Karney_2024_Testset.txt") as datei:
+    with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
         lines = datei.readlines()
     examples = []
     for i in range(num):

GHA_triaxial/panou.py

@@ -245,7 +245,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
     return p1, alpha1
 
 
-def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]:
+def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 16) -> Tuple[NDArray, float]:
     if s > np.pi / maxPartCircum * ell.ax:
         s /= 2
         point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum)
@@ -286,7 +286,7 @@ def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell:
     alpha_para = arctan2(p_para @ sigma_ell, q_para @ sigma_ell)
     sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
 
-    if np.linalg.norm(sigma_para - sigma_ell) > 1e-12:
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-9:
         raise Exception("Alpha Umrechnung fehlgeschlagen")
 
     return u, v, alpha_para
@@ -335,7 +335,7 @@ if __name__ == "__main__":
     # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
     # print(diffs_panou)
 
-    ell = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
+    ell: EllipsoidTriaxial = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
     diffs_karney = []
     # examples_karney = ne_karney.get_examples((30499, 30500, 40500))
     examples_karney = ne_karney.get_random_examples(20)
@@ -343,12 +343,12 @@ if __name__ == "__main__":
         beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
         P0 = ell.ell2cart(beta0, lamb0)
 
-        P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 5000)
+        P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 10000)
         beta1_num, lamb1_num = ell.cart2ell(P1_num)
 
         try:
             _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
-            P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 30, maxPartCircum=16)
+            P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 45, maxPartCircum=32)
             beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
         except:
             beta1_ana, lamb1_ana = np.inf, np.inf

algorithms_test.py

@@ -0,0 +1,267 @@
+import time
+import pickle
+import numpy as np
+from numpy import nan
+import winkelumrechnungen as wu
+import os
+from contextlib import contextmanager, redirect_stdout, redirect_stderr
+
+from ellipsoide import EllipsoidTriaxial
+from GHA_triaxial.panou import alpha_ell2para, alpha_para2ell
+
+from GHA_triaxial.panou import gha1_num, gha1_ana
+from GHA_triaxial.approx_gha1 import gha1_approx
+
+from GHA_triaxial.panou_2013_2GHA_num import gha2_num
+from GHA_triaxial.ES_gha2 import gha2_ES
+from GHA_triaxial.approx_gha2 import gha2_approx
+
+from GHA_triaxial.numeric_examples_panou import get_random_examples as get_examples_panou
+from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney
+
+@contextmanager
+def suppress_print():
+    with open(os.devnull, 'w') as fnull:
+        with redirect_stdout(fnull), redirect_stderr(fnull):
+            yield
+
+# steps_gha1_num = [2000, 5000, 10000, 20000]
+# maxM_gha1_ana = [20, 40, 60]
+# parts_gha1_ana = [4, 8, 16]
+# dsPart_gha1_approx = [600, 1250, 6000, 60000]   # entspricht bei der Erde ca. 10000, 5000, 1000, 100
+#
+# steps_gha2_num = [2000, 5000, 10000, 20000]
+# dsPart_gha2_ES = [600, 1250, 6000]   # entspricht bei der Erde ca. 10000, 5000, 1000
+# dsPart_gha2_approx = [600, 1250, 6000, 60000]   # entspricht bei der Erde ca. 10000, 5000, 1000, 100
+
+steps_gha1_num = [2000, 5000]
+maxM_gha1_ana = [20, 40]
+parts_gha1_ana = [4, 8]
+dsPart_gha1_approx = [600, 1250]
+
+steps_gha2_num = [2000, 5000]
+dsPart_gha2_ES = [20]
+dsPart_gha2_approx = [600, 1250]
+
+ell_karney: EllipsoidTriaxial = EllipsoidTriaxial.init_name("KarneyTest2024")
+ell_panou: EllipsoidTriaxial = EllipsoidTriaxial.init_name("BursaSima1980round")
+
+results_karney = {}
+results_panou = {}
+
+examples_karney = get_examples_karney(2, 42)
+examples_panou = get_examples_panou(2, 42)
+
+for example in examples_karney:
+    example_results = {}
+
+    beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
+    P0 = ell_karney.ell2cart(beta0, lamb0)
+    P1 = ell_karney.ell2cart(beta1, lamb1)
+    _, _, alpha0_para = alpha_ell2para(ell_karney, beta0, lamb0, alpha0_ell)
+
+    for steps in steps_gha1_num:
+        start = time.perf_counter()
+        try:
+            P1_num, alpha1_num_1 = gha1_num(ell_karney, P0, alpha0_ell, s, num=steps)
+            end = time.perf_counter()
+            beta1_num, lamb1_num = ell_karney.cart2ell(P1_num)
+            d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600
+            d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600
+            d_time = end - start
+            example_results[f"GHA1_num_{steps}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA1_num_{steps}"] = (nan, nan, nan, nan)
+
+    for maxM in maxM_gha1_ana:
+        for parts in parts_gha1_ana:
+            start = time.perf_counter()
+            try:
+                P1_ana, alpha1_ana_para = gha1_ana(ell_karney, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)
+                end = time.perf_counter()
+                beta1_ana, lamb1_ana = ell_karney.cart2ell(P1_ana)
+                _, _, alpha1_ana_ell = alpha_para2ell(ell_karney, beta1_ana, lamb1_ana, alpha1_ana_para)
+                d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600
+                d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600
+                d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600
+                d_time = end - start
+                example_results[f"GHA1_ana_{maxM}_{parts}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+            except Exception as e:
+                print(e)
+                example_results[f"GHA1_ana_{maxM}_{parts}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha1_approx:
+        ds = ell_karney.ax/dsPart
+        start = time.perf_counter()
+        try:
+            P1_approx, alpha1_approx = gha1_approx(ell_karney, P0, alpha0_ell, s, ds=ds)
+            end = time.perf_counter()
+            beta1_approx, lamb1_approx = ell_karney.cart2ell(P1_approx)
+            d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600
+            d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600
+            d_time = end - start
+            example_results[f"GHA1_approx_{ds:.3f}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA1_approx_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    for steps in steps_gha2_num:
+        start = time.perf_counter()
+        try:
+            alpha0_num, alpha1_num_2, s_num = gha2_num(ell_karney, beta0, lamb0, beta1, lamb1, n=steps)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600
+            d_s = abs(s_num - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_num_{steps}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_num_{steps}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha2_ES:
+        ds = ell_karney.ax/dsPart
+        start = time.perf_counter()
+        try:
+            with suppress_print():
+                alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell_karney, P0, P1, stepLenTarget=ds)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600
+            d_s = abs(s_ES - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_ES_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_ES_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha2_approx:
+        ds = ell_karney.ax/dsPart
+        start = time.perf_counter()
+        try:
+            alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell_karney, P0, P1, ds=ds)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600
+            d_s = abs(s_approx - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_approx_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_approx_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    results_karney[f"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}"] = example_results
+
+for example in examples_panou:
+    example_results = {}
+
+    beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example
+    P0 = ell_panou.ell2cart(beta0, lamb0)
+    P1 = ell_panou.ell2cart(beta1, lamb1)
+    _, _, alpha0_para = alpha_ell2para(ell_panou, beta0, lamb0, alpha0_ell)
+
+    for steps in steps_gha1_num:
+        start = time.perf_counter()
+        try:
+            P1_num, alpha1_num_1 = gha1_num(ell_panou, P0, alpha0_ell, s, num=steps)
+            end = time.perf_counter()
+            beta1_num, lamb1_num = ell_panou.cart2ell(P1_num)
+            d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600
+            d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600
+            d_time = end - start
+            example_results[f"GHA1_num_{steps}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA1_num_{steps}"] = (nan, nan, nan, nan)
+
+    for maxM in maxM_gha1_ana:
+        for parts in parts_gha1_ana:
+            start = time.perf_counter()
+            try:
+                P1_ana, alpha1_ana_para = gha1_ana(ell_panou, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)
+                end = time.perf_counter()
+                beta1_ana, lamb1_ana = ell_panou.cart2ell(P1_ana)
+                _, _, alpha1_ana = alpha_para2ell(ell_panou, beta1_ana, lamb1_ana, alpha1_ana_para)
+                d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600
+                d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600
+                d_alpha1 = abs(wu.rad2deg(alpha1_ana - alpha1_ell)) / 3600
+                d_time = end - start
+                example_results[f"GHA1_ana_{maxM}_{parts}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+            except Exception as e:
+                print(e)
+                example_results[f"GHA1_ana_{maxM}_{parts}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha1_approx:
+        ds = ell_panou.ax/dsPart
+        start = time.perf_counter()
+        try:
+            P1_approx, alpha1_approx = gha1_approx(ell_panou, P0, alpha0_ell, s, ds=ds)
+            end = time.perf_counter()
+            beta1_approx, lamb1_approx = ell_panou.cart2ell(P1_approx)
+            d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600
+            d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600
+            d_time = end - start
+            example_results[f"GHA1_approx_{ds:.3f}"] = (d_beta1, d_lamb1, d_alpha1, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA1_approx_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    for steps in steps_gha2_num:
+        start = time.perf_counter()
+        try:
+            alpha0_num, alpha1_num_2, s_num = gha2_num(ell_panou, beta0, lamb0, beta1, lamb1, n=steps)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600
+            d_s = abs(s_num - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_num_{steps}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_num_{steps}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha2_ES:
+        ds = ell_panou.ax/dsPart
+        start = time.perf_counter()
+        try:
+            with suppress_print():
+                alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell_panou, P0, P1, stepLenTarget=ds)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600
+            d_s = abs(s_ES - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_ES_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_ES_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    for dsPart in dsPart_gha2_approx:
+        ds = ell_panou.ax/dsPart
+        start = time.perf_counter()
+        try:
+            alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell_panou, P0, P1, ds=ds)
+            end = time.perf_counter()
+            d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600
+            d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600
+            d_s = abs(s_approx - s) / 1000
+            d_time = end - start
+            example_results[f"GHA2_approx_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time)
+        except Exception as e:
+            print(e)
+            example_results[f"GHA2_approx_{ds:.3f}"] = (nan, nan, nan, nan)
+
+    results_panou[f"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}"] = example_results
+
+print(results_karney)
+with open("results_karney.pkl", "wb") as f:
+    pickle.dump(results_karney, f)
+
+print(results_panou)
+with open("results_panou.pkl", "wb") as f:
+    pickle.dump(results_panou, f)

ellipsoide.py

@@ -327,14 +327,45 @@ class EllipsoidTriaxial:
         x, y, z = point
         beta, lamb = self.cart2ell_panou(point)
         delta_ell = np.array([np.inf, np.inf]).T
+        tiny = 1e-30
 
         i = 0
-        while np.sum(delta_ell) > eps and i < maxI:
+        while np.linalg.norm(delta_ell) > eps and i < maxI:
+            if abs(y) < eps:
+                delta_y = 1e-4
+                best_delta = np.inf
+                while True:
+                    try:
+                        y1 = y - delta_y
+                        beta1, lamb1 = self.cart2ell(np.array([x, y1, z]))
+                        point1 = self.ell2cart(beta1, lamb1)
+
+                        y2 = y + delta_y
+                        beta2, lamb2 = self.cart2ell(np.array([x, y2, z]))
+                        point2 = self.ell2cart(beta2, lamb2)
+
+                        pointM = (point1 + point2) / 2
+
+                        actual_delta = np.linalg.norm(point-pointM)
+                    except:
+                        actual_delta = np.inf
+
+                    if actual_delta < best_delta:
+                        best_delta = actual_delta
+                        delta_y /= 10
+                    else:
+                        delta_y *= 10
+
+                        y1 = y - delta_y
+                        beta1, lamb1 = self.cart2ell(np.array([x, y1, z]))
+
+                        return beta1, lamb1
+
             x0, y0, z0 = self.ell2cart(beta, lamb)
             delta_l = np.array([x-x0, y-y0, z-z0]).T
 
-            B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2
+            B = max(self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2, tiny)
-            L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2
+            L = max(self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2, tiny)
 
             J = np.array([[(-self.ax * self.Ey ** 2) / (2 * self.Ex) * sin(2 * beta) / sqrt(B) * cos(lamb),
                            -self.ax / self.Ex * sqrt(B) * sin(lamb)],
@@ -345,23 +376,21 @@ class EllipsoidTriaxial:
 
             N = J.T @ J
             det = N[0, 0] * N[1, 1] - N[0, 1] * N[1, 0]
-            if abs(det) < eps:
-                det = eps
             N_inv = 1 / det * np.array([[N[1, 1], -N[0, 1]], [-N[1, 0], N[0, 0]]])
             delta_ell = N_inv @ J.T @ delta_l
+
             beta += delta_ell[0]
             lamb += delta_ell[1]
             i += 1
 
         if i == maxI:
-            raise Exception("Umrechung ist nicht konvergiert")
+            raise Exception("Umrechnung ist nicht konvergiert")
 
         point_n = self.ell2cart(beta, lamb)
         delta_r = np.linalg.norm(point - point_n, axis=-1)
 
-        if delta_r > 1e-4:
+        if delta_r > 1e-3:
-            # raise Exception("Fehler in der Umrechnung cart2ell")
+            raise Exception("Fehler in der Umrechnung cart2ell")
-            print(f"Fehler in der Umrechnung cart2ell, deltaR = {delta_r}m")
 
         return beta, lamb
 
@@ -395,9 +424,9 @@ class EllipsoidTriaxial:
         t1, t2 = self.func_t12(point)
 
         num_beta = max(t1 - self.b ** 2, 0)
-        den_beta = max(self.ay ** 2 - t1, 0)
+        den_beta = max(self.ay ** 2 - t1, 1e-30)
         num_lamb = max(t2 - self.ay ** 2, 0)
-        den_lamb = max(self.ax ** 2 - t2, 0)
+        den_lamb = max(self.ax ** 2 - t2, 1e-30)
 
         beta = arctan(sqrt(num_beta / den_beta))
         lamb = arctan(sqrt(num_lamb / den_lamb))
@@ -675,7 +704,7 @@ class EllipsoidTriaxial:
 
 
 if __name__ == "__main__":
-    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    ell = EllipsoidTriaxial.init_name("BursaSima1980")
     diff_list = []
     diffs_para = []
     diffs_ell = []
@@ -711,6 +740,7 @@ if __name__ == "__main__":
     diffs_geod = np.array(diffs_geod)
 
     pass
+
     points = np.array(points)
     fig = plt.figure()
     ax = fig.add_subplot(projection='3d')