105 lines
3.7 KiB
Python
105 lines
3.7 KiB
Python
import numpy as np
|
|
from ellipsoide import EllipsoidTriaxial
|
|
from GHA_triaxial.panou_2013_2GHA_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
|
|
|
|
|
|
def gha2_approx(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]:
|
|
"""
|
|
Numerische Approximation für die zweite Hauptaufgabe
|
|
:param ell: Ellipsoid
|
|
:param p0: Startpunkt
|
|
:param p1: Endpunkt
|
|
:param ds: maximales Streckenelement
|
|
:param all_points: Alle Punkte ausgeben?
|
|
:return:
|
|
"""
|
|
points = np.array([p0, p1])
|
|
while True:
|
|
new_points = []
|
|
|
|
for i in range(len(points)-1):
|
|
new_points.append(points[i])
|
|
|
|
pi = points[i] + 1/2 * (points[i+1] - points[i])
|
|
pi = ell.point_onto_ellipsoid(pi)
|
|
|
|
new_points.append(pi)
|
|
|
|
new_points.append(points[-1])
|
|
points = np.array(new_points)
|
|
|
|
elements = np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)])
|
|
|
|
if np.average(elements) < ds:
|
|
break
|
|
|
|
p0i = ell.point_onto_ellipsoid(p0 + ds / 100 * (points[1] - p0) / np.linalg.norm(points[1] - p0))
|
|
sigma0 = (p0i - p0) / np.linalg.norm(p0i - p0)
|
|
alpha0 = sigma2alpha(ell, sigma0, p0)
|
|
|
|
p1i = ell.point_onto_ellipsoid(p1 - ds / 100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2]))
|
|
sigma1 = (p1 - p1i) / np.linalg.norm(p1 - p1i)
|
|
alpha1 = sigma2alpha(ell, sigma1, p1)
|
|
|
|
s = np.sum(np.array([np.linalg.norm(points[i] - points[i+1]) for i in range(len(points)-1)]))
|
|
|
|
if all_points:
|
|
return alpha0, alpha1, s, np.array(points)
|
|
else:
|
|
return alpha0, alpha1, s
|
|
|
|
def show_points(points: NDArray, points_app: NDArray, p0: NDArray, p1: NDArray):
|
|
"""
|
|
Anzeigen der Punkte
|
|
:param points: wahre Punkte der Linie
|
|
:param points_app: approximierte Punkte der Linie
|
|
:param p0: wahrer 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="green", width=3), name="Analytisch")
|
|
fig.add_scatter3d(x=points_app[:, 0], y=points_app[:, 1], z=points_app[:, 2],
|
|
mode='lines', line=dict(color="red", width=3), name="Approximiert")
|
|
fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
|
|
mode='markers', marker=dict(color="black"), name="P0")
|
|
fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
|
|
mode='markers', marker=dict(color="black"), name="P1")
|
|
|
|
fig.update_layout(
|
|
scene=dict(xaxis_title='X [km]',
|
|
yaxis_title='Y [km]',
|
|
zaxis_title='Z [km]',
|
|
aspectmode='data'))
|
|
|
|
fig.show()
|
|
|
|
|
|
if __name__ == '__main__':
|
|
ell = EllipsoidTriaxial.init_name("Bursa1970")
|
|
|
|
beta0, lamb0 = (0.1, 0.1)
|
|
P0 = ell.ell2cart(beta0, lamb0)
|
|
beta1, lamb1 = (0.3, np.pi)
|
|
P1 = ell.ell2cart(beta1, lamb1)
|
|
|
|
alpha0_app, alpha1_app, s_app, points = gha2_approx(ell, P0, P1, ds=100, all_points=True)
|
|
print("done")
|
|
alpha0, alpha1, s, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=10000, all_points=True)
|
|
points_ana = []
|
|
for beta, lamb in zip(betas, lambs):
|
|
points_ana.append(ell.ell2cart(beta, lamb))
|
|
points_ana = np.array(points_ana)
|
|
|
|
show_points(points_ana, points, P0, P1)
|
|
print(f"Differenz s: {s_app - s} m")
|
|
print(f"Differenz alpha0: {wu.rad2deg(alpha0_app - alpha0)}°")
|
|
print(f"Differenz alpha1: {wu.rad2deg(alpha1_app - alpha1)}°")
|