From 4e85eef5d7770c87fb569c910b588dec1166cded Mon Sep 17 00:00:00 2001 From: Hendrik Date: Tue, 4 Nov 2025 15:47:44 +0100 Subject: [PATCH] analytisch funktioniert, p_q in ellisoid --- GHA_triaxial/panou.py | 79 +++++++++++++------------------------------ ellipsoide.py | 57 ++++++++++++++++++++++++++----- 2 files changed, 73 insertions(+), 63 deletions(-) diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index dcf3bea..2a6925d 100644 --- a/GHA_triaxial/panou.py +++ b/GHA_triaxial/panou.py @@ -9,42 +9,9 @@ from math import comb # Panou, Korakitits 2019 -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))]) - - beta, lamb, u = ell.cart2ell(x, y, z) - 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) - 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) - # 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.array([n[1]*p[2]-n[2]*p[1], - n[2]*p[0]-n[0]*p[2], - n[1]*p[1]-n[1]*p[0]]) - - 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) + values = ell.p_q(x, y, z) H = values["H"] p = values["p"] q = values["q"] @@ -75,7 +42,7 @@ def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points): z = point[5] dzds = point[6] - values = p_q(ell, x, y, z) + values = ell.p_q(x, y, z) p = values["p"] q = values["q"] t1 = values["t1"] @@ -107,7 +74,7 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): z_m = [z] # erste Ableitungen (7-8) - sqrtH = np.sqrt(ell.H(x, y, z)) + sqrtH = np.sqrt(ell.p_q(x, y, z)["H"]) n = np.array([x / sqrtH, y / ((1-ell.ee**2) * sqrtH), z / ((1-ell.ex**2) * sqrtH)]) @@ -118,7 +85,7 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): 1/G * np.sqrt(1-ell.ex**2) * np.cos(u)]) p = np.array([q[1]*n[2] - q[2]*n[1], q[2]*n[0] - q[0]*n[2], - q[1]*n[1] - q[1]*n[0]]) + q[0]*n[1] - q[1]*n[0]]) x_m.append(p[0] * np.sin(alpha0) + q[0] * np.cos(alpha0)) y_m.append(p[1] * np.sin(alpha0) + q[1] * np.cos(alpha0)) z_m.append(p[2] * np.sin(alpha0) + q[2] * np.cos(alpha0)) @@ -130,40 +97,41 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM): # h Ableitungen (7) h_ = lambda q: np.sum([comb(q, j) * (x_m[q-j+1] * x_m[j+1] + - 1 / (1 - ell.ee ** 2) ** 2 * y_m[q-j+1] * y_m[j+1] + - 1 / (1 - ell.ex ** 2) ** 2 * z_m[q-j+1] * z_m[j+1]) for j in range(0, q + 1)]) + 1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] + + 1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1]) for j in range(0, q+1)]) # h/H Ableitungen (6) hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) * - np.sum([H_(t+1-l)/(fact(t+1-l)*fact(l-1))*hH_t[l-1] for l in range(1, t+1)])) + np.sum([H_(t+1-l) / (fact(t+1-l) * fact(l-1)) * hH_t[l-1] for l in range(1, t+1)])) # xm, ym, zm Ableitungen (6) x_ = lambda m: -np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)]) - y_ = lambda m: -1/(1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)]) - z_ = lambda m: -1/(1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)]) + y_ = lambda m: -1 / (1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)]) + z_ = lambda m: -1 / (1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)]) hH_t = [] a_m = [] b_m = [] c_m = [] - for m in range(2, maxM+1): - hH_t.append(hH_(m-2)) - x_m.append(x_(m)) + for m in range(0, maxM+1): + if m >= 2: + hH_t.append(hH_(m-2)) + x_m.append(x_(m)) + y_m.append(y_(m)) + z_m.append(z_(m)) a_m.append(x_m[m] / fact(m)) - y_m.append(y_(m)) b_m.append(y_m[m] / fact(m)) - z_m.append(z_(m)) c_m.append(z_m[m] / fact(m)) # am, bm, cm (6) x_s = 0 - for a in a_m: + for a in reversed(a_m): x_s = x_s * s + a y_s = 0 - for b in b_m: + for b in reversed(b_m): y_s = y_s * s + b z_s = 0 - for c in c_m: + for c in reversed(c_m): z_s = z_s * s + c return x_s, y_s, z_s @@ -178,15 +146,16 @@ if __name__ == "__main__": y0 = 2698193.7242382686 z0 = 1103177.6450055107 alpha0 = wu.gms2rad([20, 0, 0]) - s = 100 + 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(np.sqrt((x0-werteTri[-1][1])**2+(y0-werteTri[-1][3])**2+(z0-werteTri[-1][5])**2)) - checkLiouville(ell, werteTri) + 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(np.sqrt((x0-werteBi[0])**2+(y0-werteBi[1])**2+(z0-werteBi[2])**2)) + 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)) \ No newline at end of file + 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 diff --git a/ellipsoide.py b/ellipsoide.py index cc93d6c..1cf1a7d 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -2,6 +2,7 @@ import numpy as np import winkelumrechnungen as wu import ausgaben as aus import jacobian_Ligas +from GHA_triaxial.panou import p_q class EllipsoidBiaxial: @@ -282,19 +283,57 @@ class EllipsoidTriaxial: v = np.pi/2 - 2 * np.arctan(v_check2 / (v_check1 + np.sqrt(x**2*(1-self.ee**2)+y**2))) return u, v - def H(self, x, y, z): - return x**2 + y**2/(1-self.ee**2)**2 + z**2/(1-self.ex**2)**2 + def p_q(self, x, y, z): + 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))]) + + beta, lamb, u = self.cart2ell(x, y, z) + carts = self.ell2cart(beta, lamb, u) + 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 + + 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 + t2e = self.ax ** 2 * np.sin(lamb) ** 2 + self.ay ** 2 * np.cos(lamb) ** 2 + t1e = self.ay ** 2 * np.sin(beta) ** 2 + self.b ** 2 * np.cos(beta) ** 2 + + F = self.Ey ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(lamb) ** 2 + p1 = -np.sqrt(L / (F * t2)) * self.ax / self.Ex * np.sqrt(B) * np.sin(lamb) + p2 = np.sqrt(L / (F * t2)) * self.ay * np.cos(beta) * np.cos(lamb) + p3 = 1 / np.sqrt(F * t2) * (self.b * self.Ee ** 2) / (2 * self.Ex) * np.sin(beta) * np.sin(2 * lamb) + # p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * self.ax / (self.Ex * self.Ee) * np.sqrt(B) * np.sqrt(t2 - self.ay ** 2) + # p2 = np.sign(x) * np.sqrt(L / (F * t2)) * self.ay / (self.Ey * self.Ee) * np.sqrt((ell.ay ** 2 - t1) * (self.ax ** 2 - t2)) + # p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * self.b / (self.Ex * self.Ey) * np.sqrt( + # (t1 - self.b ** 2) * (t2 - self.ay ** 2) * (self.ax ** 2 - t2)) + p = np.array([p1, p2, p3]) + q = np.array([n[1] * p[2] - n[2] * p[1], + n[2] * p[0] - n[0] * p[2], + n[1] * p[1] - n[1] * p[0]]) + + 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 louvilleConstant(ell, x, y, z): + values = p_q(ell, x, y, z) + pass + if __name__ == "__main__": ellips = 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)) + # 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)) # stellen = 20 # geod1 = ellips.cart2geod("ligas1", 5712200, 2663400, 1106000) @@ -312,4 +351,6 @@ if __name__ == "__main__": # 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