From 6cc7245b0ff1562c32c179d7f07c9bd6209d031f Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sun, 11 Jan 2026 16:05:15 +0100 Subject: [PATCH] =?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])