analytisch funktioniert noch nicht

GHA_triaxial/panou.py

@@ -4,6 +4,8 @@ import Numerische_Integration.num_int_runge_kutta as rk
 import winkelumrechnungen as wu
 import ausgaben as aus
 import GHA.rk as ghark
+from scipy.special import factorial as fact
+from math import comb
 
 # Panou, Korakitits 2019
 
@@ -35,7 +37,9 @@ def p_q(ell: ellipsoide.EllipsoidTriaxial, x, y, z):
     # 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)
+    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}
@@ -98,30 +102,71 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, x, y, z, alpha0, s, maxM):
     :param maxM:
     :return:
     """
-    Hsp = []
+    x_m = [x]
-    # H Ableitungen (7)
+    y_m = [y]
-    hsq = []
+    z_m = [z]
-    # h Ableitungen (7)
+
-    hHst = []
-    # h/H Ableitungen (6)
-    xsm = [x]
-    ysm = [y]
-    zsm = [z]
     # erste Ableitungen (7-8)
+    sqrtH = np.sqrt(ell.H(x, y, z))
+    n = np.array([x / sqrtH,
+                  y / ((1-ell.ee**2) * sqrtH),
+                  z / ((1-ell.ex**2) * sqrtH)])
+    u, v = ell.cart2para(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),
+                  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]])
+    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))
+
+    # H Ableitungen (7)
+    H_ = lambda p: np.sum([comb(p, i) * (x_m[p - i] * x_m[i] +
+                                         1 / (1-ell.ee**2) ** 2 * y_m[p-i] * y_m[i] +
+                                         1 / (1-ell.ex**2) ** 2 * z_m[p-i] * z_m[i]) for i in range(0, p+1)])
+
+    # 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)])
+
+    # 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)]))
+
     # xm, ym, zm Ableitungen (6)
-    am = []
+    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)])
-    bm = []
+    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)])
-    cm = []
+    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))
+        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 am:
+    for a in a_m:
         x_s = x_s * s + a
     y_s = 0
-    for b in bm:
+    for b in b_m:
         y_s = y_s * s + b
     z_s = 0
-    for c in cm:
+    for c in c_m:
         z_s = z_s * s + c
+
+    return x_s, y_s, z_s
     pass
 
 
@@ -133,12 +178,15 @@ if __name__ == "__main__":
     y0 = 2698193.7242382686
     z0 = 1103177.6450055107
     alpha0 = wu.gms2rad([20, 0, 0])
-    s = 10000
+    s = 100
-    num = 10000
+    num = 100
-    # werteTri = gha1_num(ellbi, x0, y0, z0, alpha0, s, num)
+    werteTri = gha1_num(ellbi, x0, y0, z0, alpha0, s, num)
-    # print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8))
+    print(aus.xyz(werteTri[-1][1], werteTri[-1][3], werteTri[-1][5], 8))
-    # checkLiouville(ell, werteTri)
+    print(np.sqrt((x0-werteTri[-1][1])**2+(y0-werteTri[-1][3])**2+(z0-werteTri[-1][5])**2))
-    # werteBi = ghark.gha1(re, x0, y0, z0, alpha0, s, num)
+    checkLiouville(ell, werteTri)
-    # print(aus.xyz(werteBi[0], werteBi[1], werteBi[2], 8))
+    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))
 
-    gha1_ana(ell, x0, y0, z0, alpha0, s, 7)
+    werteAna = gha1_ana(ell, x0, y0, z0, alpha0, s, 7)
+    print(aus.xyz(werteAna[0], werteAna[1], werteAna[2], 8))

ellipsoide.py

@@ -267,15 +267,23 @@ class EllipsoidTriaxial:
         :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_check1 = z*np.sqrt(1 - self.ee**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))
+        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)
         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))
+            u = np.pi/2 * np.arctan(u_check2 / u_check1)
 
-        if y <= x*np.sqrt(1-self.ee**2):
+        v_check1 = y
-            v = 2*np.arctan(y/(x*np.sqrt(1-self.ee**2) + np.sqrt(x**2*(1-self.ee**2)+y**2)))
+        v_check2 = x*np.sqrt(1-self.ee**2)
+        if v_check1 <= v_check2:
+            v = 2 * np.arctan(v_check1 / (v_check2 + 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)))
+            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
 
 
 if __name__ == "__main__":

test.py

@@ -1,4 +1,6 @@
 import numpy as np
+from scipy.special import factorial as fact
+from math import comb
 
 J = np.array([
     [2, 3, 0],
@@ -22,3 +24,6 @@ print("Spaltenvektor:")
 print(res_col[0,0])
 print("Zeilenvektor:")
 print(res_row)
+t = 5
+l = 2
+print(fact(t+1-l) / (fact(t+1-l) * fact(l-1)), comb(l-1, t+1-l))