import numpy as np import ellipsoide 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 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) 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) h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2 f1 = lambda s, x, dxds, y, dyds, z, dzds: dxds f2 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * x f3 = lambda s, x, dxds, y, dyds, z, dzds: dyds f4 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * y/(1-ell.ee**2) f5 = lambda s, x, dxds, y, dyds, z, dzds: dzds f6 = lambda s, x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / H * z/(1-ell.ex**2) 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: """ 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) 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 a_m: x_s = x_s * s + a y_s = 0 for b in b_m: y_s = y_s * s + b z_s = 0 for c in c_m: z_s = z_s * s + c return x_s, y_s, z_s pass if __name__ == "__main__": ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial") re = ellipsoide.EllipsoidBiaxial.init_name("Bessel") x0 = 5672455.1954766 y0 = 2698193.7242382686 z0 = 1103177.6450055107 alpha0 = wu.gms2rad([20, 0, 0]) 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)) werteAna = gha1_ana(ell, x0, y0, z0, alpha0, s, 7) print(aus.xyz(werteAna[0], werteAna[1], werteAna[2], 8))