From 07212dcc97b9de6f8d299e97b2eaa41e0cf5b5b2 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sat, 17 Jan 2026 18:51:47 +0100 Subject: [PATCH] Algorithmen Test --- GHA_triaxial/ES_gha2.py | 6 +- GHA_triaxial/approx_gha1.py | 23 +- GHA_triaxial/approx_gha1_2.py | 120 +++++++++++ GHA_triaxial/approx_gha2.py | 8 +- GHA_triaxial/numeric_examples_karney.py | 2 +- GHA_triaxial/panou.py | 10 +- algorithms_test.py | 267 ++++++++++++++++++++++++ ellipsoide.py | 54 +++-- 8 files changed, 457 insertions(+), 33 deletions(-) create mode 100644 GHA_triaxial/approx_gha1_2.py create mode 100644 algorithms_test.py diff --git a/GHA_triaxial/ES_gha2.py b/GHA_triaxial/ES_gha2.py index 0c851e3..90c9b7c 100644 --- a/GHA_triaxial/ES_gha2.py +++ b/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) diff --git a/GHA_triaxial/approx_gha1.py b/GHA_triaxial/approx_gha1.py index 28fae4c..c571fa5 100644 --- a/GHA_triaxial/approx_gha1.py +++ b/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) - alpha0 = wu.deg2rad(90) - s = 1000000 - 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=5000, all_points=True) + P0 = ell.para2cart(0.2, 0.3) + alpha0 = wu.deg2rad(35) + s = 13000000 + P1_app, alpha1_app, points = gha1_approx(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)) diff --git a/GHA_triaxial/approx_gha1_2.py b/GHA_triaxial/approx_gha1_2.py new file mode 100644 index 0000000..2070f1f --- /dev/null +++ b/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)) diff --git a/GHA_triaxial/approx_gha2.py b/GHA_triaxial/approx_gha2.py index 70879a1..d8b969c 100644 --- a/GHA_triaxial/approx_gha2.py +++ b/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): diff --git a/GHA_triaxial/numeric_examples_karney.py b/GHA_triaxial/numeric_examples_karney.py index 173683d..e17ed57 100644 --- a/GHA_triaxial/numeric_examples_karney.py +++ b/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): diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 433b780..d4c6ac4 100644 --- a/GHA_triaxial/panou.py +++ b/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 diff --git a/algorithms_test.py b/algorithms_test.py new file mode 100644 index 0000000..7320bd8 --- /dev/null +++ b/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) \ No newline at end of file diff --git a/ellipsoide.py b/ellipsoide.py index 96b59e7..1c07198 100644 --- a/ellipsoide.py +++ b/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 - L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2 + B = max(self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2, tiny) + 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: - # raise Exception("Fehler in der Umrechnung cart2ell") - print(f"Fehler in der Umrechnung cart2ell, deltaR = {delta_r}m") + if delta_r > 1e-3: + raise Exception("Fehler in der Umrechnung cart2ell") 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')