From 797afdfd6fc8c6b6a3c841ad99fb5e8b4a214217 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Fri, 9 Jan 2026 17:49:49 +0100 Subject: [PATCH 1/3] =?UTF-8?q?Projekt=20aufger=C3=A4umt,=20gha1=20geteste?= =?UTF-8?q?t,=20Runge-Kutta=20angepasst=20(gha2=5Fnum=20sollte=20jetzt=20d?= =?UTF-8?q?eutlich=20schneller=20sein)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA/rk.py | 17 - {GHA => GHA_biaxial}/__init__.py | 0 {GHA => GHA_biaxial}/bessel.py | 0 {GHA => GHA_biaxial}/gauss.py | 0 GHA_biaxial/rk.py | 21 + GHA_triaxial/numeric_examples_karney.py | 43 ++ GHA_triaxial/numeric_examples_panou.py | 18 +- GHA_triaxial/panou.py | 300 ++++++--- GHA_triaxial/panou_2013_2GHA_num.py | 211 ++++-- Numerische_Integration/num_int_runge_kutta.py | 84 --- .../num_int_trapezformel.py | 46 -- Numerische_Integration/rk_DGLS_1Ordnung.py | 7 - Numerische_Integration/rk_DGL_1Ordnung.py | 6 - Numerische_Integration/rk_DGL_2Ordnung.py | 10 - Numerische_Integration/rk_DGL_2Ordnung_2.py | 7 - dashboard.py | 1 + ellipsoide.py | 608 ++++++++++-------- phi_ellipse.py | 55 -- runge_kutta.py | 43 ++ s_ellipse.py | 78 --- show_constant_lines.py | 30 - test.py | 31 +- test_algorithms.py | 84 --- 23 files changed, 832 insertions(+), 868 deletions(-) delete mode 100644 GHA/rk.py rename {GHA => GHA_biaxial}/__init__.py (100%) rename {GHA => GHA_biaxial}/bessel.py (100%) rename {GHA => GHA_biaxial}/gauss.py (100%) create mode 100644 GHA_biaxial/rk.py create mode 100644 GHA_triaxial/numeric_examples_karney.py delete mode 100644 Numerische_Integration/num_int_runge_kutta.py delete mode 100644 Numerische_Integration/num_int_trapezformel.py delete mode 100644 Numerische_Integration/rk_DGLS_1Ordnung.py delete mode 100644 Numerische_Integration/rk_DGL_1Ordnung.py delete mode 100644 Numerische_Integration/rk_DGL_2Ordnung.py delete mode 100644 Numerische_Integration/rk_DGL_2Ordnung_2.py delete mode 100644 phi_ellipse.py create mode 100644 runge_kutta.py delete mode 100644 s_ellipse.py delete mode 100644 show_constant_lines.py delete mode 100644 test_algorithms.py diff --git a/GHA/rk.py b/GHA/rk.py deleted file mode 100644 index 4c8913e..0000000 --- a/GHA/rk.py +++ /dev/null @@ -1,17 +0,0 @@ -import Numerische_Integration.num_int_runge_kutta as rk -from numpy import sin, cos, tan -import winkelumrechnungen as wu -from ellipsoide import EllipsoidBiaxial - -def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num): - phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0) - - f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c - f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c) - f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c - - funktionswerte = rk.verfahren([f_phi, f_lam, f_A], - [0, phi0, lamb0, A0], - s, num) - coords = re.ell2cart(funktionswerte[-1][1], funktionswerte[-1][2], h0) - return coords diff --git a/GHA/__init__.py b/GHA_biaxial/__init__.py similarity index 100% rename from GHA/__init__.py rename to GHA_biaxial/__init__.py diff --git a/GHA/bessel.py b/GHA_biaxial/bessel.py similarity index 100% rename from GHA/bessel.py rename to GHA_biaxial/bessel.py diff --git a/GHA/gauss.py b/GHA_biaxial/gauss.py similarity index 100% rename from GHA/gauss.py rename to GHA_biaxial/gauss.py diff --git a/GHA_biaxial/rk.py b/GHA_biaxial/rk.py new file mode 100644 index 0000000..85101c8 --- /dev/null +++ b/GHA_biaxial/rk.py @@ -0,0 +1,21 @@ +import runge_kutta as rk +from numpy import sin, cos, tan +import winkelumrechnungen as wu +from ellipsoide import EllipsoidBiaxial +import numpy as np + +def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num): + phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0) + + def buildODE(): + def ODE(s, v): + phi, lam, A = v + dphi = cos(A) * re.V(phi) ** 3 / re.c + dlam = sin(A) * re.V(phi) / (cos(phi) * re.c) + dA = tan(phi) * sin(A) * re.V(phi) / re.c + return np.array([dphi, dlam, dA]) + return ODE + + _, funktionswerte = rk.rk4(buildODE(), 0, np.array([phi0, lamb0, A0]), s, num) + coords = re.ell2cart(funktionswerte[-1][0], funktionswerte[-1][1], h0) + return coords diff --git a/GHA_triaxial/numeric_examples_karney.py b/GHA_triaxial/numeric_examples_karney.py new file mode 100644 index 0000000..7eecfa5 --- /dev/null +++ b/GHA_triaxial/numeric_examples_karney.py @@ -0,0 +1,43 @@ +import random +import winkelumrechnungen as wu + +def line2example(line): + split = line.split() + example = [float(value) for value in split[:7]] + for i, value in enumerate(example): + if i < 6: + example[i] = wu.deg2rad(value) + # example[i] = value + return example + +def get_random_examples(num): + """ + beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12 + :param num: + :return: + """ + random.seed(42) + with open("Karney_2024_Testset.txt") as datei: + lines = datei.readlines() + examples = [] + for i in range(num): + example = line2example(lines[random.randint(0, len(lines) - 1)]) + examples.append(example) + return examples + +def get_examples(l_i): + """ + beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12 + :param num: + :return: + """ + with open("Karney_2024_Testset.txt") as datei: + lines = datei.readlines() + examples = [] + for i in l_i: + example = line2example(lines[i]) + examples.append(example) + return examples + +if __name__ == "__main__": + get_random_examples(10) \ No newline at end of file diff --git a/GHA_triaxial/numeric_examples_panou.py b/GHA_triaxial/numeric_examples_panou.py index 393d03a..5f025cc 100644 --- a/GHA_triaxial/numeric_examples_panou.py +++ b/GHA_triaxial/numeric_examples_panou.py @@ -1,4 +1,5 @@ import winkelumrechnungen as wu +import random table1 = [ (wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(90), 1.00000000000, @@ -105,7 +106,22 @@ def get_example(table, example): def get_tables(): return tables +def get_random_examples(num): + random.seed(42) + examples = [] + for i in range(num): + table = random.randint(1, 4) + if table == 4: + example = random.randint(1, 8) + else: + example = random.randint(1, 7) + example = get_example(table, example) + examples.append(example) + return examples + + if __name__ == "__main__": - test = get_example(1, 4) + # test = get_example(1, 4) + examples = get_random_examples(5) pass diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 5b265ff..5592366 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -1,42 +1,61 @@ import numpy as np from numpy import sin, cos, sqrt, arctan2 import ellipsoide -import Numerische_Integration.num_int_runge_kutta as rk +import runge_kutta as rk import winkelumrechnungen as wu -import ausgaben as aus -import GHA.rk as ghark from scipy.special import factorial as fact from math import comb -import GHA_triaxial.numeric_examples_panou as nep - -# Panou, Korakitits 2019 +import GHA_triaxial.numeric_examples_panou as ne_panou +import GHA_triaxial.numeric_examples_karney as ne_karney +from ellipsoide import EllipsoidTriaxial +from typing import Callable -def gha1_num_old(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num): - phi, lamb, h = ell.cart2geod(point, "ligas3") - x, y, z = ell.geod2cart(phi, lamb, 0) - p, q = ell.p_q(x, y, z) +def pq_ell(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]: + """ + Berechnung von p und q in elliptischen Koordinaten + Panou, Korakitits 2019 + :param x: x + :param y: y + :param z: z + :return: p und q + """ + n = ell.func_n(x, y, z) - 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) + beta, lamb = ell.cart2ell(np.array([x, y, z])) + B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2 + L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2 - h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2 + c1 = x ** 2 + y ** 2 + z ** 2 - (ell.ax ** 2 + ell.ay ** 2 + ell.b ** 2) + c0 = (ell.ax ** 2 * ell.ay ** 2 + ell.ax ** 2 * ell.b ** 2 + ell.ay ** 2 * ell.b ** 2 - + (ell.ay ** 2 + ell.b ** 2) * x ** 2 - (ell.ax ** 2 + ell.b ** 2) * y ** 2 - ( + ell.ax ** 2 + ell.ay ** 2) * z ** 2) + t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2 - f1 = lambda x, dxds, y, dyds, z, dzds: dxds - f2 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * x - f3 = lambda x, dxds, y, dyds, z, dzds: dyds - f4 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * y/(1-ell.ee**2) - f5 = lambda x, dxds, y, dyds, z, dzds: dzds - f6 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * z/(1-ell.ex**2) + F = ell.Ey ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(lamb) ** 2 + p1 = -sqrt(L / (F * t2)) * ell.ax / ell.Ex * sqrt(B) * sin(lamb) + p2 = sqrt(L / (F * t2)) * ell.ay * cos(beta) * cos(lamb) + p3 = 1 / sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * sin(beta) * sin(2 * lamb) + p = np.array([p1, p2, p3]) + q = np.array([n[1] * p[2] - n[2] * p[1], + n[2] * p[0] - n[0] * p[2], + n[0] * p[1] - n[1] * p[0]]) - funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [x, dxds0, y, dyds0, z, dzds0], s, num, fein=False) - P2 = funktionswerte[-1] - P2 = (P2[0], P2[2], P2[4]) - return P2 + return p, q -def buildODE(ell): - def ODE(v): +def buildODE(ell: EllipsoidTriaxial) -> Callable: + """ + Aufbau des DGL-Systems + :param ell: Ellipsoid + :return: DGL-System + """ + def ODE(s: float, v: np.ndarray) -> np.ndarray: + """ + DGL-System + :param s: unabhängige Variable + :param v: abhängige Variablen + :return: Ableitungen der abhängigen Variablen + """ x, dxds, y, dyds, z, dzds = v H = ell.func_H(x, y, z) @@ -46,54 +65,78 @@ def buildODE(ell): ddy = -(h/H)*y/(1-ell.ee**2) ddz = -(h/H)*z/(1-ell.ex**2) - return [dxds, ddx, dyds, ddy, dzds, ddz] + return np.array([dxds, ddx, dyds, ddy, dzds, ddz]) return ODE -def gha1_num(ell, point, alpha0, s, num): +def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float]: + """ + Panou, Korakitits 2019 + :param ell: + :param point: + :param alpha0: + :param s: + :param num: + :return: + """ phi, lam, _ = ell.cart2geod(point, "ligas3") x0, y0, z0 = ell.geod2cart(phi, lam, 0) - p, q = ell.p_q(x0, y0, z0) + p, q = pq_ell(ell, x0, y0, z0) 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 = [x0, dxds0, y0, dyds0, z0, dzds0] + v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0]) - F = buildODE(ell) + ode = buildODE(ell) - werte = rk.rk_chat(F, v_init, s, num) - x1, _, y1, _, z1, _ = werte[-1] + _, werte = rk.rk4(ode, 0, v_init, s, num) + x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1] - return x1, y1, z1, werte + p1, q1 = pq_ell(ell, x1, y1, z1) + sigma = np.array([dx1ds, dy1ds, dz1ds]) + P = p1 @ sigma + Q = q1 @ sigma + alpha1 = arctan2(P, Q) -def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points): - constantValues = [] - for point in points: - x = point[1] - dxds = point[2] - y = point[3] - dyds = point[4] - z = point[5] - dzds = point[6] + if alpha1 < 0: + alpha1 += 2 * np.pi - values = ell.p_q(x, y, z) - p = values["p"] - q = values["q"] - t1 = values["t1"] - t2 = values["t2"] + return np.array([x1, y1, z1]), alpha1 - P = p[0]*dxds + p[1]*dyds + p[2]*dzds - Q = q[0]*dxds + q[1]*dyds + q[2]*dzds - alpha = arctan2(P, Q) +# --------------------------------------------------------------------------------------------------------------------- - c = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2) - constantValues.append(c) - pass +def pq_para(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]: + """ + Berechnung von p und q in parametrischen Koordinaten + Panou, Korakitits 2020 + :param x: x + :param y: y + :param z: z + :return: p und q + """ + n = ell.func_n(x, y, z) + u, v = ell.cart2para(np.array([x, y, z])) + # 41-47 + G = sqrt(1 - ell.ex ** 2 * cos(u) ** 2 - ell.ee ** 2 * sin(u) ** 2 * sin(v) ** 2) + q = np.array([-1 / G * sin(u) * cos(v), + -1 / G * sqrt(1 - ell.ee ** 2) * sin(u) * sin(v), + 1 / G * sqrt(1 - ell.ex ** 2) * cos(u)]) + p = np.array([q[1] * n[2] - q[2] * n[1], + q[2] * n[0] - q[0] * n[2], + q[0] * n[1] - q[1] * n[0]]) -def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): + t1 = np.dot(n, q) + t2 = np.dot(n, p) + t3 = np.dot(p, q) + if not (t1<1e-10 or t1>1-1e-10) and not (t2<1e-10 or t2>1-1e-10) and not (t3<1e-10 or t3>1-1e-10): + raise Exception("Fehler in den normierten Vektoren") + + return p, q + +def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): """ Panou, Korakitits 2020, 5ff. :param ell: @@ -104,44 +147,36 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): :return: """ x, y, z = point + + # S. 6 x_m = [x] y_m = [y] z_m = [z] - # erste Ableitungen (7-8) - H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2 - sqrtH = sqrt(H) - n = np.array([x / sqrtH, - y / ((1-ell.ee**2) * sqrtH), - z / ((1-ell.ex**2) * sqrtH)]) - u, v = ell.cart2para(np.array([x, y, z])) - G = sqrt(1 - ell.ex**2 * cos(u)**2 - ell.ee**2 * sin(u)**2 * sin(v)**2) - q = np.array([-1/G * sin(u) * cos(v), - -1/G * sqrt(1-ell.ee**2) * sin(u) * sin(v), - 1/G * sqrt(1-ell.ex**2) * cos(u)]) - p = np.array([q[1]*n[2] - q[2]*n[1], - q[2]*n[0] - q[0]*n[2], - q[0]*n[1] - q[1]*n[0]]) + p, q = pq_para(ell, x, y, z) + + # 48-50 x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0)) y_m.append(p[1] * sin(alpha0) + q[1] * cos(alpha0)) z_m.append(p[2] * sin(alpha0) + q[2] * cos(alpha0)) - # H Ableitungen (7) - H_ = lambda p: np.sum([comb(p, i) * (x_m[p - i] * x_m[i] + - 1 / (1-ell.ee**2) ** 2 * y_m[p-i] * y_m[i] + - 1 / (1-ell.ex**2) ** 2 * z_m[p-i] * z_m[i]) for i in range(0, p+1)]) + # 34 + H_ = lambda p: np.sum([comb(p, p - i) * (x_m[p - i] * x_m[i] + + 1 / (1 - ell.ee ** 2) ** 2 * y_m[p - i] * y_m[i] + + 1 / (1 - ell.ex ** 2) ** 2 * z_m[p - i] * z_m[i]) + for i in range(0, p + 1)]) - # h Ableitungen (7) - h_ = lambda q: np.sum([comb(q, j) * (x_m[q-j+1] * x_m[j+1] + - 1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] + - 1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1]) for j in range(0, q+1)]) + # 35 + h_ = lambda q: np.sum([comb(q, q-j) * (x_m[q-j+1] * x_m[j+1] + + 1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] + + 1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1]) + for j in range(0, q+1)]) - # h/H Ableitungen (6) - hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) * - np.sum([H_(t+1-l) / (fact(t+1-l) * fact(l-1)) * hH_t[l-1] for l in range(1, t+1)])) + # 31 + hH_ = lambda t: 1/H_(0) * (h_(t) - np.sum([comb(t, l-1) * H_(t+1-l) * hH_t[l-1] for l in range(1, t+1)])) - # xm, ym, zm Ableitungen (6) - x_ = lambda m: -np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)]) + # 28-30 + x_ = lambda m: - np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)]) y_ = lambda m: -1 / (1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)]) z_ = lambda m: -1 / (1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)]) @@ -155,11 +190,14 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): x_m.append(x_(m)) y_m.append(y_(m)) z_m.append(z_(m)) - a_m.append(x_m[m] / fact(m)) - b_m.append(y_m[m] / fact(m)) - c_m.append(z_m[m] / fact(m)) + fact_m = fact(m) - # am, bm, cm (6) + # 22-24 + a_m.append(x_m[m] / fact_m) + b_m.append(y_m[m] / fact_m) + c_m.append(z_m[m] / fact_m) + + # 19-21 x_s = 0 for a in reversed(a_m): x_s = x_s * s + a @@ -170,25 +208,89 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): for c in reversed(c_m): z_s = z_s * s + c - return x_s, y_s, z_s - pass + p_s, q_s = pq_para(ell, x_s, y_s, z_s) + + # 57-59 + dx_s = 0 + for i, a in reversed(list(enumerate(a_m[1:], start=1))): + dx_s = dx_s * s + i * a + + dy_s = 0 + for i, b in reversed(list(enumerate(b_m[1:], start=1))): + dy_s = dy_s * s + i * b + + dz_s = 0 + for i, c in reversed(list(enumerate(c_m[1:], start=1))): + dz_s = dz_s * s + i * c + + # 52-53 + sigma = np.array([dx_s, dy_s, dz_s]) + P = p_s @ sigma + Q = q_s @ sigma + + # 51 + alpha1 = arctan2(P, Q) + + if alpha1 < 0: + alpha1 += 2 * np.pi + + 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): + if s > np.pi / maxPartCircum * ell.ax: + s /= 2 + point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum) + point_end, alpha_end = gha1_ana(ell, point_step, alpha_step, s, maxM, maxPartCircum) + else: + point_end, alpha_end = gha1_ana_step(ell, point, alpha0, s, maxM) + + _, _, h = ell.cart2geod(point_end, "ligas3") + if h > 1e-5: + raise Exception("Analyitsche Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellpsoid") + + return point_end, alpha_end if __name__ == "__main__": - # ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") - # ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial") - re = ellipsoide.EllipsoidBiaxial.init_name("Bessel") - # Panou 2013, 7, Table 1, beta0=60° - beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = nep.get_example(table=1, example=3) - P0 = ell.ell2cart(beta0, lamb0) - P1 = ell.ell2cart(beta1, lamb1) + 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 = gha1_num(ell, P0, alpha0, s, 1000) - P1_num = gha1_num(ell, P0, alpha0, s, 10000) - P1_ana = gha1_ana(ell, P0, alpha0, s, 30) + 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) - beta, lamb = ellipsoide.EllipsoidTriaxial.cart2ell(ell, P1_num) + 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) + for example in examples_karney: + beta0, lamb0, alpha0, beta1, lamb1, 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) + try: + P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 40) + beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) + except: + beta1_ana, lamb1_ana = np.inf, np.inf + + diffs_karney.append((wu.rad2deg(abs(beta1-beta1_num)), wu.rad2deg(abs(lamb1-lamb1_num)), wu.rad2deg(abs(beta1-beta1_ana)), wu.rad2deg(abs(lamb1-lamb1_ana)))) + diffs_karney = np.array(diffs_karney) + 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 diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index ce2012c..568167e 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -1,7 +1,9 @@ import numpy as np from ellipsoide import EllipsoidTriaxial -import Numerische_Integration.num_int_runge_kutta as rk -import ausgaben as aus +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 # Panou 2013 def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30): @@ -110,27 +112,45 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep q_33, q_22, q_11, q_00) if lamb_1 != lamb_2: - def functions(): - def f_beta(lamb, beta, beta_p, X3, X4): - return beta_p + # def functions(): + # def f_beta(lamb, beta, beta_p, X3, X4): + # 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 - 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 + \ + 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 [f_beta, f_beta_p, f_X3, f_X4] + return np.array([dbeta, dbeta_p, dX3, dX4]) + + return ODE N = n @@ -144,15 +164,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep converged = False iterations = 0 - funcs = functions() + # funcs = functions() + ode_lamb = buildODElamb() for i in range(iter_max): iterations = i + 1 - startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0] + # startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0] + startwerte = np.array([beta_1, beta_0, 0.0, 1.0]) - werte = rk.verfahren(funcs, startwerte, dlamb, N) - lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1] + # werte = rk.verfahren(funcs, startwerte, dlamb, N) + lamb_list, werte = rk.rk4(ode_lamb, lamb_1, startwerte, dlamb, N, False) + # lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1] + lamb_end = lamb_list[-1] + beta_end, beta_p_end, X3_end, X4_end = werte[-1] d_beta_end_d_beta0 = X3_end delta = beta_end - beta_2 @@ -174,16 +199,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep raise RuntimeError("konvergiert nicht.") # Z - werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N) + # werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N, False) + lamb_list, werte = rk.rk4(ode_lamb, lamb_1, np.array([beta_1, beta_0, 0.0, 1.0]), dlamb, N, False) beta_arr = np.zeros(N + 1) - lamb_arr = np.zeros(N + 1) + # lamb_arr = np.zeros(N + 1) + lamb_arr = np.array(lamb_list) beta_p_arr = np.zeros(N + 1) for i, state in enumerate(werte): - lamb_arr[i] = state[0] - beta_arr[i] = state[1] - beta_p_arr[i] = state[2] + # lamb_arr[i] = state[0] + # beta_arr[i] = state[1] + # beta_p_arr[i] = state[2] + beta_arr[i] = state[0] + beta_p_arr[i] = state[1] (_, _, E1, G1, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) @@ -230,37 +259,59 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep converged = False iterations = 0 - def functions_beta(): - def g_lamb(beta, lamb, lamb_p, Y3, Y4): - return lamb_p + # 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 - 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 + \ + 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 [g_lamb, g_lamb_p, g_Y3, g_Y4] + return np.array([dlamb, dlamb_p, dY3, dY4]) + return ODE - funcs_beta = functions_beta() + # funcs_beta = functions_beta() + ode_beta = buildODEbeta() for i in range(iter_max): iterations = i + 1 - startwerte = [beta_1, lamb_1, lamb_0, 0.0, 1.0] + startwerte = [lamb_1, lamb_0, 0.0, 1.0] - werte = rk.verfahren(funcs_beta, startwerte, dbeta, N) - beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1] + # 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 @@ -279,16 +330,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep lamb_0 = lamb_0 - step - werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N) + # 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) + # beta_arr = np.zeros(N + 1) + beta_arr = np.array(beta_list) lamb_arr = np.zeros(N + 1) lambda_p_arr = np.zeros(N + 1) for i, state in enumerate(werte): - beta_arr[i] = state[0] - lamb_arr[i] = state[1] - lambda_p_arr[i] = state[2] + # beta_arr[i] = state[0] + # lamb_arr[i] = state[1] + # lambda_p_arr[i] = state[2] + lamb_arr[i] = state[0] + lambda_p_arr[i] = state[1] # Azimute (BETA1, LAMBDA1, E1, G1, @@ -318,22 +373,54 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - # beta1 = np.deg2rad(75) - # lamb1 = np.deg2rad(-90) - # beta2 = np.deg2rad(75) - # lamb2 = np.deg2rad(66) - # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2) - # print(aus.gms("a1", a1, 4)) - # print(aus.gms("a2", a2, 4)) + # ell = EllipsoidTriaxial.init_name("Fiction") + # # beta1 = np.deg2rad(75) + # # lamb1 = np.deg2rad(-90) + # # beta2 = np.deg2rad(75) + # # lamb2 = np.deg2rad(66) + # # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2) + # # print(aus.gms("a1", a1, 4)) + # # print(aus.gms("a2", a2, 4)) + # # print(s) + # cart1 = ell.para2cart(0, 0) + # cart2 = ell.para2cart(0.4, 1.4) + # beta1, lamb1 = ell.cart2ell(cart1) + # beta2, lamb2 = ell.cart2ell(cart2) + # + # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000) # print(s) - cart1 = ell.para2cart(0, 0) - cart2 = ell.para2cart(0.4, 0.4) - beta1, lamb1 = ell.cart2ell(cart1) - beta2, lamb2 = ell.cart2ell(cart2) - a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=2500) - 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("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/Numerische_Integration/num_int_runge_kutta.py b/Numerische_Integration/num_int_runge_kutta.py deleted file mode 100644 index ca12ab8..0000000 --- a/Numerische_Integration/num_int_runge_kutta.py +++ /dev/null @@ -1,84 +0,0 @@ -def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int, fein: bool = True) -> list: - """ - Runge-Kutta-Verfahren für ein beliebiges DGLS - :param funktionen: Liste mit allen Funktionen - :param startwerte: Liste mit allen Startwerten der Variablen - :param weite: gesamte Weite über die integriert werden soll - :param schritte: Anzahl der Schritte über die gesamte Weite - :return: Liste mit Listen für alle Wertepaare - """ - h = weite / schritte - werte = [startwerte] - for i in range(schritte): - - zuschlaege_grob = zuschlaege(funktionen, werte[-1], h) - werte_grob = [werte[-1][j] if j == 0 else werte[-1][j] + zuschlaege_grob[j - 1] - for j in range(len(startwerte))] - if fein: - zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2) - werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1] - for j in range(len(startwerte))] - - zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2) - werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1] - for j in range(len(startwerte))] - - werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j]) - for j in range(len(startwerte))] - - werte.append(werte_korr) - else: - werte.append(werte_grob) - return werte - - -def zuschlaege(funktionen: list, startwerte: list, h: float) -> list: - """ - Berechnung der Zuschläge eines einzelnen Schritts - :param funktionen: Liste mit allen Funktionen - :param startwerte: Liste mit allen Startwerten der Variablen - :param h: Schrittweite - :return: Liste mit Zuschlägen für die einzelnen Variablen - """ - werte = [wert for wert in startwerte] - - k1 = [h * funktion(*werte) for funktion in funktionen] - - werte = [startwerte[i] + (h / 2 if i == 0 else k1[i - 1] / 2) - for i in range(len(startwerte))] - - k2 = [h * funktion(*werte) for funktion in funktionen] - - werte = [startwerte[i] + (h / 2 if i == 0 else k2[i - 1] / 2) - for i in range(len(startwerte))] - - k3 = [h * funktion(*werte) for funktion in funktionen] - - werte = [startwerte[i] + (h if i == 0 else k3[i - 1]) - for i in range(len(startwerte))] - - k4 = [h * funktion(*werte) for funktion in funktionen] - - k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))] - - return k_ - -def rk_chat(F, v0: list, weite: float, schritte: int): - h = weite/schritte - v = v0 - werte = [v] - - for _ in range(schritte): - k1 = F(v) - k2 = F([v[i] + 0.5 * h * k1[i] for i in range(6)]) - k3 = F([v[i] + 0.5 * h * k2[i] for i in range(6)]) - k4 = F([v[i] + h * k3[i] for i in range(6)]) - - v = [ - v[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) - for i in range(6) - ] - - werte.append(v) - - return werte \ No newline at end of file diff --git a/Numerische_Integration/num_int_trapezformel.py b/Numerische_Integration/num_int_trapezformel.py deleted file mode 100644 index 4502d22..0000000 --- a/Numerische_Integration/num_int_trapezformel.py +++ /dev/null @@ -1,46 +0,0 @@ -def f(xi, yi): - ys = xi + yi - return ys - - -# Gegeben: -x0 = 0 -y0 = 0 -h = 0.2 -xmax = 0.4 -Ey = 0.0001 - -x = [x0] -y = [[y0]] - -n = 0 -while x[-1] < xmax: - x.append(x[-1]+h) - fn = f(x[n], y[n][-1]) - if n == 0: - y_neu = y[n][-1] + h * fn - else: - y_neu = y[n-1][-1] + 2*h * fn - y.append([y_neu]) - dy = 1 - while dy > Ey: - y_neu = y[n][-1] + h/2 * (fn + f(x[n+1], y[n+1][-1])) - y[-1].append(y_neu) - dy = abs(y[-1][-2]-y[-1][-1]) - n += 1 - -print(x) -print(y) - -werte = [] -for i in range(len(x)): - werte.append((x[i], y[i][-1])) - -for paar in werte: - print(f"({round(paar[0], 5)}, {round(paar[1], 5)})") - -integral = 0 -for i in range(len(werte)-1): - integral += (werte[i+1][1]+werte[i][1])/2 * (werte[i+1][0]-werte[i][0]) - -print(f"Integral = {round(integral, 5)}") diff --git a/Numerische_Integration/rk_DGLS_1Ordnung.py b/Numerische_Integration/rk_DGLS_1Ordnung.py deleted file mode 100644 index e9b8643..0000000 --- a/Numerische_Integration/rk_DGLS_1Ordnung.py +++ /dev/null @@ -1,7 +0,0 @@ -import num_int_runge_kutta as rk - -f = lambda ti, xi, yi: yi - ti -g = lambda ti, xi, yi: xi + yi - -funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 3) -print(funktionswerte) diff --git a/Numerische_Integration/rk_DGL_1Ordnung.py b/Numerische_Integration/rk_DGL_1Ordnung.py deleted file mode 100644 index 253cb85..0000000 --- a/Numerische_Integration/rk_DGL_1Ordnung.py +++ /dev/null @@ -1,6 +0,0 @@ -import num_int_runge_kutta as rk - -f = lambda xi, yi: xi + yi - -funktionswerte = rk.verfahren([f], [0, 0], 1, 5) -print(funktionswerte) diff --git a/Numerische_Integration/rk_DGL_2Ordnung.py b/Numerische_Integration/rk_DGL_2Ordnung.py deleted file mode 100644 index 37a86e4..0000000 --- a/Numerische_Integration/rk_DGL_2Ordnung.py +++ /dev/null @@ -1,10 +0,0 @@ -import numpy as np -import num_int_runge_kutta as rk - -f = lambda ti, ui, phii: -4 * np.sin(phii) -g = lambda ti, ui, phii: ui - -funktionswerte = rk.verfahren([f, g], [0, 1, 0], 0.6, 3) - -for wert in funktionswerte: - print(f"t = {round(wert[0],1)}s -> phi = {round(wert[2],5)}, phip = {round(wert[1],5)}, v = {round(2.45 * wert[1],5)}") diff --git a/Numerische_Integration/rk_DGL_2Ordnung_2.py b/Numerische_Integration/rk_DGL_2Ordnung_2.py deleted file mode 100644 index 4d8f4ab..0000000 --- a/Numerische_Integration/rk_DGL_2Ordnung_2.py +++ /dev/null @@ -1,7 +0,0 @@ -import num_int_runge_kutta as rk - -f = lambda xi, yi, ui: ui -g = lambda xi, yi, ui: 4 * ui - 4 * yi - -funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 2) -print(funktionswerte) diff --git a/dashboard.py b/dashboard.py index 9f354ea..e8841a9 100644 --- a/dashboard.py +++ b/dashboard.py @@ -169,6 +169,7 @@ app.layout = html.Div( {"label": "Bursa1970", "value": "Bursa1970"}, {"label": "BesselBiaxial", "value": "BesselBiaxial"}, {"label": "Fiction", "value": "Fiction"}, + {"label": "KarneyTest2024", "value": "KarneyTest2024"}, #{"label": "Ei", "value": "Ei"}, ], value="", diff --git a/ellipsoide.py b/ellipsoide.py index 399f12c..6cb82bb 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -1,8 +1,9 @@ import numpy as np -from numpy import sin, cos, arctan, arctan2, sqrt +from numpy import sin, cos, arctan, arctan2, sqrt, pi, arccos import winkelumrechnungen as wu import ausgaben as aus import jacobian_Ligas +import matplotlib.pyplot as plt class EllipsoidBiaxial: @@ -10,8 +11,8 @@ class EllipsoidBiaxial: self.a = a self.b = b self.c = a ** 2 / b - self.e = np.sqrt(a ** 2 - b ** 2) / a - self.e_ = np.sqrt(a ** 2 - b ** 2) / b + self.e = sqrt(a ** 2 - b ** 2) / a + self.e_ = sqrt(a ** 2 - b ** 2) / b @classmethod def init_name(cls, name: str): @@ -40,28 +41,28 @@ class EllipsoidBiaxial: b = a - a * f return cls(a, b) - V = lambda self, phi: np.sqrt(1 + self.e_ ** 2 * np.cos(phi) ** 2) + V = lambda self, phi: sqrt(1 + self.e_ ** 2 * cos(phi) ** 2) M = lambda self, phi: self.c / self.V(phi) ** 3 N = lambda self, phi: self.c / self.V(phi) - beta2psi = lambda self, beta: np.arctan(self.a / self.b * np.tan(beta)) - beta2phi = lambda self, beta: np.arctan(self.a ** 2 / self.b ** 2 * np.tan(beta)) + beta2psi = lambda self, beta: arctan(self.a / self.b * np.tan(beta)) + beta2phi = lambda self, beta: arctan(self.a ** 2 / self.b ** 2 * np.tan(beta)) - psi2beta = lambda self, psi: np.arctan(self.b / self.a * np.tan(psi)) - psi2phi = lambda self, psi: np.arctan(self.a / self.b * np.tan(psi)) + psi2beta = lambda self, psi: arctan(self.b / self.a * np.tan(psi)) + psi2phi = lambda self, psi: arctan(self.a / self.b * np.tan(psi)) - phi2beta = lambda self, phi: np.arctan(self.b ** 2 / self.a ** 2 * np.tan(phi)) - phi2psi = lambda self, phi: np.arctan(self.b / self.a * np.tan(phi)) + phi2beta = lambda self, phi: arctan(self.b ** 2 / self.a ** 2 * np.tan(phi)) + phi2psi = lambda self, phi: arctan(self.b / self.a * np.tan(phi)) - phi2p = lambda self, phi: self.N(phi) * np.cos(phi) + phi2p = lambda self, phi: self.N(phi) * cos(phi) def cart2ell(self, Eh, Ephi, x, y, z): - p = np.sqrt(x**2+y**2) + p = sqrt(x**2+y**2) # print(f"p = {round(p, 5)} m") - lamb = np.arctan(y/x) + lamb = arctan(y/x) - phi_null = np.arctan(z/p*(1-self.e**2)**-1) + phi_null = arctan(z/p*(1-self.e**2)**-1) hi = [0] phii = [phi_null] @@ -69,9 +70,9 @@ class EllipsoidBiaxial: i = 0 while True: - N = self.a*(1-self.e**2*np.sin(phii[i])**2)**(-1/2) - h = p/np.cos(phii[i])-N - phi = np.arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1)) + N = self.a*(1-self.e**2*sin(phii[i])**2)**(-1/2) + h = p/cos(phii[i])-N + phi = arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1)) hi.append(h) phii.append(phi) dh = abs(hi[i]-h) @@ -86,11 +87,11 @@ class EllipsoidBiaxial: return phi, lamb, h def ell2cart(self, phi, lamb, h): - W = np.sqrt(1 - self.e**2 * np.sin(phi)**2) + W = sqrt(1 - self.e**2 * sin(phi)**2) N = self.a / W - x = (N+h) * np.cos(phi) * np.cos(lamb) - y = (N+h) * np.cos(phi) * np.sin(lamb) - z = (N * (1-self.e**2) + h) * np.sin(lamb) + x = (N+h) * cos(phi) * cos(lamb) + y = (N+h) * cos(phi) * sin(lamb) + z = (N * (1-self.e**2) + h) * sin(lamb) return x, y, z class EllipsoidTriaxial: @@ -98,15 +99,15 @@ class EllipsoidTriaxial: self.ax = ax self.ay = ay self.b = b - self.ex = np.sqrt((self.ax**2 - self.b**2) / self.ax**2) - self.ey = np.sqrt((self.ay**2 - self.b**2) / self.ay**2) - self.ee = np.sqrt((self.ax**2 - self.ay**2) / self.ax**2) - self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2) - self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2) - self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2) - self.Ex = np.sqrt(self.ax**2 - self.b**2) - self.Ey = np.sqrt(self.ay**2 - self.b**2) - self.Ee = np.sqrt(self.ax**2 - self.ay**2) + self.ex = sqrt((self.ax**2 - self.b**2) / self.ax**2) + self.ey = sqrt((self.ay**2 - self.b**2) / self.ay**2) + self.ee = sqrt((self.ax**2 - self.ay**2) / self.ax**2) + self.ex_ = sqrt((self.ax**2 - self.b**2) / self.b**2) + self.ey_ = sqrt((self.ay**2 - self.b**2) / self.b**2) + self.ee_ = sqrt((self.ax**2 - self.ay**2) / self.ay**2) + self.Ex = sqrt(self.ax**2 - self.b**2) + self.Ey = sqrt(self.ay**2 - self.b**2) + self.Ee = sqrt(self.ax**2 - self.ay**2) @classmethod def init_name(cls, name: str): @@ -156,22 +157,21 @@ class EllipsoidTriaxial: b = 4000000 return cls(ax, ay, b) elif name == "KarneyTest2024": - ax = np.sqrt(2) + ax = sqrt(2) ay = 1 - b = 1 / np.sqrt(2) + b = 1 / sqrt(2) return cls(ax, ay, b) - def point_on(self, point: np.ndarray) -> bool: - """ - Test, ob ein Punkt auf dem Ellipsoid liegt. - :param point: kartesische 3D-Koordinaten - :return: Punkt auf dem Ellispoid? - """ - value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2 - if abs(1-value) < 0.000001: - return True - else: - return False + def func_H(self, x: float, y: float, z: float): + return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 + + def func_n(self, x: float, y: float, z: float, H: float = None): + if H is None: + H = self.func_H(x, y, z) + sqrtH = sqrt(H) + return np.array([x / sqrtH, + y / ((1 - self.ee ** 2) * sqrtH), + z / ((1 - self.ex ** 2) * sqrtH)]) def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: """ @@ -184,60 +184,12 @@ class EllipsoidTriaxial: :param u: Größe entlang der z-Achse :return: Punkt in kartesischen Koordinaten """ - # s1 = u**2 - self.b**2 - # s2 = -self.ay**2 * np.sin(beta)**2 - self.b**2 * np.cos(beta)**2 - # s3 = -self.ax**2 * np.sin(lamb)**2 - self.ay**2 * np.cos(lamb)**2 - # print(s1, s2, s3) - - # xe = np.sqrt(((self.ax**2+s1) * (self.ax**2+s2) * (self.ax**2+s3)) / - # ((self.ax**2-self.ay**2) * (self.ax**2-self.b**2))) - # ye = np.sqrt(((self.ay**2+s1) * (self.ay**2+s2) * (self.ay**2+s3)) / - # ((self.ay**2-self.ax**2) * (self.ay**2-self.b**2))) - # ze = np.sqrt(((self.b**2+s1) * (self.b**2+s2) * (self.b**2+s3)) / - # ((self.b**2-self.ax**2) * (self.b**2-self.ay**2))) - - x = np.sqrt(u**2 + self.Ex**2) * np.sqrt(np.cos(beta)**2 + self.Ee**2/self.Ex**2 * np.sin(beta)**2) * np.cos(lamb) - y = np.sqrt(u**2 + self.Ey**2) * np.cos(beta) * np.sin(lamb) - z = u * np.sin(beta) * np.sqrt(1 - self.Ee**2/self.Ex**2 * np.cos(lamb)**2) + x = sqrt(u**2 + self.Ex**2) * sqrt(cos(beta)**2 + self.Ee**2/self.Ex**2 * sin(beta)**2) * cos(lamb) + y = sqrt(u**2 + self.Ey**2) * cos(beta) * sin(lamb) + z = u * sin(beta) * sqrt(1 - self.Ee**2/self.Ex**2 * cos(lamb)**2) return np.array([x, y, z]) - def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: - """ - Panou, Korakitis 2019 2 - :param beta: elliptische Breite [rad] - :param lamb: elliptische Länge [rad] - :return: Punkt in kartesischen Koordinaten - """ - beta = np.asarray(beta, dtype=float) - lamb = np.asarray(lamb, dtype=float) - - beta, lamb = np.broadcast_arrays(beta, lamb) - - B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2 - L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2 - - x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb) - y = self.ay * np.cos(beta) * np.sin(lamb) - z = self.b / self.Ex * np.sin(beta) * np.sqrt(L) - - xyz = np.stack((x, y, z), axis=-1) - - # Pole - mask_south = beta == -np.pi / 2 - mask_north = beta == np.pi / 2 - xyz[mask_south] = np.array([0, 0, -self.b]) - xyz[mask_north] = np.array([0, 0, self.b]) - - # Äquator - mask_eq = beta == 0 - xyz[mask_eq & (lamb == -np.pi / 2)] = np.array([0, -self.ay, 0]) - xyz[mask_eq & (lamb == np.pi / 2)] = np.array([0, self.ay, 0]) - xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0]) - xyz[mask_eq & (lamb == np.pi)] = np.array([-self.ax, 0, 0]) - - return xyz - def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]: """ Panou 2014 15ff. @@ -253,25 +205,137 @@ class EllipsoidTriaxial: p = (c2**2 - 3*c1) / 9 q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54 - omega = np.arccos(q / np.sqrt(p**3)) + omega = arccos(q / sqrt(p**3)) - s1 = 2 * np.sqrt(p) * np.cos(omega/3) - c2/3 - s2 = 2 * np.sqrt(p) * np.cos(omega/3 - 2*np.pi/3) - c2/3 - s3 = 2 * np.sqrt(p) * np.cos(omega/3 - 4*np.pi/3) - c2/3 + s1 = 2 * sqrt(p) * cos(omega/3) - c2/3 + s2 = 2 * sqrt(p) * cos(omega/3 - 2*pi/3) - c2/3 + s3 = 2 * sqrt(p) * cos(omega/3 - 4*pi/3) - c2/3 # print(s1, s2, s3) - beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2))) + beta = arctan(sqrt((-self.b**2 - s2) / (self.ay**2 + s2))) if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7: - lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3))) + lamb = arctan(sqrt((-self.ay**2 - s3) / (self.ax**2 + s3))) else: lamb = 0 - u = np.sqrt(self.b**2 + s1) + u = sqrt(self.b**2 + s1) return beta, lamb, u - def cart2ell(self, point: np.ndarray) -> tuple[float, float]: + def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: """ - Panou, Korakitis 2019 2f. + Panou, Korakitis 2019 2 + :param beta: elliptische Breite [rad] + :param lamb: elliptische Länge [rad] + :return: Punkt in kartesischen Koordinaten + """ + beta = np.asarray(beta, dtype=float) + lamb = np.asarray(lamb, dtype=float) + + beta, lamb = np.broadcast_arrays(beta, lamb) + + B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2 + L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2 + + x = self.ax / self.Ex * sqrt(B) * cos(lamb) + y = self.ay * cos(beta) * sin(lamb) + z = self.b / self.Ex * sin(beta) * sqrt(L) + + xyz = np.stack((x, y, z), axis=-1) + + # Pole + mask_south = beta == -pi / 2 + mask_north = beta == pi / 2 + xyz[mask_south] = np.array([0, 0, -self.b]) + xyz[mask_north] = np.array([0, 0, self.b]) + + # Äquator + mask_eq = beta == 0 + xyz[mask_eq & (lamb == -pi / 2)] = np.array([0, -self.ay, 0]) + xyz[mask_eq & (lamb == pi / 2)] = np.array([0, self.ay, 0]) + xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0]) + xyz[mask_eq & (lamb == pi)] = np.array([-self.ax, 0, 0]) + + return xyz + + def ell2cart_bektas(self, beta: float | np.ndarray, omega: float | np.ndarray): + """ + Bektas 2015 + :param beta: + :param omega: + :return: + """ + x = self.ax * cos(omega) * sqrt((self.ax**2 - self.ay**2 * sin(beta)**2 - self.b**2 * cos(beta)**2) / (self.ax**2 - self.b**2)) + y = self.ay * cos(beta) * sin(omega) + z = self.b * sin(beta) * sqrt((self.ax**2 * sin(omega)**2 + self.ay**2 * cos(omega)**2 - self.b**2) / (self.ax**2 - self.b**2)) + + return np.array([x, y, z]) + + def ell2cart_karney(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: + """ + Karney 2025 Geographic Lib + :param beta: + :param lamb: + :return: + """ + e = sqrt(self.ax**2 - self.b**2) / self.ay + k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2) + k_ = sqrt(self.ax**2 - self.ay**2) / sqrt(self.ax**2 - self.b**2) + X = self.ax * cos(lamb) * sqrt(k**2*cos(beta)**2+k_**2) + Y = self.ay * cos(beta) * sin(lamb) + Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2) + return np.array([X, Y, Z]) + + def cart2ell(self, point, eps=1e-12, maxI = 100) -> tuple[float, float]: + """ + Panou, Korakitis 2019 3f. (num) + :param point: + :param eps: + :return: + """ + x, y, z = point + beta, lamb = self.cart2ell_panou(point) + delta_ell = np.array([np.inf, np.inf]).T + + i = 0 + while np.sum(delta_ell) > eps and i < maxI: + 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 + + 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)], + [-self.ay * sin(beta) * sin(lamb), + self.ay * cos(beta) * cos(lamb)], + [self.b / self.Ex * cos(beta) * sqrt(L), + (self.b * self.Ee ** 2) / (2 * self.Ex) * sin(beta) * sin(2 * lamb) / sqrt(L)]]) + + 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") + + 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") + + return beta, lamb + + def cart2ell_panou(self, point: np.ndarray) -> tuple[float, float]: + """ + Panou, Korakitis 2019 2f. (analytisch -> Näherung) :param point: Punkt in kartesischen Koordinaten :return: elliptische Breite, elliptische Länge """ @@ -280,18 +344,18 @@ class EllipsoidTriaxial: eps = 1e-9 if abs(x) < eps and abs(y) < eps: # Punkt in der z-Achse - beta = np.pi / 2 if z > 0 else -np.pi / 2 + beta = pi / 2 if z > 0 else -pi / 2 lamb = 0.0 return beta, lamb elif abs(x) < eps and abs(z) < eps: # Punkt in der y-Achse beta = 0.0 - lamb = np.pi / 2 if y > 0 else -np.pi / 2 + lamb = pi / 2 if y > 0 else -pi / 2 return beta, lamb elif abs(y) < eps and abs(z) < eps: # Punkt in der x-Achse beta = 0.0 - lamb = 0.0 if x > 0 else np.pi + lamb = 0.0 if x > 0 else pi return beta, lamb # ---- Allgemeiner Fall ----- @@ -300,15 +364,18 @@ class EllipsoidTriaxial: 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 + np.sqrt(c1 ** 2 - 4 * c0)) / 2 + t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2 + if t2 == 0: + t2 = 1e-14 t1 = c0 / t2 num_beta = max(t1 - self.b ** 2, 0) den_beta = max(self.ay ** 2 - t1, 0) num_lamb = max(t2 - self.ay ** 2, 0) den_lamb = max(self.ax ** 2 - t2, 0) - beta = np.arctan(np.sqrt(num_beta / den_beta)) - lamb = np.arctan(np.sqrt(num_lamb / den_lamb)) + + beta = arctan(sqrt(num_beta / den_beta)) + lamb = arctan(sqrt(num_lamb / den_lamb)) if z < 0: beta = -beta @@ -317,19 +384,62 @@ class EllipsoidTriaxial: lamb = -lamb if x < 0: - lamb = np.pi - lamb + lamb = pi - lamb if abs(x) < eps: - lamb = -np.pi/2 if y < 0 else np.pi/2 + lamb = -pi/2 if y < 0 else pi/2 elif abs(y) < eps: - lamb = 0 if x > 0 else np.pi + lamb = 0 if x > 0 else pi elif abs(z) < eps: beta = 0 return beta, lamb + def cart2ell_bektas(self, point, eps=1e-12, maxI = 100) -> tuple[float, float]: + """ + Bektas 2015 + :param point: + :param eps: + :param maxI: + :return: + """ + x, y, z = point + phi, lamb = self.cart2para(point) + p = sqrt((self.ax**2 - self.ay**2) / (self.ax**2 - self.b**2)) + d_phi = np.inf + d_lamb = np.inf + + i = 0 + while d_phi > eps and d_lamb > eps and i < maxI: + lamb_new = arctan2(self.ax * y * sqrt((p**2-1) * sin(phi)**2 + 1), self.ay * x * cos(phi)) + phi_new = arctan2(self.ay * z * sin(lamb), self.b * y * sqrt(1 - p**2 * cos(lamb)**2)) + d_phi = abs(phi_new - phi) + phi = phi_new + d_lamb = abs(lamb_new - lamb) + lamb = lamb_new + i += 1 + + if i == maxI: + raise Exception("Umrechung ist nicht konvergiert") + + return phi, lamb + + def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray: + """ + Ligas 2012, 250 + :param phi: geodätische Breite [rad] + :param lamb: geodätische Länge [rad] + :param h: Höhe über dem Ellipsoid + :return: kartesische Koordinaten + """ + v = self.ax / sqrt(1 - self.ex**2*sin(phi)**2-self.ee**2*cos(phi)**2*sin(lamb)**2) + xG = (v + h) * cos(phi) * cos(lamb) + yG = (v * (1-self.ee**2) + h) * cos(phi) * sin(lamb) + zG = (v * (1-self.ex**2) + h) * sin(phi) + return np.array([xG, yG, zG]) + def cart2geod(self, point: np.ndarray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: """ Ligas 2012 @@ -344,24 +454,24 @@ class EllipsoidTriaxial: eps = 1e-9 if abs(xG) < eps and abs(yG) < eps: # Punkt in der z-Achse - phi = np.pi / 2 if zG > 0 else -np.pi / 2 + phi = pi / 2 if zG > 0 else -pi / 2 lamb = 0.0 h = abs(zG) - self.b return phi, lamb, h elif abs(xG) < eps and abs(zG) < eps: # Punkt in der y-Achse phi = 0.0 - lamb = np.pi / 2 if yG > 0 else -np.pi / 2 + lamb = pi / 2 if yG > 0 else -pi / 2 h = abs(yG) - self.ay return phi, lamb, h elif abs(yG) < eps and abs(zG) < eps: # Punkt in der x-Achse phi = 0.0 - lamb = 0.0 if xG > 0 else np.pi + lamb = 0.0 if xG > 0 else pi h = abs(xG) - self.ax return phi, lamb, h - rG = np.sqrt(xG ** 2 + yG ** 2 + zG ** 2) + rG = sqrt(xG ** 2 + yG ** 2 + zG ** 2) pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) E = 1 / self.ax**2 @@ -380,38 +490,98 @@ class EllipsoidTriaxial: invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE) pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1) pEi = pEi.reshape(1, -1).flatten() - loa = np.sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2) + loa = sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2) pE = pEi i += 1 - phi = np.arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / np.sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2)) - lamb = np.arctan(1/(1-self.ee**2) * pE[1]/pE[0]) - h = np.sign(zG - pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2) + phi = arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2)) + lamb = arctan(1/(1-self.ee**2) * pE[1]/pE[0]) + h = np.sign(zG - pE[2]) * np.sign(pE[2]) * sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2) if xG < 0 and yG < 0: - lamb = -np.pi + lamb + lamb = -pi + lamb elif xG < 0: - lamb = np.pi + lamb + lamb = pi + lamb if abs(zG) < eps: phi = 0 return phi, lamb, h - def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray: + def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray: """ - Ligas 2012, 250 - :param phi: geodätische Breite [rad] - :param lamb: geodätische Länge [rad] - :param h: Höhe über dem Ellipsoid - :return: kartesische Koordinaten + Panou, Korakitits 2020, 4 + :param u: Parameter u + :param v: Parameter v + :return: Punkt in kartesischen Koordinaten """ - v = self.ax / np.sqrt(1 - self.ex**2*np.sin(phi)**2-self.ee**2*np.cos(phi)**2*np.sin(lamb)**2) - xG = (v + h) * np.cos(phi) * np.cos(lamb) - yG = (v * (1-self.ee**2) + h) * np.cos(phi) * np.sin(lamb) - zG = (v * (1-self.ex**2) + h) * np.sin(phi) - return np.array([xG, yG, zG]) + x = self.ax * cos(u) * cos(v) + y = self.ay * cos(u) * sin(v) + z = self.b * sin(u) + z = np.broadcast_to(z, np.shape(x)) + return np.array([x, y, z]) + + def cart2para(self, point: np.ndarray) -> tuple[float, float]: + """ + Panou, Korakitits 2020, 4 + :param point: Punkt in kartesischen Koordinaten + :return: parametrische Koordinaten + """ + x, y, z = point + + u_check1 = z*sqrt(1 - self.ee**2) + u_check2 = sqrt(x**2 * (1-self.ee**2) + y**2) * sqrt(1-self.ex**2) + if u_check1 <= u_check2: + u = arctan2(u_check1, u_check2) + else: + u = pi/2 - arctan2(u_check2, u_check1) + + v_check1 = y + v_check2 = x*sqrt(1-self.ee**2) + v_factor = sqrt(x**2*(1-self.ee**2)+y**2) + if v_check1 <= v_check2: + v = 2 * arctan2(v_check1, v_check2 + v_factor) + else: + v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor) + + return u, v + + def ell2para(self, beta: float, lamb: float) -> tuple[float, float]: + cart = self.ell2cart(beta, lamb) + return self.cart2para(cart) + + def para2ell(self, u: float, v: float) -> tuple[float, float]: + cart = self.para2cart(u, v) + return self.cart2ell(cart) + + def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + cart = self.para2cart(u, v) + return self.cart2geod(cart, mode, maxIter, maxLoa) + + def geod2para(self, phi: float, lamb: float, h: float) -> tuple[float, float]: + cart = self.geod2cart(phi, lamb, h) + return self.cart2para(cart) + + def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + cart = self.ell2cart(beta, lamb) + return self.cart2geod(cart, mode, maxIter, maxLoa) + + def geod2ell(self, phi: float, lamb: float, h: float) -> tuple[float, float]: + cart = self.geod2cart(phi, lamb, h) + return self.cart2ell(cart) + + def point_on(self, point: np.ndarray) -> bool: + """ + Test, ob ein Punkt auf dem Ellipsoid liegt. + :param point: kartesische 3D-Koordinaten + :return: Punkt auf dem Ellispoid? + """ + value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2 + if abs(1-value) < 1e-6: + return True + else: + return False def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]: """ @@ -434,144 +604,64 @@ class EllipsoidTriaxial: pointH = self. geod2cart(phi, lamb, h) return pointH - def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray: - """ - Panou, Korakitits 2020, 4 - :param u: Parameter u - :param v: Parameter v - :return: Punkt in kartesischen Koordinaten - """ - x = self.ax * np.cos(u) * np.cos(v) - y = self.ay * np.cos(u) * np.sin(v) - z = self.b * np.sin(u) - z = np.broadcast_to(z, np.shape(x)) - return np.array([x, y, z]) - - def cart2para(self, point: np.ndarray) -> tuple[float, float]: - """ - Panou, Korakitits 2020, 4 - :param point: Punkt in kartesischen Koordinaten - :return: parametrische Koordinaten - """ - x, y, z = point - - u_check1 = z*np.sqrt(1 - self.ee**2) - u_check2 = np.sqrt(x**2 * (1-self.ee**2) + y**2) * np.sqrt(1-self.ex**2) - if u_check1 <= u_check2: - u = np.arctan2(u_check1, u_check2) - else: - u = np.pi/2 - np.arctan2(u_check2, u_check1) - - v_check1 = y - v_check2 = x*np.sqrt(1-self.ee**2) - v_factor = np.sqrt(x**2*(1-self.ee**2)+y**2) - if v_check1 <= v_check2: - v = 2 * np.arctan2(v_check1, v_check2 + v_factor) - else: - v = np.pi/2 - 2 * np.arctan2(v_check2, v_check1 + v_factor) - - return u, v - - def ell2para(self, beta, lamb) -> tuple[float, float]: - cart = self.ell2cart(beta, lamb) - return self.cart2para(cart) - - def para2ell(self, u, v) -> tuple[float, float]: - cart = self.para2cart(u, v) - return self.cart2ell(cart) - - def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: - cart = self.para2cart(u, v) - return self.cart2geod(cart, mode, maxIter, maxLoa) - - def geod2para(self, phi, lamb, h) -> tuple[float, float]: - cart = self.geod2cart(phi, lamb, h) - return self.cart2para(cart) - - def ell2geod(self, beta, lamb, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: - cart = self.ell2cart(beta, lamb) - return self.cart2geod(cart, mode, maxIter, maxLoa) - - def func_H(self, x, y, z): - return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 - - def func_n(self, x, y, z, H=None): - if H is None: - H = self.func_H(x, y, z) - return np.array([x / sqrt(H), - y / ((1 - self.ee ** 2) * sqrt(H)), - z / ((1 - self.ex ** 2) * sqrt(H))]) - - def p_q(self, x, y, z) -> tuple[np.ndarray, np.ndarray]: - """ - Berechnung sämtlicher Größen - :param x: x - :param y: y - :param z: z - :return: Dictionary sämtlicher Größen - """ - n = self.func_n(x, y, z) - - beta, lamb = self.cart2ell(np.array([x, y, z])) - B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2 - L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2 - - 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 + np.sqrt(c1 ** 2 - 4 * c0)) / 2 - t1 = c0 / t2 - t2e = self.ax ** 2 * np.sin(lamb) ** 2 + self.ay ** 2 * np.cos(lamb) ** 2 - t1e = self.ay ** 2 * np.sin(beta) ** 2 + self.b ** 2 * np.cos(beta) ** 2 - - F = self.Ey ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(lamb) ** 2 - p1 = -np.sqrt(L / (F * t2)) * self.ax / self.Ex * np.sqrt(B) * np.sin(lamb) - p2 = np.sqrt(L / (F * t2)) * self.ay * np.cos(beta) * np.cos(lamb) - p3 = 1 / np.sqrt(F * t2) * (self.b * self.Ee ** 2) / (2 * self.Ex) * np.sin(beta) * np.sin(2 * lamb) - # p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * self.ax / (self.Ex * self.Ee) * np.sqrt(B) * np.sqrt(t2 - self.ay ** 2) - # p2 = np.sign(x) * np.sqrt(L / (F * t2)) * self.ay / (self.Ey * self.Ee) * np.sqrt((ell.ay ** 2 - t1) * (self.ax ** 2 - t2)) - # p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * self.b / (self.Ex * self.Ey) * np.sqrt( - # (t1 - self.b ** 2) * (t2 - self.ay ** 2) * (self.ax ** 2 - t2)) - p = np.array([p1, p2, p3]) - q = np.array([n[1] * p[2] - n[2] * p[1], - n[2] * p[0] - n[0] * p[2], - n[0] * p[1] - n[1] * p[0]]) - - return p, q - if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("Eitschberger1978") + ell = EllipsoidTriaxial.init_name("BursaSima1980round") diff_list = [] - for beta_deg in range(-150, 210, 30): - for lamb_deg in range(-150, 210, 30): - beta = wu.deg2rad(beta_deg) - lamb = wu.deg2rad(lamb_deg) - point = ell.ell2cart(beta, lamb) + diffs_para = [] + diffs_ell = [] + diffs_geod = [] + points = [] + for v_deg in range(-180, 181, 5): + for u_deg in range(-90, 91, 5): + v = wu.deg2rad(v_deg) + u = wu.deg2rad(u_deg) + point = ell.para2cart(u, v) + points.append(point) elli = ell.cart2ell(point) cart_elli = ell.ell2cart(elli[0], elli[1]) - diff_ell = np.sum(np.abs(point-cart_elli)) + diff_ell = np.linalg.norm(point - cart_elli, axis=-1) para = ell.cart2para(point) cart_para = ell.para2cart(para[0], para[1]) - diff_para = np.sum(np.abs(point-cart_para)) - - # geod = ell.cart2geod(point, "ligas1") - # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - # diff_geod1 = np.sum(np.abs(point-cart_geod)) - # - # geod = ell.cart2geod(point, "ligas2") - # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - # diff_geod2 = np.sum(np.abs(point-cart_geod)) + diff_para = np.linalg.norm(point - cart_para, axis=-1) geod = ell.cart2geod(point, "ligas3") cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - diff_geod3 = np.sum(np.abs(point-cart_geod)) + diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1) - diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod3]) - diff_list.append([diff_ell]) + diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3]) # + diffs_ell.append([diff_ell]) + diffs_para.append([diff_para]) + diffs_geod.append([diff_geod3]) diff_list = np.array(diff_list) - pass \ No newline at end of file + diffs_ell = np.array(diffs_ell) + diffs_para = np.array(diffs_para) + diffs_geod = np.array(diffs_geod) + + pass + points = np.array(points) + fig = plt.figure() + ax = fig.add_subplot(projection='3d') + + sc = ax.scatter( + points[:, 0], + points[:, 1], + points[:, 2], + c=diffs_ell, # Farbcode = diff + cmap='viridis', # Colormap + s=10 + 20 * diffs_ell, # optional: Größe abhängig vom diff + alpha=0.8 + ) + + # Farbskala + cbar = plt.colorbar(sc) + cbar.set_label("diff") + + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Z") + + plt.show() \ No newline at end of file diff --git a/phi_ellipse.py b/phi_ellipse.py deleted file mode 100644 index 29b1a37..0000000 --- a/phi_ellipse.py +++ /dev/null @@ -1,55 +0,0 @@ -import winkelumrechnungen as wu - - -def polyapp_tscheby_hayford(s: float) -> float: - """ - Berechnung der ellipsoidisch geodätischen Breite. - Polynomapproximation mittels Tschebyscheff-Polynomen. - Auf dem Hayford-Ellipsoid. - :param s: Strecke auf einer Ellipse vom Äquator aus - :type s: float - :return: ellipsoidisch geodätische Breite - :rtype: float - """ - c0 = 1 - c1 = -0.00837809325 - c2 = 0.00428127367 - c3 = -0.00114523986 - c4 = 0.00023219707 - c5 = -0.00004421222 - c6 = 0.00000570244 - alpha = wu.gms2rad([0, 0, 325643.97199]) - s90 = 10002288.2990 - - xi = s/s90 - - phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) + - c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6)) - return phi - - -def polyapp_tscheby_bessel(s: float) -> float: - """ - Berechnung der ellipsoidisch geodätischen Breite. - Polynomapproximation mittels Tschebyscheff-Polynomen. - Auf dem Bessel-Ellipsoid. - :param s: Strecke auf einer Ellipse vom Äquator aus - :type s: float - :return: ellipsoidisch geodätische Breite - :rtype: float - """ - c0 = 1 - c1 = -0.00831729565 - c2 = 0.00424914906 - c3 = -0.00113566119 - c4 = 0.00022976983 - c5 = -0.00004363980 - c6 = 0.00000562025 - alpha = wu.gms2rad([0, 0, 325632.08677]) - s90 = 10000855.7644 - - xi = s/s90 - - phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) + - c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6)) - return phi diff --git a/runge_kutta.py b/runge_kutta.py new file mode 100644 index 0000000..4ab408e --- /dev/null +++ b/runge_kutta.py @@ -0,0 +1,43 @@ +import numpy as np + + +def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]: + """ + Standard Runge-Kutta Verfahren 4. Ordnung + :param ode: ODE-System als Funktion + :param t0: Startwert der unabhängigen Variable + :param v0: Startwerte + :param weite: Integrationsweite + :param schritte: Schrittzahl + :param fein: + :return: Variable und Funktionswerte an jedem Stützpunkt + """ + h = weite/schritte + + t_list = [t0] + werte = [v0] + + for _ in range(schritte): + t = t_list[-1] + v = werte[-1] + + if not fein: + v_next = rk4_step(ode, t, v, h) + + else: + v_grob = rk4_step(ode, t, v, h) + v_half = rk4_step(ode, t, v, 0.5 * h) + v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h) + v_next = v_fein + (v_fein - v_grob) / 15.0 + + t_list.append(t + h) + werte.append(v_next) + + return t_list, werte + +def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray: + k1 = ode(t, v) + k2 = ode(t + 0.5 * h, v + 0.5 * h * k1) + k3 = ode(t + 0.5 * h, v + 0.5 * h * k2) + k4 = ode(t + h, v + h * k3) + return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4) \ No newline at end of file diff --git a/s_ellipse.py b/s_ellipse.py deleted file mode 100644 index b9ab052..0000000 --- a/s_ellipse.py +++ /dev/null @@ -1,78 +0,0 @@ -import numpy as np - - -def reihenentwicklung(es: float, c: float, phi: float) -> float: - """ - Berechnung der Strecke auf einer Ellipse. - Reihenentwicklung. - :param es: zweite numerische Exzentrizität - :type es: float - :param c: Polkrümmungshalbmesser - :type c: float - :param phi: ellipsoidisch geodästische Breite in Radiant - :type phi: float - :return: Strecke auf einer Ellipse vom Äquator aus - :rtype: float - """ - Ass = 1 - 3/4*es**2 + 45/64*es**4 - 175/256*es**6 + 11025/16384*es**8 - Bss = - 3/4*es**2 + 15/16*es**4 - 525/512*es**6 + 2205/2048*es**8 - Css = 15/64*es**4 - 105/256*es**6 + 2205/4096*es**8 - Dss = - 35/512*es**6 + 315/2048*es**8 - print(f"A'' = {round(Ass, 10):.10f}\nB'' = {round(Bss, 10):.10f}\nC'' = {round(Css, 10):.10f}\nD'' = {round(Dss, 10):.10f}") - - s = c * (Ass*phi + 1/2*Bss * np.sin(2*phi) + 1/4*Css * np.sin(4*phi) + 1/6*Dss * np.sin(6*phi)) - return s - - -def polyapp_tscheby_hayford(phi: float) -> float: - """ - Berechnung der Strecke auf einer Ellipse. - Polynomapproximation mittels Tschebyscheff-Polynomen. - Auf dem Hayford-Ellipsoid. - :param phi: ellipsoidisch geodästische Breite in Radiant - :type phi: float - :return: Strecke auf einer Ellipse vom Äquator aus - :rtype: float - """ - c1s = 0.00829376218 - c2s = -0.00398963425 - c3s = 0.00084200710 - c4s = -0.0000648906 - c5s = -0.00001075680 - c6s = 0.00000396474 - c7s = -0.00000046347 - alpha = 9951793.0123 - - xi = 2 / np.pi * phi - xis = xi ^ 2 - - s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3 - + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7) - return s - - -def polyapp_tscheby_bessel(phi: float) -> float: - """ - Berechnung der Strecke auf einer Ellipse. - Polynomapproximation mittels Tschebyscheff-Polynomen. - Auf dem Bessel-Ellipsoid. - :param phi: ellipsoidisch geodästische Breite in Radiant - :type phi: float - :return: Strecke auf einer Ellipse vom Äquator aus - :rtype: float - """ - c1s = 0.00823417717 - c2s = -0.00396170744 - c3s = 0.00083680249 - c4s = -0.00006488462 - c5s = -0.00001053242 - c6s = 0.00000390854 - c7s = -0.00000045768 - alpha = 9950730.8876 - - xi = 2 / np.pi * phi - xis = xi ** 2 - - s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3 - + c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7) - return s diff --git a/show_constant_lines.py b/show_constant_lines.py deleted file mode 100644 index d20c6e6..0000000 --- a/show_constant_lines.py +++ /dev/null @@ -1,30 +0,0 @@ -import numpy as np -import plotly.graph_objects as go -from ellipsoide import EllipsoidTriaxial -import winkelumrechnungen as wu -from dashboard import ellipsoid_figure - -u = np.linspace(0, 2*np.pi, 51) -v = np.linspace(0, np.pi, 51) -ell = EllipsoidTriaxial.init_name("BursaSima1980round") -points = [] -lines = [] -for u_i, u_value in enumerate(u): - for v_i, v_value in enumerate(v): - cart = ell.ell2cart(u_value, v_value) - if u_i != 0 and v_i != 0: - lines.append((points[-1], cart, "red")) - points.append(cart) -points = [] -for v_i, v_value in enumerate(v): - for u_i, u_value in enumerate(u): - cart = ell.ell2cart(u_value, v_value) - if u_i != 0 and v_i != 0: - lines.append((points[-1], cart, "blue")) - points.append(cart) -ax = ell.ax -ay = ell.ay -b = ell.b - -figu = ellipsoid_figure(ax, ay, b, lines=lines) -figu.show() \ No newline at end of file diff --git a/test.py b/test.py index 01cc87c..4d8bd3b 100644 --- a/test.py +++ b/test.py @@ -1,29 +1,14 @@ import numpy as np from scipy.special import factorial as fact from math import comb +import matplotlib.pyplot as plt +import ellipsoide +from GHA_triaxial.panou import pq -J = np.array([ - [2, 3, 0], - [0, 3, 0], - [6, 0, 4] -]) +ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") -xi = np.array([1, 2, 3]) -xi_col = xi.reshape(-1, 1) -print(xi_col) -xi_row = xi_col.reshape(1, -1).flatten() -print(xi_row) +x, y, z = ell.para2cart(0.5, 0.7) -# Spaltenvektor-Variante -res_col = xi[:, None] - J @ xi[:, None] - -# Zeilenvektor-Variante -res_row = xi[None, :] - xi[None, :] @ J - -print("Spaltenvektor:") -print(res_col[0,0]) -print("Zeilenvektor:") -print(res_row) -t = 5 -l = 2 -print(fact(t+1-l) / (fact(t+1-l) * fact(l-1)), comb(l-1, t+1-l)) \ No newline at end of file +x1 = pq(ell, x, y, z) +x2 = ell.p_q(x, y, z) +pass \ No newline at end of file diff --git a/test_algorithms.py b/test_algorithms.py deleted file mode 100644 index f48298b..0000000 --- a/test_algorithms.py +++ /dev/null @@ -1,84 +0,0 @@ -import GHA_triaxial.numeric_examples_panou as nep -import ellipsoide -from GHA_triaxial.panou_2013_2GHA_num import gha2_num -from GHA_triaxial.panou import gha1_ana, gha1_num -import numpy as np -import time - -def test(): - ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") - - tables = nep.get_tables() - - diffs_gha1_num = [] - diffs_gha1_ana = [] - diffs_gha2_num = [] - times_gha1_num = [] - times_gha1_ana = [] - times_gha2_num = [] - - for table in tables: - diffs_gha1_num.append([]) - diffs_gha1_ana.append([]) - diffs_gha2_num.append([]) - times_gha1_num.append([]) - times_gha1_ana.append([]) - times_gha2_num.append([]) - - for example in table: - beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = example - P0 = ell.ell2cart(beta0, lamb0) - P1 = ell.ell2cart(beta1, lamb1) - - start = time.perf_counter() - try: - P1_num = gha1_num(ell, P0, alpha0, s, 10000) - end = time.perf_counter() - diff_P1_num = np.linalg.norm(P1 - P1_num) - except: - end = time.perf_counter() - diff_P1_num = None - time_gha1_num = end - start - - start = time.perf_counter() - try: - P1_ana = gha1_ana(ell, P0, alpha0, s, 50) - end = time.perf_counter() - diff_P1_ana = np.linalg.norm(P1 - P1_ana) - except: - end = time.perf_counter() - diff_P1_ana = None - time_gha1_ana = end - start - - start = time.perf_counter() - try: - alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000) - end = time.perf_counter() - diff_s_num = abs(s - s_num) - except: - end = time.perf_counter() - diff_s_num = None - time_gha2_num = None - time_gha2_num = end - start - - diffs_gha1_num[-1].append(diff_P1_num) - diffs_gha1_ana[-1].append(diff_P1_ana) - diffs_gha2_num[-1].append(diff_s_num) - times_gha1_num[-1].append(time_gha1_num) - times_gha1_ana[-1].append(time_gha1_ana) - times_gha2_num[-1].append(time_gha2_num) - print(diffs_gha1_num, diffs_gha1_ana, diffs_gha2_num) - print(times_gha1_num, times_gha1_ana, times_gha2_num) - -def display(): - diffs = [[{'gha1_num': np.float64(3.410763124264611e-05), 'gha1_ana': np.float64(3.393273802112796e-05), 'gha2_num': np.float64(3.3931806683540344e-05)}, {'gha1_num': np.float64(0.0008736425000530604), 'gha1_ana': np.float64(0.0008736458415010259), 'gha2_num': None}, {'gha1_num': np.float64(0.0007739730058338136), 'gha1_ana': np.float64(0.0007739621469802854), 'gha2_num': np.float64(1.5832483768463135e-07)}, {'gha1_num': np.float64(0.00010554956741100295), 'gha1_ana': np.float64(8.814246009944831), 'gha2_num': np.float64(4.864111542701721e-05)}, {'gha1_num': np.float64(0.0002135908394614854), 'gha1_ana': np.float64(0.0002138610897967267), 'gha2_num': np.float64(5.0179407158866525)}, {'gha1_num': np.float64(0.00032727226891456654), 'gha1_ana': np.float64(0.00032734569198545905), 'gha2_num': np.float64(9.735533967614174e-05)}, {'gha1_num': np.float64(0.0005195973303787956), 'gha1_ana': np.float64(0.0005197766935509641), 'gha2_num': None}], [{'gha1_num': np.float64(1.780250537652368e-05), 'gha1_ana': np.float64(1.996805145339501e-05), 'gha2_num': np.float64(1.8164515495300293e-05)}, {'gha1_num': np.float64(4.8607540473363564e-05), 'gha1_ana': np.float64(2205539.954949392), 'gha2_num': None}, {'gha1_num': np.float64(0.00017376854985685854), 'gha1_ana': np.float64(328124.1513636429), 'gha2_num': np.float64(0.17443156614899635)}, {'gha1_num': np.float64(5.83429352558999e-05), 'gha1_ana': np.float64(0.01891628037258558), 'gha2_num': np.float64(1.4207654744386673)}, {'gha1_num': np.float64(0.0006421087024666934), 'gha1_ana': np.float64(0.0006420400127297228), 'gha2_num': np.float64(0.12751091085374355)}, {'gha1_num': np.float64(0.0004456207867164434), 'gha1_ana': np.float64(0.0004455649707698245), 'gha2_num': np.float64(0.00922046648338437)}, {'gha1_num': np.float64(0.0002340879908275419), 'gha1_ana': np.float64(0.00023422217242111216), 'gha2_num': np.float64(0.001307751052081585)}], [{'gha1_num': np.float64(976.6580096633622), 'gha1_ana': np.float64(976.6580096562798), 'gha2_num': np.float64(6.96033239364624e-05)}, {'gha1_num': np.float64(2825.2936643258527), 'gha1_ana': np.float64(2794.954866417055), 'gha2_num': np.float64(1.3615936040878296e-05)}, {'gha1_num': np.float64(1248.8942058074501), 'gha1_ana': np.float64(538.5550561841195), 'gha2_num': np.float64(3.722589462995529e-05)}, {'gha1_num': np.float64(2201.1793359793814), 'gha1_ana': np.float64(3735.376499414938), 'gha2_num': np.float64(1.4525838196277618e-05)}, {'gha1_num': np.float64(2262.134819997246), 'gha1_ana': np.float64(25549.567793410763), 'gha2_num': np.float64(9.328126907348633e-06)}, {'gha1_num': np.float64(2673.219788119847), 'gha1_ana': np.float64(21760.866677295206), 'gha2_num': np.float64(8.635222911834717e-06)}, {'gha1_num': np.float64(1708.758419275875), 'gha1_ana': np.float64(3792.1128807063437), 'gha2_num': np.float64(2.4085864424705505e-05)}], [{'gha1_num': np.float64(0.7854659044152204), 'gha1_ana': np.float64(0.785466068424286), 'gha2_num': np.float64(0.785466069355607)}, {'gha1_num': np.float64(237.79878717216718), 'gha1_ana': np.float64(1905080.064324282), 'gha2_num': None}, {'gha1_num': np.float64(55204.601699830164), 'gha1_ana': np.float64(55204.60175211949), 'gha2_num': None}, {'gha1_num': np.float64(12766.348063015519), 'gha1_ana': np.float64(12766.376619517901), 'gha2_num': np.float64(12582.786206113175)}, {'gha1_num': np.float64(29703.049988324146), 'gha1_ana': np.float64(29703.056427749252), 'gha2_num': np.float64(28933.668131249025)}, {'gha1_num': np.float64(43912.03007182513), 'gha1_ana': np.float64(43912.03007528712), 'gha2_num': None}, {'gha1_num': np.float64(28522.29828970693), 'gha1_ana': np.float64(28522.29830145182), 'gha2_num': None}, {'gha1_num': np.float64(17769.115549537233), 'gha1_ana': np.float64(17769.115549483362), 'gha2_num': np.float64(17769.121286311187)}]] - arr = [] - for table in diffs: - for example in table: - arr.append([example['gha1_num'], example['gha1_ana'], example['gha2_num']]) - arr = np.array(arr) - pass - -if __name__ == "__main__": - # test() - display() From 4d5b6fcc3ef980cc1c4fb1060cc7e420daa44ff1 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sun, 11 Jan 2026 13:14:43 +0100 Subject: [PATCH 2/3] =?UTF-8?q?Umrechnung=20alpha,=20N=C3=A4herungsl=C3=B6?= =?UTF-8?q?sung=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) From 6cc7245b0ff1562c32c179d7f07c9bd6209d031f Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sun, 11 Jan 2026 16:05:15 +0100 Subject: [PATCH 3/3] =?UTF-8?q?N=C3=A4herungsl=C3=B6sung=20GHA=202?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA_triaxial/approx_gha1.py | 19 ++--- GHA_triaxial/approx_gha2.py | 111 ++++++++++++++++++++++++++++ GHA_triaxial/approx_louville.py | 108 --------------------------- GHA_triaxial/panou.py | 103 +++++++++++++------------- GHA_triaxial/panou_2013_2GHA_num.py | 31 +++++--- dashboard.py | 32 +++----- ellipsoide.py | 91 ++++++++++++++--------- 7 files changed, 260 insertions(+), 235 deletions(-) create mode 100644 GHA_triaxial/approx_gha2.py delete mode 100644 GHA_triaxial/approx_louville.py diff --git a/GHA_triaxial/approx_gha1.py b/GHA_triaxial/approx_gha1.py index 67fdd2c..f87180d 100644 --- a/GHA_triaxial/approx_gha1.py +++ b/GHA_triaxial/approx_gha1.py @@ -1,12 +1,10 @@ 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): +def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]: l0 = louville_constant(ell, p0, alpha0) points = [p0] alphas = [alpha0] @@ -17,10 +15,9 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: in break p1 = points[-1] alpha1 = alphas[-1] - x1, y1, z1 = p1 - sigma = func_sigma_ell(ell, x1, y1, z1, alpha1) + sigma = func_sigma_ell(ell, p1, alpha1) p2 = p1 + ds_step * sigma - p2, _, _, _ = ell.cartonell(p2) + p2 = ell.cartonell(p2) ds_step = np.linalg.norm(p2 - p1) points.append(p2) @@ -30,7 +27,11 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: in alpha2 = alpha1 + (l0 - l2) / dl_dalpha alphas.append(alpha2) s_curr += ds_step - return points[-1], alphas[-1], np.array(points) + + if all_points: + return points[-1], alphas[-1], np.array(points) + else: + return points[-1], alphas[-1] def show_points(points, p0, p1): fig = go.Figure() @@ -57,7 +58,7 @@ if __name__ == '__main__': 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) + P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=32) + P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, ds=5000, all_points=True) 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 new file mode 100644 index 0000000..1fe8b40 --- /dev/null +++ b/GHA_triaxial/approx_gha2.py @@ -0,0 +1,111 @@ +import numpy as np +from numpy import arctan2 +from ellipsoide import EllipsoidTriaxial +from panou import pq_ell +from panou_2013_2GHA_num import gha2_num +import plotly.graph_objects as go +import winkelumrechnungen as wu +from numpy.typing import NDArray +from typing import Tuple + +def sigma2alpha(sigma: NDArray, point: NDArray) -> float: + """ + Berechnung des Richtungswinkels an einem Punkt anhand der Ableitung zu den kartesischen Koordinaten + :param sigma: Ableitungsvektor ver kartesischen Koordinaten + :param point: Punkt + :return: Richtungswinkel + """ + "" + p, q = pq_ell(ell, point) + P = float(p @ sigma) + Q = float(q @ sigma) + + alpha = arctan2(P, Q) + return alpha + +def gha2(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 + :param p0: Startpunkt + :param p1: Endpunkt + :param ds: maximales Streckenelement + :param all_points: Alle Punkte ausgeben? + :return: + """ + points = np.array([p0, p1]) + while True: + new_points = [] + + for i in range(len(points)-1): + new_points.append(points[i]) + + pi = points[i] + 1/2 * (points[i+1] - points[i]) + pi = ell.cartonell(pi) + + new_points.append(pi) + + new_points.append(points[-1]) + points = np.array(new_points) + + elements = np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)]) + + if np.average(elements) < ds: + break + + p0i = ell.cartonell(p0 + ds/100 * (points[1] - p0) / np.linalg.norm(points[1] - p0)) + sigma0 = (p0i - p0) / np.linalg.norm(p0i - p0) + alpha0 = sigma2alpha(sigma0, p0) + + p1i = ell.cartonell(p1 - ds/100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2])) + sigma1 = (p1 - p1i) / np.linalg.norm(p1 - p1i) + alpha1 = sigma2alpha(sigma1, p1) + + s = np.sum(np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)])) + + if all_points: + return alpha0, alpha1, s, np.array(points) + else: + return alpha0, alpha1, s + +def show_points(points: NDArray, points_app: NDArray, p0: NDArray, p1: NDArray): + fig = go.Figure() + + fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2], + mode='lines', line=dict(color="green", width=3), name="Analytisch") + fig.add_scatter3d(x=points_app[:, 0], y=points_app[:, 1], z=points_app[:, 2], + mode='lines', line=dict(color="red", width=3), name="Approximiert") + fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]], + mode='markers', marker=dict(color="black"), name="P0") + fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]], + mode='markers', marker=dict(color="black"), name="P1") + + fig.update_layout( + scene=dict(xaxis_title='X [km]', + yaxis_title='Y [km]', + zaxis_title='Z [km]', + aspectmode='data')) + + fig.show() + + +if __name__ == '__main__': + ell = EllipsoidTriaxial.init_name("KarneyTest2024") + + 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, 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): + points_ana.append(ell.ell2cart(beta, lamb)) + points_ana = np.array(points_ana) + + show_points(points_ana, points, P0, P1) + print(f"Differenz s: {s_app - s} m") + print(f"Differenz alpha0: {wu.rad2deg(alpha0_app - alpha0)}°") + print(f"Differenz alpha1: {wu.rad2deg(alpha1_app - alpha1)}°") diff --git a/GHA_triaxial/approx_louville.py b/GHA_triaxial/approx_louville.py deleted file mode 100644 index f97835c..0000000 --- a/GHA_triaxial/approx_louville.py +++ /dev/null @@ -1,108 +0,0 @@ -import numpy as np -from numpy import sin, cos, arcsin, arccos, arctan2 -from ellipsoide import EllipsoidTriaxial -import matplotlib.pyplot as plt - -dbeta_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ey**2 * sin(alpha)**2 * cos(beta) * sin(beta) -dlamb_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ee**2 * cos(alpha)**2 * sin(lamb) * cos(lamb) -dalpha_dc = lambda ell, beta, lamb, alpha: (2 * sin(alpha) * cos(alpha) * - (ell.Ey**2 * cos(beta)**2 + ell.Ee**2 * sin(lamb)**2)) - -lamb2_sphere = lambda r, phi1, lamb1, a12, s: lamb1 + arctan2(sin(s/r) * sin(a12), - cos(phi1) * cos(s/r) - sin(s/r) * cos(a12)) -phi2_sphere = lambda r, phi1, lamb1, a12, s: arcsin(sin(phi1) * cos(s/r) + cos(phi1) * sin(s/r) * cos(a12)) - -a12_sphere = lambda phi1, lamb1, phi2, lamb2: arctan2(cos(phi2) * sin(lamb2 - lamb1), - cos(phi1) * cos(phi2) - - sin(phi1) * cos(phi2) * cos(lamb2 - lamb1)) -s_sphere = lambda r, phi1, lamb1, phi2, lamb2: r * arccos(sin(phi1) * sin(phi2) + - cos(phi1) * cos(phi2) * cos(lamb2 - lamb1)) - -louville = lambda beta, lamb, alpha: (ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 - - ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2) - -def points_approx_gha2(r: float, phi1: np.ndarray, lamb1: np.ndarray, phi2: np.ndarray, lamb2: np.ndarray, num: int = None, step_size: float = 10000): - s_approx = s_sphere(r, phi1, lamb1, phi2, lamb2) - if num is not None: - step_size = s_approx / (num+1) - a_approx = a12_sphere(phi1, lamb1, phi2, lamb2) - - points = [np.array([phi1, lamb1, a_approx])] - current_s = step_size - while current_s < s_approx: - phi_n = phi2_sphere(r, phi1, lamb1, a_approx, current_s) - lamb_n = lamb2_sphere(r, phi1, lamb1, a_approx, current_s) - points.append(np.array([phi_n, lamb_n, a_approx])) - current_s += step_size - points.append(np.array([phi2, lamb2, a_approx])) - return points - - -def num_update(ell: EllipsoidTriaxial, points, diffs): - for i, (beta, lamb, alpha) in enumerate(points): - dalpha = dalpha_dc(ell, beta, lamb, alpha) - if i == 0 or i == len(points) - 1: - grad = np.array([0, 0, dalpha]) - else: - dbeta = dbeta_dc(ell, beta, lamb, alpha) - dlamb = dlamb_dc(ell, beta, lamb, alpha) - grad = np.array([dbeta, dlamb, dalpha]) - - delta = -diffs[i] * grad / np.dot(grad, grad) - points[i] += delta - return points - - -def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int): - beta1, lamb1 = ell.cart2ell(p1) - beta2, lamb2 = ell.cart2ell(p2) - points = points_approx_gha2(ell.ax, beta1, lamb1, beta2, lamb2, 5) - - for j in range(maxI): - constants = [louville(point[0], point[1], point[2]) for point in points] - mean_constant = np.mean(constants) - diffs = constants - mean_constant - if np.mean(np.abs(diffs)) > 10: - points = num_update(ell, points, diffs) - else: - break - for k in range(maxI): - - - last_diff_alpha = points[-2][-1] - points[-3][-1] - alpha_extrap = points[-2][-1] + last_diff_alpha - if abs(alpha_extrap - points[-1][-1]) > 0.0005: - pass - else: - break - pass - pass - return points - -def show_points(ell: EllipsoidTriaxial, points): - points_cart = [] - for point in points: - points_cart.append(ell.ell2cart(point[0], point[1])) - points_cart = np.array(points_cart) - - fig = plt.figure() - ax = fig.add_subplot(111, projection='3d') - - ax.plot(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2]) - ax.scatter(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2]) - - ax.set_xlabel('X') - ax.set_ylabel('Y') - ax.set_zlabel('Z') - - plt.show() - - -if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("Eitschberger1978") - p1 = np.array([4189000, 812000, 4735000]) - p2 = np.array([4090000, 868000, 4808000]) - p1, phi1, lamb1, h1 = ell.cartonell(p1) - p2, phi2, lamb2, h2 = ell.cartonell(p2) - points = gha2(ell, p1, p2, 10) - show_points(ell, points) diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 3579763..433b780 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -8,21 +8,22 @@ from math import comb import GHA_triaxial.numeric_examples_panou as ne_panou import GHA_triaxial.numeric_examples_karney as ne_karney from ellipsoide import EllipsoidTriaxial -from typing import Callable +from typing import Callable, Tuple, List +from numpy.typing import NDArray -def pq_ell(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]: +def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: """ Berechnung von p und q in elliptischen Koordinaten Panou, Korakitits 2019 - :param x: x - :param y: y - :param z: z + :param ell: Ellipsoid + :param point: Punkt :return: p und q """ - n = ell.func_n(x, y, z) + x, y, z = point + n = ell.func_n(point) - beta, lamb = ell.cart2ell(np.array([x, y, z])) + beta, lamb = ell.cart2ell(point) B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2 L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2 @@ -49,7 +50,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable: :param ell: Ellipsoid :return: DGL-System """ - def ODE(s: float, v: np.ndarray) -> np.ndarray: + def ODE(s: float, v: NDArray) -> NDArray: """ DGL-System :param s: unabhängige Variable @@ -58,7 +59,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable: """ x, dxds, y, dyds, z, dzds = v - H = ell.func_H(x, y, z) + H = ell.func_H(np.array([x, y, z])) h = dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2 ddx = -(h/H)*x @@ -68,7 +69,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, list]: +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: @@ -76,12 +77,14 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, :param alpha0: :param s: :param num: + :param all_points: :return: """ phi, lam, _ = ell.cart2geod(point, "ligas3") - x0, y0, z0 = ell.geod2cart(phi, lam, 0) + p0 = ell.geod2cart(phi, lam, 0) + x0, y0, z0 = p0 - p, q = pq_ell(ell, x0, y0, z0) + 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) @@ -93,7 +96,9 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, _, werte = rk.rk4(ode, 0, v_init, s, num) x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1] - p1, q1 = pq_ell(ell, x1, y1, z1) + point1 = np.array([x1, y1, z1]) + + p1, q1 = pq_ell(ell, point1) sigma = np.array([dx1ds, dy1ds, dz1ds]) P = float(p1 @ sigma) Q = float(q1 @ sigma) @@ -103,21 +108,23 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, if alpha1 < 0: alpha1 += 2 * np.pi - return np.array([x1, y1, z1]), alpha1, werte + if all_points: + return point1, alpha1, werte + else: + return point1, alpha1 # --------------------------------------------------------------------------------------------------------------------- -def pq_para(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]: +def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: """ Berechnung von p und q in parametrischen Koordinaten Panou, Korakitits 2020 - :param x: x - :param y: y - :param z: z + :param ell: Ellipsoid + :param point: Punkt :return: p und q """ - n = ell.func_n(x, y, z) - u, v = ell.cart2para(np.array([x, y, z])) + n = ell.func_n(point) + u, v = ell.cart2para(point) # 41-47 G = sqrt(1 - ell.ex ** 2 * cos(u) ** 2 - ell.ee ** 2 * sin(u) ** 2 * sin(v) ** 2) @@ -131,12 +138,12 @@ def pq_para(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.nd t1 = np.dot(n, q) t2 = np.dot(n, p) t3 = np.dot(p, q) - if not (t1<1e-10 or t1>1-1e-10) and not (t2<1e-10 or t2>1-1e-10) and not (t3<1e-10 or t3>1-1e-10): + if not (t1 < 1e-10 or t1 > 1-1e-10) and not (t2 < 1e-10 or t2 > 1-1e-10) and not (t3 < 1e-10 or t3 > 1-1e-10): raise Exception("Fehler in den normierten Vektoren") return p, q -def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): +def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]: """ Panou, Korakitits 2020, 5ff. :param ell: @@ -153,7 +160,7 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): y_m = [y] z_m = [z] - p, q = pq_para(ell, x, y, z) + p, q = pq_para(ell, point) # 48-50 x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0)) @@ -208,7 +215,8 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): for c in reversed(c_m): z_s = z_s * s + c - p_s, q_s = pq_para(ell, x_s, y_s, z_s) + p1 = np.array([x_s, y_s, z_s]) + p_s, q_s = pq_para(ell, p1) # 57-59 dx_s = 0 @@ -234,10 +242,10 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): if alpha1 < 0: alpha1 += 2 * np.pi - return np.array([x_s, y_s, z_s]), alpha1 + return p1, alpha1 -def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> tuple[np.ndarray, float]: +def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]: if s > np.pi / maxPartCircum * ell.ax: s /= 2 point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum) @@ -251,14 +259,14 @@ 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) +def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]: + point = ell.para2cart(u, v) beta, lamb = ell.para2ell(u, v) - p_para, q_para = pq_para(ell, x, y, z) + p_para, q_para = pq_para(ell, point) sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para) - p_ell, q_ell = pq_ell(ell, x, y, z) + p_ell, q_ell = pq_ell(ell, point) alpha_ell = arctan2(p_ell @ sigma_para, q_ell @ sigma_para) sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) @@ -267,49 +275,42 @@ def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float 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) +def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]: + point = ell.ell2cart(beta, lamb) u, v = ell.ell2para(beta, lamb) - p_ell, q_ell = pq_ell(ell, x, y, z) + p_ell, q_ell = pq_ell(ell, point) sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) - p_para, q_para = pq_para(ell, x, y, z) + p_para, q_para = pq_para(ell, point) 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) +def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha: float) -> NDArray: + p, q = pq_ell(ell, point) 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) +def func_sigma_para(ell: EllipsoidTriaxial, point: NDArray, alpha: float) -> NDArray: + p, q = pq_para(ell, point) sigma = p * sin(alpha) + q * cos(alpha) return sigma -def louville_constant(ell: EllipsoidTriaxial, p0: np.ndarray, alpha: float) -> float: +def louville_constant(ell: EllipsoidTriaxial, p0: 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): +def louville_l2c(ell: EllipsoidTriaxial, l: float) -> float: return sqrt((l + ell.Ee**2) / ell.Ex**2) -def louville_c2l(ell, c): +def louville_c2l(ell: EllipsoidTriaxial, c: float) -> float: return ell.Ex**2 * c**2 - ell.Ee**2 @@ -322,7 +323,7 @@ if __name__ == "__main__": # 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) + # 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) @@ -342,7 +343,7 @@ 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, 5000) beta1_num, lamb1_num = ell.cart2ell(P1_num) try: @@ -357,5 +358,3 @@ 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 diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index 713c394..694b300 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -4,9 +4,13 @@ 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 -def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30): +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 @@ -17,6 +21,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep :param n: Anzahl Schritte :param epsilon: :param iter_max: Maximale Anzhal Iterationen + :param all_points: :return: """ @@ -25,7 +30,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep def arccot(x): return np.arctan2(1.0, x) - 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) @@ -60,7 +64,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep 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, @@ -244,15 +247,21 @@ 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, beta_arr, lamb_arr + if all_points: + return alpha_1, alpha_2, s, beta_arr, lamb_arr + else: + return alpha_1, alpha_2, s - if lamb_1 == lamb_2: + else: # lamb_1 == lamb_2 N = n dbeta = beta_2 - beta_1 if abs(dbeta) < 10**-15: - return 0, 0, 0, np.array([]), np.array([]) + if all_points: + return 0, 0, 0, np.array([]), np.array([]) + else: + return 0, 0, 0 lamb_0 = 0 @@ -369,7 +378,10 @@ 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, beta_arr, lamb_arr + if all_points: + return alpha_1, alpha_2, s, beta_arr, lamb_arr + else: + return alpha_1, alpha_2, s if __name__ == "__main__": @@ -390,7 +402,6 @@ if __name__ == "__main__": # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000) # print(s) - # ell = EllipsoidTriaxial.init_name("BursaSima1980round") # diffs_panou = [] # examples_panou = ne_panou.get_random_examples(4) @@ -421,8 +432,4 @@ if __name__ == "__main__": # diffs_karney = np.array(diffs_karney) # print(diffs_karney) - pass - - - diff --git a/dashboard.py b/dashboard.py index bed22e7..68a2999 100644 --- a/dashboard.py +++ b/dashboard.py @@ -150,7 +150,6 @@ def figure_lines(fig, line, color): return fig - app.layout = html.Div( style={"fontFamily": "Arial", "padding": "5px", "width": "70%", "margin-left": "auto"}, children=[ @@ -170,7 +169,7 @@ app.layout = html.Div( {"label": "BesselBiaxial", "value": "BesselBiaxial"}, {"label": "Fiction", "value": "Fiction"}, {"label": "KarneyTest2024", "value": "KarneyTest2024"}, - #{"label": "Ei", "value": "Ei"}, + # {"label": "Ei", "value": "Ei"}, ], value="", style={"width": "300px", "marginBottom": "20px"}, @@ -254,7 +253,6 @@ def fill_inputs_from_dropdown(selected_ell): if not selected_ell: return None, None, None - ell = EllipsoidTriaxial.init_name(selected_ell) ax = ell.ax ay = ell.ay @@ -408,16 +406,16 @@ def calc_and_plot(n1, n2, alpha_rad = wu.deg2rad(float(a_deg)) s_val = float(s) - p1 = tuple(map(float, ell.ell2cart(beta_rad, lamb_rad))) + p1 = ell.ell2cart(beta_rad, lamb_rad) out1 = [] if "analytisch" in method1: # ana - (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]) + p2_ana, alpha2 = gha1_ana(ell, p1, alpha_rad, s_val, 70) + x2, y2, z2 = p2_ana + beta2, lamb2 = ell.cart2ell(p2_ana) - #out1 += f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)}," + # out1 += f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)}," out1.append( html.Div([ html.Strong("Analytisch: "), @@ -430,8 +428,7 @@ def calc_and_plot(n1, n2, if "numerisch" in method1: # num - (x1, y1, z1), alpha1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000) - p2_num = x1, y1, z1 + p2_num, alpha1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000, all_points=True) beta2_num, lamb2_num = ell.cart2ell(p2_num) out1.append( @@ -448,7 +445,6 @@ def calc_and_plot(n1, n2, for x1, _, y1, _, z1, _ in werte: geo_line_num1.append([x1, y1, z1]) - if "stochastisch" in method1: # stoch p2_stoch = "noch nicht implementiert.." @@ -464,15 +460,15 @@ def calc_and_plot(n1, n2, return html.Span("Bitte Berechnungsverfahren auswählen!", style={"color": "red"}), "", go.Figure() fig = ellipsoid_figure(ell, title="Erste Hauptaufgabe - analystisch") - #fig = figure_constant_lines(fig, ell, "geod") + # fig = figure_constant_lines(fig, ell, "geod") fig = figure_constant_lines(fig, ell, "ell") - #fig = figure_constant_lines(fig, ell, "para") + # fig = figure_constant_lines(fig, ell, "para") if "analytisch" in method1: fig = figure_points(fig, [("P1", p1, "black"), ("P2", p2_ana, "red")]) if "numerisch" in method1: fig = figure_lines(fig, geo_line_num1, "#ff8c00") - #out1 = f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)}, {p2_num}" + # out1 = f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)}, {p2_num}" return out1, "", fig if dash.ctx.triggered_id == "button-calc-gha2": @@ -488,15 +484,14 @@ def calc_and_plot(n1, n2, alpha_1, alpha_2, s12, beta_arr, lamb_arr = gha2_num( ell, np.deg2rad(float(beta21)), np.deg2rad(float(lamb21)), - np.deg2rad(float(beta22)), np.deg2rad(float(lamb22)) + np.deg2rad(float(beta22)), np.deg2rad(float(lamb22)), + all_points=True ) geo_line_num = [] for beta, lamb in zip(beta_arr, lamb_arr): point = ell.ell2cart(beta, lamb) geo_line_num.append(point) - - out2.append( html.Div([ html.Strong("Numerisch: "), @@ -504,7 +499,6 @@ def calc_and_plot(n1, n2, ]) ) - if "stochastisch" in method2: # stoch a_stoch = "noch nicht implementiert.." @@ -524,8 +518,6 @@ def calc_and_plot(n1, n2, if "numerisch" in method2: fig = figure_lines(fig, geo_line_num, "#ff8c00") fig = figure_points(fig, [("P1", p1, "black"), ("P2", p2, "red")]) - - return "", out2, fig diff --git a/ellipsoide.py b/ellipsoide.py index 6a79e92..76bbc24 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -1,9 +1,10 @@ import numpy as np from numpy import sin, cos, arctan, arctan2, sqrt, pi, arccos import winkelumrechnungen as wu -import ausgaben as aus import jacobian_Ligas import matplotlib.pyplot as plt +from typing import Tuple +from numpy.typing import NDArray class EllipsoidBiaxial: @@ -162,18 +163,40 @@ class EllipsoidTriaxial: b = 1 / sqrt(2) return cls(ax, ay, b) - def func_H(self, x: float, y: float, z: float): + def func_H(self, point: NDArray) -> float: + """ + Berechnung H + Panou, Korakitis 2019 [43] + :param point: Punkt + :return: H + """ + x, y, z = point return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 - def func_n(self, x: float, y: float, z: float, H: float = None): + def func_n(self, point: NDArray, H: float = None) -> NDArray: + """ + Berechnung normalen Vektor + Panou, Korakitis 2019 [9-12] + :param point: Punkt + :param H: + :return: + """ if H is None: - H = self.func_H(x, y, z) + H = self.func_H(point) sqrtH = sqrt(H) + x, y, z = point return np.array([x / sqrtH, y / ((1 - self.ee ** 2) * sqrtH), z / ((1 - self.ex ** 2) * sqrtH)]) - def func_t12(self, x, y, z): + def func_t12(self, point: NDArray) -> Tuple[float, float]: + """ + Berechnung Wurzeln + Panou, Korakitis 2019 [9-12] + :param point: Punkt + :return: Wurzeln t1, t2 + """ + x, y, z = point 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 - ( @@ -184,7 +207,7 @@ class EllipsoidTriaxial: t1 = c0 / t2 return t1, t2 - def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: + def ellu2cart(self, beta: float, lamb: float, u: float) -> NDArray: """ Panou 2014 12ff. Elliptische Breite+Länge sind nicht gleich der geodätischen @@ -201,7 +224,7 @@ class EllipsoidTriaxial: return np.array([x, y, z]) - def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]: + def cart2ellu(self, point: NDArray) -> Tuple[float, float, float]: """ Panou 2014 15ff. :param point: Punkt in kartesischen Koordinaten @@ -232,7 +255,7 @@ class EllipsoidTriaxial: return beta, lamb, u - def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: + def ell2cart(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray: """ Panou, Korakitis 2019 2 :param beta: elliptische Breite [rad] @@ -268,7 +291,7 @@ class EllipsoidTriaxial: return xyz - def ell2cart_bektas(self, beta: float | np.ndarray, omega: float | np.ndarray): + def ell2cart_bektas(self, beta: float | NDArray, omega: float | NDArray) -> NDArray: """ Bektas 2015 :param beta: @@ -281,14 +304,13 @@ class EllipsoidTriaxial: return np.array([x, y, z]) - def ell2cart_karney(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray: + def ell2cart_karney(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray: """ Karney 2025 Geographic Lib :param beta: :param lamb: :return: """ - e = sqrt(self.ax**2 - self.b**2) / self.ay k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2) k_ = sqrt(self.ax**2 - self.ay**2) / sqrt(self.ax**2 - self.b**2) X = self.ax * cos(lamb) * sqrt(k**2*cos(beta)**2+k_**2) @@ -296,11 +318,12 @@ class EllipsoidTriaxial: Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2) return np.array([X, Y, Z]) - def cart2ell(self, point, eps=1e-12, maxI = 100) -> tuple[float, float]: + def cart2ell(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]: """ Panou, Korakitis 2019 3f. (num) :param point: :param eps: + :param maxI: :return: """ x, y, z = point @@ -323,10 +346,10 @@ class EllipsoidTriaxial: (self.b * self.Ee ** 2) / (2 * self.Ex) * sin(beta) * sin(2 * lamb) / sqrt(L)]]) N = J.T @ J - det = N[0,0] * N[1,1] - N[0,1] * N[1,0] + 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]]]) + 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] @@ -344,7 +367,7 @@ class EllipsoidTriaxial: return beta, lamb - def cart2ell_panou(self, point: np.ndarray) -> tuple[float, float]: + def cart2ell_panou(self, point: NDArray) -> Tuple[float, float]: """ Panou, Korakitis 2019 2f. (analytisch -> Näherung) :param point: Punkt in kartesischen Koordinaten @@ -371,7 +394,7 @@ class EllipsoidTriaxial: # ---- Allgemeiner Fall ----- - t1, t2 = self.func_t12(x, y, z) + t1, t2 = self.func_t12(point) num_beta = max(t1 - self.b ** 2, 0) den_beta = max(self.ay ** 2 - t1, 0) @@ -401,7 +424,7 @@ class EllipsoidTriaxial: return beta, lamb - def cart2ell_bektas(self, point, eps=1e-12, maxI = 100) -> tuple[float, float]: + def cart2ell_bektas(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]: """ Bektas 2015 :param point: @@ -430,7 +453,7 @@ class EllipsoidTriaxial: return phi, lamb - def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray: + def geod2cart(self, phi: float | NDArray, lamb: float | NDArray, h: float) -> NDArray: """ Ligas 2012, 250 :param phi: geodätische Breite [rad] @@ -444,7 +467,7 @@ class EllipsoidTriaxial: zG = (v * (1-self.ex**2) + h) * sin(phi) return np.array([xG, yG, zG]) - def cart2geod(self, point: np.ndarray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + def cart2geod(self, point: NDArray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: """ Ligas 2012 :param mode: ligas1, ligas2, oder ligas3 @@ -513,7 +536,7 @@ class EllipsoidTriaxial: return phi, lamb, h - def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray: + def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray: """ Panou, Korakitits 2020, 4 :param u: Parameter u @@ -526,7 +549,7 @@ class EllipsoidTriaxial: z = np.broadcast_to(z, np.shape(x)) return np.array([x, y, z]) - def cart2para(self, point: np.ndarray) -> tuple[float, float]: + def cart2para(self, point: NDArray) -> Tuple[float, float]: """ Panou, Korakitits 2020, 4 :param point: Punkt in kartesischen Koordinaten @@ -551,31 +574,31 @@ class EllipsoidTriaxial: return u, v - def ell2para(self, beta: float, lamb: float) -> tuple[float, float]: + def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]: cart = self.ell2cart(beta, lamb) return self.cart2para(cart) - def para2ell(self, u: float, v: float) -> tuple[float, float]: + def para2ell(self, u: float, v: float) -> Tuple[float, float]: cart = self.para2cart(u, v) return self.cart2ell(cart) - def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: cart = self.para2cart(u, v) return self.cart2geod(cart, mode, maxIter, maxLoa) - def geod2para(self, phi: float, lamb: float, h: float) -> tuple[float, float]: + def geod2para(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: cart = self.geod2cart(phi, lamb, h) return self.cart2para(cart) - def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: + def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: cart = self.ell2cart(beta, lamb) return self.cart2geod(cart, mode, maxIter, maxLoa) - def geod2ell(self, phi: float, lamb: float, h: float) -> tuple[float, float]: + def geod2ell(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: cart = self.geod2cart(phi, lamb, h) return self.cart2ell(cart) - def point_on(self, point: np.ndarray) -> bool: + def point_on(self, point: NDArray) -> bool: """ Test, ob ein Punkt auf dem Ellipsoid liegt. :param point: kartesische 3D-Koordinaten @@ -587,17 +610,17 @@ class EllipsoidTriaxial: else: return False - def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]: + def cartonell(self, point: NDArray) -> NDArray: """ Berechnung des Lotpunktes auf einem Ellipsoiden :param point: Punkt in kartesischen Koordinaten, der gelotet werden soll - :return: Lotpunkt in kartesischen Koordinaten, geodätische Koordinaten des Punktes + :return: Lotpunkt in kartesischen Koordinaten """ phi, lamb, h = self.cart2geod(point, "ligas3") - x, y, z = self. geod2cart(phi, lamb, 0) - return np.array([x, y, z]), phi, lamb, h + p = self. geod2cart(phi, lamb, 0) + return p - def cartellh(self, point: np.ndarray, h: float) -> np.ndarray: + def cartellh(self, point: NDArray, h: float) -> NDArray: """ Punkt auf Ellipsoid hoch loten :param point: Punkt auf dem Ellipsoid @@ -635,7 +658,7 @@ if __name__ == "__main__": cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1) - diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3]) # + diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3]) diffs_ell.append([diff_ell]) diffs_para.append([diff_para]) diffs_geod.append([diff_geod3])