From b7d437e0ca4be2190d13c028ea200c618042865d Mon Sep 17 00:00:00 2001 From: Hendrik Date: Wed, 22 Oct 2025 18:07:17 +0200 Subject: [PATCH] Umrechnung geod cart --- ausgaben.py | 4 +-- ellipsoide.py | 72 ++++++++++++++++++++++++++++++++++++++++++++--- jacobian_Ligas.py | 69 +++++++++++++++++++++++++++++++++++++++++++++ test.py | 24 ++++++++++++++++ 4 files changed, 162 insertions(+), 7 deletions(-) create mode 100644 jacobian_Ligas.py create mode 100644 test.py diff --git a/ausgaben.py b/ausgaben.py index 8c3e55f..05a328c 100644 --- a/ausgaben.py +++ b/ausgaben.py @@ -10,9 +10,7 @@ def xyz(x: float, y: float, z: float, stellen: int) -> str: :param stellen: Anzahl Nachkommastellen :return: String zur Ausgabe der Koordinaten """ - return f"""x = {(round(x,stellen))} m -y = {(round(y,stellen))} m -z = {(round(z,stellen))} m""" + return f"""x = {(round(x,stellen))} m y = {(round(y,stellen))} m z = {(round(z,stellen))} m""" def gms(name: str, rad: float, stellen: int) -> str: diff --git a/ellipsoide.py b/ellipsoide.py index e77c434..ad732bb 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -1,6 +1,7 @@ import numpy as np import winkelumrechnungen as wu import ausgaben as aus +import jacobian_Ligas class EllipsoidBiaxial: @@ -87,9 +88,12 @@ class EllipsoidTriaxial: self.ax = ax self.ay = ay self.b = b - 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 = 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) @classmethod def init_name(cls, name: str): @@ -170,6 +174,53 @@ class EllipsoidTriaxial: return beta, lamb, u + def cart2geod(self, mode: str, xG, yG, zG, maxIter=30, maxLoa=0.005): + """ + Ligas 2012 + :param mode: + :param xG: + :param yG: + :param zG: + :param maxIter: + :param maxLoa: + :return: + """ + 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 + F = 1 / self.ay**2 + G = 1 / self.b**2 + + i = 0 + loa = np.inf + + while i < maxIter and loa > maxLoa: + if mode == "ligas1": + invJ, fxE = jacobian_Ligas.case1(E, F, G, np.array([xG, yG, zG]), pE) + elif mode == "ligas2": + invJ, fxE = jacobian_Ligas.case2(E, F, G, np.array([xG, yG, zG]), pE) + elif mode == "ligas3": + 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) + 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) + + return phi, lamb, h + + def geod2cart(self, phi, lamb, h): + 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 + if __name__ == "__main__": ellips = EllipsoidTriaxial.init_name("Eitschberger1978") @@ -177,6 +228,19 @@ if __name__ == "__main__": # carts = ellips.ell2cart(10, 30, 6378172) # ells = ellips.cart2ell(carts[0], carts[1], carts[2]) - carts = ellips.ell2cart(90, 0, 6356754.4) + # carts = ellips.ell2cart(10, 25, 6378293.435) # print(aus.gms("beta", ells[0], 3), aus.gms("lambda", ells[1], 3), "u =", ells[2]) + 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", 5712216.95426783, 2663487.024865021, 1106098.8415910944) + print(aus.gms("phi", geod2[0], stellen), aus.gms("lambda", geod2[1], stellen), "h =", geod2[2]) + geod3 = ellips.cart2geod("ligas3", 5712216.95426783, 2663487.024865021, 1106098.8415910944) + 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)) pass diff --git a/jacobian_Ligas.py b/jacobian_Ligas.py new file mode 100644 index 0000000..293b660 --- /dev/null +++ b/jacobian_Ligas.py @@ -0,0 +1,69 @@ +import numpy as np + +def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): + j11 = 2 * E * pE[0] + j12 = 2 * F * pE[1] + j13 = 2 * G * pE[2] + j21 = F * pE[1] - (pE[1] - pG[1]) * E + j22 = (pE[0] - pG[0]) * F - E * pE[0] + j23 = 0 + j31 = G * pE[2] - (pE[2] - pG[2]) * E + j32 = 0 + j33 = (pE[0] - pG[0]) * G - E * pE[0] + + detJ = j11 * j22 * j33 - j21 * j12 * j33 - j31 * j13 * j22 + + invJ = 1/detJ * np.array([[j22*j33, -j12*j33, -j13*j22], + [-j21*j33, j11*j33-j13*j31, j13*j21], + [-j22*j31, j12*j31, 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[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0]]) + + return invJ, fxE + +def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): + j11 = 2 * E * pE[0] + j12 = 2 * F * pE[1] + j13 = 2 * G * pE[2] + j21 = F * pE[1] - (pE[1] - pG[1]) * E + j22 = (pE[0] - pG[0]) * F - E * pE[0] + j23 = 0 + j31 = 0 + j32 = G * pE[2] - (pE[2] - pG[2]) * F + 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]]) + + 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 + +def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): + j11 = 2 * E * pE[0] + j12 = 2 * F * pE[1] + j13 = 2 * G * pE[2] + j21 = G * pE[2] - (pE[2] - pG[2]) * E + j22 = 0 + j23 = (pE[0] - pG[0]) * G - E * pE[0] + j31 = 0 + j32 = G * pE[2] - (pE[2] - pG[2]) * F + j33 = (pE[1] - pG[1]) * G - F * pE[1] + + detJ = -j11 * j23 * j32 - j21 * j12 * j33 + j21 * j13 * j32 + + 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]]) + + return invJ, fxE diff --git a/test.py b/test.py new file mode 100644 index 0000000..be76988 --- /dev/null +++ b/test.py @@ -0,0 +1,24 @@ +import numpy as np + +J = np.array([ + [2, 3, 0], + [0, 3, 0], + [6, 0, 4] +]) + +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) + +# 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) \ No newline at end of file