import numpy as np from numpy import sin, cos, sqrt, arctan2 import ellipsoide import runge_kutta as rk import winkelumrechnungen as wu from scipy.special import factorial as fact from math import comb import GHA_triaxial.numeric_examples_panou as ne_panou import GHA_triaxial.numeric_examples_karney as ne_karney from ellipsoide import EllipsoidTriaxial from typing import Callable, Tuple, List from numpy.typing import NDArray def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: """ Berechnung von p und q in elliptischen Koordinaten Panou, Korakitits 2019 :param ell: Ellipsoid :param point: Punkt :return: p und q """ x, y, z = point n = ell.func_n(point) beta, lamb = ell.cart2ell(point) B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2 L = ell.Ex ** 2 - ell.Ee ** 2 * 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 + sqrt(c1 ** 2 - 4 * c0)) / 2 F = ell.Ey ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(lamb) ** 2 p1 = -sqrt(L / (F * t2)) * ell.ax / ell.Ex * sqrt(B) * sin(lamb) p2 = sqrt(L / (F * t2)) * ell.ay * cos(beta) * cos(lamb) p3 = 1 / sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * sin(beta) * sin(2 * lamb) 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[0] * p[1] - n[1] * p[0]]) return p, q def buildODE(ell: EllipsoidTriaxial) -> Callable: """ Aufbau des DGL-Systems :param ell: Ellipsoid :return: DGL-System """ def ODE(s: float, v: NDArray) -> NDArray: """ DGL-System :param s: unabhängige Variable :param v: abhängige Variablen :return: Ableitungen der abhängigen Variablen """ x, dxds, y, dyds, z, dzds = v H = ell.func_H(np.array([x, y, z])) h = dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2 ddx = -(h/H)*x ddy = -(h/H)*y/(1-ell.ee**2) ddz = -(h/H)*z/(1-ell.ex**2) return np.array([dxds, ddx, dyds, ddy, dzds, ddz]) return ODE def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, num: int, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, List]: """ Panou, Korakitits 2019 :param ell: :param point: :param alpha0: :param s: :param num: :param all_points: :return: """ phi, lam, _ = ell.cart2geod(point, "ligas3") p0 = ell.geod2cart(phi, lam, 0) x0, y0, z0 = p0 p, q = pq_ell(ell, p0) dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0) dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0) dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0) v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0]) ode = buildODE(ell) _, werte = rk.rk4(ode, 0, v_init, s, num) x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1] point1 = np.array([x1, y1, z1]) p1, q1 = pq_ell(ell, point1) sigma = np.array([dx1ds, dy1ds, dz1ds]) P = float(p1 @ sigma) Q = float(q1 @ sigma) alpha1 = arctan2(P, Q) if alpha1 < 0: alpha1 += 2 * np.pi if all_points: return point1, alpha1, werte else: return point1, alpha1 # --------------------------------------------------------------------------------------------------------------------- def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: """ Berechnung von p und q in parametrischen Koordinaten Panou, Korakitits 2020 :param ell: Ellipsoid :param point: Punkt :return: p und q """ n = ell.func_n(point) u, v = ell.cart2para(point) # 41-47 G = sqrt(1 - ell.ex ** 2 * cos(u) ** 2 - ell.ee ** 2 * sin(u) ** 2 * sin(v) ** 2) q = np.array([-1 / G * sin(u) * cos(v), -1 / G * sqrt(1 - ell.ee ** 2) * sin(u) * sin(v), 1 / G * sqrt(1 - ell.ex ** 2) * cos(u)]) p = np.array([q[1] * n[2] - q[2] * n[1], q[2] * n[0] - q[0] * n[2], q[0] * n[1] - q[1] * n[0]]) t1 = np.dot(n, q) t2 = np.dot(n, p) t3 = np.dot(p, q) if not (t1 < 1e-10 or t1 > 1-1e-10) and not (t2 < 1e-10 or t2 > 1-1e-10) and not (t3 < 1e-10 or t3 > 1-1e-10): raise Exception("Fehler in den normierten Vektoren") return p, q def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]: """ Panou, Korakitits 2020, 5ff. :param ell: :param point: :param alpha0: :param s: :param maxM: :return: """ x, y, z = point # S. 6 x_m = [x] y_m = [y] z_m = [z] p, q = pq_para(ell, point) # 48-50 x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0)) y_m.append(p[1] * sin(alpha0) + q[1] * cos(alpha0)) z_m.append(p[2] * sin(alpha0) + q[2] * cos(alpha0)) # 34 H_ = lambda p: np.sum([comb(p, 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)]) # 35 h_ = lambda q: np.sum([comb(q, q-j) * (x_m[q-j+1] * x_m[j+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)]) # 31 hH_ = lambda t: 1/H_(0) * (h_(t) - np.sum([comb(t, l-1) * H_(t+1-l) * hH_t[l-1] for l in range(1, t+1)])) # 28-30 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(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)) fact_m = fact(m) # 22-24 a_m.append(x_m[m] / fact_m) b_m.append(y_m[m] / fact_m) c_m.append(z_m[m] / fact_m) # 19-21 x_s = 0 for a in reversed(a_m): x_s = x_s * s + a y_s = 0 for b in reversed(b_m): y_s = y_s * s + b z_s = 0 for c in reversed(c_m): z_s = z_s * s + c p1 = np.array([x_s, y_s, z_s]) p_s, q_s = pq_para(ell, p1) # 57-59 dx_s = 0 for i, a in reversed(list(enumerate(a_m[1:], start=1))): dx_s = dx_s * s + i * a dy_s = 0 for i, b in reversed(list(enumerate(b_m[1:], start=1))): dy_s = dy_s * s + i * b dz_s = 0 for i, c in reversed(list(enumerate(c_m[1:], start=1))): dz_s = dz_s * s + i * c # 52-53 sigma = np.array([dx_s, dy_s, dz_s]) P = float(p_s @ sigma) Q = float(q_s @ sigma) # 51 alpha1 = arctan2(P, Q) if alpha1 < 0: alpha1 += 2 * np.pi return p1, alpha1 def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 16) -> Tuple[NDArray, float]: if s > np.pi / maxPartCircum * ell.ax: s /= 2 point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum) point_end, alpha_end = gha1_ana(ell, point_step, alpha_step, s, maxM, maxPartCircum) else: point_end, alpha_end = gha1_ana_step(ell, point, alpha0, s, maxM) _, _, h = ell.cart2geod(point_end, "ligas3") if h > 1e-5: raise Exception("Analyitsche Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellpsoid") return point_end, alpha_end def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]: point = ell.para2cart(u, v) beta, lamb = ell.para2ell(u, v) p_para, q_para = pq_para(ell, point) sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para) p_ell, q_ell = pq_ell(ell, point) alpha_ell = arctan2(p_ell @ sigma_para, q_ell @ sigma_para) sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) if np.linalg.norm(sigma_para - sigma_ell) > 1e-12: raise Exception("Alpha Umrechnung fehlgeschlagen") return beta, lamb, alpha_ell def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]: point = ell.ell2cart(beta, lamb) u, v = ell.ell2para(beta, lamb) p_ell, q_ell = pq_ell(ell, point) sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell) p_para, q_para = pq_para(ell, point) alpha_para = arctan2(p_para @ sigma_ell, q_para @ sigma_ell) sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para) if np.linalg.norm(sigma_para - sigma_ell) > 1e-9: raise Exception("Alpha Umrechnung fehlgeschlagen") return u, v, alpha_para def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha: float) -> NDArray: p, q = pq_ell(ell, point) sigma = p * sin(alpha) + q * cos(alpha) return sigma def func_sigma_para(ell: EllipsoidTriaxial, point: NDArray, alpha: float) -> NDArray: p, q = pq_para(ell, point) sigma = p * sin(alpha) + q * cos(alpha) return sigma def louville_constant(ell: EllipsoidTriaxial, p0: NDArray, alpha: float) -> float: beta, lamb = ell.cart2ell(p0) l = ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 - ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2 return l def louville_l2c(ell: EllipsoidTriaxial, l: float) -> float: return sqrt((l + ell.Ee**2) / ell.Ex**2) def louville_c2l(ell: EllipsoidTriaxial, c: float) -> float: return ell.Ex**2 * c**2 - ell.Ee**2 if __name__ == "__main__": # ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round") # diffs_panou = [] # examples_panou = ne_panou.get_random_examples(5) # for example in examples_panou: # beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example # P0 = ell.ell2cart(beta0, lamb0) # # P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 100) # beta1_num, lamb1_num = ell.cart2ell(P1_num) # # _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell) # P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 60) # beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) # diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana))) # diffs_panou = np.array(diffs_panou) # mask_360 = (diffs_panou > 359) & (diffs_panou < 361) # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360) # print(diffs_panou) ell: EllipsoidTriaxial = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024") diffs_karney = [] # examples_karney = ne_karney.get_examples((30499, 30500, 40500)) examples_karney = ne_karney.get_random_examples(20) for example in examples_karney: beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example P0 = ell.ell2cart(beta0, lamb0) P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 10000) beta1_num, lamb1_num = ell.cart2ell(P1_num) try: _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell) P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 45, maxPartCircum=32) beta1_ana, lamb1_ana = ell.cart2ell(P1_ana) except: beta1_ana, lamb1_ana = np.inf, np.inf diffs_karney.append((wu.rad2deg(abs(beta1-beta1_num)), wu.rad2deg(abs(lamb1-lamb1_num)), wu.rad2deg(abs(beta1-beta1_ana)), wu.rad2deg(abs(lamb1-lamb1_ana)))) diffs_karney = np.array(diffs_karney) mask_360 = (diffs_karney > 359) & (diffs_karney < 361) diffs_karney[mask_360] = np.abs(diffs_karney[mask_360] - 360) print(diffs_karney)