Umbenennung, Umstrukturierung, Doc-Strings
This commit is contained in:
@@ -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
136
GHA_triaxial/gha1_ana.py
Normal 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
|
||||
@@ -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
129
GHA_triaxial/gha1_num.py
Normal 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)
|
||||
@@ -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
|
||||
@@ -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]:
|
||||
@@ -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
183
GHA_triaxial/utils.py
Normal 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
|
||||
Reference in New Issue
Block a user