analytisch funktioniert, p_q in ellisoid

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,8 +97,8 @@ 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.ee ** 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.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) *
@@ -146,24 +113,25 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM):
     a_m = []
     b_m = []
     c_m = []
-    for m in range(2, maxM+1):
+    for m in range(0, maxM+1):
+        if m >= 2:
             hH_t.append(hH_(m-2))
             x_m.append(x_(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))
+        a_m.append(x_m[m] / fact(m))
+        b_m.append(y_m[m] / fact(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))
+    print("Distanz Triaxial Numerisch", np.sqrt((x0-werteTri[-1][1])**2+(y0-werteTri[-1][3])**2+(z0-werteTri[-1][5])**2))
-    checkLiouville(ell, werteTri)
+    # 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))
+    print("Distanz Triaxial Analytisch", np.sqrt((x0-werteAna[0])**2+(y0-werteAna[1])**2+(z0-werteAna[2])**2))

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):
+    def p_q(self, x, y, z):
-        return x**2 + y**2/(1-self.ee**2)**2 + z**2/(1-self.ex**2)**2
+        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)
+    # carts = ellips.ell2cart(wu.deg2rad(10), wu.deg2rad(30), 6378172)
-    ells = ellips.cart2ell(carts[0], carts[1], carts[2])
+    # 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])
+    # print(aus.gms("beta", ells[0], 3), aus.gms("lambda", ells[1], 3), "u =", ells[2])
-    ells2 = ellips.cart2ell(5712200, 2663400, 1106000)
+    # ells2 = ellips.cart2ell(5712200, 2663400, 1106000)
-    carts2 = ellips.ell2cart(ells2[0], ells2[1], ells2[2])
+    # carts2 = ellips.ell2cart(ells2[0], ells2[1], ells2[2])
-    print(aus.xyz(carts2[0], carts2[1], carts2[2], 10))
+    # 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