From 1fbfb555a4f3b0376a49992a4cae695a797a3d17 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Tue, 10 Feb 2026 21:10:11 +0100 Subject: [PATCH] wraps --- GHA_triaxial/gha1_ES.py | 18 +-- GHA_triaxial/gha1_ana.py | 7 +- GHA_triaxial/gha1_approx.py | 30 +++-- GHA_triaxial/gha1_num.py | 4 +- GHA_triaxial/gha2_num.py | 16 +-- GHA_triaxial/numeric_examples_karney.py | 4 +- GHA_triaxial/utils.py | 29 ++--- ellipsoide.py | 147 +++++++++++++----------- utils_angle.py | 17 ++- 9 files changed, 158 insertions(+), 114 deletions(-) diff --git a/GHA_triaxial/gha1_ES.py b/GHA_triaxial/gha1_ES.py index dba0b6f..8296b95 100644 --- a/GHA_triaxial/gha1_ES.py +++ b/GHA_triaxial/gha1_ES.py @@ -7,7 +7,7 @@ from ellipsoide import EllipsoidTriaxial from GHA_triaxial.gha1_ana import gha1_ana from GHA_triaxial.gha1_approx import gha1_approx from Hansen_ES_CMA import escma -from utils_angle import wrap_to_pi +from utils_angle import wrap_mpi_pi from numpy.typing import NDArray import winkelumrechnungen as wu @@ -40,7 +40,7 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \ R (XYZ) = Punkt in XYZ """ # Berechnungshilfen - omega = wrap_to_pi(omega) + omega = wrap_mpi_pi(omega) cb = np.cos(beta) sb = np.sin(beta) co = np.cos(omega) @@ -121,7 +121,7 @@ def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_ sE = float(np.dot(vT_hat, E_hat_curr)) sN = float(np.dot(vT_hat, N_hat_curr)) - return wrap_to_pi(float(np.arctan2(sE, sN))) + return wrap_mpi_pi(float(np.arctan2(sE, sN))) def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float, gamma0: float, @@ -158,7 +158,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float #d_beta = ds * float(np.cos(alpha_i)) / Nn_i #d_omega = ds * float(np.sin(alpha_i)) / En_i beta_pred = beta_i + d_beta - omega_pred = wrap_to_pi(omega_i + d_omega) + omega_pred = wrap_mpi_pi(omega_i + d_omega) xmean = np.array([beta_pred, omega_pred], dtype=float) @@ -175,7 +175,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float :return: Fitnesswert (f) """ beta = x[0] - omega = wrap_to_pi(x[1]) + omega = wrap_mpi_pi(x[1]) P = ell.ell2cart_karney(beta, omega) # in kartesischer Koordinaten d = float(np.linalg.norm(P - P_i)) # Distanz zwischen @@ -201,7 +201,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES beta_best = xb[0] - omega_best = wrap_to_pi(xb[1]) + omega_best = wrap_mpi_pi(xb[1]) P_best = ell.ell2cart_karney(beta_best, omega_best) E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta_best, omega_best, ell) alpha_end = azimuth_at_ESpoint(P_i, P_best, E_j, N_j, U_j) @@ -223,8 +223,8 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, :return: Zielpunkt Pk, Azimut am Zielpunkt und Punktliste """ beta = float(beta0) - omega = wrap_to_pi(float(omega0)) - alpha = wrap_to_pi(float(alpha0)) + omega = wrap_mpi_pi(float(omega0)) + alpha = wrap_mpi_pi(float(alpha0)) gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0 @@ -243,7 +243,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, ell=ell, maxSegLen=maxSegLen) s_acc += ds P_all.append(P) - alpha_end.append(alpha) + alpha_end.append(wrap_mpi_pi(alpha)) if step > nsteps_est + 50: raise RuntimeError("GHA1_ES: Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.") Pk = P_all[-1] diff --git a/GHA_triaxial/gha1_ana.py b/GHA_triaxial/gha1_ana.py index 1ac772e..69d0ca0 100644 --- a/GHA_triaxial/gha1_ana.py +++ b/GHA_triaxial/gha1_ana.py @@ -4,8 +4,9 @@ from typing import Tuple import numpy as np from numpy import sin, cos, arctan2 -from numpy._typing import NDArray +from numpy.typing import NDArray import winkelumrechnungen as wu +from utils_angle import wrap_0_2pi from ellipsoide import EllipsoidTriaxial from GHA_triaxial.utils import pq_para @@ -110,7 +111,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa if alpha1 < 0: alpha1 += 2 * np.pi - return p1, alpha1 + return p1, wrap_0_2pi(alpha1) def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]: @@ -134,7 +135,7 @@ def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, ma if h > 1e-5: raise Exception("GHA1_ana: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid") - return point_end, alpha_end + return point_end, wrap_0_2pi(alpha_end) if __name__ == "__main__": diff --git a/GHA_triaxial/gha1_approx.py b/GHA_triaxial/gha1_approx.py index a84edfe..e81c043 100644 --- a/GHA_triaxial/gha1_approx.py +++ b/GHA_triaxial/gha1_approx.py @@ -6,6 +6,7 @@ from GHA_triaxial.gha1_ana import gha1_ana from GHA_triaxial.utils import func_sigma_ell, louville_constant, pq_ell import plotly.graph_objects as go import winkelumrechnungen as wu +from utils_angle import wrap_0_2pi, wrap_mhalfpi_halfpi, wrap_mpi_pi def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]: """ @@ -37,19 +38,24 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, if last_p is not None and np.dot(p, last_p) < 0: p = -p q = -q + last_p = p sigma = p * sin(alpha1) + q * cos(alpha1) if last_sigma is not None and np.dot(sigma, last_sigma) < 0: sigma = -sigma + alpha1 += np.pi + alpha1 = wrap_0_2pi(alpha1) p2 = p1 + ds_step * sigma p2 = ell.point_onto_ellipsoid(p2) - dalpha = 1e-6 + dalpha = 1e-9 l2 = louville_constant(ell, p2, alpha1) dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha - alpha2 = alpha1 + (l0 - l2) / dl_dalpha - + if abs(dl_dalpha) < 1e-20: + alpha2 = alpha1 + 0 + else: + alpha2 = alpha1 + (l0 - l2) / dl_dalpha points.append(p2) - alphas.append(alpha2) + alphas.append(wrap_0_2pi(alpha2)) ds_step = np.linalg.norm(p2 - p1) s_curr += ds_step @@ -88,11 +94,11 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray): if __name__ == '__main__': - ell = EllipsoidTriaxial.init_name("BursaSima1980round") - P0 = ell.ell2cart(wu.deg2rad(89), wu.deg2rad(1)) - alpha0 = wu.deg2rad(2) - s = 200000 - P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=100, all_points=True) - P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=20, maxPartCircum=2) - print(np.linalg.norm(P1_app - P1_ana)) - show_points(points, P0, P1_ana) + ell = EllipsoidTriaxial.init_name("KarneyTest2024") + P0 = ell.ell2cart(wu.deg2rad(15), wu.deg2rad(15)) + alpha0 = wu.deg2rad(270) + s = 1 + P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=0.1, all_points=True) + # P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=40, maxPartCircum=32) + # print(np.linalg.norm(P1_app - P1_ana)) + # show_points(points, P0, P0) diff --git a/GHA_triaxial/gha1_num.py b/GHA_triaxial/gha1_num.py index 1b7e101..25dcca3 100644 --- a/GHA_triaxial/gha1_num.py +++ b/GHA_triaxial/gha1_num.py @@ -10,6 +10,7 @@ from typing import Callable, Tuple, List from numpy.typing import NDArray from GHA_triaxial.utils import alpha_ell2para, pq_ell +from utils_angle import wrap_0_2pi def buildODE(ell: EllipsoidTriaxial) -> Callable: @@ -75,8 +76,7 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu alpha1 = arctan2(P, Q) - if alpha1 < 0: - alpha1 += 2 * np.pi + alpha1 = wrap_0_2pi(alpha1) _, _, h = ell.cart2geod(point1, "ligas3") if h > 1e-5: diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index eb54c5d..6a920e8 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -7,7 +7,7 @@ import winkelumrechnungen as wu from typing import Tuple from numpy.typing import NDArray import ausgaben as aus -from utils_angle import cot, arccot, wrap_to_pi +from utils_angle import cot, arccot, wrap_mpi_pi, wrap_0_2pi def norm_a(a: float) -> float: @@ -22,7 +22,7 @@ def azimut(E: float, G: float, dbeta_du: float, dlamb_du: float) -> float: def sph_azimuth(beta1, lam1, beta2, lam2): - dlam = wrap_to_pi(lam2 - lam1) + dlam = wrap_mpi_pi(lam2 - lam1) y = np.sin(dlam) * np.cos(beta2) x = np.cos(beta1) * np.sin(beta2) - np.sin(beta1) * np.cos(beta2) * np.cos(dlam) a = np.arctan2(y, x) @@ -347,8 +347,8 @@ def gha2_num( return best[0], best[1], sgn, dbeta, ode_beta - lamb0 = float(wrap_to_pi(lamb_0)) - lamb1 = float(wrap_to_pi(lamb_1)) + lamb0 = float(wrap_mpi_pi(lamb_0)) + lamb1 = float(wrap_mpi_pi(lamb_1)) beta0 = float(beta_0) beta1 = float(beta_1) @@ -491,7 +491,7 @@ def gha2_num( else: s = np.trapz(integrand, dx=h) - return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr + return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N_full, integrand_lambda) beta_end, beta_p_end, _, _ = y_end @@ -502,7 +502,7 @@ def gha2_num( (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), float(lamb0 + dlamb)) alpha_1 = azimut(E_end, G_end, dbeta_du=float(beta_p_end) * sgn, dlamb_du=1.0 * sgn) - return float(alpha_0), float(alpha_1), float(s) + return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s) # Fall 2 (lambda_0 == lambda_1) N = int(n) @@ -574,7 +574,7 @@ def gha2_num( else: s = np.trapz(integrand, dx=h) - return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr + return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr _, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta) lamb_end, lamb_p_end, _, _ = y_end @@ -585,7 +585,7 @@ def gha2_num( (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, float(lamb_end)) alpha_1 = azimut(E_end, G_end, dbeta_du=1.0 * sgn, dlamb_du=float(lamb_p_end) * sgn) - return float(alpha_0), float(alpha_1), float(s) + return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s) if __name__ == "__main__": diff --git a/GHA_triaxial/numeric_examples_karney.py b/GHA_triaxial/numeric_examples_karney.py index 32a59fe..207fb31 100644 --- a/GHA_triaxial/numeric_examples_karney.py +++ b/GHA_triaxial/numeric_examples_karney.py @@ -113,7 +113,7 @@ def get_random_examples_gamma(group: str, num: int, seed: int = None, length: st beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example gamma = jacobi_konstante(beta0, lamb0, alpha0_ell, ell) - if group not in ["a", "b", "c", "d", "e"]: + if group not in ["a", "b", "c", "d", "e", "de"]: break elif group == "a" and not 1 >= gamma >= 0.01: continue @@ -125,6 +125,8 @@ def get_random_examples_gamma(group: str, num: int, seed: int = None, length: st continue elif group == "e" and not -1e-17 >= gamma >= -1: continue + elif group == "de" and not -eps > gamma > -1: + continue if length == "short": if example[6] < long_short: diff --git a/GHA_triaxial/utils.py b/GHA_triaxial/utils.py index 6bfc5c0..8eed03d 100644 --- a/GHA_triaxial/utils.py +++ b/GHA_triaxial/utils.py @@ -2,8 +2,8 @@ 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 utils_angle import wrap_mpi_pi, wrap_0_2pi, wrap_mhalfpi_halfpi from ellipsoide import EllipsoidTriaxial @@ -21,7 +21,7 @@ def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float Q = float(q @ sigma) alpha = arctan2(P, Q) - return alpha + return wrap_0_2pi(alpha) def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]: @@ -43,10 +43,10 @@ def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float 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-9: + if np.linalg.norm(sigma_para - sigma_ell) > 1e-7: raise Exception("alpha_para2ell: Differenz in den Richtungsableitungen") - return beta, lamb, alpha_ell + return beta, lamb, wrap_0_2pi(alpha_ell) def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]: @@ -68,10 +68,10 @@ def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: 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: + if np.linalg.norm(sigma_para - sigma_ell) > 1e-7: raise Exception("alpha_ell2para: Differenz in den Richtungsableitungen") - return u, v, alpha_para + return u, v, wrap_0_2pi(alpha_para) def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha_ell: float) -> NDArray: @@ -124,11 +124,10 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: :param point: Punkt :return: p und q """ - x, y, z = point n = ell.func_n(point) beta, lamb = ell.cart2ell(point) - if abs(cos(beta)) < 1e-12 and abs(np.sin(lamb)) < 1e-12: + if abs(cos(beta)) < 1e-15 and abs(np.sin(lamb)) < 1e-15: if beta > 0: p = np.array([0, -1, 0]) else: @@ -137,11 +136,7 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: 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 + _, t2 = ell.func_t12(point) 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) @@ -181,3 +176,11 @@ def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]: q = q / np.linalg.norm(q) return p, q + + +if __name__ == "__main__": + ell = EllipsoidTriaxial.init_name("KarneyTest2024") + alpha_para = 0 + u, v = ell.ell2para(np.pi/2, 0) + alpha_ell = alpha_para2ell(ell, u, v, alpha_para) + pass \ No newline at end of file diff --git a/ellipsoide.py b/ellipsoide.py index 2208699..f5d9f8e 100644 --- a/ellipsoide.py +++ b/ellipsoide.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt from typing import Tuple from numpy.typing import NDArray import math +from utils_angle import wrap_mpi_pi, wrap_0_2pi, wrap_mhalfpi_halfpi class EllipsoidBiaxial: @@ -218,8 +219,11 @@ class EllipsoidTriaxial: c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 - (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - ( self.ax ** 2 + self.ay ** 2) * z ** 2) - if c1 ** 2 - 4 * c0 < 0: + if c1 ** 2 - 4 * c0 < -1e-9: t2 = np.nan + raise Exception("t1, t2: Negativer Wurzelterm") + elif c1 ** 2 - 4 * c0 < 0: + t2 = 0 else: t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2 if t2 == 0: @@ -284,6 +288,11 @@ class EllipsoidTriaxial: beta, lamb = np.broadcast_arrays(beta, lamb) + beta = np.where( + np.isclose(np.abs(beta), np.pi / 2, atol=1e-15), + beta * 8999999999999999 / 9000000000000000, + beta + ) B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2 L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2 @@ -419,7 +428,7 @@ class EllipsoidTriaxial: if delta_r > 1e-6: raise Exception("Umrechnung cart2ell: Punktdifferenz") - return beta, lamb + return wrap_mhalfpi_halfpi(beta), wrap_mpi_pi(lamb) except Exception as e: # Wenn die Berechnung fehlschlägt auf Grund von sehr kleinem y, solange anpassen, bis Umrechnung ohne Fehler @@ -605,8 +614,8 @@ class EllipsoidTriaxial: if abs(zG) < eps: phi = 0 - - return phi, lamb, h + wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb) + return wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb), h def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray: """ @@ -643,8 +652,8 @@ class EllipsoidTriaxial: v = 2 * arctan2(v_check1, v_check2 + v_factor) else: v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor) - - return u, v + wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v) + return wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v) def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]: """ @@ -749,63 +758,71 @@ class EllipsoidTriaxial: if __name__ == "__main__": - ell = EllipsoidTriaxial.init_name("BursaSima1980") - diff_list = [] - diffs_para = [] - diffs_ell = [] - diffs_geod = [] - points = [] - for v_deg in range(-180, 181, 5): - for u_deg in range(-90, 91, 5): - v = wu.deg2rad(v_deg) - u = wu.deg2rad(u_deg) - point = ell.para2cart(u, v) - points.append(point) + ell = EllipsoidTriaxial.init_name("KarneyTest2024") + # cart = ell.ell2cart(np.pi/2, 0) + # print(cart) + # cart = ell.ell2cart(np.pi/2*8999999999999999/9000000000000000, 0) + # print(cart) + elli = ell.cart2ell([0, 0.0, 1/np.sqrt(2)]) + print(elli) - elli = ell.cart2ell(point) - cart_elli = ell.ell2cart(elli[0], elli[1]) - diff_ell = np.linalg.norm(point - cart_elli, axis=-1) - - para = ell.cart2para(point) - cart_para = ell.para2cart(para[0], para[1]) - diff_para = np.linalg.norm(point - cart_para, axis=-1) - - geod = ell.cart2geod(point, "ligas3") - cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) - diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1) - - diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3]) - diffs_ell.append([diff_ell]) - diffs_para.append([diff_para]) - diffs_geod.append([diff_geod3]) - - diff_list = np.array(diff_list) - diffs_ell = np.array(diffs_ell) - diffs_para = np.array(diffs_para) - diffs_geod = np.array(diffs_geod) - - pass - - points = np.array(points) - fig = plt.figure() - ax = fig.add_subplot(projection='3d') - - sc = ax.scatter( - points[:, 0], - points[:, 1], - points[:, 2], - c=diffs_ell, # Farbcode = diff - cmap='viridis', # Colormap - s=10 + 20 * diffs_ell, # optional: Größe abhängig vom diff - alpha=0.8 - ) - - # Farbskala - cbar = plt.colorbar(sc) - cbar.set_label("diff") - - ax.set_xlabel("X") - ax.set_ylabel("Y") - ax.set_zlabel("Z") - - plt.show() \ No newline at end of file + # ell = EllipsoidTriaxial.init_name("BursaSima1980") + # diff_list = [] + # diffs_para = [] + # diffs_ell = [] + # diffs_geod = [] + # points = [] + # for v_deg in range(-180, 181, 5): + # for u_deg in range(-90, 91, 5): + # v = wu.deg2rad(v_deg) + # u = wu.deg2rad(u_deg) + # point = ell.para2cart(u, v) + # points.append(point) + # + # elli = ell.cart2ell(point) + # cart_elli = ell.ell2cart(elli[0], elli[1]) + # diff_ell = np.linalg.norm(point - cart_elli, axis=-1) + # + # para = ell.cart2para(point) + # cart_para = ell.para2cart(para[0], para[1]) + # diff_para = np.linalg.norm(point - cart_para, axis=-1) + # + # geod = ell.cart2geod(point, "ligas3") + # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2]) + # diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1) + # + # diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3]) + # diffs_ell.append([diff_ell]) + # diffs_para.append([diff_para]) + # diffs_geod.append([diff_geod3]) + # + # diff_list = np.array(diff_list) + # diffs_ell = np.array(diffs_ell) + # diffs_para = np.array(diffs_para) + # diffs_geod = np.array(diffs_geod) + # + # pass + # + # points = np.array(points) + # fig = plt.figure() + # ax = fig.add_subplot(projection='3d') + # + # sc = ax.scatter( + # points[:, 0], + # points[:, 1], + # points[:, 2], + # c=diffs_ell, # Farbcode = diff + # cmap='viridis', # Colormap + # s=10 + 20 * diffs_ell, # optional: Größe abhängig vom diff + # alpha=0.8 + # ) + # + # # Farbskala + # cbar = plt.colorbar(sc) + # cbar.set_label("diff") + # + # ax.set_xlabel("X") + # ax.set_ylabel("Y") + # ax.set_zlabel("Z") + # + # plt.show() \ No newline at end of file diff --git a/utils_angle.py b/utils_angle.py index e560739..abbd95c 100644 --- a/utils_angle.py +++ b/utils_angle.py @@ -1,4 +1,5 @@ import numpy as np +import winkelumrechnungen as wu def arccot(x): @@ -9,5 +10,19 @@ def cot(a): return np.cos(a) / np.sin(a) -def wrap_to_pi(x): +def wrap_mpi_pi(x): return (x + np.pi) % (2 * np.pi) - np.pi + + +def wrap_mhalfpi_halfpi(x): + return (x + np.pi / 2) % np.pi - np.pi / 2 + + +def wrap_0_2pi(x): + return x % (2 * np.pi) + + +if __name__ == "__main__": + print(wu.rad2deg(wrap_mhalfpi_halfpi(wu.deg2rad(181)))) + print(wu.rad2deg(wrap_0_2pi(wu.deg2rad(181)))) + print(wu.rad2deg(wrap_mpi_pi(wu.deg2rad(181))))