From d76859d17b966bb7f463134a2b4fb321b4db8b6f Mon Sep 17 00:00:00 2001 From: Hendrik Date: Wed, 26 Nov 2025 11:05:18 +0100 Subject: [PATCH] =?UTF-8?q?Koordinatenumrechnungen=20funktionieren=20inkl.?= =?UTF-8?q?=20Randf=C3=A4lle,=20GHA2=5Fnum=20funktioniert=20mit=20Standard?= =?UTF-8?q?=20Ellipsoid?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA_triaxial/approx_louville.py | 6 +- GHA_triaxial/panou.py | 12 +- GHA_triaxial/panou_2013_2GHA_num.py | 67 +--------- ellipsoide.py | 190 ++++++++++++++++++++++------ jacobian_Ligas.py | 31 +++-- 5 files changed, 181 insertions(+), 125 deletions(-) diff --git a/GHA_triaxial/approx_louville.py b/GHA_triaxial/approx_louville.py index 41f4da1..f97835c 100644 --- a/GHA_triaxial/approx_louville.py +++ b/GHA_triaxial/approx_louville.py @@ -54,8 +54,8 @@ def num_update(ell: EllipsoidTriaxial, points, diffs): def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int): - beta1, lamb1 = ell.cart2ell2(p1) - beta2, lamb2 = ell.cart2ell2(p2) + 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): @@ -82,7 +82,7 @@ def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int): def show_points(ell: EllipsoidTriaxial, points): points_cart = [] for point in points: - points_cart.append(ell.ell2cart2(point[0], point[1])) + points_cart.append(ell.ell2cart(point[0], point[1])) points_cart = np.array(points_cart) fig = plt.figure() diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index b3452fd..7d3b9c7 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -152,16 +152,16 @@ if __name__ == "__main__": lamb1 = wu.deg2rad(0) beta2 = wu.deg2rad(60) lamb2 = wu.deg2rad(175) - P1 = ell.ell2cart2(wu.deg2rad(60), wu.deg2rad(0)) - P2 = ell.ell2cart2(wu.deg2rad(60), wu.deg2rad(175)) + P1 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(0)) + P2 = ell.ell2cart(wu.deg2rad(60), wu.deg2rad(175)) para1 = ell.cart2para(P1) para2 = ell.cart2para(P2) cart1 = ell.para2cart(para1[0], para1[1]) cart2 = ell.para2cart(para2[0], para2[1]) - ell11 = ell.cart2ell2(P1) - ell21 = ell.cart2ell2(P2) - ell1 = ell.cart2ell2(cart1) - ell2 = ell.cart2ell2(cart2) + ell11 = ell.cart2ell(P1) + ell21 = ell.cart2ell(P2) + ell1 = ell.cart2ell(cart1) + ell2 = ell.cart2ell(cart2) c = 0.06207487624 alpha0 = wu.gms2rad([2, 52, 26.2393]) diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/panou_2013_2GHA_num.py index e26e264..680e335 100644 --- a/GHA_triaxial/panou_2013_2GHA_num.py +++ b/GHA_triaxial/panou_2013_2GHA_num.py @@ -1,68 +1,9 @@ import numpy as np -#from ellipsoide import EllipsoidTriaxial +from ellipsoide import EllipsoidTriaxial import Numerische_Integration.num_int_runge_kutta as rk - +import ausgaben as aus # Panou 2013 - -class EllipsoidTriaxial: - def __init__(self, ax: float, ay: float, b: float): - 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) - - @classmethod - def init_name(cls, name: str): - if name == "BursaFialova1993": - ax = 6378171.36 - ay = 6378101.61 - b = 6356751.84 - return cls(ax, ay, b) - elif name == "BursaSima1980": - ax = 6378172 - ay = 6378102.7 - b = 6356752.6 - return cls(ax, ay, b) - elif name == "BursaSima1980round": - # Panou 2013 - ax = 6378172 - ay = 6378103 - b = 6356753 - return cls(ax, ay, b) - elif name == "Eitschberger1978": - ax = 6378173.435 - ay = 6378103.9 - b = 6356754.4 - return cls(ax, ay, b) - elif name == "Bursa1972": - ax = 6378173 - ay = 6378104 - b = 6356754 - return cls(ax, ay, b) - elif name == "Bursa1970": - ax = 6378173 - ay = 6378105 - b = 6356754 - return cls(ax, ay, b) - elif name == "Bessel-biaxial": - ax = 6377397.15509 - ay = 6377397.15508 - b = 6356078.96290 - return cls(ax, ay, b) - - - - - def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30): """ @@ -383,8 +324,8 @@ if __name__ == "__main__": beta2 = np.deg2rad(75) lamb2 = np.deg2rad(66) a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2) - print(np.rad2deg(a1)) - print(np.rad2deg(a2)) + print(aus.gms("a1", a1, 4)) + print(aus.gms("a2", a2, 4)) print(s) diff --git a/ellipsoide.py b/ellipsoide.py index dd3f451..606e912 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -119,6 +119,12 @@ class EllipsoidTriaxial: b = 6356751.84 return cls(ax, ay, b) elif name == "BursaSima1980": + ax = 6378172 + ay = 6378102.7 + b = 6356752.6 + return cls(ax, ay, b) + elif name == "BursaSima1980round": + # Panou 2013 ax = 6378172 ay = 6378103 b = 6356753 @@ -161,7 +167,7 @@ class EllipsoidTriaxial: else: return False - def ell2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: + def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: """ Panou 2014 12ff. Ellipsoidische Breite+Länge sind nicht gleich der geodätischen @@ -190,21 +196,34 @@ class EllipsoidTriaxial: return np.array([x, y, z]) - def ell2cart2(self, beta: float, lamb: float) -> np.ndarray: + def ell2cart(self, beta: float, lamb: float) -> np.ndarray: """ Panou, Korakitis 2019 2 :param beta: ellipsoidische Breite [rad] :param lamb: ellipsoidische Länge [rad] :return: Punkt in kartesischen Koordinaten """ - 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) - return np.array([x, y, z]) + if beta == -np.pi/2: + return np.array([0, 0, -self.b]) + elif beta == np.pi/2: + return np.array([0, 0, self.b]) + elif beta == 0 and lamb == -np.pi/2: + return np.array([0, -self.ay, 0]) + elif beta == 0 and lamb == np.pi/2: + return np.array([0, self.ay, 0]) + elif beta == 0 and lamb == 0: + return np.array([self.ax, 0, 0]) + elif beta == 0 and lamb == np.pi: + return np.array([-self.ax, 0, 0]) + else: + 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) + return np.array([x, y, z]) - def cart2ell(self, point: np.ndarray) -> tuple[float, float, float]: + def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]: """ Panou 2014 15ff. :param point: Punkt in kartesischen Koordinaten @@ -232,21 +251,65 @@ class EllipsoidTriaxial: return beta, lamb, u - def cart2ell2(self, point: np.ndarray) -> tuple[float, float]: + def cart2ell(self, point: np.ndarray) -> tuple[float, float]: """ Panou, Korakitis 2019 2f. :param point: Punkt in kartesischen Koordinaten :return: ellipsoidische Breite, ellipsoidische Länge """ x, y, z = point + + 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 + 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 + 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 + return beta, lamb + + # ---- 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 + np.sqrt(c1 ** 2 - 4 * c0)) / 2 t1 = c0 / t2 - beta = np.arctan(np.sqrt((t1 - self.b**2) / (self.ay**2 - t1))) - lamb = np.arctan(np.sqrt((t2 - self.ay**2) / (self.ax**2 - 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)) + + if z < 0: + beta = -beta + + if y < 0: + lamb = -lamb + + if x < 0: + lamb = np.pi - lamb + + if abs(x) < eps: + lamb = -np.pi/2 if y < 0 else np.pi/2 + + elif abs(y) < eps: + lamb = 0 if x > 0 else np.pi + + elif abs(z) < eps: + beta = 0 + return beta, lamb def cart2geod(self, mode: str, point: np.ndarray, maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: @@ -259,6 +322,45 @@ class EllipsoidTriaxial: :return: phi, lambda, h """ xG, yG, zG = point + + 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 + lamb = 0.0 + h = abs(zG) - ell.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 + h = abs(yG) - ell.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 + h = abs(xG) - ell.ax + return phi, lamb, h + + # elif abs(zG) < eps: # Punkt in der xy-Ebene + # phi = 0 + # lamb = np.arctan2(yG / ell.ay**2, xG / ell.ax**2) + # rG = np.sqrt(xG ** 2 + yG ** 2) + # pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) + # rE = np.sqrt(pE[0] ** 2 + pE[1] ** 2) + # h = rG - rE + # return phi, lamb, h + # + # elif abs(yG) < eps: # Punkt in der xz-Ebene + # phi = np.arctan2(zG / ell.b**2, xG / ell.ax**2) + # lamb = 0 if xG > 0 else np.pi + # rG = np.sqrt(xG ** 2 + zG ** 2) + # pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64) + # rE = np.sqrt(pE[0] ** 2 + pE[2] ** 2) + # h = rG - rE + # return phi, lamb, h + rG = np.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) @@ -286,6 +388,15 @@ class EllipsoidTriaxial: 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) + if xG < 0 and yG < 0: + lamb = -np.pi + lamb + + elif xG < 0: + lamb = np.pi + lamb + + if abs(zG) < eps: + phi = 0 + return phi, lamb, h def geod2cart(self, phi: float, lamb: float, h: float) -> np.ndarray: @@ -342,12 +453,13 @@ class EllipsoidTriaxial: :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) + u = np.pi/2 - np.arctan2(u_check2, u_check1) v_check1 = y v_check2 = x*np.sqrt(1-self.ee**2) @@ -356,6 +468,7 @@ class EllipsoidTriaxial: 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 p_q(self, x, y, z) -> dict: @@ -370,7 +483,7 @@ class EllipsoidTriaxial: n = np.array([x / np.sqrt(H), y / ((1 - self.ee ** 2) * np.sqrt(H)), z / ((1 - self.ex ** 2) * np.sqrt(H))]) - beta, lamb, u = self.cart2ell(np.array([x, y, z])) + beta, lamb, u = self.cart2ellu(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 @@ -403,37 +516,34 @@ class EllipsoidTriaxial: if __name__ == "__main__": ell = EllipsoidTriaxial.init_name("Eitschberger1978") + 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) - point = np.array([5672455, 2698193, 1103177]) - pointell, _, _, _ = ell.cartonell(point) - values1 = ell.cart2ell(point) - values2 = ell.cart2ell2(pointell) - carts1 = ell.ell2cart(values2[0], values2[1], ell.b) - carts2 = ell.ell2cart2(values2[0], values2[1]) + elli = ell.cart2ell(point) + cart_elli = ell.ell2cart(elli[0], elli[1]) + diff_ell = np.sum(np.abs(point-cart_elli)) + para = ell.cart2para(point) + cart_para = ell.para2cart(para[0], para[1]) + diff_para = np.sum(np.abs(point-cart_para)) - pass + geod = ell.cart2geod("ligas1", point) + cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + diff_geod1 = np.sum(np.abs(point-cart_geod)) - # cart_x, cart_y, cart_z = np.array([5672455, 2698193, 1103177]) - # cart_x, cart_y, cart_z = np.array([3415031.1337320395, 3414993.9029805865, 588647.548260936]) - # cart_x, cart_y, cart_z = np.array([6378173.435, 0, 0]) - # cart_x, cart_y, cart_z = np.array([0, 6378103.9, 0]) - # cart_x, cart_y, cart_z = np.array([0, 0, 6356754.4]) - cart_x, cart_y, cart_z = np.array([1000, 1000, 1000]) - print(aus.xyz(cart_x, cart_y, cart_z, 10), " Startwerte") + geod = ell.cart2geod("ligas2", point) + cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + diff_geod2 = np.sum(np.abs(point-cart_geod)) - beta, lamb_ell, u = ell.cart2ell(np.array([cart_x, cart_y, cart_z])) - phi, lamb_geod, h = ell.cart2geod("ligas3", np.array([cart_x, cart_y, cart_z])) - u_para, v_para = ell.cart2para(np.array([cart_x, cart_y, cart_z])) + geod = ell.cart2geod("ligas3", point) + cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + diff_geod3 = np.sum(np.abs(point-cart_geod)) - e_x, e_y, e_z = ell.ell2cart(beta, lamb_ell, u) + diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod1, diff_geod2, diff_geod3]) - print(aus.xyz(e_x, e_y, e_z, 10), f", Distanz ellipsoidisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([e_x, e_y, e_z]))} m") - g_x, g_y, g_z = ell.geod2cart(phi, lamb_geod, h) - print(aus.xyz(g_x, g_y, g_z, 10), f", Distanz geodätisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([g_x, g_y, g_z]))} m") - p_x, p_y, p_z = ell.para2cart(u_para, v_para) - print(aus.xyz(p_x, p_y, p_z, 10), f", Distanz parametrisch: {np.linalg.norm(np.array([cart_x, cart_y, cart_z]) - np.array([p_x, p_y, p_z]))} m") - pass - - carts = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(20), 6356754.4) + diff_list = np.array(diff_list) pass \ No newline at end of file diff --git a/jacobian_Ligas.py b/jacobian_Ligas.py index 293b660..7f49a5f 100644 --- a/jacobian_Ligas.py +++ b/jacobian_Ligas.py @@ -35,13 +35,16 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): j33 = (pE[1] - pG[1]) * G - F * pE[1] detJ = j11 * j22 * j33 - j21 * j12 * j33 + j21 * j13 * j32 - invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22], - [-j21*j33, j11*j33, j13*j21], - [j21*j32, -j11*j32, j11*j22-j12*j21]]) + if detJ == 0: + invJ, fxE = case3(E, F, G, pG, pE) + else: + invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22], + [-j21*j33, j11*j33, j13*j21], + [j21*j32, -j11*j32, j11*j22-j12*j21]]) - fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1, - (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0], - (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]]) + fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1, + (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0], + (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]]) return invJ, fxE @@ -57,13 +60,15 @@ def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): j33 = (pE[1] - pG[1]) * G - F * pE[1] detJ = -j11 * j23 * j32 - j21 * j12 * j33 + j21 * j13 * j32 + if detJ == 0: + invJ, fxE = case2(E, F, G, pG, pE) + else: + invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23], + [-j21*j33, j11*j33, -(j11*j23-j13*j21)], + [j21*j32, -j11*j32, -j12*j21]]) - invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23], - [-j21*j33, j11*j33, -(j11*j23-j13*j21)], - [j21*j32, -j11*j32, -j12*j21]]) - - fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1, - (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0], - (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]]) + fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1, + (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0], + (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]]) return invJ, fxE