diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py index 71bca73..dcf3bea 100644 --- a/GHA_triaxial/panou.py +++ b/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 = [] - # H Ableitungen (7) - hsq = [] - # h Ableitungen (7) - hHst = [] - # h/H Ableitungen (6) - xsm = [x] - ysm = [y] - zsm = [z] + x_m = [x] + y_m = [y] + z_m = [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 = [] - bm = [] - cm = [] + 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)]) + + 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 - num = 10000 - # 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)) + s = 100 + 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) + 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) \ No newline at end of file + 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 diff --git a/ellipsoide.py b/ellipsoide.py index 98ee63c..cc93d6c 100644 --- a/ellipsoide.py +++ b/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 = 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_check1 = z*np.sqrt(1 - self.ee**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 = 2*np.arctan(y/(x*np.sqrt(1-self.ee**2) + np.sqrt(x**2*(1-self.ee**2)+y**2))) + v_check1 = y + 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__": diff --git a/test.py b/test.py index be76988..01cc87c 100644 --- a/test.py +++ b/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], @@ -21,4 +23,7 @@ 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 +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)) \ No newline at end of file