Merge remote-tracking branch 'origin/main'

This commit is contained in:
Tammo.Weber
2026-02-05 10:59:45 +01:00
24 changed files with 3287 additions and 1683 deletions

View File

@@ -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)

View File

@@ -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))

136
GHA_triaxial/gha1_ana.py Normal file
View File

@@ -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

View File

@@ -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

129
GHA_triaxial/gha1_num.py Normal file
View File

@@ -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)

View File

@@ -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

View File

@@ -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]:

View File

@@ -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)

183
GHA_triaxial/utils.py Normal file
View File

@@ -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

2655
Tests/algorithms_test.ipynb Normal file

File diff suppressed because it is too large Load Diff

View File

@@ -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",

34
Tests/test_biaxial.py Normal file
View File

@@ -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

View File

@@ -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
}

View File

@@ -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)

View File

@@ -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)

View File

@@ -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?
"""

View File

@@ -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]

View File

@@ -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_

View File

@@ -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": {

13
test.py
View File

@@ -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
cart = ell.para2cart(0, np.pi/2)
print(cart)

View File

@@ -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