From 610ea3dc28966faeed136fa2153ab07d8e5819df Mon Sep 17 00:00:00 2001 From: Hendrik Date: Sat, 1 Nov 2025 17:59:57 +0100 Subject: [PATCH] =?UTF-8?q?Aufger=C3=A4umt,=20Grundkonstrukt=20analytische?= =?UTF-8?q?=20L=C3=B6sung?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GHA_triaxial/panou.py | 114 ++++++++++++++++++++++++++++++++++++------ Richtungsreduktion.py | 14 ------ ellipsoide.py | 101 ++++++++++++++++++++++++++----------- hausarbeit.py | 78 ----------------------------- 4 files changed, 169 insertions(+), 138 deletions(-) delete mode 100644 Richtungsreduktion.py delete mode 100644 hausarbeit.py diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 2fe98de..71bca73 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -7,28 +7,44 @@ import GHA.rk as ghark # Panou, Korakitits 2019 -def gha1(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num): - H = x**2 + y**2 / (1-ell.ee**2)**2 + z**2/(1-ell.ex**2)**2 +def p_q(ell: ellipsoide.EllipsoidTriaxial, x, y, z): + H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2 - n = np.array([x/np.sqrt(H), y/((1-ell.ee**2)*np.sqrt(H)), z/((1-ell.ex**2)*np.sqrt(H))]) + n = np.array([x / np.sqrt(H), y / ((1 - ell.ee ** 2) * np.sqrt(H)), z / ((1 - ell.ex ** 2) * np.sqrt(H))]) beta, lamb, u = ell.cart2ell(x, y, z) - B = ell.Ex**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(beta)**2 - L = ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2 + carts = ell.ell2cart(beta, lamb, u) + B = ell.Ex ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(beta) ** 2 + L = ell.Ex ** 2 - ell.Ee ** 2 * np.cos(lamb) ** 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) + 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 + np.sqrt(c1**2 - 4*c0)) / 2 t1 = c0 / t2 + t2e = ell.ax ** 2 * np.sin(lamb) ** 2 + ell.ay ** 2 * np.cos(lamb) ** 2 + t1e = ell.ay ** 2 * np.sin(beta) ** 2 + ell.b ** 2 * np.cos(beta) ** 2 - F = ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2 - p1 = -np.sqrt(L/(F*t2)) * ell.ax/ell.Ex * np.sqrt(B) * np.sin(lamb) - p2 = np.sqrt(L/(F*t2)) * ell.ay * np.cos(beta) * np.cos(lamb) - p3 = 1 / np.sqrt(F*t2) * (ell.b*ell.Ee**2)/(2*ell.Ex) * np.sin(beta) * np.sin(2*lamb) + F = ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2 + p1 = -np.sqrt(L / (F * t2)) * ell.ax / ell.Ex * np.sqrt(B) * np.sin(lamb) + p2 = np.sqrt(L / (F * t2)) * ell.ay * np.cos(beta) * np.cos(lamb) + p3 = 1 / np.sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * np.sin(beta) * np.sin(2 * lamb) + # p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * ell.ax / (ell.Ex * ell.Ee) * np.sqrt(B) * np.sqrt(t2 - ell.ay ** 2) + # p2 = np.sign(x) * np.sqrt(L / (F * t2)) * ell.ay / (ell.Ey * ell.Ee) * np.sqrt((ell.ay ** 2 - t1) * (ell.ax ** 2 - t2)) + # p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * ell.b / (ell.Ex * ell.Ey) * np.sqrt( + # (t1 - ell.b ** 2) * (t2 - ell.ay ** 2) * (ell.ax ** 2 - t2)) p = np.array([p1, p2, p3]) q = np.cross(n, p) + return {"H": H, "n": n, "beta": beta, "lamb": lamb, "u": u, "B": B, "L": L, "c1": c1, "c0": c0, "t1": t1, "t2": t2, + "F": F, "p": p, "q": q} +def gha1_num(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num): + values = p_q(ell, x, y, z) + H = values["H"] + p = values["p"] + q = values["q"] + dxds0 = p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0) dyds0 = p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0) dzds0 = p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0) @@ -45,6 +61,69 @@ def gha1(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, num): funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [0, x, dxds0, y, dyds0, z, dzds0], s, num) return funktionswerte +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] + + values = p_q(ell, x, y, z) + p = values["p"] + q = values["q"] + t1 = values["t1"] + t2 = values["t2"] + + P = p[0]*dxds + p[1]*dyds + p[2]*dzds + Q = q[0]*dxds + q[1]*dyds + q[2]*dzds + alpha = np.arctan(P/Q) + + c = ell.ay**2 - (t1 * np.sin(alpha)**2 + t2 * np.cos(alpha)**2) + constantValues.append(c) + pass + + +def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): + """ + Panou, Korakitits 2020, 5ff. + :param ell: + :param x: + :param y: + :param z: + :param alpha0: + :param s: + :param maxM: + :return: + """ + Hsp = [] + # H Ableitungen (7) + hsq = [] + # h Ableitungen (7) + hHst = [] + # h/H Ableitungen (6) + xsm = [x] + ysm = [y] + zsm = [z] + # erste Ableitungen (7-8) + # xm, ym, zm Ableitungen (6) + am = [] + bm = [] + cm = [] + # am, bm, cm (6) + x_s = 0 + for a in am: + x_s = x_s * s + a + y_s = 0 + for b in bm: + y_s = y_s * s + b + z_s = 0 + for c in cm: + z_s = z_s * s + c + pass + if __name__ == "__main__": ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") @@ -56,7 +135,10 @@ if __name__ == "__main__": alpha0 = wu.gms2rad([20, 0, 0]) s = 10000 num = 10000 - werteTri = gha1(ellbi, x0, y0, z0, alpha0, s, num) - print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8)) - werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num) - print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8)) \ No newline at end of file + # werteTri = gha1_num(ellbi, x0, y0, z0, alpha0, s, num) + # print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8)) + # checkLiouville(ell, werteTri) + # werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num) + # print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8)) + + gha1_ana(ell, x0, y0, z0, alpha0, s, 7) \ No newline at end of file diff --git a/Richtungsreduktion.py b/Richtungsreduktion.py deleted file mode 100644 index 688f5eb..0000000 --- a/Richtungsreduktion.py +++ /dev/null @@ -1,14 +0,0 @@ -from numpy import sqrt, cos, sin -import ellipsoide -import ausgaben as aus -import winkelumrechnungen as wu - - -re = ellipsoide.EllipsoidBiaxial.init_name("Bessel") -eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2) - -A = wu.gms2rad([327,0,0]) -phi = wu.gms2rad([35,0,0]) -h = 3500 -dA = eta(phi)**2 * h * sin(A)*cos(A) / re.N(phi) -print(aus.gms("dA", dA, 10)) \ No newline at end of file diff --git a/ellipsoide.py b/ellipsoide.py index abc3220..98ee63c 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -135,7 +135,7 @@ class EllipsoidTriaxial: b = 6356754 return cls(ax, ay, b) elif name == "Bessel-biaxial": - ax = 6377397.155085 + ax = 6377397.15509 ay = 6377397.15508 b = 6356078.96290 return cls(ax, ay, b) @@ -148,17 +148,20 @@ class EllipsoidTriaxial: :param u: Höhe :return: kartesische Koordinaten """ - beta = wu.deg2rad(beta) - lamb = wu.deg2rad(lamb) - - # 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))) - 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 = 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 @@ -232,35 +235,73 @@ class EllipsoidTriaxial: return phi, lamb, h def geod2cart(self, phi, lamb, h): + """ + Ligas 2012, 250 + :param phi: + :param lamb: + :param h: + :return: + """ 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 + def para2cart(self, u, v): + """ + Panou, Korakitits 2020, 4 + :param u: + :param v: + :return: + """ + x = self.ax * np.cos(u) * np.cos(v) + y = self.ay * np.cos(u) * np.cos(v) + z = self.b * np.sin(u) + + def cart2para(self, x, y, z): + """ + Panou, Korakitits 2020, 4 + :param x: + :param y: + :param z: + :return: + """ + if z*np.sqrt(1-self.ee**2) <= np.sqrt(x**2 * (1-self.ee**2)+y**2) * np.sqrt(1-self.ex**2): + u = np.arctan(z*np.sqrt(1-self.ee**2) / np.sqrt(x**2 * (1-self.ee**2)+y**2) * np.sqrt(1-self.ex**2)) + else: + u = np.arctan(np.sqrt(x**2 * (1-self.ee**2)+y**2) * np.sqrt(1-self.ex**2) / z*np.sqrt(1-self.ee**2)) + + if y <= x*np.sqrt(1-self.ee**2): + v = 2*np.arctan(y/(x*np.sqrt(1-self.ee**2) + np.sqrt(x**2*(1-self.ee**2)+y**2))) + else: + v = np.pi/2 - 2*np.arctan(x*np.sqrt(1-self.ee**2) / (y + np.sqrt(x**2*(1-self.ee**2)+y**2))) + if __name__ == "__main__": ellips = EllipsoidTriaxial.init_name("Eitschberger1978") - # ellips = EllipsoidTriaxial.init_name("Bursa1972") - # carts = ellips.ell2cart(10, 30, 6378172) - # ells = ellips.cart2ell(carts[0], carts[1], carts[2]) - # 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", 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)) + 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)) - test_cart = ellips.geod2cart(0.175, 0.444, 100) - print(aus.xyz(test_cart[0], test_cart[1], test_cart[2], 10)) + # 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)) pass diff --git a/hausarbeit.py b/hausarbeit.py deleted file mode 100644 index 6d30043..0000000 --- a/hausarbeit.py +++ /dev/null @@ -1,78 +0,0 @@ -from numpy import cos, sin, tan -import winkelumrechnungen as wu -import s_ellipse as s_ell -import ausgaben as aus -from ellipsoide import EllipsoidBiaxial -import Numerische_Integration.num_int_runge_kutta as rk -import GHA.gauss as gauss -import GHA.bessel as bessel - -matrikelnummer = "6044051" -print(f"Matrikelnummer: {matrikelnummer}") - -m4 = matrikelnummer[-4] -m3 = matrikelnummer[-3] -m2 = matrikelnummer[-2] -m1 = matrikelnummer[-1] -print(f"m1={m1}\tm2={m2}\tm3={m3}\tm4={m4}") - -re = EllipsoidBiaxial.init_name("Bessel") -nks = 3 - -print(f"\na = {re.a} m\nb = {re.b} m\nc = {re.c} m\ne' = {re.e_}") - -print("\n\nAufgabe 1 (via Reihenentwicklung)") -phi1 = re.beta2phi(wu.gms2rad([int("1"+m1), int("1"+m2), int("1"+m3)])) -print(f"{aus.gms('phi_P1', phi1, 5)}") -s = s_ell.reihenentwicklung(re.e_, re.c, phi1) -print(f'\ns_P1 = {round(s,nks)} m') -# s_poly = s_ell.polyapp_tscheby_bessel(phi1) -# print(f'Meridianbogenlänge: {round(s_poly,nks)} m (Polynomapproximation)') - -print("\n\nAufgabe 2 (via Gauß´schen Mittelbreitenformeln)") - -lambda1 = wu.gms2rad([int("1"+m3), int("1"+m2), int("1"+m4)]) -A12 = wu.gms2rad([90, 0, 0]) -s12 = float("1"+m4+m3+m2+"."+m1+"0") - -# print("\nvia Runge-Kutta-Verfahren") -# re = Ellipsoid.init_name("Bessel") -# 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, phi1, lambda1, A12], -# s12, 1) -# -# for s, phi, lam, A in funktionswerte: -# print(f"{s} m, {aus.gms('phi', phi, nks)}, {aus.gms('lambda', lam, nks)}, {aus.gms('A', A, nks)}") - -# print("via Gauß´schen Mittelbreitenformeln") -phi_p2, lambda_p2, A_p2 = gauss.gha1(EllipsoidBiaxial.init_name("Bessel"), - phi_p1=phi1, - lambda_p1=lambda1, - A_p1=A12, - s=s12, eps=wu.gms2rad([1*10**-8, 0, 00])) - -print(f"\nP2: {aus.gms('phi', phi_p2[-1], nks)}\t\t{aus.gms('lambda', lambda_p2[-1], nks)}\t{aus.gms('A', A_p2[-1], nks)}") - -# print("\nvia Verfahren nach Bessel") -# phi_p2, lambda_p2, A_p2 = bessel.gha1(Ellipsoid.init_name("Bessel"), -# phi_p1=phi1, -# lambda_p1=lambda1, -# A_p1=A12, -# s=s12) -# print(f"P2: {aus.gms('phi', phi_p2, nks)}, {aus.gms('lambda', lambda_p2, nks)}, {aus.gms('A', A_p2, nks)}") - -print("\n\nAufgabe 3") -p = re.phi2p(phi_p2[-1]) -print(f"p = {round(p,2)} m") - - -print("\n\nAufgabe 4") -x = float("4308"+m3+"94.556") -y = float("1214"+m2+"88.242") -z = float("4529"+m4+"03.878") -phi_p3, lambda_p3, h_p3 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x, y, z) -print(f"\nP3: {aus.gms('phi', phi_p3, nks)}, {aus.gms('lambda', lambda_p3, 5)}, h = {round(h_p3,nks)} m")