From 4d5b6fcc3ef980cc1c4fb1060cc7e420daa44ff1 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sun, 11 Jan 2026 13:14:43 +0100 Subject: [PATCH] =?UTF-8?q?Umrechnung=20alpha,=20N=C3=A4herungsl=C3=B6sung?= =?UTF-8?q?=20GHA=201?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA_triaxial/approx_gha1.py | 63 ++++++++++++ GHA_triaxial/numeric_examples_karney.py | 2 +- GHA_triaxial/panou.py | 125 ++++++++++++++++++------ GHA_triaxial/panou_2013_2GHA_num.py | 62 ++++++------ dashboard.py | 4 +- ellipsoide.py | 20 ++-- 6 files changed, 205 insertions(+), 71 deletions(-) create mode 100644 GHA_triaxial/approx_gha1.py diff --git a/GHA_triaxial/approx_gha1.py b/GHA_triaxial/approx_gha1.py new file mode 100644 index 0000000..67fdd2c --- /dev/null +++ b/GHA_triaxial/approx_gha1.py @@ -0,0 +1,63 @@ +import numpy as np +from numpy import sin, cos, arcsin, arccos, arctan2 +from ellipsoide import EllipsoidTriaxial +import matplotlib.pyplot as plt +from panou import louville_constant, func_sigma_ell, gha1_ana +import plotly.graph_objects as go +import winkelumrechnungen as wu + +def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: int): + l0 = louville_constant(ell, p0, alpha0) + 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] + x1, y1, z1 = p1 + sigma = func_sigma_ell(ell, x1, y1, z1, alpha1) + p2 = p1 + ds_step * sigma + p2, _, _, _ = ell.cartonell(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 + alphas.append(alpha2) + s_curr += ds_step + return points[-1], alphas[-1], np.array(points) + +def show_points(points, p0, p1): + 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, 0) + alpha0 = wu.deg2rad(90) + s = 1000000 + P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 60, maxPartCircum=32) + P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, 5000) + show_points(points, P0, P1_ana) + print(np.linalg.norm(P1_app - P1_ana)) diff --git a/GHA_triaxial/numeric_examples_karney.py b/GHA_triaxial/numeric_examples_karney.py index 7eecfa5..5303b10 100644 --- a/GHA_triaxial/numeric_examples_karney.py +++ b/GHA_triaxial/numeric_examples_karney.py @@ -16,7 +16,7 @@ def get_random_examples(num): :param num: :return: """ - random.seed(42) + # random.seed(42) with open("Karney_2024_Testset.txt") as datei: lines = datei.readlines() examples = [] diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 5592366..3579763 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -68,7 +68,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable: return np.array([dxds, ddx, dyds, ddy, dzds, ddz]) return ODE -def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float]: +def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float, list]: """ Panou, Korakitits 2019 :param ell: @@ -95,15 +95,15 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, p1, q1 = pq_ell(ell, x1, y1, z1) sigma = np.array([dx1ds, dy1ds, dz1ds]) - P = p1 @ sigma - Q = q1 @ sigma + P = float(p1 @ sigma) + Q = float(q1 @ sigma) alpha1 = arctan2(P, Q) if alpha1 < 0: alpha1 += 2 * np.pi - return np.array([x1, y1, z1]), alpha1 + return np.array([x1, y1, z1]), alpha1, werte # --------------------------------------------------------------------------------------------------------------------- @@ -225,8 +225,8 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): # 52-53 sigma = np.array([dx_s, dy_s, dz_s]) - P = p_s @ sigma - Q = q_s @ sigma + P = float(p_s @ sigma) + Q = float(q_s @ sigma) # 51 alpha1 = arctan2(P, Q) @@ -237,7 +237,7 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): return np.array([x_s, y_s, z_s]), alpha1 -def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4): +def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> tuple[np.ndarray, float]: if s > np.pi / maxPartCircum * ell.ax: s /= 2 point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum) @@ -251,39 +251,103 @@ def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, return point_end, alpha_end +def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> tuple[float, float, float]: + x, y, z = ell.para2cart(u, v) + beta, lamb = ell.para2ell(u, v) + + p_para, q_para = pq_para(ell, x, y, z) + sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para) + + p_ell, q_ell = pq_ell(ell, x, y, z) + alpha_ell = arctan2(p_ell @ sigma_para, q_ell @ sigma_para) + sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) + + if np.linalg.norm(sigma_para - sigma_ell) > 1e-12: + raise Exception("Alpha Umrechnung fehlgeschlagen") + + return beta, lamb, alpha_ell + +def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> tuple[float, float, float]: + x, y, z = ell.ell2cart(beta, lamb) + u, v = ell.ell2para(beta, lamb) + + p_ell, q_ell = pq_ell(ell, x, y, z) + sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) + + p_para, q_para = pq_para(ell, x, y, z) + 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: + raise Exception("Alpha Umrechnung fehlgeschlagen") + print("Alpha Umrechnung fehlgeschlagen:", np.linalg.norm(sigma_para - sigma_ell)) + + return u, v, alpha_para + +def func_sigma_ell(ell, x, y, z, alpha): + p, q = pq_ell(ell, x, y, z) + sigma = p * sin(alpha) + q * cos(alpha) + return sigma + +def func_sigma_para(ell, x, y, z, alpha): + p, q = pq_para(ell, x, y, z) + sigma = p * sin(alpha) + q * cos(alpha) + return sigma + + +def louville_constant(ell: EllipsoidTriaxial, p0: np.ndarray, alpha: float) -> float: + beta, lamb = ell.cart2ell(p0) + l = ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 - ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2 + # x, y, z = p0 + # t1, t2 = ell.func_t12(x, y, z) + # l_cart = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2) + # if abs(l - l_cart) > 1e-12: + # # raise Exception("Louville constant fehlgeschlagen") + # print("Diff zwischen constant:", abs(l - l_cart)) + return l + +def louville_l2c(ell, l): + return sqrt((l + ell.Ee**2) / ell.Ex**2) + +def louville_c2l(ell, c): + return ell.Ex**2 * c**2 - ell.Ee**2 + if __name__ == "__main__": - ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") - diffs_panou = [] - examples_panou = ne_panou.get_random_examples(5) - for example in examples_panou: - beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example - P0 = ell.ell2cart(beta0, lamb0) - - P1_num, alpha1_num = gha1_num(ell, P0, alpha0, s, 100) - beta1_num, lamb1_num = ell.cart2ell(P1_num) - P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 60) - beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) - diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana))) - diffs_panou = np.array(diffs_panou) - mask_360 = (diffs_panou > 359) & (diffs_panou < 361) - diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360) - print(diffs_panou) + # ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") + # diffs_panou = [] + # examples_panou = ne_panou.get_random_examples(5) + # for example in examples_panou: + # beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example + # P0 = ell.ell2cart(beta0, lamb0) + # + # P1_num, alpha1_num, _ = gha1_num(ell, P0, alpha0_ell, s, 100) + # beta1_num, lamb1_num = ell.cart2ell(P1_num) + # + # _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell) + # P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 60) + # beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) + # diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana))) + # diffs_panou = np.array(diffs_panou) + # mask_360 = (diffs_panou > 359) & (diffs_panou < 361) + # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360) + # print(diffs_panou) ell = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024") - diffs_karney = [] - examples_karney = ne_karney.get_examples((30499, 30500, 40500)) - # examples_karney = ne_karney.get_random_examples(5) + # examples_karney = ne_karney.get_examples((30499, 30500, 40500)) + examples_karney = ne_karney.get_random_examples(20) for example in examples_karney: - beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example + beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example P0 = ell.ell2cart(beta0, lamb0) - P1_num, alpha1_num = gha1_num(ell, P0, alpha0, s, 100) + P1_num, alpha1_num, _ = gha1_num(ell, P0, alpha0_ell, s, 5000) beta1_num, lamb1_num = ell.cart2ell(P1_num) + try: - P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 40) + _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell) + P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 30, maxPartCircum=16) beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) except: beta1_ana, lamb1_ana = np.inf, np.inf @@ -293,4 +357,5 @@ if __name__ == "__main__": mask_360 = (diffs_karney > 359) & (diffs_karney < 361) diffs_karney[mask_360] = np.abs(diffs_karney[mask_360] - 360) print(diffs_karney) - pass \ No newline at end of file + + pass diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index 568167e..713c394 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -244,7 +244,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep + (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2 ) - return alpha_1, alpha_2, s + return alpha_1, alpha_2, s, beta_arr, lamb_arr if lamb_1 == lamb_2: @@ -252,7 +252,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep dbeta = beta_2 - beta_1 if abs(dbeta) < 10**-15: - return 0, 0, 0 + return 0, 0, 0, np.array([]), np.array([]) lamb_0 = 0 @@ -369,7 +369,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep else: s = np.trapz(integrand, dx=h) - return alpha_1, alpha_2, s + return alpha_1, alpha_2, s, beta_arr, lamb_arr if __name__ == "__main__": @@ -391,35 +391,37 @@ if __name__ == "__main__": # print(s) - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - diffs_panou = [] - examples_panou = ne_panou.get_random_examples(4) - for example in examples_panou: - beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example - P0 = ell.ell2cart(beta0, lamb0) - try: - alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10) - diffs_panou.append( - (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num))) - except: - print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}") - diffs_panou = np.array(diffs_panou) - print(diffs_panou) + # ell = EllipsoidTriaxial.init_name("BursaSima1980round") + # diffs_panou = [] + # examples_panou = ne_panou.get_random_examples(4) + # for example in examples_panou: + # beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example + # P0 = ell.ell2cart(beta0, lamb0) + # try: + # alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10) + # diffs_panou.append( + # (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num))) + # except: + # print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}") + # diffs_panou = np.array(diffs_panou) + # print(diffs_panou) + # + # ell = EllipsoidTriaxial.init_name("KarneyTest2024") + # diffs_karney = [] + # # examples_karney = ne_karney.get_examples((30500, 40500)) + # examples_karney = ne_karney.get_random_examples(2) + # for example in examples_karney: + # beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example + # + # try: + # alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10) + # diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num))) + # except: + # print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}") + # diffs_karney = np.array(diffs_karney) + # print(diffs_karney) - ell = EllipsoidTriaxial.init_name("KarneyTest2024") - diffs_karney = [] - # examples_karney = ne_karney.get_examples((30500, 40500)) - examples_karney = ne_karney.get_random_examples(2) - for example in examples_karney: - beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example - try: - alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10) - diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num))) - except: - print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}") - diffs_karney = np.array(diffs_karney) - print(diffs_karney) pass diff --git a/dashboard.py b/dashboard.py index e8841a9..bed22e7 100644 --- a/dashboard.py +++ b/dashboard.py @@ -413,7 +413,7 @@ def calc_and_plot(n1, n2, if "analytisch" in method1: # ana - x2, y2, z2 = gha1_ana(ell, p1, alpha_rad, s_val, 70) + (x2, y2, z2), alpha2 = gha1_ana(ell, p1, alpha_rad, s_val, 70) p2_ana = (float(x2), float(y2), float(z2)) beta2, lamb2 = ell.cart2ell([x2, y2, z2]) @@ -430,7 +430,7 @@ def calc_and_plot(n1, n2, if "numerisch" in method1: # num - x1, y1, z1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000) + (x1, y1, z1), alpha1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000) p2_num = x1, y1, z1 beta2_num, lamb2_num = ell.cart2ell(p2_num) diff --git a/ellipsoide.py b/ellipsoide.py index 6cb82bb..6a79e92 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -173,6 +173,17 @@ class EllipsoidTriaxial: y / ((1 - self.ee ** 2) * sqrtH), z / ((1 - self.ex ** 2) * sqrtH)]) + def func_t12(self, x, y, z): + c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2) + c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 - + (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - ( + self.ax ** 2 + self.ay ** 2) * z ** 2) + t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2 + if t2 == 0: + t2 = 1e-18 + t1 = c0 / t2 + return t1, t2 + def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: """ Panou 2014 12ff. @@ -360,14 +371,7 @@ class EllipsoidTriaxial: # ---- Allgemeiner Fall ----- - c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2) - c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 - - (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - ( - self.ax ** 2 + self.ay ** 2) * z ** 2) - t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2 - if t2 == 0: - t2 = 1e-14 - t1 = c0 / t2 + t1, t2 = self.func_t12(x, y, z) num_beta = max(t1 - self.b ** 2, 0) den_beta = max(self.ay ** 2 - t1, 0)