Näherungslösung GHA 2
This commit is contained in:
@@ -1,12 +1,10 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arcsin, arccos, arctan2
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import matplotlib.pyplot as plt
|
||||
from panou import louville_constant, func_sigma_ell, gha1_ana
|
||||
import plotly.graph_objects as go
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: int):
|
||||
def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]:
|
||||
l0 = louville_constant(ell, p0, alpha0)
|
||||
points = [p0]
|
||||
alphas = [alpha0]
|
||||
@@ -17,10 +15,9 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: in
|
||||
break
|
||||
p1 = points[-1]
|
||||
alpha1 = alphas[-1]
|
||||
x1, y1, z1 = p1
|
||||
sigma = func_sigma_ell(ell, x1, y1, z1, alpha1)
|
||||
sigma = func_sigma_ell(ell, p1, alpha1)
|
||||
p2 = p1 + ds_step * sigma
|
||||
p2, _, _, _ = ell.cartonell(p2)
|
||||
p2 = ell.cartonell(p2)
|
||||
ds_step = np.linalg.norm(p2 - p1)
|
||||
|
||||
points.append(p2)
|
||||
@@ -30,7 +27,11 @@ def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: in
|
||||
alpha2 = alpha1 + (l0 - l2) / dl_dalpha
|
||||
alphas.append(alpha2)
|
||||
s_curr += ds_step
|
||||
return points[-1], alphas[-1], np.array(points)
|
||||
|
||||
if all_points:
|
||||
return points[-1], alphas[-1], np.array(points)
|
||||
else:
|
||||
return points[-1], alphas[-1]
|
||||
|
||||
def show_points(points, p0, p1):
|
||||
fig = go.Figure()
|
||||
@@ -57,7 +58,7 @@ if __name__ == '__main__':
|
||||
P0 = ell.para2cart(0, 0)
|
||||
alpha0 = wu.deg2rad(90)
|
||||
s = 1000000
|
||||
P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 60, maxPartCircum=32)
|
||||
P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, 5000)
|
||||
P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=32)
|
||||
P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, ds=5000, all_points=True)
|
||||
show_points(points, P0, P1_ana)
|
||||
print(np.linalg.norm(P1_app - P1_ana))
|
||||
|
||||
111
GHA_triaxial/approx_gha2.py
Normal file
111
GHA_triaxial/approx_gha2.py
Normal file
@@ -0,0 +1,111 @@
|
||||
import numpy as np
|
||||
from numpy import arctan2
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from panou import pq_ell
|
||||
from 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
|
||||
|
||||
def sigma2alpha(sigma: NDArray, point: NDArray) -> float:
|
||||
"""
|
||||
Berechnung des Richtungswinkels an einem Punkt anhand der Ableitung zu den kartesischen Koordinaten
|
||||
: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
|
||||
|
||||
def gha2(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.cartonell(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.cartonell(p0 + ds/100 * (points[1] - p0) / np.linalg.norm(points[1] - p0))
|
||||
sigma0 = (p0i - p0) / np.linalg.norm(p0i - p0)
|
||||
alpha0 = sigma2alpha(sigma0, p0)
|
||||
|
||||
p1i = ell.cartonell(p1 - ds/100 * (p1 - points[-2]) / np.linalg.norm(p1 - points[-2]))
|
||||
sigma1 = (p1 - p1i) / np.linalg.norm(p1 - p1i)
|
||||
alpha1 = sigma2alpha(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):
|
||||
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("KarneyTest2024")
|
||||
|
||||
beta0, lamb0 = (0.2, 0.1)
|
||||
P0 = ell.ell2cart(beta0, lamb0)
|
||||
beta1, lamb1 = (0.7, 0.3)
|
||||
P1 = ell.ell2cart(beta1, lamb1)
|
||||
|
||||
alpha0_app, alpha1_app, s_app, points = gha2(ell, P0, P1, ds=1e-4, all_points=True)
|
||||
|
||||
alpha0, alpha1, s, betas, lambs = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=5000, 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)}°")
|
||||
@@ -1,108 +0,0 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arcsin, arccos, arctan2
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
dbeta_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ey**2 * sin(alpha)**2 * cos(beta) * sin(beta)
|
||||
dlamb_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ee**2 * cos(alpha)**2 * sin(lamb) * cos(lamb)
|
||||
dalpha_dc = lambda ell, beta, lamb, alpha: (2 * sin(alpha) * cos(alpha) *
|
||||
(ell.Ey**2 * cos(beta)**2 + ell.Ee**2 * sin(lamb)**2))
|
||||
|
||||
lamb2_sphere = lambda r, phi1, lamb1, a12, s: lamb1 + arctan2(sin(s/r) * sin(a12),
|
||||
cos(phi1) * cos(s/r) - sin(s/r) * cos(a12))
|
||||
phi2_sphere = lambda r, phi1, lamb1, a12, s: arcsin(sin(phi1) * cos(s/r) + cos(phi1) * sin(s/r) * cos(a12))
|
||||
|
||||
a12_sphere = lambda phi1, lamb1, phi2, lamb2: arctan2(cos(phi2) * sin(lamb2 - lamb1),
|
||||
cos(phi1) * cos(phi2) -
|
||||
sin(phi1) * cos(phi2) * cos(lamb2 - lamb1))
|
||||
s_sphere = lambda r, phi1, lamb1, phi2, lamb2: r * arccos(sin(phi1) * sin(phi2) +
|
||||
cos(phi1) * cos(phi2) * cos(lamb2 - lamb1))
|
||||
|
||||
louville = lambda beta, lamb, alpha: (ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 -
|
||||
ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2)
|
||||
|
||||
def points_approx_gha2(r: float, phi1: np.ndarray, lamb1: np.ndarray, phi2: np.ndarray, lamb2: np.ndarray, num: int = None, step_size: float = 10000):
|
||||
s_approx = s_sphere(r, phi1, lamb1, phi2, lamb2)
|
||||
if num is not None:
|
||||
step_size = s_approx / (num+1)
|
||||
a_approx = a12_sphere(phi1, lamb1, phi2, lamb2)
|
||||
|
||||
points = [np.array([phi1, lamb1, a_approx])]
|
||||
current_s = step_size
|
||||
while current_s < s_approx:
|
||||
phi_n = phi2_sphere(r, phi1, lamb1, a_approx, current_s)
|
||||
lamb_n = lamb2_sphere(r, phi1, lamb1, a_approx, current_s)
|
||||
points.append(np.array([phi_n, lamb_n, a_approx]))
|
||||
current_s += step_size
|
||||
points.append(np.array([phi2, lamb2, a_approx]))
|
||||
return points
|
||||
|
||||
|
||||
def num_update(ell: EllipsoidTriaxial, points, diffs):
|
||||
for i, (beta, lamb, alpha) in enumerate(points):
|
||||
dalpha = dalpha_dc(ell, beta, lamb, alpha)
|
||||
if i == 0 or i == len(points) - 1:
|
||||
grad = np.array([0, 0, dalpha])
|
||||
else:
|
||||
dbeta = dbeta_dc(ell, beta, lamb, alpha)
|
||||
dlamb = dlamb_dc(ell, beta, lamb, alpha)
|
||||
grad = np.array([dbeta, dlamb, dalpha])
|
||||
|
||||
delta = -diffs[i] * grad / np.dot(grad, grad)
|
||||
points[i] += delta
|
||||
return points
|
||||
|
||||
|
||||
def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int):
|
||||
beta1, lamb1 = ell.cart2ell(p1)
|
||||
beta2, lamb2 = ell.cart2ell(p2)
|
||||
points = points_approx_gha2(ell.ax, beta1, lamb1, beta2, lamb2, 5)
|
||||
|
||||
for j in range(maxI):
|
||||
constants = [louville(point[0], point[1], point[2]) for point in points]
|
||||
mean_constant = np.mean(constants)
|
||||
diffs = constants - mean_constant
|
||||
if np.mean(np.abs(diffs)) > 10:
|
||||
points = num_update(ell, points, diffs)
|
||||
else:
|
||||
break
|
||||
for k in range(maxI):
|
||||
|
||||
|
||||
last_diff_alpha = points[-2][-1] - points[-3][-1]
|
||||
alpha_extrap = points[-2][-1] + last_diff_alpha
|
||||
if abs(alpha_extrap - points[-1][-1]) > 0.0005:
|
||||
pass
|
||||
else:
|
||||
break
|
||||
pass
|
||||
pass
|
||||
return points
|
||||
|
||||
def show_points(ell: EllipsoidTriaxial, points):
|
||||
points_cart = []
|
||||
for point in points:
|
||||
points_cart.append(ell.ell2cart(point[0], point[1]))
|
||||
points_cart = np.array(points_cart)
|
||||
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111, projection='3d')
|
||||
|
||||
ax.plot(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2])
|
||||
ax.scatter(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2])
|
||||
|
||||
ax.set_xlabel('X')
|
||||
ax.set_ylabel('Y')
|
||||
ax.set_zlabel('Z')
|
||||
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("Eitschberger1978")
|
||||
p1 = np.array([4189000, 812000, 4735000])
|
||||
p2 = np.array([4090000, 868000, 4808000])
|
||||
p1, phi1, lamb1, h1 = ell.cartonell(p1)
|
||||
p2, phi2, lamb2, h2 = ell.cartonell(p2)
|
||||
points = gha2(ell, p1, p2, 10)
|
||||
show_points(ell, points)
|
||||
@@ -8,21 +8,22 @@ 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
|
||||
from typing import Callable, Tuple, List
|
||||
from numpy.typing import NDArray
|
||||
|
||||
|
||||
def pq_ell(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]:
|
||||
def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
"""
|
||||
Berechnung von p und q in elliptischen Koordinaten
|
||||
Panou, Korakitits 2019
|
||||
:param x: x
|
||||
:param y: y
|
||||
:param z: z
|
||||
:param ell: Ellipsoid
|
||||
:param point: Punkt
|
||||
:return: p und q
|
||||
"""
|
||||
n = ell.func_n(x, y, z)
|
||||
x, y, z = point
|
||||
n = ell.func_n(point)
|
||||
|
||||
beta, lamb = ell.cart2ell(np.array([x, y, z]))
|
||||
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
|
||||
|
||||
@@ -49,7 +50,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
:param ell: Ellipsoid
|
||||
:return: DGL-System
|
||||
"""
|
||||
def ODE(s: float, v: np.ndarray) -> np.ndarray:
|
||||
def ODE(s: float, v: NDArray) -> NDArray:
|
||||
"""
|
||||
DGL-System
|
||||
:param s: unabhängige Variable
|
||||
@@ -58,7 +59,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
"""
|
||||
x, dxds, y, dyds, z, dzds = v
|
||||
|
||||
H = ell.func_H(x, y, z)
|
||||
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
|
||||
@@ -68,7 +69,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
|
||||
return ODE
|
||||
|
||||
def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float, list]:
|
||||
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:
|
||||
@@ -76,12 +77,14 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,
|
||||
:param alpha0:
|
||||
:param s:
|
||||
:param num:
|
||||
:param all_points:
|
||||
:return:
|
||||
"""
|
||||
phi, lam, _ = ell.cart2geod(point, "ligas3")
|
||||
x0, y0, z0 = ell.geod2cart(phi, lam, 0)
|
||||
p0 = ell.geod2cart(phi, lam, 0)
|
||||
x0, y0, z0 = p0
|
||||
|
||||
p, q = pq_ell(ell, x0, y0, z0)
|
||||
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)
|
||||
@@ -93,7 +96,9 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,
|
||||
_, werte = rk.rk4(ode, 0, v_init, s, num)
|
||||
x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1]
|
||||
|
||||
p1, q1 = pq_ell(ell, x1, y1, z1)
|
||||
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)
|
||||
@@ -103,21 +108,23 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
|
||||
return np.array([x1, y1, z1]), alpha1, werte
|
||||
if all_points:
|
||||
return point1, alpha1, werte
|
||||
else:
|
||||
return point1, alpha1
|
||||
|
||||
# ---------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
def pq_para(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.ndarray, np.ndarray]:
|
||||
def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
"""
|
||||
Berechnung von p und q in parametrischen Koordinaten
|
||||
Panou, Korakitits 2020
|
||||
:param x: x
|
||||
:param y: y
|
||||
:param z: z
|
||||
:param ell: Ellipsoid
|
||||
:param point: Punkt
|
||||
:return: p und q
|
||||
"""
|
||||
n = ell.func_n(x, y, z)
|
||||
u, v = ell.cart2para(np.array([x, y, z]))
|
||||
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)
|
||||
@@ -131,12 +138,12 @@ def pq_para(ell: EllipsoidTriaxial, x: float, y: float, z: float) -> tuple[np.nd
|
||||
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):
|
||||
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")
|
||||
|
||||
return p, q
|
||||
|
||||
def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]:
|
||||
"""
|
||||
Panou, Korakitits 2020, 5ff.
|
||||
:param ell:
|
||||
@@ -153,7 +160,7 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
y_m = [y]
|
||||
z_m = [z]
|
||||
|
||||
p, q = pq_para(ell, x, y, z)
|
||||
p, q = pq_para(ell, point)
|
||||
|
||||
# 48-50
|
||||
x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0))
|
||||
@@ -208,7 +215,8 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
for c in reversed(c_m):
|
||||
z_s = z_s * s + c
|
||||
|
||||
p_s, q_s = pq_para(ell, x_s, y_s, z_s)
|
||||
p1 = np.array([x_s, y_s, z_s])
|
||||
p_s, q_s = pq_para(ell, p1)
|
||||
|
||||
# 57-59
|
||||
dx_s = 0
|
||||
@@ -234,10 +242,10 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
|
||||
return np.array([x_s, y_s, z_s]), alpha1
|
||||
return p1, alpha1
|
||||
|
||||
|
||||
def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> tuple[np.ndarray, float]:
|
||||
def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]:
|
||||
if s > np.pi / maxPartCircum * ell.ax:
|
||||
s /= 2
|
||||
point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum)
|
||||
@@ -251,14 +259,14 @@ def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,
|
||||
|
||||
return point_end, alpha_end
|
||||
|
||||
def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> tuple[float, float, float]:
|
||||
x, y, z = ell.para2cart(u, v)
|
||||
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, x, y, z)
|
||||
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, x, y, z)
|
||||
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)
|
||||
|
||||
@@ -267,49 +275,42 @@ def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float
|
||||
|
||||
return beta, lamb, alpha_ell
|
||||
|
||||
def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> tuple[float, float, float]:
|
||||
x, y, z = ell.ell2cart(beta, lamb)
|
||||
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, x, y, z)
|
||||
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, x, y, z)
|
||||
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-12:
|
||||
raise Exception("Alpha Umrechnung fehlgeschlagen")
|
||||
print("Alpha Umrechnung fehlgeschlagen:", np.linalg.norm(sigma_para - sigma_ell))
|
||||
|
||||
return u, v, alpha_para
|
||||
|
||||
def func_sigma_ell(ell, x, y, z, alpha):
|
||||
p, q = pq_ell(ell, x, y, z)
|
||||
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, x, y, z, alpha):
|
||||
p, q = pq_para(ell, x, y, z)
|
||||
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: np.ndarray, alpha: float) -> float:
|
||||
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
|
||||
# x, y, z = p0
|
||||
# t1, t2 = ell.func_t12(x, y, z)
|
||||
# l_cart = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2)
|
||||
# if abs(l - l_cart) > 1e-12:
|
||||
# # raise Exception("Louville constant fehlgeschlagen")
|
||||
# print("Diff zwischen constant:", abs(l - l_cart))
|
||||
return l
|
||||
|
||||
def louville_l2c(ell, l):
|
||||
def louville_l2c(ell: EllipsoidTriaxial, l: float) -> float:
|
||||
return sqrt((l + ell.Ee**2) / ell.Ex**2)
|
||||
|
||||
def louville_c2l(ell, c):
|
||||
def louville_c2l(ell: EllipsoidTriaxial, c: float) -> float:
|
||||
return ell.Ex**2 * c**2 - ell.Ee**2
|
||||
|
||||
|
||||
@@ -322,7 +323,7 @@ if __name__ == "__main__":
|
||||
# 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)
|
||||
# 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)
|
||||
@@ -342,7 +343,7 @@ if __name__ == "__main__":
|
||||
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, 5000)
|
||||
P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 5000)
|
||||
beta1_num, lamb1_num = ell.cart2ell(P1_num)
|
||||
|
||||
try:
|
||||
@@ -357,5 +358,3 @@ if __name__ == "__main__":
|
||||
mask_360 = (diffs_karney > 359) & (diffs_karney < 361)
|
||||
diffs_karney[mask_360] = np.abs(diffs_karney[mask_360] - 360)
|
||||
print(diffs_karney)
|
||||
|
||||
pass
|
||||
|
||||
@@ -4,9 +4,13 @@ import runge_kutta as rk
|
||||
import GHA_triaxial.numeric_examples_karney as ne_karney
|
||||
import GHA_triaxial.numeric_examples_panou as ne_panou
|
||||
import winkelumrechnungen as wu
|
||||
from typing import Tuple
|
||||
from numpy.typing import NDArray
|
||||
|
||||
# Panou 2013
|
||||
def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30):
|
||||
def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float,
|
||||
n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False
|
||||
) -> Tuple[float, float, float]| Tuple[float, float, float, NDArray, NDArray]:
|
||||
"""
|
||||
|
||||
:param ell: triaxiales Ellipsoid
|
||||
@@ -17,6 +21,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
:param n: Anzahl Schritte
|
||||
:param epsilon:
|
||||
:param iter_max: Maximale Anzhal Iterationen
|
||||
:param all_points:
|
||||
:return:
|
||||
"""
|
||||
|
||||
@@ -25,7 +30,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
def arccot(x):
|
||||
return np.arctan2(1.0, x)
|
||||
|
||||
|
||||
def BETA_LAMBDA(beta, lamb):
|
||||
|
||||
BETA = (ell.ay**2 * np.sin(beta)**2 + ell.b**2 * np.cos(beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)
|
||||
@@ -60,7 +64,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
G_beta_lamb = - LAMBDA_ * ell.Ey**2 * np.sin(2*beta)
|
||||
G_lamb_lamb = LAMBDA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + 2 * LAMBDA_ * ell.Ee**2 * np.sin(2*lamb) + 2 * LAMBDA * ell.Ee**2 * np.cos(2*lamb)
|
||||
|
||||
|
||||
return (BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
@@ -244,15 +247,21 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
+ (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2
|
||||
)
|
||||
|
||||
return alpha_1, alpha_2, s, beta_arr, lamb_arr
|
||||
if all_points:
|
||||
return alpha_1, alpha_2, s, beta_arr, lamb_arr
|
||||
else:
|
||||
return alpha_1, alpha_2, s
|
||||
|
||||
if lamb_1 == lamb_2:
|
||||
else: # lamb_1 == lamb_2
|
||||
|
||||
N = n
|
||||
dbeta = beta_2 - beta_1
|
||||
|
||||
if abs(dbeta) < 10**-15:
|
||||
return 0, 0, 0, np.array([]), np.array([])
|
||||
if all_points:
|
||||
return 0, 0, 0, np.array([]), np.array([])
|
||||
else:
|
||||
return 0, 0, 0
|
||||
|
||||
lamb_0 = 0
|
||||
|
||||
@@ -369,7 +378,10 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
else:
|
||||
s = np.trapz(integrand, dx=h)
|
||||
|
||||
return alpha_1, alpha_2, s, beta_arr, lamb_arr
|
||||
if all_points:
|
||||
return alpha_1, alpha_2, s, beta_arr, lamb_arr
|
||||
else:
|
||||
return alpha_1, alpha_2, s
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
@@ -390,7 +402,6 @@ if __name__ == "__main__":
|
||||
# a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
|
||||
# print(s)
|
||||
|
||||
|
||||
# ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
# diffs_panou = []
|
||||
# examples_panou = ne_panou.get_random_examples(4)
|
||||
@@ -421,8 +432,4 @@ if __name__ == "__main__":
|
||||
# diffs_karney = np.array(diffs_karney)
|
||||
# print(diffs_karney)
|
||||
|
||||
|
||||
pass
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user