diff --git a/1,lambda-ES-CMA.py b/1,lambda-ES-CMA.py deleted file mode 100644 index 7887b35..0000000 --- a/1,lambda-ES-CMA.py +++ /dev/null @@ -1,364 +0,0 @@ -import numpy as np -import plotly.graph_objects as go -from ellipsoide import EllipsoidTriaxial - - -def ell2cart(ell: EllipsoidTriaxial, beta, lamb): - x = ell.ax * np.cos(beta) * np.cos(lamb) - y = ell.ay * np.cos(beta) * np.sin(lamb) - z = ell.b * np.sin(beta) - return np.array([x, y, z], dtype=float) - - -def liouville(ell: EllipsoidTriaxial, beta, lamb, alpha): - k = (ell.Ee**2) / (ell.Ex**2) - c_sq = (np.cos(beta)**2 + k * np.sin(beta)**2) * np.sin(alpha)**2 \ - + k * np.cos(lamb)**2 * np.cos(alpha)**2 - return c_sq - - -def normalize_angles(beta, lamb): - beta = np.clip(beta, -np.pi / 2.0, np.pi / 2.0) - lamb = (lamb + np.pi) % (2 * np.pi) - np.pi - return beta, lamb - - -def compute_azimuth(beta0: float, lamb0: float, - beta1: float, lamb1: float) -> float: - dlam = lamb1 - lamb0 - y = np.cos(beta1) * np.sin(dlam) - x = np.cos(beta0) * np.sin(beta1) - np.sin(beta0) * np.cos(beta1) * np.cos(dlam) - alpha = np.arctan2(y, x) - return alpha - - -def sphere_forward_step(beta0: float, - lamb0: float, - alpha0: float, - s: float, - R: float) -> tuple[float, float]: - - delta = s / R - sin_beta0 = np.sin(beta0) - cos_beta0 = np.cos(beta0) - - sin_beta2 = sin_beta0 * np.cos(delta) + cos_beta0 * np.sin(delta) * np.cos(alpha0) - beta2 = np.arcsin(sin_beta2) - - y = np.sin(alpha0) * np.sin(delta) * cos_beta0 - x = np.cos(delta) - sin_beta0 * sin_beta2 - dlam = np.arctan2(y, x) - lamb2 = lamb0 + dlam - - beta2, lamb2 = normalize_angles(beta2, lamb2) - return beta2, lamb2 - - -def local_step_objective(candidate: np.ndarray, - beta_start: float, - lamb_start: float, - alpha_prev: float, - c_target: float, - step_length: float, - ell: EllipsoidTriaxial, - beta_pred: float, - lamb_pred: float, - w_L: float = 5.0, - w_d: float = 5.0, - w_p: float = 2.0, - w_a: float = 2.0) -> float: - - beta1, lamb1 = candidate - beta1, lamb1 = normalize_angles(beta1, lamb1) - - P0 = ell2cart(ell, beta_start, lamb_start) - P1 = ell2cart(ell, beta1, lamb1) - d = np.linalg.norm(P1 - P0) - - alpha1 = compute_azimuth(beta_start, lamb_start, beta1, lamb1) - - c1 = liouville(ell, beta1, lamb1, alpha1) - J_L = (c1 - c_target) ** 2 - - J_d = (d - step_length) ** 2 - - d_beta = beta1 - beta_pred - d_lamb = lamb1 - lamb_pred - d_ang2 = d_beta**2 + (np.cos(beta_pred) * d_lamb)**2 - J_p = d_ang2 - - d_alpha = np.arctan2(np.sin(alpha1 - alpha_prev), - np.cos(alpha1 - alpha_prev)) - J_a = d_alpha**2 - - return w_L * J_L + w_d * J_d + w_p * J_p + w_a * J_a - - - - -def ES_CMA_step(beta_start: float, - lamb_start: float, - alpha_prev: float, - c_target: float, - step_length: float, - ell: EllipsoidTriaxial, - beta_pred: float, - lamb_pred: float, - sigma0: float, - stopfitness: float = 1e-18) -> tuple[float, float]: - - N = 2 - xmean = np.array([beta_pred, lamb_pred], dtype=float) - sigma = sigma0 - stopeval = int(400 * N**2) - - lamb = 30 - mu = 1 - weights = np.array([1.0]) - mueff = 1.0 - - cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N) - cs = (mueff + 2) / (N + mueff + 5) - c1 = 2 / ((N + 1.3)**2 + mueff) - cmu = min(1 - c1, - 2 * (mueff - 2 + 1/mueff) / ((N + 2)**2 + mueff)) - damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs - - pc = np.zeros(N) - ps = np.zeros(N) - B = np.eye(N) - D = np.eye(N) - C = B @ D @ (B @ D).T - eigeneval = 0 - chiN = np.sqrt(N) * (1 - 1 / (4 * N) + 1 / (21 * N**2)) - - counteval = 0 - arx = np.zeros((N, lamb)) - arz = np.zeros((N, lamb)) - arfitness = np.zeros(lamb) - - while counteval < stopeval: - for k in range(lamb): - arz[:, k] = np.random.randn(N) - arx[:, k] = xmean + sigma * (B @ D @ arz[:, k]) - - arfitness[k] = local_step_objective( - arx[:, k], - beta_start, lamb_start, - alpha_prev, - c_target, - step_length, - ell, - beta_pred, lamb_pred - ) - counteval += 1 - - idx = np.argsort(arfitness) - arfitness = arfitness[idx] - arindex = idx - - xold = xmean.copy() - xmean = arx[:, arindex[:mu]] @ weights - zmean = arz[:, arindex[:mu]] @ weights - - ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * (B @ zmean) - norm_ps = np.linalg.norm(ps) - hsig = norm_ps / np.sqrt(1 - (1 - cs)**(2 * counteval / lamb)) / chiN < 1.4 + 2 / (N + 1) - hsig = 1.0 if hsig else 0.0 - pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (B @ D @ zmean) - - BDz = B @ D @ arz[:, arindex[:mu]] - C = (1 - c1 - cmu) * C \ - + c1 * (np.outer(pc, pc) + (1 - hsig) * cc * (2 - cc) * C) \ - + cmu * BDz @ np.diag(weights) @ BDz.T - - sigma = sigma * np.exp((cs / damps) * (norm_ps / chiN - 1)) - - if counteval - eigeneval > lamb / ((c1 + cmu) * N * 10): - eigeneval = counteval - C = (C + C.T) / 2.0 - eigvals, B = np.linalg.eigh(C) - D = np.diag(np.sqrt(np.maximum(eigvals, 1e-20))) - - if arfitness[0] <= stopfitness: - break - - xmin = arx[:, arindex[0]] - beta1, lamb1 = normalize_angles(xmin[0], xmin[1]) - return beta1, lamb1 - - -def march_geodesic(beta1: float, - lamb1: float, - alpha1: float, - S_total: float, - step_length: float, - ell: EllipsoidTriaxial): - - beta_curr = beta1 - lamb_curr = lamb1 - alpha_curr = alpha1 - - betas = [beta_curr] - lambs = [lamb_curr] - alphas = [alpha_curr] - - c_target = liouville(ell, beta_curr, lamb_curr, alpha_curr) - - total_distance = 0.0 - R_sphere = ell.ax - sigma0 = 1e-10 - - while total_distance < S_total - 1e-6: - remaining = S_total - total_distance - this_step = min(step_length, remaining) - - beta_pred, lamb_pred = sphere_forward_step( - beta_curr, lamb_curr, alpha_curr, - this_step, R_sphere - ) - - beta_next, lamb_next = ES_CMA_step( - beta_curr, lamb_curr, alpha_curr, - c_target, - this_step, - ell, - beta_pred, - lamb_pred, - sigma0=sigma0 - ) - - P_curr = ell2cart(ell, beta_curr, lamb_curr) - P_next = ell2cart(ell, beta_next, lamb_next) - d_step = np.linalg.norm(P_next - P_curr) - total_distance += d_step - - alpha_next = compute_azimuth(beta_curr, lamb_curr, beta_next, lamb_next) - - beta_curr, lamb_curr, alpha_curr = beta_next, lamb_next, alpha_next - betas.append(beta_curr) - lambs.append(lamb_curr) - alphas.append(alpha_curr) - - return np.array(betas), np.array(lambs), np.array(alphas), total_distance - - - -if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - - beta1 = np.deg2rad(10.0) - lamb1 = np.deg2rad(10.0) - alpha1 = 0.7494562507041596 # Ergebnis 20°, 20° - - STEP_LENGTH = 500.0 - S_total = 1542703.1458877102 - - betas, lambs, alphas, S_real = march_geodesic( - beta1, lamb1, alpha1, - S_total, - STEP_LENGTH, - ell - ) - - print("Anzahl Schritte:", len(betas) - 1) - print("Resultierende Gesamtstrecke (Chord, Ellipsoid):", S_real, "m") - print("Letzter Punkt (beta, lambda) in Grad:", - np.rad2deg(betas[-1]), np.rad2deg(lambs[-1])) - -def plot_geodesic_3d(ell: EllipsoidTriaxial, - betas: np.ndarray, - lambs: np.ndarray, - n_beta: int = 60, - n_lamb: int = 120): - - beta_grid = np.linspace(-np.pi/2, np.pi/2, n_beta) - lamb_grid = np.linspace(-np.pi, np.pi, n_lamb) - B, L = np.meshgrid(beta_grid, lamb_grid, indexing="ij") - - Xs = ell.ax * np.cos(B) * np.cos(L) - Ys = ell.ay * np.cos(B) * np.sin(L) - Zs = ell.b * np.sin(B) - - Xp = ell.ax * np.cos(betas) * np.cos(lambs) - Yp = ell.ay * np.cos(betas) * np.sin(lambs) - Zp = ell.b * np.sin(betas) - - fig = go.Figure() - - fig.add_trace( - go.Surface( - x=Xs, - y=Ys, - z=Zs, - opacity=0.6, - showscale=False, - colorscale="Viridis", - name="Ellipsoid" - ) - ) - - fig.add_trace( - go.Scatter3d( - x=Xp, - y=Yp, - z=Zp, - mode="lines+markers", - line=dict(width=5, color="red"), - marker=dict(size=4, color="black"), - name="ES-CMA Schritte" - ) - ) - - fig.add_trace( - go.Scatter3d( - x=[Xp[0]], y=[Yp[0]], z=[Zp[0]], - mode="markers", - marker=dict(size=6, color="green"), - name="Start" - ) - ) - fig.add_trace( - go.Scatter3d( - x=[Xp[-1]], y=[Yp[-1]], z=[Zp[-1]], - mode="markers", - marker=dict(size=6, color="blue"), - name="Ende" - ) - ) - - fig.update_layout( - title="ES-CMA", - scene=dict( - xaxis_title="X [m]", - yaxis_title="Y [m]", - zaxis_title="Z [m]", - aspectmode="data" - ) - ) - - fig.show() - -plot_geodesic_3d(ell, betas, lambs) - - - - - - - - - - - - - - - - - - - - - - diff --git a/GHA_triaxial/approx_gha1_2.py b/GHA_triaxial/approx_gha1_2.py deleted file mode 100644 index 2070f1f..0000000 --- a/GHA_triaxial/approx_gha1_2.py +++ /dev/null @@ -1,120 +0,0 @@ -import numpy as np -from ellipsoide import EllipsoidTriaxial -from panou import louville_constant, func_sigma_ell, gha1_ana -import plotly.graph_objects as go -import winkelumrechnungen as wu -from numpy import sin, cos, arccos - -def Bogenlaenge(P1: NDArray, P2: NDArray) -> float: - """ - Berechnung der mittleren Bogenlänge zwischen zwei kartesischen Punkten - :param P1: kartesische Koordinate Punkt 1 - :param P2: kartesische Koordinate Punkt 2 - :return: Bogenlänge s - """ - R1 = np.linalg.norm(P1) - R2 = np.linalg.norm(P2) - R = 0.5 * (R1 + R2) - if P1 @ P2 / (R1 * R2) > 1: - s = np.linalg.norm(P1 - P2) - else: - theta = arccos(P1 @ P2 / (R1 * R2)) - s = float(R * theta) - return s - -def gha1_approx2(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]: - """ - Berechung einer Näherungslösung der ersten Hauptaufgabe - :param ell: Ellipsoid - :param p0: Anfangspunkt - :param alpha0: Azimut im Anfangspunkt - :param s: Strecke bis zum Endpunkt - :param ds: Länge einzelner Streckenelemente - :param all_points: Ausgabe aller Punkte als Array? - :return: Endpunkt, Azimut im Endpunkt, optional alle Punkte - """ - l0 = louville_constant(ell, p0, alpha0) - points = [p0] - alphas = [alpha0] - s_curr = 0.0 - - while s_curr < s: - ds_target = min(ds, s - s_curr) - if ds_target < 1e-8: - break - - p1 = points[-1] - alpha1 = alphas[-1] - alpha1_mid = alphas[-1] - p2 = points[-1] - alpha2 = alphas[-1] - - i = 0 - while i < 2: - i += 1 - - sigma = func_sigma_ell(ell, p1, alpha1_mid) - p2_new = p1 + ds_target * sigma - p2_new = ell.point_onto_ellipsoid(p2_new) - p2 = p2_new - - j = 0 - while j < 2: - j += 1 - - dalpha = 1e-6 - l2 = louville_constant(ell, p2, alpha2) - dl_dalpha = (louville_constant(ell, p2, alpha2 + dalpha) - l2) / dalpha - alpha2_new = alpha2 + (l0 - l2) / dl_dalpha - alpha2 = alpha2_new - - alpha1_mid = (alpha1 + alpha2) / 2 - - points.append(p2) - alphas.append(alpha2) - - ds_actual = np.linalg.norm(p2 - p1) - s_curr += ds_actual - if s_curr > 10000000: - pass - - if all_points: - return points[-1], alphas[-1], np.array(points) - else: - return points[-1], alphas[-1] - -def show_points(points: NDArray, p0: NDArray, p1: NDArray): - """ - Anzeigen der Punkte - :param points: Array aller approximierten Punkte - :param p0: Startpunkt - :param p1: wahrer Endpunkt - """ - fig = go.Figure() - - fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2], - mode='lines', line=dict(color="red", width=3), name="Approx") - fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]], - mode='markers', marker=dict(color="green"), name="P0") - fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]], - mode='markers', marker=dict(color="green"), name="P1") - - fig.update_layout( - scene=dict(xaxis_title='X [km]', - yaxis_title='Y [km]', - zaxis_title='Z [km]', - aspectmode='data'), - title="CHAMP") - - fig.show() - - -if __name__ == '__main__': - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - P0 = ell.para2cart(0.2, 0.3) - alpha0 = wu.deg2rad(35) - s = 13000000 - P1_app, alpha1_app, points = gha1_approx2(ell, P0, alpha0, s, ds=10000, all_points=True) - P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16) - show_points(points, P0, P1_ana) - print(np.linalg.norm(P1_app - P1_ana)) diff --git a/GHA_triaxial/gha1_ana.py b/GHA_triaxial/gha1_ana.py new file mode 100644 index 0000000..c9cb35b --- /dev/null +++ b/GHA_triaxial/gha1_ana.py @@ -0,0 +1,136 @@ +from math import comb +from typing import Tuple + +import numpy as np +from numpy import sin, cos, arctan2 +from numpy._typing import NDArray +from scipy.special import factorial as fact + +from ellipsoide import EllipsoidTriaxial +from GHA_triaxial.utils import pq_para + + +def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]: + """ + Panou, Korakitits 2020, 5ff. + :param ell: Ellipsoid + :param point: Punkt in kartesischen Koordinaten + :param alpha0: Azimut im Startpunkt + :param s: Strecke + :param maxM: maximale Ordnung + :return: Zwischenpunkt, Azimut im Zwischenpunkt + """ + 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]: + """ + :param ell: Ellipsoid + :param point: Punkt in kartesischen Koordinaten + :param alpha0: Azimut im Startpunkt + :param s: Strecke + :param maxM: maximale Ordnung + :param maxPartCircum: maximale Aufteilung (1/x halber Ellipsoidumfang) + :return: Zielpunkt, Azimut im Zielpunkt + """ + 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("Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid") + + return point_end, alpha_end diff --git a/GHA_triaxial/approx_gha1.py b/GHA_triaxial/gha1_approx.py similarity index 96% rename from GHA_triaxial/approx_gha1.py rename to GHA_triaxial/gha1_approx.py index c571fa5..042f54a 100644 --- a/GHA_triaxial/approx_gha1.py +++ b/GHA_triaxial/gha1_approx.py @@ -1,6 +1,7 @@ import numpy as np from ellipsoide import EllipsoidTriaxial -from GHA_triaxial.panou import louville_constant, func_sigma_ell, gha1_ana +from GHA_triaxial.gha1_ana import gha1_ana +from GHA_triaxial.utils import func_sigma_ell, louville_constant import plotly.graph_objects as go import winkelumrechnungen as wu diff --git a/GHA_triaxial/gha1_num.py b/GHA_triaxial/gha1_num.py new file mode 100644 index 0000000..c527216 --- /dev/null +++ b/GHA_triaxial/gha1_num.py @@ -0,0 +1,129 @@ +import numpy as np +from numpy import sin, cos, arctan2 +import ellipsoide +import runge_kutta as rk +import winkelumrechnungen as wu +import GHA_triaxial.numeric_examples_karney as ne_karney +from GHA_triaxial.gha1_ana import gha1_ana +from ellipsoide import EllipsoidTriaxial +from typing import Callable, Tuple, List +from numpy.typing import NDArray + +from GHA_triaxial.utils import alpha_ell2para, pq_ell + + +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: Ellipsoid + :param point: Punkt in kartesischen Koordinaten + :param alpha0: Azimut im Startpunkt + :param s: Strecke + :param num: Anzahl Zwischenpunkte + :param all_points: Ausgabe aller Punkte? + :return: Zielpunkt, Azimut im Zielpunkt (, alle Punkte) + """ + 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]) + + 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 + + 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 + +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) diff --git a/GHA_triaxial/ES_gha2.py b/GHA_triaxial/gha2_ES.py similarity index 98% rename from GHA_triaxial/ES_gha2.py rename to GHA_triaxial/gha2_ES.py index 2bc4713..d403b02 100644 --- a/GHA_triaxial/ES_gha2.py +++ b/GHA_triaxial/gha2_ES.py @@ -1,11 +1,10 @@ import numpy as np -from numpy import arccos from Hansen_ES_CMA import escma from ellipsoide import EllipsoidTriaxial from numpy.typing import NDArray import plotly.graph_objects as go -from GHA_triaxial.panou_2013_2GHA_num import gha2_num -from utils import sigma2alpha +from GHA_triaxial.gha2_num import gha2_num +from GHA_triaxial.utils import sigma2alpha ell_ES: EllipsoidTriaxial = None P_left: NDArray = None diff --git a/GHA_triaxial/approx_gha2.py b/GHA_triaxial/gha2_approx.py similarity index 97% rename from GHA_triaxial/approx_gha2.py rename to GHA_triaxial/gha2_approx.py index 07567e2..dd5316d 100644 --- a/GHA_triaxial/approx_gha2.py +++ b/GHA_triaxial/gha2_approx.py @@ -1,12 +1,12 @@ import numpy as np from ellipsoide import EllipsoidTriaxial -from GHA_triaxial.panou_2013_2GHA_num import gha2_num +from GHA_triaxial.gha2_num import gha2_num import plotly.graph_objects as go import winkelumrechnungen as wu from numpy.typing import NDArray from typing import Tuple -from utils import sigma2alpha +from GHA_triaxial.utils import sigma2alpha def gha2_approx(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]: diff --git a/GHA_triaxial/panou_2013_2GHA_num.py b/GHA_triaxial/gha2_num.py similarity index 100% rename from GHA_triaxial/panou_2013_2GHA_num.py rename to GHA_triaxial/gha2_num.py diff --git a/GHA_triaxial/panou.py b/GHA_triaxial/panou.py deleted file mode 100644 index 536892e..0000000 --- a/GHA_triaxial/panou.py +++ /dev/null @@ -1,365 +0,0 @@ -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]) - p = p / np.linalg.norm(p) - 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]]) - q = q / np.linalg.norm(q) - - 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") - - p = p / np.linalg.norm(p) - q = q / np.linalg.norm(q) - - 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("Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid") - - 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) diff --git a/GHA_triaxial/utils.py b/GHA_triaxial/utils.py new file mode 100644 index 0000000..0581234 --- /dev/null +++ b/GHA_triaxial/utils.py @@ -0,0 +1,183 @@ +from typing import Tuple + +import numpy as np +from numpy import arctan2, sin, cos, sqrt +from numpy._typing import NDArray +from numpy.typing import NDArray + +from ellipsoide import EllipsoidTriaxial + + +def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float: + """ + Berechnung des Azimuts an einem Punkt anhand der Ableitung zu den kartesischen Koordinaten + :param ell: Ellipsoid + :param sigma: Ableitungsvektor ver kartesischen Koordinaten + :param point: Punkt + :return: Azimuts + """ + p, q = pq_ell(ell, point) + P = float(p @ sigma) + Q = float(q @ sigma) + + alpha = arctan2(P, Q) + return alpha + + +def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]: + """ + Umrechnung des Azimuts bezogen auf parametrische Koordinaten zu ellipsoidischen + :param ell: Ellipsoid + :param u: parametrische Breite + :param v: parametrische Länge + :param alpha_para: Azimut bezogen auf parametrische Koordinaten + :return: Azimut bezogen auf ellipsoidische Koordinaten + """ + 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]: + """ + Umrechnung des Azimuts bezogen auf ellipsoidische Koordinaten zu parametrischen + :param ell: Ellipsoid + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge + :param alpha_ell: Azimut bezogen auf ellipsoidische Koordinaten + :return: Azimut bezogen auf parametrische Koordinaten + """ + 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_ell: float) -> NDArray: + """ + Berechnung der Richtungsableitungen in kartesischen Koordinaten aus ellipsoidischem Azimut + Panou (2019) [6] + :param ell: Ellipsoid + :param point: Punkt in kartesischen Koordinaten + :param alpha_ell: ellipsoidischer Azimut + :return: Richtungsableitungen in kartesischen Koordinaten + """ + p, q = pq_ell(ell, point) + sigma = p * sin(alpha_ell) + q * cos(alpha_ell) + return sigma + + +def func_sigma_para(ell: EllipsoidTriaxial, point: NDArray, alpha_para: float) -> NDArray: + """ + Berechnung der Richtungsableitungen in kartesischen Koordinaten aus parametischem Azimut + Panou, Korakitis (2019) [6] + :param ell: Ellipsoid + :param point: Punkt in kartesischen Koordinaten + :param alpha_para: parametrischer Azimut + :return: Richtungsableitungen in kartesischen Koordinaten + """ + p, q = pq_para(ell, point) + sigma = p * sin(alpha_para) + q * cos(alpha_para) + return sigma + + +def louville_constant(ell: EllipsoidTriaxial, p0: NDArray, alpha_ell: float) -> float: + """ + Berechnung der Louville Konstanten + Panou, Korakitis (2019) [6] + :param ell: Ellipsoid + :param p0: Punkt in kartesischen Koordinaten + :param alpha_ell: ellipsoidischer Azimut + :return: + """ + beta, lamb = ell.cart2ell(p0) + l = ell.Ey**2 * cos(beta)**2 * sin(alpha_ell)**2 - ell.Ee**2 * sin(lamb)**2 * cos(alpha_ell)**2 + return l + + +def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: + """ + Berechnung von p (Tangente entlang konstantem beta) und q (Tangente entlang konstantem lambda) + Panou, Korakitits (2019) [5f.] + :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]) + p = p / np.linalg.norm(p) + 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]]) + q = q / np.linalg.norm(q) + + return p, q + + +def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: + """ + Berechnung von p (Tangente entlang konstantem u) und q (Tangente entlang konstantem v) + 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") + + p = p / np.linalg.norm(p) + q = q / np.linalg.norm(q) + + return p, q diff --git a/Numerische_Integration/__init__.py b/Tests/__init__.py similarity index 100% rename from Numerische_Integration/__init__.py rename to Tests/__init__.py diff --git a/Tests/algorithms_test.ipynb b/Tests/algorithms_test.ipynb new file mode 100644 index 0000000..bc8e2b5 --- /dev/null +++ b/Tests/algorithms_test.ipynb @@ -0,0 +1,2655 @@ +{ + "cells": [ + { + "cell_type": "code", + "id": "initial_id", + "metadata": { + "collapsed": true, + "ExecuteTime": { + "end_time": "2026-02-04T17:56:25.185202Z", + "start_time": "2026-02-04T17:56:22.991833Z" + } + }, + "source": [ + "%load_ext autoreload\n", + "%autoreload 2" + ], + "outputs": [], + "execution_count": 1 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T17:56:35.920624Z", + "start_time": "2026-02-04T17:56:25.201916Z" + } + }, + "cell_type": "code", + "source": [ + "%reload_ext autoreload\n", + "%autoreload 2\n", + "import time\n", + "from numpy import nan\n", + "import winkelumrechnungen as wu\n", + "import os\n", + "from contextlib import contextmanager, redirect_stdout, redirect_stderr\n", + "import plotly.graph_objects as go\n", + "import warnings\n", + "\n", + "from ellipsoide import EllipsoidTriaxial\n", + "from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n", + "\n", + "from GHA_triaxial.gha1_num import gha1_num\n", + "from GHA_triaxial.gha1_ana import gha1_ana\n", + "from GHA_triaxial.gha1_approx import gha1_approx\n", + "\n", + "from GHA_triaxial.gha2_num import gha2_num\n", + "from GHA_triaxial.gha2_ES import gha2_ES\n", + "from GHA_triaxial.gha2_approx import gha2_approx\n", + "\n", + "from GHA_triaxial.numeric_examples_panou import get_tables as get_tables_panou\n", + "from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney" + ], + "id": "961cb22764c5bcb9", + "outputs": [], + "execution_count": 2 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T17:56:36.470974Z", + "start_time": "2026-02-04T17:56:35.937646Z" + } + }, + "cell_type": "code", + "source": [ + "@contextmanager\n", + "def suppress_print():\n", + " with open(os.devnull, 'w') as fnull:\n", + " with redirect_stdout(fnull), redirect_stderr(fnull):\n", + " yield" + ], + "id": "e4740625a366ac13", + "outputs": [], + "execution_count": 3 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T18:11:14.405810Z", + "start_time": "2026-02-04T18:11:14.169636Z" + } + }, + "cell_type": "code", + "source": [ + "# steps_gha1_num = [2000, 5000, 10000, 20000]\n", + "# maxM_gha1_ana = [20, 50]\n", + "# parts_gha1_ana = [4, 8]\n", + "# dsPart_gha1_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", + "#\n", + "# steps_gha2_num = [2000, 5000, 10000, 20000]\n", + "# dsPart_gha2_ES = [600, 1250, 6000] # entspricht bei der Erde ca. 10000, 5000, 1000\n", + "# dsPart_gha2_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", + "\n", + "steps_gha1_num = [2000, 5000, 10000]\n", + "maxM_gha1_ana = [20, 50]\n", + "parts_gha1_ana = [2, 8, 32]\n", + "dsPart_gha1_approx = [600, 1250, 6000, 60000]\n", + "\n", + "steps_gha2_num = [100, 1000]\n", + "dsPart_gha2_ES = [20, 100]\n", + "dsPart_gha2_approx = [600, 1250, 6000, 60000]" + ], + "id": "96093cdde03f8d57", + "outputs": [], + "execution_count": 16 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T18:11:15.988345Z", + "start_time": "2026-02-04T18:11:15.774321Z" + } + }, + "cell_type": "code", + "source": [ + "test = \"Karney\"\n", + "# test = \"Panou\"\n", + "if test == \"Karney\":\n", + " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n", + " examples = get_examples_karney(10, 42)\n", + "elif test == \"Panou\":\n", + " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n", + " tables = get_tables_panou()\n", + " table_indices = []\n", + " examples = []\n", + " for i, table in enumerate(tables):\n", + " for example in table:\n", + " table_indices.append(i+1)\n", + " examples.append(example)" + ], + "id": "6e384cc01c2dbe", + "outputs": [], + "execution_count": 17 + }, + { + "metadata": { + "jupyter": { + "is_executing": true + }, + "ExecuteTime": { + "start_time": "2026-02-04T18:11:18.785718Z" + } + }, + "cell_type": "code", + "source": [ + "results = {}\n", + "for i, example in enumerate(examples):\n", + " print(f\"----- Beispiel {i+1}/{len(examples)}\")\n", + " example_results = {}\n", + "\n", + " beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n", + " P0 = ell.ell2cart(beta0, lamb0)\n", + " P1 = ell.ell2cart(beta1, lamb1)\n", + " _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n", + "\n", + " for steps in steps_gha1_num:\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n", + " end = time.perf_counter()\n", + " beta1_num, lamb1_num = ell.cart2ell(P1_num)\n", + " d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600\n", + " d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600\n", + " d_time = end - start\n", + " example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n", + "\n", + " for maxM in maxM_gha1_ana:\n", + " for parts in parts_gha1_ana:\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n", + " end = time.perf_counter()\n", + " beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", + " _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", + " d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n", + " d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n", + " d_time = end - start\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha1_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n", + " end = time.perf_counter()\n", + " beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n", + " d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n", + " d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", + " d_time = end - start\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " for steps in steps_gha2_num:\n", + " start = time.perf_counter()\n", + " try:\n", + " with warnings.catch_warnings():\n", + " warnings.simplefilter(\"ignore\", RuntimeWarning)\n", + " alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n", + " d_s = abs(s_num - s) / 1000\n", + " d_time = end - start\n", + " example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha2_ES:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " with suppress_print():\n", + " alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n", + " d_s = abs(s_ES - s) / 1000\n", + " d_time = end - start\n", + " example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " for dsPart in dsPart_gha2_approx:\n", + " ds = ell.ax/dsPart\n", + " start = time.perf_counter()\n", + " try:\n", + " alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n", + " end = time.perf_counter()\n", + " d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n", + " d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", + " d_s = abs(s_approx - s) / 1000\n", + " d_time = end - start\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", + " except Exception as e:\n", + " print(e)\n", + " example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n", + "\n", + " results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results" + ], + "id": "ef8849908ed3231e", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "----- Beispiel 1/10\n", + "Keine Multi-Start-Variante konvergiert.\n", + "Keine Multi-Start-Variante konvergiert.\n", + "----- Beispiel 2/10\n", + "Keine Multi-Start-Variante konvergiert.\n", + "Keine Multi-Start-Variante konvergiert.\n", + "----- Beispiel 3/10\n", + "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", + "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", + "Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid\n", + "Fehler in der Umrechnung cart2ell\n", + "Keine Multi-Start-Variante konvergiert.\n", + "Keine Multi-Start-Variante konvergiert.\n" + ] + } + ], + "execution_count": null + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T14:35:09.900827Z", + "start_time": "2026-02-04T14:35:09.690955Z" + } + }, + "cell_type": "code", + "source": [ + "# with open(f\"gha_results{test}.pkl\", \"wb\") as f:\n", + "# pickle.dump(results, f)" + ], + "id": "664105ea35d50a7b", + "outputs": [], + "execution_count": 7 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T14:35:10.113762Z", + "start_time": "2026-02-04T14:35:09.910687Z" + } + }, + "cell_type": "code", + "source": [ + "# with open(f\"gha_results{test}.pkl\", \"rb\") as f:\n", + "# results = pickle.load(f)" + ], + "id": "ad8aa9dcc8af4a05", + "outputs": [], + "execution_count": 8 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T18:01:30.114488Z", + "start_time": "2026-02-04T18:01:25.723109Z" + } + }, + "cell_type": "code", + "source": [ + "def format_max(values, is_angle=False):\n", + " arr = np.array(values, dtype=float)\n", + " if arr.size==0:\n", + " return np.nan\n", + " maxi = np.nanmax(np.abs(arr))\n", + " if maxi is None or (isinstance(maxi,float) and (math.isnan(maxi))):\n", + " return \"nan\"\n", + " # i = 2\n", + " # while maxi > np.nanmean(np.abs(arr)):\n", + " # maxi = np.sort(np.abs(arr[~np.isnan(arr)]))[-i]\n", + " # i += 1\n", + " if is_angle:\n", + " maxi = wu.rad2deg(maxi)*3600\n", + " if f\"{maxi:.3g}\" == 0:\n", + " pass\n", + " return f\"{maxi:.3g}\"\n", + "\n", + "\n", + "def build_max_table(gha_prefix, title, group_value = None):\n", + " # metrics\n", + " if gha_prefix==\"GHA1\":\n", + " metrics = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n", + " angle_mask = [True, True, True, False]\n", + " else:\n", + " metrics = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n", + " angle_mask = [True, True, False, False]\n", + "\n", + " # collect keys in this group\n", + " if group_value is None:\n", + " example_keys = [example_key for example_key in results.keys()]\n", + " else:\n", + " example_keys = [example_key for example_key, group_index in zip(results.keys(), table_indices) if group_index==group_value]\n", + " # all variant keys (inner dict keys) matching prefix\n", + " algorithms = sorted({algorithm for example_key in example_keys for algorithm in results[example_key].keys() if algorithm.startswith(gha_prefix)})\n", + "\n", + " header = [\"Algorithmus\", \"Parameter\", \"NaN\"] + list(metrics)\n", + " cells = [[] for i in range(len(metrics) + 3)]\n", + " for algorithm in algorithms:\n", + "\n", + " ghaNr, variant, params = algorithm.split(\"_\", 2)\n", + " cells[0].append(variant)\n", + " cells[1].append(params)\n", + " nan_values = []\n", + " for example_key in example_keys:\n", + " nan_values.append(results[example_key][algorithm][0])\n", + " cells[2].append(np.sum(np.isnan(nan_values)))\n", + " for i, metric in enumerate(metrics):\n", + " values = []\n", + " for example_key in example_keys:\n", + " values.append(results[example_key][algorithm][i])\n", + " cells[i+3].append(format_max(values, is_angle=angle_mask[i]))\n", + "\n", + " header = dict(\n", + " values=header,\n", + " align=\"center\",\n", + " fill_color=\"lightgrey\",\n", + " font=dict(size=13)\n", + " )\n", + " cells = dict(\n", + " values=cells,\n", + " align=\"center\"\n", + " )\n", + "\n", + " fig = go.Figure(data=[go.Table(header=header, cells=cells)])\n", + " fig.update_layout(title=title,\n", + " template=\"simple_white\",\n", + " width=800,\n", + " height=280,\n", + " margin=dict(l=20, r=20, t=60, b=20))\n", + " return fig\n", + "\n", + "figs = []\n", + "if test == \"Panou\":\n", + " for table_index in sorted(set(table_indices)):\n", + " fig1 = build_max_table(\"GHA1\", f\"{test} - Gruppe {table_index} - GHA1\", table_index)\n", + " fig2 = build_max_table(\"GHA2\", f\"{test} - Gruppe {table_index} - GHA2\", table_index)\n", + " figs.append(fig1)\n", + " figs.append(fig2)\n", + "elif test == \"Karney\":\n", + " fig1 = build_max_table(\"GHA1\", f\"{test} - GHA1\")\n", + " fig2 = build_max_table(\"GHA2\", f\"{test} - GHA2\")\n", + " figs.append(fig1)\n", + " figs.append(fig2)\n", + "\n", + "for fig in figs:\n", + " fig.show()" + ], + "id": "b46d57fc0d794e28", + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\moell\\AppData\\Local\\Temp\\ipykernel_25588\\2480185505.py:5: RuntimeWarning:\n", + "\n", + "All-NaN slice encountered\n", + "\n" + ] + }, + { + "data": { + "application/vnd.plotly.v1+json": { + "data": [ + { + "cells": { + "align": "center", + "values": [ + [ + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "ana", + "approx", + "approx", + "num", + "num" + ], + [ + "20_16", + "20_2", + "20_24", + "20_4", + "20_8", + "40_16", + "40_2", + "40_24", + "40_4", + "40_8", + "60_16", + "60_2", + "60_24", + "60_4", + "60_8", + "1250", + "600", + "2000", + "5000" + ], + [ + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0, + 0 + ], + [ + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "5.1e+03", + "7.16e+03", + "7.15e+03", + "3.67e+03", + "3.67e+03" + ], + [ + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.06e+04", + "1.67e+04", + "1.67e+04", + "3.99e+03", + "3.99e+03" + ], + [ + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "1.55e+04", + "2.47e+03", + "2.46e+03", + "1.82e+03", + "1.82e+03" + ], + [ + "0.02", + "0.00746", + "0.0353", + "0.0123", + "0.0103", + "0.0868", + "0.0558", + "0.201", + "0.0561", + "0.0435", + "0.251", + "0.0799", + "0.495", + "0.0816", + "0.126", + "0.731", + "0.294", + "0.0847", + "0.231" + ] + ] + }, + "header": { + "align": "center", + "fill": { + "color": "lightgrey" + }, + "font": { + "size": 13 + }, + "values": [ + "Algorithmus", + "Parameter", + "NaN", + "dBeta [\"]", + "dLambda [\"]", + "dAlpha1 [\"]", + "time [s]" + ] + }, + "type": "table" + } + ], + "layout": { + "template": { + "data": { + "barpolar": [ + { + "marker": { + "line": { + "color": "white", + "width": 0.5 + }, + "pattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + } + }, + "type": "barpolar" + } + ], + "bar": [ + { + "error_x": { + "color": "rgb(36,36,36)" + }, + "error_y": { + "color": "rgb(36,36,36)" + }, + "marker": { + "line": { + "color": "white", + "width": 0.5 + }, + "pattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + } + }, + "type": "bar" + } + ], + "carpet": [ + { + "aaxis": { + "endlinecolor": "rgb(36,36,36)", + "gridcolor": "white", + "linecolor": "white", + "minorgridcolor": "white", + "startlinecolor": "rgb(36,36,36)" + }, + "baxis": { + "endlinecolor": "rgb(36,36,36)", + "gridcolor": "white", + "linecolor": "white", + "minorgridcolor": "white", + "startlinecolor": "rgb(36,36,36)" + }, + "type": "carpet" + } + ], + "choropleth": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "choropleth" + } + ], + "contourcarpet": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "contourcarpet" + } + ], + "contour": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "contour" + } + ], + "heatmap": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "heatmap" + } + ], + "histogram2dcontour": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "histogram2dcontour" + } + ], + "histogram2d": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "histogram2d" + } + ], + "histogram": [ + { + "marker": { + "line": { + "color": "white", + "width": 0.6 + } + }, + "type": "histogram" + } + ], + "mesh3d": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "mesh3d" + } + ], + "parcoords": [ + { + "line": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "parcoords" + } + ], + "pie": [ + { + "automargin": true, + "type": "pie" + } + ], + "scatter3d": [ + { + "line": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatter3d" + } + ], + "scattercarpet": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattercarpet" + } + ], + "scattergeo": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattergeo" + } + ], + "scattergl": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattergl" + } + ], + "scattermapbox": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattermapbox" + } + ], + "scattermap": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattermap" + } + ], + "scatterpolargl": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterpolargl" + } + ], + "scatterpolar": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterpolar" + } + ], + "scatter": [ + { + "fillpattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + }, + "type": "scatter" + } + ], + "scatterternary": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterternary" + } + ], + "surface": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "surface" + } + ], + "table": [ + { + "cells": { + "fill": { + "color": "rgb(237,237,237)" + }, + "line": { + "color": "white" + } + }, + "header": { + "fill": { + "color": "rgb(217,217,217)" + }, + "line": { + "color": "white" + } + }, + "type": "table" + } + ] + }, + "layout": { + "annotationdefaults": { + "arrowhead": 0, + "arrowwidth": 1 + }, + "autotypenumbers": "strict", + "coloraxis": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "colorscale": { + "diverging": [ + [ + 0.0, + "rgb(103,0,31)" + ], + [ + 0.1, + "rgb(178,24,43)" + ], + [ + 0.2, + "rgb(214,96,77)" + ], + [ + 0.3, + "rgb(244,165,130)" + ], + [ + 0.4, + "rgb(253,219,199)" + ], + [ + 0.5, + "rgb(247,247,247)" + ], + [ + 0.6, + "rgb(209,229,240)" + ], + [ + 0.7, + "rgb(146,197,222)" + ], + [ + 0.8, + "rgb(67,147,195)" + ], + [ + 0.9, + "rgb(33,102,172)" + ], + [ + 1.0, + "rgb(5,48,97)" + ] + ], + "sequential": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "sequentialminus": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ] + }, + "colorway": [ + "#1F77B4", + "#FF7F0E", + "#2CA02C", + "#D62728", + "#9467BD", + "#8C564B", + "#E377C2", + "#7F7F7F", + "#BCBD22", + "#17BECF" + ], + "font": { + "color": "rgb(36,36,36)" + }, + "geo": { + "bgcolor": "white", + "lakecolor": "white", + "landcolor": "white", + "showlakes": true, + "showland": true, + "subunitcolor": "white" + }, + "hoverlabel": { + "align": "left" + }, + "hovermode": "closest", + "mapbox": { + "style": "light" + }, + "paper_bgcolor": "white", + "plot_bgcolor": "white", + "polar": { + "angularaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "bgcolor": "white", + "radialaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + } + }, + "scene": { + "xaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "yaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "zaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + } + }, + "shapedefaults": { + "fillcolor": "black", + "line": { + "width": 0 + }, + "opacity": 0.3 + }, + "ternary": { + "aaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "baxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "bgcolor": "white", + "caxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + } + }, + "title": { + "x": 0.05 + }, + "xaxis": { + "automargin": true, + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside", + "title": { + "standoff": 15 + }, + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "yaxis": { + "automargin": true, + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside", + "title": { + "standoff": 15 + }, + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + } + } + }, + "margin": { + "l": 20, + "r": 20, + "t": 60, + "b": 20 + }, + "title": { + "text": "Karney - GHA1" + }, + "width": 800, + "height": 280 + }, + "config": { + "plotlyServerURL": "https://plot.ly" + } + } + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + }, + { + "data": { + "application/vnd.plotly.v1+json": { + "data": [ + { + "cells": { + "align": "center", + "values": [ + [ + "ES", + "approx", + "num" + ], + [ + "20", + "600", + "2000" + ], + [ + 0, + 0, + 2 + ], + [ + "1.07e+04", + "9.95e+03", + "nan" + ], + [ + "4.35e+03", + "2.77e+03", + "nan" + ], + [ + "0.000886", + "0.000888", + "nan" + ], + [ + "1.47", + "0.108", + "nan" + ] + ] + }, + "header": { + "align": "center", + "fill": { + "color": "lightgrey" + }, + "font": { + "size": 13 + }, + "values": [ + "Algorithmus", + "Parameter", + "NaN", + "dAlpha0 [\"]", + "dAlpha1 [\"]", + "dStrecke [m]", + "time [s]" + ] + }, + "type": "table" + } + ], + "layout": { + "template": { + "data": { + "barpolar": [ + { + "marker": { + "line": { + "color": "white", + "width": 0.5 + }, + "pattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + } + }, + "type": "barpolar" + } + ], + "bar": [ + { + "error_x": { + "color": "rgb(36,36,36)" + }, + "error_y": { + "color": "rgb(36,36,36)" + }, + "marker": { + "line": { + "color": "white", + "width": 0.5 + }, + "pattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + } + }, + "type": "bar" + } + ], + "carpet": [ + { + "aaxis": { + "endlinecolor": "rgb(36,36,36)", + "gridcolor": "white", + "linecolor": "white", + "minorgridcolor": "white", + "startlinecolor": "rgb(36,36,36)" + }, + "baxis": { + "endlinecolor": "rgb(36,36,36)", + "gridcolor": "white", + "linecolor": "white", + "minorgridcolor": "white", + "startlinecolor": "rgb(36,36,36)" + }, + "type": "carpet" + } + ], + "choropleth": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "choropleth" + } + ], + "contourcarpet": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "contourcarpet" + } + ], + "contour": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "contour" + } + ], + "heatmap": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "heatmap" + } + ], + "histogram2dcontour": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "histogram2dcontour" + } + ], + "histogram2d": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "histogram2d" + } + ], + "histogram": [ + { + "marker": { + "line": { + "color": "white", + "width": 0.6 + } + }, + "type": "histogram" + } + ], + "mesh3d": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "type": "mesh3d" + } + ], + "parcoords": [ + { + "line": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "parcoords" + } + ], + "pie": [ + { + "automargin": true, + "type": "pie" + } + ], + "scatter3d": [ + { + "line": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatter3d" + } + ], + "scattercarpet": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattercarpet" + } + ], + "scattergeo": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattergeo" + } + ], + "scattergl": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattergl" + } + ], + "scattermapbox": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattermapbox" + } + ], + "scattermap": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scattermap" + } + ], + "scatterpolargl": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterpolargl" + } + ], + "scatterpolar": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterpolar" + } + ], + "scatter": [ + { + "fillpattern": { + "fillmode": "overlay", + "size": 10, + "solidity": 0.2 + }, + "type": "scatter" + } + ], + "scatterternary": [ + { + "marker": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "type": "scatterternary" + } + ], + "surface": [ + { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + }, + "colorscale": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "type": "surface" + } + ], + "table": [ + { + "cells": { + "fill": { + "color": "rgb(237,237,237)" + }, + "line": { + "color": "white" + } + }, + "header": { + "fill": { + "color": "rgb(217,217,217)" + }, + "line": { + "color": "white" + } + }, + "type": "table" + } + ] + }, + "layout": { + "annotationdefaults": { + "arrowhead": 0, + "arrowwidth": 1 + }, + "autotypenumbers": "strict", + "coloraxis": { + "colorbar": { + "outlinewidth": 1, + "tickcolor": "rgb(36,36,36)", + "ticks": "outside" + } + }, + "colorscale": { + "diverging": [ + [ + 0.0, + "rgb(103,0,31)" + ], + [ + 0.1, + "rgb(178,24,43)" + ], + [ + 0.2, + "rgb(214,96,77)" + ], + [ + 0.3, + "rgb(244,165,130)" + ], + [ + 0.4, + "rgb(253,219,199)" + ], + [ + 0.5, + "rgb(247,247,247)" + ], + [ + 0.6, + "rgb(209,229,240)" + ], + [ + 0.7, + "rgb(146,197,222)" + ], + [ + 0.8, + "rgb(67,147,195)" + ], + [ + 0.9, + "rgb(33,102,172)" + ], + [ + 1.0, + "rgb(5,48,97)" + ] + ], + "sequential": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ], + "sequentialminus": [ + [ + 0.0, + "#440154" + ], + [ + 0.1111111111111111, + "#482878" + ], + [ + 0.2222222222222222, + "#3e4989" + ], + [ + 0.3333333333333333, + "#31688e" + ], + [ + 0.4444444444444444, + "#26828e" + ], + [ + 0.5555555555555556, + "#1f9e89" + ], + [ + 0.6666666666666666, + "#35b779" + ], + [ + 0.7777777777777778, + "#6ece58" + ], + [ + 0.8888888888888888, + "#b5de2b" + ], + [ + 1.0, + "#fde725" + ] + ] + }, + "colorway": [ + "#1F77B4", + "#FF7F0E", + "#2CA02C", + "#D62728", + "#9467BD", + "#8C564B", + "#E377C2", + "#7F7F7F", + "#BCBD22", + "#17BECF" + ], + "font": { + "color": "rgb(36,36,36)" + }, + "geo": { + "bgcolor": "white", + "lakecolor": "white", + "landcolor": "white", + "showlakes": true, + "showland": true, + "subunitcolor": "white" + }, + "hoverlabel": { + "align": "left" + }, + "hovermode": "closest", + "mapbox": { + "style": "light" + }, + "paper_bgcolor": "white", + "plot_bgcolor": "white", + "polar": { + "angularaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "bgcolor": "white", + "radialaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + } + }, + "scene": { + "xaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "yaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "zaxis": { + "backgroundcolor": "white", + "gridcolor": "rgb(232,232,232)", + "gridwidth": 2, + "linecolor": "rgb(36,36,36)", + "showbackground": true, + "showgrid": false, + "showline": true, + "ticks": "outside", + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + } + }, + "shapedefaults": { + "fillcolor": "black", + "line": { + "width": 0 + }, + "opacity": 0.3 + }, + "ternary": { + "aaxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "baxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + }, + "bgcolor": "white", + "caxis": { + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside" + } + }, + "title": { + "x": 0.05 + }, + "xaxis": { + "automargin": true, + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside", + "title": { + "standoff": 15 + }, + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + }, + "yaxis": { + "automargin": true, + "gridcolor": "rgb(232,232,232)", + "linecolor": "rgb(36,36,36)", + "showgrid": false, + "showline": true, + "ticks": "outside", + "title": { + "standoff": 15 + }, + "zeroline": false, + "zerolinecolor": "rgb(36,36,36)" + } + } + }, + "margin": { + "l": 20, + "r": 20, + "t": 60, + "b": 20 + }, + "title": { + "text": "Karney - GHA2" + }, + "width": 800, + "height": 280 + }, + "config": { + "plotlyServerURL": "https://plot.ly" + } + } + }, + "metadata": {}, + "output_type": "display_data", + "jetTransient": { + "display_id": null + } + } + ], + "execution_count": 10 + }, + { + "metadata": { + "ExecuteTime": { + "end_time": "2026-02-04T18:08:09.603687Z", + "start_time": "2026-02-04T18:08:09.269698Z" + } + }, + "cell_type": "code", + "source": [ + "import numpy as np\n", + "import math\n", + "from collections import defaultdict\n", + "\n", + "def to_latex_sci(x, sig=3):\n", + " \"\"\"\n", + " Formatiert eine Zahl als LaTeX, z.B. 8.02e-08 -> $8.02\\\\cdot10^{-8}$.\n", + " Für normale Zahlen -> $0.216$.\n", + " Für nan -> nan (ohne $...$).\n", + " \"\"\"\n", + " if x is None:\n", + " return \"nan\"\n", + " try:\n", + " xf = float(x)\n", + " except Exception:\n", + " return str(x)\n", + "\n", + " if math.isnan(xf):\n", + " return \"nan\"\n", + "\n", + " # nahe 0 -> 0\n", + " if xf == 0.0:\n", + " return \"$0$\"\n", + "\n", + " s = f\"{xf:.{sig}g}\" # z.B. '8.02e-08' oder '0.216'\n", + " if \"e\" in s or \"E\" in s:\n", + " mant, exp = s.lower().split(\"e\")\n", + " exp = int(exp)\n", + " return f\"${mant}\\\\cdot10^{{{exp}}}$\"\n", + " else:\n", + " return f\"${s}$\"\n", + "\n", + "\n", + "def nan_count_for_algorithm(results, example_keys, algorithm):\n", + " vals = []\n", + " for k in example_keys:\n", + " vals.extend(results[k][algorithm])\n", + " return int(np.sum(np.isnan(np.array(vals, dtype=float))))\n", + "\n", + "\n", + "def max_abs_metric(results, example_keys, algorithm, metric_index, is_angle=False, wu=None):\n", + " \"\"\"\n", + " Echter Max(|...|) über alle example_keys.\n", + " Wenn is_angle=True: Werte werden als rad angenommen und in Bogensekunden umgerechnet.\n", + " \"\"\"\n", + " arr = []\n", + " for k in example_keys:\n", + " arr.append(results[k][algorithm][metric_index])\n", + " arr = np.array(arr, dtype=float)\n", + "\n", + " if arr.size == 0:\n", + " return np.nan\n", + " m = np.nanmax(np.abs(arr))\n", + " if is_angle:\n", + " if wu is None:\n", + " raise ValueError(\"wu wird benötigt für rad2deg, wenn is_angle=True\")\n", + " m = wu.rad2deg(m) * 3600.0\n", + " return m\n", + "\n", + "\n", + "def parse_variant_params(algorithm_key):\n", + " \"\"\"\n", + " Erwartet: GHA1_ana_20_4 -> variant='ana', params='20_4'\n", + " Gibt zusätzlich eine hübsche Parameterdarstellung zurück.\n", + " \"\"\"\n", + " ghaNr, variant, params = algorithm_key.split(\"_\", 2)\n", + "\n", + " # Standard: params 그대로\n", + " pretty = params\n", + "\n", + " # Für deine Tabellen: ana hat z.B. '20_4' -> '4, 20' (Segmentgröße, Ordnung)\n", + " if variant == \"ana\":\n", + " # bei dir: params = '20_4' oder '50_8' -> Ordnung_Segment\n", + " # du willst aber: Segment, Ordnung -> '4, 20'\n", + " a, b = params.split(\"_\")\n", + " order = a\n", + " seg = b\n", + " pretty = f\"{seg}, {order}\"\n", + " else:\n", + " # num: '2000' bleibt '2000'\n", + " # approx/ES: oft eine Zahl mit Punkt/Unterstrich? -> 그대로\n", + " pretty = params.replace(\"_\", \", \")\n", + "\n", + " return variant, params, pretty\n", + "\n", + "\n", + "def build_latex_table_from_results(\n", + " results,\n", + " gha_prefix=\"GHA1\",\n", + " example_keys=None,\n", + " caption=\"\",\n", + " label=\"tab:results_algorithms\",\n", + " include_nan_col=False,\n", + " wu=None\n", + "):\n", + " \"\"\"\n", + " Erzeugt LaTeX Tabular (inkl. table-Umgebung) im gewünschten Stil.\n", + " \"\"\"\n", + "\n", + " if example_keys is None:\n", + " example_keys = list(results.keys())\n", + "\n", + " # Metriken & Winkel-Maske wie bei dir\n", + " if gha_prefix == \"GHA1\":\n", + " metric_headers = [\n", + " r\"$\\max(|\\Delta \\beta|)$ [$''$]\",\n", + " r\"$\\max(|\\Delta \\lambda|)$ [$''$]\",\n", + " r\"$\\max(|\\Delta \\alpha_1|)$ [$''$]\",\n", + " r\"time [s]\"\n", + " ]\n", + " angle_mask = [True, True, True, False]\n", + " else:\n", + " metric_headers = [\n", + " r\"$\\max(|\\Delta \\alpha_0|)$ [$''$]\",\n", + " r\"$\\max(|\\Delta \\alpha_1|)$ [$''$]\",\n", + " r\"$\\max(|\\Delta s|)$ [m]\",\n", + " r\"time [s]\"\n", + " ]\n", + " angle_mask = [True, True, False, False]\n", + "\n", + " # Alle Algorithmen sammeln\n", + " algorithms = sorted({\n", + " alg for k in example_keys\n", + " for alg in results[k].keys()\n", + " if alg.startswith(gha_prefix)\n", + " })\n", + "\n", + " # Gruppieren nach variant (ana, num, approx, ES, ...)\n", + " grouped = defaultdict(list)\n", + " for alg in algorithms:\n", + " variant, raw_params, pretty = parse_variant_params(alg)\n", + " grouped[variant].append((alg, pretty))\n", + "\n", + " # gewünschte Reihenfolge (falls vorhanden)\n", + " variant_order = [\"ana\", \"num\", \"approx\", \"ES\"]\n", + " ordered_variants = [v for v in variant_order if v in grouped] + [v for v in grouped.keys() if v not in variant_order]\n", + "\n", + " # LaTeX Header\n", + " cols = 2 + (1 if include_nan_col else 0) + len(metric_headers)\n", + " colspec = \"|\" + \"|\".join([\"c\"] * cols) + \"|\"\n", + "\n", + " header_cells = [\"Methode\", \"Parameterwerte\"]\n", + " if include_nan_col:\n", + " header_cells.append(\"NaN\")\n", + " header_cells += metric_headers\n", + "\n", + " lines = []\n", + " lines.append(r\"\\begin{table}[H]\")\n", + " lines.append(r\"\\centering\")\n", + " if caption:\n", + " lines.append(rf\"\\caption{{{caption}}}\")\n", + " lines.append(rf\"\\label{{{label}}}\")\n", + " lines.append(\"\")\n", + " lines.append(rf\"\\begin{{tabular}}{{{colspec}}}\")\n", + " lines.append(r\"\\hline\")\n", + " lines.append(\" & \".join(header_cells) + r\" \\\\\")\n", + " lines.append(r\"\\Xhline{1.5pt}\")\n", + "\n", + " # Zeilen bauen\n", + " for variant in ordered_variants:\n", + " rows = grouped[variant]\n", + " # für stable output: sort by pretty param (numerisch)\n", + " # (du kannst das ändern, wenn du eine andere Reihenfolge willst)\n", + " def sort_key(t):\n", + " _, pretty = t\n", + " # versuche numerisch zu sortieren\n", + " try:\n", + " parts = [float(p.strip()) for p in pretty.split(\",\")]\n", + " return parts\n", + " except Exception:\n", + " return [pretty]\n", + " rows = sorted(rows, key=sort_key)\n", + "\n", + " multi_n = len(rows)\n", + " method_name = {\n", + " \"ana\": \"analytisch\",\n", + " \"num\": \"numerisch\",\n", + " \"approx\": \"approximiert\"\n", + " }.get(variant, variant)\n", + "\n", + " for idx, (alg, pretty_param) in enumerate(rows):\n", + " row_cells = []\n", + "\n", + " if multi_n > 1:\n", + " if idx == 0:\n", + " row_cells.append(rf\"\\multirow{{{multi_n}}}{{*}}{{{method_name}}}\")\n", + " else:\n", + " row_cells.append(\"\") # Multirow fortsetzen\n", + " else:\n", + " row_cells.append(method_name)\n", + "\n", + " row_cells.append(pretty_param)\n", + "\n", + " if include_nan_col:\n", + " nans = nan_count_for_algorithm(results, example_keys, alg)\n", + " row_cells.append(str(nans))\n", + "\n", + " # Metriken\n", + " for mi in range(len(metric_headers)):\n", + " m = max_abs_metric(results, example_keys, alg, mi, is_angle=angle_mask[mi], wu=wu)\n", + " row_cells.append(to_latex_sci(m, sig=3))\n", + "\n", + " # Zeile + Linienlogik wie in deinem Beispiel\n", + " latex_row = \" & \".join(row_cells) + r\" \\\\\"\n", + "\n", + " # nach letzter Zeile einer Methode eine \\hline (wie bei dir)\n", + " if idx == multi_n - 1:\n", + " latex_row += r\"\\hline\"\n", + " lines.append(latex_row)\n", + "\n", + " lines.append(r\"\\end{tabular}\")\n", + " lines.append(\"\")\n", + " lines.append(r\"\\end{table}\")\n", + "\n", + " return \"\\n\".join(lines)\n", + "\n", + "\n", + "# --- Beispielaufruf ---\n", + "example_keys = list(results.keys()) # oder gefiltert auf Panou Gruppe etc.\n", + "latex = build_latex_table_from_results(\n", + " results,\n", + " gha_prefix=\"GHA2\",\n", + " example_keys=example_keys,\n", + " caption=r\"Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid\",\n", + " label=\"tab:results_algorithms\",\n", + " include_nan_col=False, # True, wenn du NaN-Spalte willst\n", + " wu=wu # wichtig, falls Winkelwerte rad sind\n", + ")\n", + "print(latex)" + ], + "id": "98de5e2a4f419483", + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\\begin{table}[H]\n", + "\\centering\n", + "\\caption{Ergebnisse der Lösungsmethoden der 1. \\gls{GHA} auf einem erdähnlichen Ellipsoid}\n", + "\\label{tab:results_algorithms}\n", + "\n", + "\\begin{tabular}{|c|c|c|c|c|c|}\n", + "\\hline\n", + "Methode & Parameterwerte & $\\max(|\\Delta \\alpha_0|)$ [$''$] & $\\max(|\\Delta \\alpha_1|)$ [$''$] & $\\max(|\\Delta s|)$ [m] & time [s] \\\\\n", + "\\Xhline{1.5pt}\n", + "numerisch & 2000 & nan & nan & nan & nan \\\\\\hline\n", + "approximiert & 600 & $9.95\\cdot10^{3}$ & $2.77\\cdot10^{3}$ & $0.000888$ & $0.108$ \\\\\\hline\n", + "ES & 20 & $1.07\\cdot10^{4}$ & $4.35\\cdot10^{3}$ & $0.000886$ & $1.47$ \\\\\\hline\n", + "\\end{tabular}\n", + "\n", + "\\end{table}\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "C:\\Users\\moell\\AppData\\Local\\Temp\\ipykernel_25588\\1918125892.py:53: RuntimeWarning:\n", + "\n", + "All-NaN slice encountered\n", + "\n" + ] + } + ], + "execution_count": 15 + }, + { + "metadata": {}, + "cell_type": "code", + "outputs": [], + "execution_count": null, + "source": "", + "id": "a431f99070cdee3c" + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/alpha_conversion_test.ipynb b/Tests/alpha_conversion_test.ipynb similarity index 98% rename from alpha_conversion_test.ipynb rename to Tests/alpha_conversion_test.ipynb index 5657bd8..62c501a 100644 --- a/alpha_conversion_test.ipynb +++ b/Tests/alpha_conversion_test.ipynb @@ -30,7 +30,7 @@ "%autoreload 2\n", "import winkelumrechnungen as wu\n", "from ellipsoide import EllipsoidTriaxial\n", - "from GHA_triaxial.panou import alpha_ell2para, alpha_para2ell\n", + "from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n", "import numpy as np" ], "id": "9ad815aea55574e3", diff --git a/conversions_test.ipynb b/Tests/conversions_test.ipynb similarity index 100% rename from conversions_test.ipynb rename to Tests/conversions_test.ipynb diff --git a/Tests/test_biaxial.py b/Tests/test_biaxial.py new file mode 100644 index 0000000..a36e45b --- /dev/null +++ b/Tests/test_biaxial.py @@ -0,0 +1,34 @@ +import numpy as np +from ellipsoide import EllipsoidBiaxial +from GHA_biaxial.bessel import gha1 as gha1_bessel +from GHA_biaxial.gauss import gha1 as gha1_gauss +from GHA_biaxial.rk import gha1 as gha1_rk +from GHA_biaxial.gauss import gha2 as gha2_gauss + +re = EllipsoidBiaxial.init_name("Bessel") + +# phi0 = 0.6 +# lamb0 = 1.2 +# alpha0 = 0.45 +# s = 123456 +# +# values_bessel = gha1_bessel(re, phi0, lamb0, alpha0, s) +# alpha1_bessel = values_bessel[-1] +# p1_bessel = re.bi_ell2cart(values_bessel[0], values_bessel[1], 0) +# +# values_gauss1 = gha1_gauss(re, phi0, lamb0, alpha0, s) +# alpha1_gauss1 = values_gauss1[-1] +# p1_gauss = re.bi_ell2cart(values_gauss1[0], values_gauss1[1], 0) +# +# values_rk = gha1_rk(re, phi0, lamb0 , alpha0, s, 10000) +# alpha1_rk = values_rk[-1] +# p1_rk = re.bi_ell2cart(values_rk[0], values_rk[1], 0) +# +# alpha0_gauss, alpha1_gauss2, s_gauss = gha2_gauss(re, phi0, lamb0, values_gauss1[0], values_gauss1[1]) + +phi0 = 0.6 +lamb0 = 1.2 + +cart = re.bi_ell2cart(phi0, lamb0, 0) +ell = re.bi_cart2ell(cart) +pass \ No newline at end of file diff --git a/algorithms_test.ipynb b/algorithms_test.ipynb deleted file mode 100644 index d5738b6..0000000 --- a/algorithms_test.ipynb +++ /dev/null @@ -1,363 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "id": "initial_id", - "metadata": { - "collapsed": true - }, - "source": [ - "%load_ext autoreload\n", - "%autoreload 2" - ], - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "%reload_ext autoreload\n", - "%autoreload 2\n", - "import time\n", - "import pickle\n", - "import numpy as np\n", - "from numpy import nan\n", - "import winkelumrechnungen as wu\n", - "import os\n", - "from contextlib import contextmanager, redirect_stdout, redirect_stderr\n", - "import pandas as pd\n", - "import plotly.graph_objects as go\n", - "import math\n", - "\n", - "from ellipsoide import EllipsoidTriaxial\n", - "from GHA_triaxial.panou import alpha_ell2para, alpha_para2ell\n", - "\n", - "from GHA_triaxial.panou import gha1_num, gha1_ana\n", - "from GHA_triaxial.approx_gha1 import gha1_approx\n", - "\n", - "from GHA_triaxial.panou_2013_2GHA_num import gha2_num\n", - "from GHA_triaxial.ES_gha2 import gha2_ES\n", - "from GHA_triaxial.approx_gha2 import gha2_approx\n", - "\n", - "from GHA_triaxial.numeric_examples_panou import get_tables as get_tables_panou\n", - "from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney" - ], - "id": "961cb22764c5bcb9", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "@contextmanager\n", - "def suppress_print():\n", - " with open(os.devnull, 'w') as fnull:\n", - " with redirect_stdout(fnull), redirect_stderr(fnull):\n", - " yield" - ], - "id": "e4740625a366ac13", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "# steps_gha1_num = [2000, 5000, 10000, 20000]\n", - "# maxM_gha1_ana = [20, 50]\n", - "# parts_gha1_ana = [4, 8]\n", - "# dsPart_gha1_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", - "#\n", - "# steps_gha2_num = [2000, 5000, 10000, 20000]\n", - "# dsPart_gha2_ES = [600, 1250, 6000] # entspricht bei der Erde ca. 10000, 5000, 1000\n", - "# dsPart_gha2_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n", - "\n", - "steps_gha1_num = [2000]\n", - "maxM_gha1_ana = [10, 20, 40, 60, 80]\n", - "parts_gha1_ana = [4, 8, 12, 16]\n", - "dsPart_gha1_approx = [600]\n", - "\n", - "steps_gha2_num = [16000]\n", - "dsPart_gha2_ES = [20]\n", - "dsPart_gha2_approx = [600]" - ], - "id": "96093cdde03f8d57", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "# test = \"Karney\"\n", - "test = \"Panou\"\n", - "if test == \"Karney\":\n", - " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n", - " examples = get_examples_karney(10, 42)\n", - "elif test == \"Panou\":\n", - " ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n", - " tables = get_tables_panou()\n", - " table_indices = []\n", - " examples = []\n", - " for i, table in enumerate(tables):\n", - " for example in table:\n", - " table_indices.append(i+1)\n", - " examples.append(example)" - ], - "id": "6e384cc01c2dbe", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "results = {}\n", - "for example in examples:\n", - " example_results = {}\n", - "\n", - " beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n", - " P0 = ell.ell2cart(beta0, lamb0)\n", - " P1 = ell.ell2cart(beta1, lamb1)\n", - " _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n", - "\n", - " # for steps in steps_gha1_num:\n", - " # start = time.perf_counter()\n", - " # try:\n", - " # P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n", - " # end = time.perf_counter()\n", - " # beta1_num, lamb1_num = ell.cart2ell(P1_num)\n", - " # d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600\n", - " # d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600\n", - " # d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600\n", - " # d_time = end - start\n", - " # example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", - " # except Exception as e:\n", - " # print(e)\n", - " # example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n", - " #\n", - " for maxM in maxM_gha1_ana:\n", - " for parts in parts_gha1_ana:\n", - " start = time.perf_counter()\n", - " try:\n", - " P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n", - " end = time.perf_counter()\n", - " beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n", - " _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n", - " d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n", - " d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n", - " d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n", - " d_time = end - start\n", - " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", - " except Exception as e:\n", - " print(e)\n", - " example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n", - " #\n", - " # for dsPart in dsPart_gha1_approx:\n", - " # ds = ell.ax/dsPart\n", - " # start = time.perf_counter()\n", - " # try:\n", - " # P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n", - " # end = time.perf_counter()\n", - " # beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n", - " # d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n", - " # d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n", - " # d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", - " # d_time = end - start\n", - " # example_results[f\"GHA1_approx_{ds:.3f}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n", - " # except Exception as e:\n", - " # print(e)\n", - " # example_results[f\"GHA1_approx_{ds:.3f}\"] = (nan, nan, nan, nan)\n", - "\n", - " # for steps in steps_gha2_num:\n", - " # start = time.perf_counter()\n", - " # try:\n", - " # alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n", - " # end = time.perf_counter()\n", - " # d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n", - " # d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n", - " # d_s = abs(s_num - s) / 1000\n", - " # d_time = end - start\n", - " # example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " # print(\"konvergiert\")\n", - " # except Exception as e:\n", - " # print(e)\n", - " # print(example)\n", - " # example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n", - "\n", - " # for dsPart in dsPart_gha2_ES:\n", - " # ds = ell.ax/dsPart\n", - " # start = time.perf_counter()\n", - " # try:\n", - " # with suppress_print():\n", - " # alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n", - " # end = time.perf_counter()\n", - " # d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n", - " # d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n", - " # d_s = abs(s_ES - s) / 1000\n", - " # d_time = end - start\n", - " # example_results[f\"GHA2_ES_{ds:.3f}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " # except Exception as e:\n", - " # print(e)\n", - " # example_results[f\"GHA2_ES_{ds:.3f}\"] = (nan, nan, nan, nan)\n", - " #\n", - " # for dsPart in dsPart_gha2_approx:\n", - " # ds = ell.ax/dsPart\n", - " # start = time.perf_counter()\n", - " # try:\n", - " # alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n", - " # end = time.perf_counter()\n", - " # d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n", - " # d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n", - " # d_s = abs(s_approx - s) / 1000\n", - " # d_time = end - start\n", - " # example_results[f\"GHA2_approx_{ds:.3f}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n", - " # except Exception as e:\n", - " # print(e)\n", - " # example_results[f\"GHA2_approx_{ds:.3f}\"] = (nan, nan, nan, nan)\n", - "\n", - " results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results" - ], - "id": "ef8849908ed3231e", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "# with open(f\"gha_results{test}.pkl\", \"wb\") as f:\n", - "# pickle.dump(results, f)" - ], - "id": "664105ea35d50a7b", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "# with open(f\"gha_results{test}.pkl\", \"rb\") as f:\n", - "# results = pickle.load(f)" - ], - "id": "ad8aa9dcc8af4a05", - "outputs": [], - "execution_count": null - }, - { - "metadata": {}, - "cell_type": "code", - "source": [ - "def format_max(values, is_angle=False):\n", - " arr = np.array(values, dtype=float)\n", - " if arr.size==0:\n", - " return np.nan\n", - " maxi = np.nanmax(np.abs(arr))\n", - " if maxi is None or (isinstance(maxi,float) and (math.isnan(maxi))):\n", - " return \"nan\"\n", - " i = 2\n", - " while maxi > np.nanmean(np.abs(arr)):\n", - " maxi = np.sort(np.abs(arr[~np.isnan(arr)]))[-i]\n", - " i += 1\n", - " if is_angle:\n", - " maxi = wu.rad2deg(maxi)*3600\n", - " if f\"{maxi:.3g}\" == 0:\n", - " pass\n", - " return f\"{maxi:.3g}\"\n", - "\n", - "\n", - "def build_max_table(gha_prefix, title, group_value = None):\n", - " # metrics\n", - " if gha_prefix==\"GHA1\":\n", - " metrics = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n", - " angle_mask = [True, True, True, False]\n", - " else:\n", - " metrics = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n", - " angle_mask = [True, True, False, False]\n", - "\n", - " # collect keys in this group\n", - " if group_value is None:\n", - " example_keys = [example_key for example_key in results.keys()]\n", - " else:\n", - " example_keys = [example_key for example_key, group_index in zip(results.keys(), table_indices) if group_index==group_value]\n", - " # all variant keys (inner dict keys) matching prefix\n", - " algorithms = sorted({algorithm for example_key in example_keys for algorithm in results[example_key].keys() if algorithm.startswith(gha_prefix)})\n", - "\n", - " header = [\"Algorithmus\", \"Parameter\"] + list(metrics)\n", - " cells = [[] for i in range(len(metrics) + 2)]\n", - " for algorithm in algorithms:\n", - " if algorithm == \"GHA1_ana_80_16\":\n", - " pass\n", - " ghaNr, variant, params = algorithm.split(\"_\", 2)\n", - " cells[0].append(variant)\n", - " cells[1].append(params)\n", - " for i, metric in enumerate(metrics):\n", - " values = []\n", - " for example_key in example_keys:\n", - " values.append(results[example_key][algorithm][i])\n", - " cells[i+2].append(format_max(values, is_angle=angle_mask[i]))\n", - "\n", - " header = dict(\n", - " values=header,\n", - " align=\"center\",\n", - " fill_color=\"lightgrey\",\n", - " font=dict(size=13)\n", - " )\n", - " cells = dict(\n", - " values=cells,\n", - " align=\"center\"\n", - " )\n", - "\n", - " fig = go.Figure(data=[go.Table(header=header, cells=cells)])\n", - " fig.update_layout(title=title,\n", - " template=\"simple_white\",\n", - " width=800,\n", - " height=280,\n", - " margin=dict(l=20, r=20, t=60, b=20))\n", - " return fig\n", - "\n", - "figs = []\n", - "if test == \"Panou\":\n", - " for table_index in sorted(set(table_indices))[:-1]:\n", - " fig1 = build_max_table(\"GHA1\", f\"{test} - Gruppe {table_index} - GHA1\", table_index)\n", - " fig2 = build_max_table(\"GHA2\", f\"{test} - Gruppe {table_index} - GHA2\", table_index)\n", - " figs.append(fig1)\n", - " figs.append(fig2)\n", - "elif test == \"Karney\":\n", - " fig1 = build_max_table(\"GHA1\", f\"{test} - GHA1\")\n", - " fig2 = build_max_table(\"GHA2\", f\"{test} - GHA2\")\n", - " figs.append(fig1)\n", - " figs.append(fig2)\n", - "\n", - "for fig in figs:\n", - " fig.show()" - ], - "id": "b46d57fc0d794e28", - "outputs": [], - "execution_count": null - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 2 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.6" - } - }, - "nbformat": 4, - "nbformat_minor": 5 -} diff --git a/algorithms_test.py b/algorithms_test.py deleted file mode 100644 index 7320bd8..0000000 --- a/algorithms_test.py +++ /dev/null @@ -1,267 +0,0 @@ -import time -import pickle -import numpy as np -from numpy import nan -import winkelumrechnungen as wu -import os -from contextlib import contextmanager, redirect_stdout, redirect_stderr - -from ellipsoide import EllipsoidTriaxial -from GHA_triaxial.panou import alpha_ell2para, alpha_para2ell - -from GHA_triaxial.panou import gha1_num, gha1_ana -from GHA_triaxial.approx_gha1 import gha1_approx - -from GHA_triaxial.panou_2013_2GHA_num import gha2_num -from GHA_triaxial.ES_gha2 import gha2_ES -from GHA_triaxial.approx_gha2 import gha2_approx - -from GHA_triaxial.numeric_examples_panou import get_random_examples as get_examples_panou -from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney - -@contextmanager -def suppress_print(): - with open(os.devnull, 'w') as fnull: - with redirect_stdout(fnull), redirect_stderr(fnull): - yield - -# steps_gha1_num = [2000, 5000, 10000, 20000] -# maxM_gha1_ana = [20, 40, 60] -# parts_gha1_ana = [4, 8, 16] -# dsPart_gha1_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100 -# -# steps_gha2_num = [2000, 5000, 10000, 20000] -# dsPart_gha2_ES = [600, 1250, 6000] # entspricht bei der Erde ca. 10000, 5000, 1000 -# dsPart_gha2_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100 - -steps_gha1_num = [2000, 5000] -maxM_gha1_ana = [20, 40] -parts_gha1_ana = [4, 8] -dsPart_gha1_approx = [600, 1250] - -steps_gha2_num = [2000, 5000] -dsPart_gha2_ES = [20] -dsPart_gha2_approx = [600, 1250] - -ell_karney: EllipsoidTriaxial = EllipsoidTriaxial.init_name("KarneyTest2024") -ell_panou: EllipsoidTriaxial = EllipsoidTriaxial.init_name("BursaSima1980round") - -results_karney = {} -results_panou = {} - -examples_karney = get_examples_karney(2, 42) -examples_panou = get_examples_panou(2, 42) - -for example in examples_karney: - example_results = {} - - beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example - P0 = ell_karney.ell2cart(beta0, lamb0) - P1 = ell_karney.ell2cart(beta1, lamb1) - _, _, alpha0_para = alpha_ell2para(ell_karney, beta0, lamb0, alpha0_ell) - - for steps in steps_gha1_num: - start = time.perf_counter() - try: - P1_num, alpha1_num_1 = gha1_num(ell_karney, P0, alpha0_ell, s, num=steps) - end = time.perf_counter() - beta1_num, lamb1_num = ell_karney.cart2ell(P1_num) - d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_num_{steps}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_num_{steps}"] = (nan, nan, nan, nan) - - for maxM in maxM_gha1_ana: - for parts in parts_gha1_ana: - start = time.perf_counter() - try: - P1_ana, alpha1_ana_para = gha1_ana(ell_karney, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts) - end = time.perf_counter() - beta1_ana, lamb1_ana = ell_karney.cart2ell(P1_ana) - _, _, alpha1_ana_ell = alpha_para2ell(ell_karney, beta1_ana, lamb1_ana, alpha1_ana_para) - d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_ana_{maxM}_{parts}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_ana_{maxM}_{parts}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha1_approx: - ds = ell_karney.ax/dsPart - start = time.perf_counter() - try: - P1_approx, alpha1_approx = gha1_approx(ell_karney, P0, alpha0_ell, s, ds=ds) - end = time.perf_counter() - beta1_approx, lamb1_approx = ell_karney.cart2ell(P1_approx) - d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_approx_{ds:.3f}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_approx_{ds:.3f}"] = (nan, nan, nan, nan) - - for steps in steps_gha2_num: - start = time.perf_counter() - try: - alpha0_num, alpha1_num_2, s_num = gha2_num(ell_karney, beta0, lamb0, beta1, lamb1, n=steps) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600 - d_s = abs(s_num - s) / 1000 - d_time = end - start - example_results[f"GHA2_num_{steps}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_num_{steps}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha2_ES: - ds = ell_karney.ax/dsPart - start = time.perf_counter() - try: - with suppress_print(): - alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell_karney, P0, P1, stepLenTarget=ds) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600 - d_s = abs(s_ES - s) / 1000 - d_time = end - start - example_results[f"GHA2_ES_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_ES_{ds:.3f}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha2_approx: - ds = ell_karney.ax/dsPart - start = time.perf_counter() - try: - alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell_karney, P0, P1, ds=ds) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600 - d_s = abs(s_approx - s) / 1000 - d_time = end - start - example_results[f"GHA2_approx_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_approx_{ds:.3f}"] = (nan, nan, nan, nan) - - results_karney[f"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}"] = example_results - -for example in examples_panou: - example_results = {} - - beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example - P0 = ell_panou.ell2cart(beta0, lamb0) - P1 = ell_panou.ell2cart(beta1, lamb1) - _, _, alpha0_para = alpha_ell2para(ell_panou, beta0, lamb0, alpha0_ell) - - for steps in steps_gha1_num: - start = time.perf_counter() - try: - P1_num, alpha1_num_1 = gha1_num(ell_panou, P0, alpha0_ell, s, num=steps) - end = time.perf_counter() - beta1_num, lamb1_num = ell_panou.cart2ell(P1_num) - d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_num_{steps}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_num_{steps}"] = (nan, nan, nan, nan) - - for maxM in maxM_gha1_ana: - for parts in parts_gha1_ana: - start = time.perf_counter() - try: - P1_ana, alpha1_ana_para = gha1_ana(ell_panou, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts) - end = time.perf_counter() - beta1_ana, lamb1_ana = ell_panou.cart2ell(P1_ana) - _, _, alpha1_ana = alpha_para2ell(ell_panou, beta1_ana, lamb1_ana, alpha1_ana_para) - d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_ana - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_ana_{maxM}_{parts}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_ana_{maxM}_{parts}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha1_approx: - ds = ell_panou.ax/dsPart - start = time.perf_counter() - try: - P1_approx, alpha1_approx = gha1_approx(ell_panou, P0, alpha0_ell, s, ds=ds) - end = time.perf_counter() - beta1_approx, lamb1_approx = ell_panou.cart2ell(P1_approx) - d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600 - d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600 - d_time = end - start - example_results[f"GHA1_approx_{ds:.3f}"] = (d_beta1, d_lamb1, d_alpha1, d_time) - except Exception as e: - print(e) - example_results[f"GHA1_approx_{ds:.3f}"] = (nan, nan, nan, nan) - - for steps in steps_gha2_num: - start = time.perf_counter() - try: - alpha0_num, alpha1_num_2, s_num = gha2_num(ell_panou, beta0, lamb0, beta1, lamb1, n=steps) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600 - d_s = abs(s_num - s) / 1000 - d_time = end - start - example_results[f"GHA2_num_{steps}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_num_{steps}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha2_ES: - ds = ell_panou.ax/dsPart - start = time.perf_counter() - try: - with suppress_print(): - alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell_panou, P0, P1, stepLenTarget=ds) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600 - d_s = abs(s_ES - s) / 1000 - d_time = end - start - example_results[f"GHA2_ES_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_ES_{ds:.3f}"] = (nan, nan, nan, nan) - - for dsPart in dsPart_gha2_approx: - ds = ell_panou.ax/dsPart - start = time.perf_counter() - try: - alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell_panou, P0, P1, ds=ds) - end = time.perf_counter() - d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600 - d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600 - d_s = abs(s_approx - s) / 1000 - d_time = end - start - example_results[f"GHA2_approx_{ds:.3f}"] = (d_alpha0, d_alpha1, d_s, d_time) - except Exception as e: - print(e) - example_results[f"GHA2_approx_{ds:.3f}"] = (nan, nan, nan, nan) - - results_panou[f"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}"] = example_results - -print(results_karney) -with open("results_karney.pkl", "wb") as f: - pickle.dump(results_karney, f) - -print(results_panou) -with open("results_panou.pkl", "wb") as f: - pickle.dump(results_panou, f) \ No newline at end of file diff --git a/conversions_test.py b/conversions_test.py deleted file mode 100644 index f545fcf..0000000 --- a/conversions_test.py +++ /dev/null @@ -1,106 +0,0 @@ -import time -import pickle -import numpy as np -from numpy import nan -import winkelumrechnungen as wu -import os -from contextlib import contextmanager, redirect_stdout, redirect_stderr -from itertools import product -import pandas as pd - -from ellipsoide import EllipsoidTriaxial - -# ellips = "KarneyTest2024" -# ellips = "BursaSima1980" -ellips = "Fiction" - -ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(ellips) - - -def deg_range(start, stop, step): - return [float(x) for x in range(start, stop + step, step)] - -def asymptotic_range(start, direction="up", max_decimals=4): - values = [] - for d in range(0, max_decimals + 1): - step = 10 ** -d - if direction == "up": - values.append(start + (1 - step)) - else: - values.append(start - (1 - step)) - return values - - -beta_5_85 = deg_range(5, 85, 5) -lambda_5_85 = deg_range(5, 85, 5) - -beta_5_90 = deg_range(5, 90, 5) -lambda_5_90 = deg_range(5, 90, 5) - -beta_0_90 = deg_range(0, 90, 5) -lambda_0_90 = deg_range(0, 90, 5) - -beta_90 = [90.0] -lambda_90 = [90.0] -beta_0 = [0.0] -lambda_0 = [0.0] - -beta_asym_89 = asymptotic_range(89.0, direction="up") -lambda_asym_0 = asymptotic_range(1.0, direction="down") - -groups = { - 1: list(product(beta_5_85, lambda_5_85)), - # 2: list(product(beta_0, lambda_0_90)), - # 3: list(product(beta_5_85, lambda_0)), - # 4: list(product(beta_90, lambda_5_90)), - # 5: list(product(beta_asym_89, lambda_asym_0)), - # 6: list(product(beta_5_85, lambda_90)), - 7: list(product(lambda_asym_0, lambda_0_90)), - # 8: list(product(beta_0_90, lambda_asym_0)), - # 9: list(product(beta_asym_89, lambda_0_90)), - # 10: list(product(beta_0_90, beta_asym_89)), -} - -for nr, points in groups.items(): - points_cart = [] - for point in points: - beta, lamb = point - cart = ell.ell2cart(wu.deg2rad(beta), wu.deg2rad(lamb)) - points_cart.append(cart) - groups[nr] = points_cart - -results = {} - -for nr, points in groups.items(): - group_results = {"ell": [], - "para": [], - "geod": []} - for point in points: - elli = ell.cart2ell(point) - cart_elli = ell.ell2cart(elli[0], elli[1]) - group_results["ell"].append(np.linalg.norm(point - cart_elli, axis=-1)) - - para = ell.cart2para(point) - cart_para = ell.para2cart(para[0], para[1]) - group_results["para"].append(np.linalg.norm(point - cart_para, axis=-1)) - - geod = ell.cart2geod(point, "ligas3") - cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - group_results["geod"].append(np.linalg.norm(point - cart_geod, axis=-1)) - - group_results["ell"] = np.array(group_results["ell"]) - group_results["para"] = np.array(group_results["para"]) - group_results["geod"] = np.array(group_results["geod"]) - results[nr] = group_results - -with open(f"conversion_results_{ellips}.pkl", "wb") as f: - pickle.dump(results, f) - -df = pd.DataFrame({ - "Gruppe": [nr for nr in results.keys()], - "max_Δr_ell": [f"{max(result["ell"]):.3g}" for result in results.values()], - "max_Δr_para": [f"{max(result["para"]):.3g}" for result in results.values()], - "max_Δr_geod": [f"{max(result["geod"]):.3g}" for result in results.values()] -}) - -print(df) diff --git a/dashboard.py b/dashboard.py index 3681ecc..71adcf4 100644 --- a/dashboard.py +++ b/dashboard.py @@ -6,13 +6,13 @@ from ellipsoide import EllipsoidTriaxial import winkelumrechnungen as wu import ausgaben as aus -from GHA_triaxial.panou import gha1_ana -from GHA_triaxial.panou import gha1_num -from GHA_triaxial.approx_gha1 import gha1_approx +from GHA_triaxial.gha1_ana import gha1_ana +from GHA_triaxial.gha1_num import gha1_num +from GHA_triaxial.gha1_approx import gha1_approx -from GHA_triaxial.panou_2013_2GHA_num import gha2_num -from GHA_triaxial.ES_gha2 import gha2_ES -from GHA_triaxial.approx_gha2 import gha2_approx +from GHA_triaxial.gha2_num import gha2_num +from GHA_triaxial.gha2_ES import gha2_ES +from GHA_triaxial.gha2_approx import gha2_approx app = Dash(__name__, suppress_callback_exceptions=True) @@ -715,7 +715,7 @@ def compute_gha2_approx(n2, beta0, lamb0, beta1, lamb1, ax, ay, b, method2): P0 = ell.ell2cart(beta0_rad, lamb0_rad) P1 = ell.ell2cart(beta1_rad, lamb1_rad) - a0_app, a1_app, s_app, points = gha2_approx(ell, P0, P1, ds=1e-4, all_points=True) + a0_app, a1_app, s_app, points = gha2_approx(ell, P0, P1, ds=1000, all_points=True) out = html.Div([ html.Strong("Approximiert: "), diff --git a/ellipsoide.py b/ellipsoide.py index bb3aa88..189b49b 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -59,6 +59,14 @@ class EllipsoidBiaxial: phi2p = lambda self, phi: self.N(phi) * cos(phi) def bi_cart2ell(self, point: NDArrayself, Eh: float = 0.001, Ephi: float = wu.gms2rad([0, 0, 0.001])) -> Tuple[float, float, float]: + """ + Umrechnung von kartesischen in ellipsoidische Koordinaten auf einem Rotationsellipsoid + # TODO: Quelle + :param point: Punkt in kartesischen Koordinaten + :param Eh: Grenzwert für die Höhe + :param Ephi: Grenzwert für die Breite + :return: ellipsoidische Breite, Länge, geodätische Höhe + """ x, y, z = point lamb = arctan2(y, x) @@ -87,6 +95,14 @@ class EllipsoidBiaxial: return phi, lamb, h def bi_ell2cart(self, phi: float, lamb: float, h: float) -> NDArray: + """ + Umrechnung von ellipsoidischen in kartesische Koordinaten auf einem Rotationsellipsoid + # TODO: Quelle + :param phi: ellipsoidische Breite + :param lamb: ellipsoidische Länge + :param h: geodätische Höhe + :return: Punkt in kartesischen Koordinaten + """ W = sqrt(1 - self.e**2 * sin(phi)**2) N = self.a / W x = (N+h) * cos(phi) * cos(lamb) @@ -112,7 +128,8 @@ class EllipsoidTriaxial: @classmethod def init_name(cls, name: str): """ - Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972, Bursa1970, BesselBiaxial + Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, BursaSima1980round, Eitschberger1978, Bursa1972, + Bursa1970, BesselBiaxial, Fiction, KarneyTest2024 Panou et al (2020) :param name: Name des dreiachsigen Ellipsoids """ @@ -213,12 +230,9 @@ class EllipsoidTriaxial: def ellu2cart(self, beta: float, lamb: float, u: float) -> NDArray: """ Panou 2014 12ff. - Elliptische Breite+Länge sind nicht gleich der geodätischen - Verhältnisse des Ellipsoids bekannt, Größe verändern bis Punkt erreicht, - dann ist u die Größe entlang der z-Achse - :param beta: ellipsoidische Breite [rad] - :param lamb: ellipsoidische Länge [rad] - :param u: Größe entlang der z-Achse + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge + :param u: radiale Koordinate entlang der kleinen Halbachse :return: Punkt in kartesischen Koordinaten """ x = sqrt(u**2 + self.Ex**2) * sqrt(cos(beta)**2 + self.Ee**2/self.Ex**2 * sin(beta)**2) * cos(lamb) @@ -231,7 +245,7 @@ class EllipsoidTriaxial: """ Panou 2014 15ff. :param point: Punkt in kartesischen Koordinaten - :return: elliptische Breite, elliptische Länge, Größe entlang der z-Achse + :return: ellipsoidische Breite, ellipsoidische Länge, radiale Koordinate entlang der kleinen Halbachse """ x, y, z = point c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2 @@ -261,8 +275,8 @@ class EllipsoidTriaxial: def ell2cart(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray: """ Panou, Korakitis 2019 2 - :param beta: elliptische Breite [rad] - :param lamb: elliptische Länge [rad] + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge :return: Punkt in kartesischen Koordinaten """ beta = np.asarray(beta, dtype=float) @@ -297,8 +311,8 @@ class EllipsoidTriaxial: def ell2cart_bektas(self, beta: float | NDArray, omega: float | NDArray) -> NDArray: """ Bektas 2015 - :param beta: elliptische Breite [rad] - :param omega: elliptische Länge [rad] + :param beta: ellipsoidische Breite + :param omega: ellipsoidische Länge :return: Punkt in kartesischen Koordinaten """ x = self.ax * cos(omega) * sqrt((self.ax**2 - self.ay**2 * sin(beta)**2 - self.b**2 * cos(beta)**2) / (self.ax**2 - self.b**2)) @@ -310,8 +324,8 @@ class EllipsoidTriaxial: def ell2cart_karney(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray: """ Karney 2025 Geographic Lib - :param beta: elliptische Breite [rad] - :param lamb: elliptische Länge [rad] + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge :return: Punkt in kartesischen Koordinaten """ k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2) @@ -321,11 +335,11 @@ class EllipsoidTriaxial: Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2) return np.array([X, Y, Z]) - def cart2ell_yFake(self, point: NDArray, start_delta) -> Tuple[float, float]: + def cart2ell_yFake(self, point: NDArray) -> Tuple[float, float]: """ - - :param point: - :return: + Bei Fehlschlagen von cart2ell + :param point: Punkt in kartesischen Koordinaten + :return: ellipsoidische Breite und Länge """ x, y, z = point delta_y = 1e-4 @@ -363,8 +377,8 @@ class EllipsoidTriaxial: :param point: Punkt in kartesischen Koordinaten :param eps: zu erreichende Genauigkeit :param maxI: maximale Anzahl Iterationen - :param noFake: - :return: elliptische Breite und Länge [rad] + :param noFake: y numerisch anpassen? + :return: ellipsoidische Breite und Länge """ x, y, z = point beta, lamb = self.cart2ell_panou(point) @@ -408,6 +422,7 @@ class EllipsoidTriaxial: return beta, lamb except Exception as e: + # Wenn die Berechnung fehlschlägt auf Grund von sehr kleinem y, solange anpassen, bis Umrechnung ohne Fehler delta_y = 10 ** math.floor(math.log10(abs(self.ay/1000))) if abs(y) < delta_y and not noFake: return self.cart2ell_yFake(point, delta_y) @@ -418,7 +433,7 @@ class EllipsoidTriaxial: """ Panou, Korakitis 2019 2f. (analytisch -> Näherung) :param point: Punkt in kartesischen Koordinaten - :return: elliptische Breite, elliptische Länge + :return: ellipsoidische Breite und Länge """ x, y, z = point @@ -477,7 +492,7 @@ class EllipsoidTriaxial: :param point: Punkt in kartesischen Koordinaten :param eps: zu erreichende Genauigkeit :param maxI: maximale Anzahl Iterationen - :return: elliptische Breite und Länge [rad] + :return: ellipsoidische Breite und Länge """ x, y, z = point phi, lamb = self.cart2para(point) @@ -503,9 +518,9 @@ class EllipsoidTriaxial: def geod2cart(self, phi: float | NDArray, lamb: float | NDArray, h: float) -> NDArray: """ Ligas 2012, 250 - :param phi: geodätische Breite [rad] - :param lamb: geodätische Länge [rad] - :param h: Höhe über dem Ellipsoid + :param phi: geodätische Breite + :param lamb: geodätische Länge + :param h: geodätische Höhe :return: kartesische Koordinaten """ v = self.ax / sqrt(1 - self.ex**2*sin(phi)**2-self.ee**2*cos(phi)**2*sin(lamb)**2) @@ -521,7 +536,7 @@ class EllipsoidTriaxial: :param point: Punkt in kartesischen Koordinaten :param maxIter: maximale Anzahl Iterationen :param maxLoa: Level of Accuracy, das erreicht werden soll - :return: phi, lambda, h + :return: geodätische Breite, Länge, Höhe """ xG, yG, zG = point @@ -596,8 +611,8 @@ class EllipsoidTriaxial: def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray: """ Panou, Korakitits 2020, 4 - :param u: Parameter u - :param v: Parameter v + :param u: parametrische Breite + :param v: parametrische Länge :return: Punkt in kartesischen Koordinaten """ x = self.ax * cos(u) * cos(v) @@ -610,7 +625,7 @@ class EllipsoidTriaxial: """ Panou, Korakitits 2020, 4 :param point: Punkt in kartesischen Koordinaten - :return: parametrische Koordinaten + :return: parametrische Breite, Länge """ x, y, z = point @@ -633,20 +648,20 @@ class EllipsoidTriaxial: def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]: """ - Umrechung von elliptischen in parametrische Koordinaten (über kartesische Koordinaten) - :param beta: elliptische Breite - :param lamb: elliptische Länge - :return: parametrische Koordinaten + Umrechung von ellipsoidischen in parametrische Koordinaten (über kartesische Koordinaten) + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge + :return: parametrische Breite, Länge """ cart = self.ell2cart(beta, lamb) return self.cart2para(cart) def para2ell(self, u: float, v: float) -> Tuple[float, float]: """ - Umrechung von parametrischen in elliptische Koordinaten (über kartesische Koordinaten) - :param u: u - :param v: v - :return: elliptische Koordinaten + Umrechung von parametrischen in ellipsoidische Koordinaten (über kartesische Koordinaten) + :param u: parametrische Breite + :param v: parametrische Länge + :return: ellipsoidische Breite, Länge """ cart = self.para2cart(u, v) return self.cart2ell(cart) @@ -654,12 +669,12 @@ class EllipsoidTriaxial: def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: """ Umrechung von parametrischen in geodätische Koordinaten (über kartesische Koordinaten) - :param u: u - :param v: v + :param u: parametrische Breite + :param v: parametrische Länge :param mode: ligas1, ligas2, oder ligas3 :param maxIter: maximale Anzahl Iterationen :param maxLoa: Level of Accuracy, das erreicht werden soll - :return: geodätische Koordinaten + :return: geodätische Breite, Länge, Höhe """ cart = self.para2cart(u, v) return self.cart2geod(cart, mode, maxIter, maxLoa) @@ -667,41 +682,41 @@ class EllipsoidTriaxial: def geod2para(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: """ Umrechung von geodätischen in parametrische Koordinaten (über kartesische Koordinaten) - :param phi: u - :param lamb: v + :param phi: geodätische Breite + :param lamb: geodätische Länge :param h: geodätische Höhe - :return: parametrische Koordinaten + :return: parametrische Breite, Länge """ cart = self.geod2cart(phi, lamb, h) return self.cart2para(cart) def ell2geod(self, beta: float, lamb: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]: """ - Umrechung von elliptischen in geodätische Koordinaten (über kartesische Koordinaten) - :param beta: elliptische Breite - :param lamb: eliptische Länge + Umrechung von ellipsoidischen in geodätische Koordinaten (über kartesische Koordinaten) + :param beta: ellipsoidische Breite + :param lamb: ellipsoidische Länge :param mode: ligas1, ligas2, oder ligas3 :param maxIter: maximale Anzahl Iterationen :param maxLoa: Level of Accuracy, das erreicht werden soll - :return: geodätische Koordinaten + :return: geodätische Breite, Länge, Höhe """ cart = self.ell2cart(beta, lamb) return self.cart2geod(cart, mode, maxIter, maxLoa) def geod2ell(self, phi: float, lamb: float, h: float) -> Tuple[float, float]: """ - Umrechung von geodätischen in elliptische Koordinaten (über kartesische Koordinaten) - :param phi: u - :param lamb: v + Umrechung von geodätischen in ellipsoidische Koordinaten (über kartesische Koordinaten) + :param phi: geodätische Breite + :param lamb: geodätische Länge :param h: geodätische Höhe - :return: elliptische Koordinaten + :return: ellipsoidische Breite, Länge """ cart = self.geod2cart(phi, lamb, h) return self.cart2ell(cart) def point_on(self, point: NDArray) -> bool: """ - Test, ob ein Punkt auf dem Ellipsoid liegt. + Test, ob ein Punkt auf dem Ellipsoid liegt :param point: kartesische 3D-Koordinaten :return: Punkt auf dem Ellispoid? """ diff --git a/jacobian_Ligas.py b/jacobian_Ligas.py index 7f49a5f..0e143b1 100644 --- a/jacobian_Ligas.py +++ b/jacobian_Ligas.py @@ -1,6 +1,17 @@ import numpy as np +from numpy.typing import NDArray +from typing import Tuple -def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): +def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]: + """ + Aufstellen des Gleichungssystem für den ersten Fall + :param E: Konstante E + :param F: Konstante F + :param G: Konstante G + :param pG: Punkt über dem Ellipsoid + :param pE: Punkt auf dem Ellipsoid + :return: inverse Jacobi-Matrix, Gleichungssystem + """ j11 = 2 * E * pE[0] j12 = 2 * F * pE[1] j13 = 2 * G * pE[2] @@ -23,7 +34,16 @@ def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): return invJ, fxE -def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): +def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, NDArray]: + """ + Aufstellen des Gleichungssystem für den zweiten Fall + :param E: Konstante E + :param F: Konstante F + :param G: Konstante G + :param pG: Punkt über dem Ellipsoid + :param pE: Punkt auf dem Ellipsoid + :return: inverse Jacobi-Matrix, Gleichungssystem + """ j11 = 2 * E * pE[0] j12 = 2 * F * pE[1] j13 = 2 * G * pE[2] @@ -48,7 +68,16 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): return invJ, fxE -def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray): +def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, NDArray]: + """ + Aufstellen des Gleichungssystem für den dritten Fall + :param E: Konstante E + :param F: Konstante F + :param G: Konstante G + :param pG: Punkt über dem Ellipsoid + :param pE: Punkt auf dem Ellipsoid + :return: inverse Jacobi-Matrix, Gleichungssystem + """ j11 = 2 * E * pE[0] j12 = 2 * F * pE[1] j13 = 2 * G * pE[2] diff --git a/kugel.py b/kugel.py index 70427bb..0abca94 100644 --- a/kugel.py +++ b/kugel.py @@ -6,6 +6,12 @@ import winkelumrechnungen as wu def cart2sph(point: NDArray) -> Tuple[float, float, float]: + """ + Umrechnung von kartesischen in sphärische Koordinaten + # TODO: Quelle + :param point: Punkt in kartesischen Koordinaten + :return: Radius, Breite, Länge + """ x, y, z = point r = sqrt(x**2 + y**2 + z**2) phi = arctan2(z, sqrt(x**2 + y**2)) @@ -15,6 +21,14 @@ def cart2sph(point: NDArray) -> Tuple[float, float, float]: def sph2cart(r: float, phi: float, lamb: float) -> NDArray: + """ + Umrechnung von sphärischen in kartesische Koordinaten + # TODO: Quelle + :param r: Radius + :param phi: Breite + :param lamb: Länge + :return: Punkt in kartesischen Koordinaten + """ x = r * cos(phi) * cos(lamb) y = r * cos(phi) * sin(lamb) z = r * sin(phi) @@ -22,7 +36,17 @@ def sph2cart(r: float, phi: float, lamb: float) -> NDArray: return np.array([x, y, z]) -def gha1(R, phi0, lamb0, s, alpha0): +def gha1(R: float, phi0: float, lamb0: float, s: float, alpha0: float) -> Tuple[float, float]: + """ + Berechnung der 1. GHA auf der Kugel + # TODO: Quelle + :param R: Radius + :param phi0: Breite des Startpunktes + :param lamb0: Länge des Startpunktes + :param s: Strecke + :param alpha0: Azimut + :return: Breite, Länge des Zielpunktes + """ s_ = s / R lamb1 = lamb0 + arctan2(sin(s_) * sin(alpha0), @@ -33,7 +57,17 @@ def gha1(R, phi0, lamb0, s, alpha0): return phi1, lamb1 -def gha2(R, phi0, lamb0, phi1, lamb1): +def gha2(R: float, phi0: float, lamb0: float, phi1: float, lamb1: float) -> Tuple[float, float, float]: + """ + Berechnung der 2. GHA auf der Kugel + # TODO: Quelle + :param R: Radius + :param phi0: Breite des Startpunktes + :param lamb0: Länge des Startpunktes + :param phi1: Breite des Zielpunktes + :param lamb1: Länge des Zielpunktes + :return: Azimut im Startpunkt, Azimut im Zielpunkt, Strecke + """ s_ = arccos(sin(phi0) * sin(phi1) + cos(phi0) * cos(phi1) * cos(lamb1 - lamb0)) s = R * s_ diff --git a/plot.ipynb b/plots.ipynb similarity index 99% rename from plot.ipynb rename to plots.ipynb index fec66c6..86ce718 100644 --- a/plot.ipynb +++ b/plots.ipynb @@ -110,9 +110,11 @@ }, { "metadata": { + "jupyter": { + "is_executing": true + }, "ExecuteTime": { - "end_time": "2026-02-04T08:51:31.229813Z", - "start_time": "2026-02-04T08:51:31.209338Z" + "start_time": "2026-02-04T10:52:53.553864Z" } }, "cell_type": "code", @@ -196,7 +198,7 @@ ], "id": "2d061b2f6ef6a78e", "outputs": [], - "execution_count": 20 + "execution_count": null }, { "metadata": { diff --git a/test.py b/test.py index 4d8bd3b..863672a 100644 --- a/test.py +++ b/test.py @@ -1,14 +1,7 @@ import numpy as np -from scipy.special import factorial as fact -from math import comb -import matplotlib.pyplot as plt import ellipsoide -from GHA_triaxial.panou import pq -ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978") +ell = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024") -x, y, z = ell.para2cart(0.5, 0.7) - -x1 = pq(ell, x, y, z) -x2 = ell.p_q(x, y, z) -pass \ No newline at end of file +cart = ell.para2cart(0, np.pi/2) +print(cart) \ No newline at end of file diff --git a/utils.py b/utils.py deleted file mode 100644 index 78403d3..0000000 --- a/utils.py +++ /dev/null @@ -1,21 +0,0 @@ -from numpy import arctan2 -from numpy.typing import NDArray - -from GHA_triaxial.panou import pq_ell -from ellipsoide import EllipsoidTriaxial - - -def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float: - """ - Berechnung des Richtungswinkels an einem Punkt anhand der Ableitung zu den kartesischen Koordinaten - :param ell: Ellipsoid - :param sigma: Ableitungsvektor ver kartesischen Koordinaten - :param point: Punkt - :return: Richtungswinkel - """ - p, q = pq_ell(ell, point) - P = float(p @ sigma) - Q = float(q @ sigma) - - alpha = arctan2(P, Q) - return alpha