From 9031a123128e46a4c7ff177ec4f331136d9bca5b Mon Sep 17 00:00:00 2001 From: Hendrik Date: Tue, 25 Nov 2025 14:49:06 +0100 Subject: [PATCH] =?UTF-8?q?diverse=20=C3=84nderungen,=20Versuch=20der=20L?= =?UTF-8?q?=C3=B6sung=20der=20zweiten=20GHA=20mit=20Louville?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA/rk.py | 2 +- GHA_triaxial/approx_louville.py | 108 +++++++++++++ GHA_triaxial/panou.py | 53 ++++--- T.py | 1 - ellipsoide.py | 265 ++++++++++++++++++++------------ 5 files changed, 308 insertions(+), 121 deletions(-) create mode 100644 GHA_triaxial/approx_louville.py delete mode 100644 T.py diff --git a/GHA/rk.py b/GHA/rk.py index 6d3574c..4c8913e 100644 --- a/GHA/rk.py +++ b/GHA/rk.py @@ -3,7 +3,7 @@ from numpy import sin, cos, tan import winkelumrechnungen as wu from ellipsoide import EllipsoidBiaxial -def gha1(re, x0, y0, z0, A0, s, num): +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 diff --git a/GHA_triaxial/approx_louville.py b/GHA_triaxial/approx_louville.py new file mode 100644 index 0000000..41f4da1 --- /dev/null +++ b/GHA_triaxial/approx_louville.py @@ -0,0 +1,108 @@ +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.cart2ell2(p1) + beta2, lamb2 = ell.cart2ell2(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.ell2cart2(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 2a6925d..b3452fd 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -10,7 +10,9 @@ from math import comb # Panou, Korakitits 2019 -def gha1_num(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num): +def gha1_num(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num): + phi, lamb, h = ell.cart2geod("ligas3", point) + x, y, z = ell.geod2cart(phi, lamb, 0) values = ell.p_q(x, y, z) H = values["H"] p = values["p"] @@ -57,7 +59,7 @@ def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points): pass -def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): +def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM): """ Panou, Korakitits 2020, 5ff. :param ell: @@ -69,6 +71,7 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): :param maxM: :return: """ + x, y, z = point x_m = [x] y_m = [y] z_m = [z] @@ -78,7 +81,7 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): n = np.array([x / sqrtH, y / ((1-ell.ee**2) * sqrtH), z / ((1-ell.ex**2) * sqrtH)]) - u, v = ell.cart2para(x, y, z) + u, v = ell.cart2para(np.array([x, y, z])) G = np.sqrt(1 - ell.ex**2 * np.cos(u)**2 - ell.ee**2 * np.sin(u)**2 * np.sin(v)**2) q = np.array([-1/G * np.sin(u) * np.cos(v), -1/G * np.sqrt(1-ell.ee**2) * np.sin(u) * np.sin(v), @@ -139,23 +142,33 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): if __name__ == "__main__": - ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") + # ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") + ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980") ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial") re = ellipsoide.EllipsoidBiaxial.init_name("Bessel") - x0 = 5672455.1954766 - y0 = 2698193.7242382686 - z0 = 1103177.6450055107 - alpha0 = wu.gms2rad([20, 0, 0]) - s = 500000 - num = 100 - werteTri = gha1_num(ellbi, x0, y0, z0, alpha0, s, num) - print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8)) - print("Distanz Triaxial Numerisch", np.sqrt((x0-werteTri[-1][1])**2+(y0-werteTri[-1][3])**2+(z0-werteTri[-1][5])**2)) - # checkLiouville(ell, werteTri) - werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num) - print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8)) - print("Distanz Biaxial", np.sqrt((x0-werteBi[0])**2+(y0-werteBi[1])**2+(z0-werteBi[2])**2)) - werteAna = gha1_ana(ell, x0, y0, z0, alpha0, s, 7) - print(aus.xyz(werteAna[0], werteAna[1], werteAna[2], 8)) - print("Distanz Triaxial Analytisch", np.sqrt((x0-werteAna[0])**2+(y0-werteAna[1])**2+(z0-werteAna[2])**2)) \ No newline at end of file + # Panou 2013, 7, Table 1, beta0=60° + beta1 = wu.deg2rad(60) + 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)) + 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) + + c = 0.06207487624 + alpha0 = wu.gms2rad([2, 52, 26.2393]) + alpha1 = wu.gms2rad([177, 4, 13.6373]) + s = 6705715.1610 + pass + + P2_num = gha1_num(ell, P1, alpha0, s, 1000) + P2_ana = gha1_ana(ell, P1, alpha0, s, 70) + pass \ No newline at end of file diff --git a/T.py b/T.py deleted file mode 100644 index cdef671..0000000 --- a/T.py +++ /dev/null @@ -1 +0,0 @@ -print("T") diff --git a/ellipsoide.py b/ellipsoide.py index a58e444..dd3f451 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -2,7 +2,6 @@ import numpy as np import winkelumrechnungen as wu import ausgaben as aus import jacobian_Ligas -from GHA_triaxial.panou import p_q class EllipsoidBiaxial: @@ -110,6 +109,10 @@ class EllipsoidTriaxial: @classmethod def init_name(cls, name: str): + """ + Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972, Bursa1970, BesselBiaxial + :param name: Name des dreiachsigen Ellipsoids + """ if name == "BursaFialova1993": ax = 6378171.36 ay = 6378101.61 @@ -117,8 +120,8 @@ class EllipsoidTriaxial: return cls(ax, ay, b) elif name == "BursaSima1980": ax = 6378172 - ay = 6378102.7 - b = 6356752.6 + ay = 6378103 + b = 6356753 return cls(ax, ay, b) elif name == "Eitschberger1978": ax = 6378173.435 @@ -135,61 +138,79 @@ class EllipsoidTriaxial: ay = 6378105 b = 6356754 return cls(ax, ay, b) - elif name == "Bessel-biaxial": + elif name == "BesselBiaxial": ax = 6377397.15509 ay = 6377397.15508 b = 6356078.96290 return cls(ax, ay, b) + elif name == "Fiction": + ax = 5500000 + ay = 4500000 + b = 4000000 + return cls(ax, ay, b) - def ell2cart(self, beta, lamb, u): + 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 ell2cart(self, beta: float, lamb: float, u: float) -> np.ndarray: """ Panou 2014 12ff. - :param beta: ellipsoidische Breite - :param lamb: ellipsoidische Länge - :param u: Höhe - :return: kartesische Koordinaten + Ellipsoidische Breite+Länge sind nicht gleich der geodätischen + Verhältnisse des Ellipsoids bekannt, Größe verändern bis Punkt erreicht, + dann ist u die Größe entlang der z-Achse + :param beta: ellipsoidische Breite [rad] + :param lamb: ellipsoidische Länge [rad] + :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 + # 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.ax**2))) + + # 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) - return x, y, z - - def ell2cart_P(self, lamb, beta): + return np.array([x, y, z]) + def ell2cart2(self, beta: float, lamb: float) -> np.ndarray: """ - Panou 2013, 241 - - :param lamb: ellipsoidische Länge - :param beta: reduzierte ellipsoidische Breite - :return: kartesische Koordinaten x y z + Panou, Korakitis 2019 2 + :param beta: ellipsoidische Breite [rad] + :param lamb: ellipsoidische Länge [rad] + :return: Punkt in kartesischen Koordinaten """ - - x = self.ax * np.sqrt(np.cos(beta) ** 2 + self.Ee ** 2 / self.Ex ** 2 * np.sin(beta) ** 2) * np.cos(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 * np.sin(beta) * np.sqrt(1 - self.Ee ** 2 / self.Ex ** 2 * np.cos(lamb) ** 2) + z = self.b / self.Ex * np.sin(beta) * np.sqrt(L) + return np.array([x, y, z]) - return x, y, z - - def cart2ell(self, x, y, z): + def cart2ell(self, point: np.ndarray) -> tuple[float, float, float]: """ Panou 2014 15ff. - :param x: - :param y: - :param z: - :return: + :param point: Punkt in kartesischen Koordinaten + :return: ellipsoidische Breite, ellipsoidische Länge, Größe entlang der z-Achse """ + x, y, z = point c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2 c1 = (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) @@ -211,18 +232,34 @@ class EllipsoidTriaxial: return beta, lamb, u - def cart2geod(self, mode: str, xG, yG, zG, maxIter=30, maxLoa=0.005): + def cart2ell2(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 + 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))) + return beta, lamb + + def cart2geod(self, mode: str, point: np.ndarray, maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]: """ Ligas 2012 - :param mode: - :param xG: - :param yG: - :param zG: - :param maxIter: - :param maxLoa: - :return: + :param mode: ligas1, ligas2, oder ligas3 + :param point: Punkt in kartesischen Koordinaten + :param maxIter: maximale Anzahl Iterationen + :param maxLoa: Level of Accuracy, das erreicht werden soll + :return: phi, lambda, h """ - rG = np.sqrt(xG**2 + yG**2 + zG**2) + xG, yG, zG = point + 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) E = 1 / self.ax**2 @@ -247,65 +284,93 @@ class EllipsoidTriaxial: 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) + h = np.sign(zG - pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2) return phi, lamb, h - def geod2cart(self, phi, lamb, h): + def geod2cart(self, phi: float, lamb: float, h: float) -> np.ndarray: """ Ligas 2012, 250 - :param phi: - :param lamb: - :param h: - :return: + :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 / 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 xG, yG, zG + return np.array([xG, yG, zG]) - def para2cart(self, u, v): + def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]: + """ + 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 + """ + phi, lamb, h = self.cart2geod("ligas3", point) + x, y, z = self. geod2cart(phi, lamb, 0) + return np.array([x, y, z]), phi, lamb, h + + def cartellh(self, point: np.ndarray, h: float) -> np.ndarray: + """ + Punkt auf Ellipsoid hoch loten + :param point: Punkt auf dem Ellipsoid + :param h: Höhe über dem Ellipsoid + :return: hochgeloteter Punkt + """ + phi, lamb, _ = self.cart2geod("ligas3", point) + pointH = self. geod2cart(phi, lamb, h) + return pointH + + def para2cart(self, u: float, v: float) -> np.ndarray: """ Panou, Korakitits 2020, 4 - :param u: - :param v: - :return: + :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.cos(v) + y = self.ay * np.cos(u) * np.sin(v) z = self.b * np.sin(u) + return np.array([x, y, z]) - def cart2para(self, x, y, z): + def cart2para(self, point: np.ndarray) -> tuple[float, float]: """ Panou, Korakitits 2020, 4 - :param x: - :param y: - :param z: - :return: + :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.arctan(u_check1 / u_check2) + u = np.arctan2(u_check1, u_check2) else: - u = np.pi/2 * np.arctan(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) + v_factor = np.sqrt(x**2*(1-self.ee**2)+y**2) if v_check1 <= v_check2: - v = 2 * np.arctan(v_check1 / (v_check2 + np.sqrt(x**2*(1-self.ee**2)+y**2))) + v = 2 * np.arctan2(v_check1, v_check2 + v_factor) else: - v = np.pi/2 - 2 * np.arctan(v_check2 / (v_check1 + np.sqrt(x**2*(1-self.ee**2)+y**2))) + v = np.pi/2 - 2 * np.arctan2(v_check2, v_check1 + v_factor) return u, v - def p_q(self, x, y, z): + def p_q(self, x, y, z) -> dict: + """ + Berechnung sämtlicher Größen + :param x: x + :param y: y + :param z: z + :return: Dictionary sämtlicher Größen + """ H = x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2 - n = np.array([x / np.sqrt(H), y / ((1 - self.ee ** 2) * np.sqrt(H)), z / ((1 - ell.ex ** 2) * np.sqrt(H))]) + 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(x, y, z) - carts = self.ell2cart(beta, lamb, u) + beta, lamb, u = 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 @@ -335,38 +400,40 @@ class EllipsoidTriaxial: "t2": t2, "F": F, "p": p, "q": q} -def louvilleConstant(ell, x, y, z): - values = p_q(ell, x, y, z) - pass - - if __name__ == "__main__": - ellips = EllipsoidTriaxial.init_name("Eitschberger1978") + ell = EllipsoidTriaxial.init_name("Eitschberger1978") - # carts = ellips.ell2cart(wu.deg2rad(10), wu.deg2rad(30), 6378172) - # ells = ellips.cart2ell(carts[0], carts[1], carts[2]) - # print(aus.gms("beta", ells[0], 3), aus.gms("lambda", ells[1], 3), "u =", ells[2]) - # ells2 = ellips.cart2ell(5712200, 2663400, 1106000) - # carts2 = ellips.ell2cart(ells2[0], ells2[1], ells2[2]) - # print(aus.xyz(carts2[0], carts2[1], carts2[2], 10)) + 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]) - # stellen = 20 - # geod1 = ellips.cart2geod("ligas1", 5712200, 2663400, 1106000) - # print(aus.gms("phi", geod1[0], stellen), aus.gms("lambda", geod1[1], stellen), "h =", geod1[2]) - # geod2 = ellips.cart2geod("ligas2", 5712200, 2663400, 1106000) - # print(aus.gms("phi", geod2[0], stellen), aus.gms("lambda", geod2[1], stellen), "h =", geod2[2]) - # geod3 = ellips.cart2geod("ligas3", 5712200, 2663400, 1106000) - # print(aus.gms("phi", geod3[0], stellen), aus.gms("lambda", geod3[1], stellen), "h =", geod3[2]) - # cart1 = ellips.geod2cart(geod1[0], geod1[1], geod1[2]) - # print(aus.xyz(cart1[0], cart1[1], cart1[2], 10)) - # cart2 = ellips.geod2cart(geod2[0], geod2[1], geod2[2]) - # print(aus.xyz(cart2[0], cart2[1], cart2[2], 10)) - # cart3 = ellips.geod2cart(geod3[0], geod3[1], geod3[2]) - # print(aus.xyz(cart3[0], cart3[1], cart3[2], 10)) - # test_cart = ellips.geod2cart(0.175, 0.444, 100) - # print(aus.xyz(test_cart[0], test_cart[1], test_cart[2], 10)) - - ellips.louvilleConstant(5712200, 2663400, 1106000) pass + + # 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") + + 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])) + + e_x, e_y, e_z = ell.ell2cart(beta, lamb_ell, u) + + 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) + pass \ No newline at end of file