Merge remote-tracking branch 'origin/main'
# Conflicts: # dashboard.py
This commit is contained in:
Tammo.Weber
17
GHA/rk.py
17
GHA/rk.py
@@ -1,17 +0,0 @@
|
||||
import Numerische_Integration.num_int_runge_kutta as rk
|
||||
from numpy import sin, cos, tan
|
||||
import winkelumrechnungen as wu
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
|
||||
def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num):
|
||||
phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0)
|
||||
|
||||
f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c
|
||||
f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c)
|
||||
f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c
|
||||
|
||||
funktionswerte = rk.verfahren([f_phi, f_lam, f_A],
|
||||
[0, phi0, lamb0, A0],
|
||||
s, num)
|
||||
coords = re.ell2cart(funktionswerte[-1][1], funktionswerte[-1][2], h0)
|
||||
return coords
|
||||
21
GHA_biaxial/rk.py
Normal file
21
GHA_biaxial/rk.py
Normal file
@@ -0,0 +1,21 @@
|
||||
import runge_kutta as rk
|
||||
from numpy import sin, cos, tan
|
||||
import winkelumrechnungen as wu
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
import numpy as np
|
||||
|
||||
def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num):
|
||||
phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0)
|
||||
|
||||
def buildODE():
|
||||
def ODE(s, v):
|
||||
phi, lam, A = v
|
||||
dphi = cos(A) * re.V(phi) ** 3 / re.c
|
||||
dlam = sin(A) * re.V(phi) / (cos(phi) * re.c)
|
||||
dA = tan(phi) * sin(A) * re.V(phi) / re.c
|
||||
return np.array([dphi, dlam, dA])
|
||||
return ODE
|
||||
|
||||
_, funktionswerte = rk.rk4(buildODE(), 0, np.array([phi0, lamb0, A0]), s, num)
|
||||
coords = re.ell2cart(funktionswerte[-1][0], funktionswerte[-1][1], h0)
|
||||
return coords
|
||||
64
GHA_triaxial/approx_gha1.py
Normal file
64
GHA_triaxial/approx_gha1.py
Normal file
@@ -0,0 +1,64 @@
|
||||
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
|
||||
|
||||
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]
|
||||
s_curr = 0.0
|
||||
while s_curr < s:
|
||||
ds_step = min(ds, s - s_curr)
|
||||
if ds_step < 1e-8:
|
||||
break
|
||||
p1 = points[-1]
|
||||
alpha1 = alphas[-1]
|
||||
sigma = func_sigma_ell(ell, p1, alpha1)
|
||||
p2 = p1 + ds_step * sigma
|
||||
p2 = ell.cartonell(p2)
|
||||
ds_step = np.linalg.norm(p2 - p1)
|
||||
|
||||
points.append(p2)
|
||||
dalpha = 1e-6
|
||||
l2 = louville_constant(ell, p2, alpha1)
|
||||
dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
|
||||
alpha2 = alpha1 + (l0 - l2) / dl_dalpha
|
||||
alphas.append(alpha2)
|
||||
s_curr += ds_step
|
||||
|
||||
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()
|
||||
|
||||
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, 0)
|
||||
alpha0 = wu.deg2rad(90)
|
||||
s = 1000000
|
||||
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)
|
||||
43
GHA_triaxial/numeric_examples_karney.py
Normal file
43
GHA_triaxial/numeric_examples_karney.py
Normal file
@@ -0,0 +1,43 @@
|
||||
import random
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
def line2example(line):
|
||||
split = line.split()
|
||||
example = [float(value) for value in split[:7]]
|
||||
for i, value in enumerate(example):
|
||||
if i < 6:
|
||||
example[i] = wu.deg2rad(value)
|
||||
# example[i] = value
|
||||
return example
|
||||
|
||||
def get_random_examples(num):
|
||||
"""
|
||||
beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12
|
||||
:param num:
|
||||
:return:
|
||||
"""
|
||||
# random.seed(42)
|
||||
with open("Karney_2024_Testset.txt") as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
for i in range(num):
|
||||
example = line2example(lines[random.randint(0, len(lines) - 1)])
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
def get_examples(l_i):
|
||||
"""
|
||||
beta0, lamb0, alpha0, beta1, lamb1, alpha1, s12
|
||||
:param num:
|
||||
:return:
|
||||
"""
|
||||
with open("Karney_2024_Testset.txt") as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
for i in l_i:
|
||||
example = line2example(lines[i])
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
if __name__ == "__main__":
|
||||
get_random_examples(10)
|
||||
@@ -1,4 +1,5 @@
|
||||
import winkelumrechnungen as wu
|
||||
import random
|
||||
|
||||
table1 = [
|
||||
(wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(0), wu.deg2rad(90), 1.00000000000,
|
||||
@@ -105,7 +106,22 @@ def get_example(table, example):
|
||||
def get_tables():
|
||||
return tables
|
||||
|
||||
def get_random_examples(num):
|
||||
random.seed(42)
|
||||
examples = []
|
||||
for i in range(num):
|
||||
table = random.randint(1, 4)
|
||||
if table == 4:
|
||||
example = random.randint(1, 8)
|
||||
else:
|
||||
example = random.randint(1, 7)
|
||||
example = get_example(table, example)
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
test = get_example(1, 4)
|
||||
# test = get_example(1, 4)
|
||||
examples = get_random_examples(5)
|
||||
pass
|
||||
|
||||
@@ -1,99 +1,149 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos, sqrt, arctan2
|
||||
import ellipsoide
|
||||
import Numerische_Integration.num_int_runge_kutta as rk
|
||||
import runge_kutta as rk
|
||||
import winkelumrechnungen as wu
|
||||
import ausgaben as aus
|
||||
import GHA.rk as ghark
|
||||
from scipy.special import factorial as fact
|
||||
from math import comb
|
||||
import GHA_triaxial.numeric_examples_panou as nep
|
||||
|
||||
# Panou, Korakitits 2019
|
||||
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 gha1_num_old(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num):
|
||||
phi, lamb, h = ell.cart2geod(point, "ligas3")
|
||||
x, y, z = ell.geod2cart(phi, lamb, 0)
|
||||
p, q = ell.p_q(x, y, z)
|
||||
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)
|
||||
|
||||
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)
|
||||
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
|
||||
|
||||
h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**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
|
||||
|
||||
f1 = lambda x, dxds, y, dyds, z, dzds: dxds
|
||||
f2 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * x
|
||||
f3 = lambda x, dxds, y, dyds, z, dzds: dyds
|
||||
f4 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * y/(1-ell.ee**2)
|
||||
f5 = lambda x, dxds, y, dyds, z, dzds: dzds
|
||||
f6 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * z/(1-ell.ex**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])
|
||||
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]])
|
||||
|
||||
funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [x, dxds0, y, dyds0, z, dzds0], s, num, fein=False)
|
||||
P2 = funktionswerte[-1]
|
||||
P2 = (P2[0], P2[2], P2[4])
|
||||
return P2
|
||||
return p, q
|
||||
|
||||
def buildODE(ell):
|
||||
def ODE(v):
|
||||
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(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
|
||||
ddy = -(h/H)*y/(1-ell.ee**2)
|
||||
ddz = -(h/H)*z/(1-ell.ex**2)
|
||||
|
||||
return [dxds, ddx, dyds, ddy, dzds, ddz]
|
||||
return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
|
||||
return ODE
|
||||
|
||||
def gha1_num(ell, point, alpha0, s, num):
|
||||
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")
|
||||
x0, y0, z0 = ell.geod2cart(phi, lam, 0)
|
||||
p0 = ell.geod2cart(phi, lam, 0)
|
||||
x0, y0, z0 = p0
|
||||
|
||||
p, q = ell.p_q(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)
|
||||
|
||||
v_init = [x0, dxds0, y0, dyds0, z0, dzds0]
|
||||
v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0])
|
||||
|
||||
F = buildODE(ell)
|
||||
ode = buildODE(ell)
|
||||
|
||||
werte = rk.rk_chat(F, v_init, s, num)
|
||||
x1, _, y1, _, z1, _ = werte[-1]
|
||||
_, werte = rk.rk4(ode, 0, v_init, s, num)
|
||||
x1, dx1ds, y1, dy1ds, z1, dz1ds = werte[-1]
|
||||
|
||||
return x1, y1, z1, werte
|
||||
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)
|
||||
|
||||
def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points):
|
||||
constantValues = []
|
||||
for point in points:
|
||||
x = point[1]
|
||||
dxds = point[2]
|
||||
y = point[3]
|
||||
dyds = point[4]
|
||||
z = point[5]
|
||||
dzds = point[6]
|
||||
alpha1 = arctan2(P, Q)
|
||||
|
||||
values = ell.p_q(x, y, z)
|
||||
p = values["p"]
|
||||
q = values["q"]
|
||||
t1 = values["t1"]
|
||||
t2 = values["t2"]
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
|
||||
P = p[0]*dxds + p[1]*dyds + p[2]*dzds
|
||||
Q = q[0]*dxds + q[1]*dyds + q[2]*dzds
|
||||
alpha = arctan2(P, Q)
|
||||
if all_points:
|
||||
return point1, alpha1, werte
|
||||
else:
|
||||
return point1, alpha1
|
||||
|
||||
c = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2)
|
||||
constantValues.append(c)
|
||||
pass
|
||||
# ---------------------------------------------------------------------------------------------------------------------
|
||||
|
||||
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)
|
||||
|
||||
def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
# 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")
|
||||
|
||||
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:
|
||||
@@ -104,44 +154,36 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
:return:
|
||||
"""
|
||||
x, y, z = point
|
||||
|
||||
# S. 6
|
||||
x_m = [x]
|
||||
y_m = [y]
|
||||
z_m = [z]
|
||||
|
||||
# erste Ableitungen (7-8)
|
||||
H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2
|
||||
sqrtH = sqrt(H)
|
||||
n = np.array([x / sqrtH,
|
||||
y / ((1-ell.ee**2) * sqrtH),
|
||||
z / ((1-ell.ex**2) * sqrtH)])
|
||||
u, v = ell.cart2para(np.array([x, y, z]))
|
||||
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]])
|
||||
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))
|
||||
|
||||
# H Ableitungen (7)
|
||||
H_ = lambda p: np.sum([comb(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)])
|
||||
# 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)])
|
||||
|
||||
# h Ableitungen (7)
|
||||
h_ = lambda q: np.sum([comb(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)])
|
||||
# 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)])
|
||||
|
||||
# h/H Ableitungen (6)
|
||||
hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) *
|
||||
np.sum([H_(t+1-l) / (fact(t+1-l) * fact(l-1)) * hH_t[l-1] for l in range(1, t+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)]))
|
||||
|
||||
# xm, ym, zm Ableitungen (6)
|
||||
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)])
|
||||
# 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)])
|
||||
|
||||
@@ -155,11 +197,14 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
x_m.append(x_(m))
|
||||
y_m.append(y_(m))
|
||||
z_m.append(z_(m))
|
||||
a_m.append(x_m[m] / fact(m))
|
||||
b_m.append(y_m[m] / fact(m))
|
||||
c_m.append(z_m[m] / fact(m))
|
||||
fact_m = fact(m)
|
||||
|
||||
# am, bm, cm (6)
|
||||
# 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
|
||||
@@ -170,25 +215,146 @@ def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
|
||||
for c in reversed(c_m):
|
||||
z_s = z_s * s + c
|
||||
|
||||
return x_s, y_s, z_s
|
||||
pass
|
||||
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 = 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)
|
||||
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("Analyitsche Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellpsoid")
|
||||
|
||||
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-12:
|
||||
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("Eitschberger1978")
|
||||
ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
# ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial")
|
||||
re = ellipsoide.EllipsoidBiaxial.init_name("Bessel")
|
||||
|
||||
# Panou 2013, 7, Table 1, beta0=60°
|
||||
beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = nep.get_example(table=1, example=3)
|
||||
P0 = ell.ell2cart(beta0, lamb0)
|
||||
P1 = ell.ell2cart(beta1, lamb1)
|
||||
# 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)
|
||||
|
||||
# P1_num = gha1_num(ell, P0, alpha0, s, 1000)
|
||||
P1_num = gha1_num(ell, P0, alpha0, s, 10000)
|
||||
P1_ana = gha1_ana(ell, P0, alpha0, s, 30)
|
||||
ell = 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)
|
||||
|
||||
beta, lamb = ellipsoide.EllipsoidTriaxial.cart2ell(ell, P1_num)
|
||||
P1_num, alpha1_num = gha1_num(ell, P0, alpha0_ell, s, 5000)
|
||||
beta1_num, lamb1_num = ell.cart2ell(P1_num)
|
||||
|
||||
pass
|
||||
try:
|
||||
_, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
|
||||
P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 30, maxPartCircum=16)
|
||||
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,10 +1,16 @@
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import Numerische_Integration.num_int_runge_kutta as rk
|
||||
import ausgaben as aus
|
||||
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
|
||||
@@ -15,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:
|
||||
"""
|
||||
|
||||
@@ -23,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)
|
||||
@@ -58,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,
|
||||
@@ -110,27 +115,45 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
q_33, q_22, q_11, q_00)
|
||||
|
||||
if lamb_1 != lamb_2:
|
||||
def functions():
|
||||
def f_beta(lamb, beta, beta_p, X3, X4):
|
||||
return beta_p
|
||||
# def functions():
|
||||
# def f_beta(lamb, beta, beta_p, X3, X4):
|
||||
# return beta_p
|
||||
#
|
||||
# def f_beta_p(lamb, beta, beta_p, X3, X4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# p_3, p_2, p_1, p_0,
|
||||
# p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
# return p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
#
|
||||
# def f_X3(lamb, beta, beta_p, X3, X4):
|
||||
# return X4
|
||||
#
|
||||
# def f_X4(lamb, beta, beta_p, X3, X4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# p_3, p_2, p_1, p_0,
|
||||
# p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
# return (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
# (3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
|
||||
#
|
||||
# return [f_beta, f_beta_p, f_X3, f_X4]
|
||||
|
||||
def buildODElamb():
|
||||
def ODE(lamb, v):
|
||||
beta, beta_p, X3, X4 = v
|
||||
|
||||
def f_beta_p(lamb, beta, beta_p, X3, X4):
|
||||
(BETA, LAMBDA, E, G,
|
||||
p_3, p_2, p_1, p_0,
|
||||
p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
return p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
|
||||
def f_X3(lamb, beta, beta_p, X3, X4):
|
||||
return X4
|
||||
|
||||
def f_X4(lamb, beta, beta_p, X3, X4):
|
||||
(BETA, LAMBDA, E, G,
|
||||
p_3, p_2, p_1, p_0,
|
||||
p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
return (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
dbeta = beta_p
|
||||
dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
dX3 = X4
|
||||
dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
(3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
|
||||
|
||||
return [f_beta, f_beta_p, f_X3, f_X4]
|
||||
return np.array([dbeta, dbeta_p, dX3, dX4])
|
||||
|
||||
return ODE
|
||||
|
||||
N = n
|
||||
|
||||
@@ -144,15 +167,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
converged = False
|
||||
iterations = 0
|
||||
|
||||
funcs = functions()
|
||||
# funcs = functions()
|
||||
ode_lamb = buildODElamb()
|
||||
|
||||
for i in range(iter_max):
|
||||
iterations = i + 1
|
||||
|
||||
startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0]
|
||||
# startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0]
|
||||
startwerte = np.array([beta_1, beta_0, 0.0, 1.0])
|
||||
|
||||
werte = rk.verfahren(funcs, startwerte, dlamb, N)
|
||||
lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1]
|
||||
# werte = rk.verfahren(funcs, startwerte, dlamb, N)
|
||||
lamb_list, werte = rk.rk4(ode_lamb, lamb_1, startwerte, dlamb, N, False)
|
||||
# lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1]
|
||||
lamb_end = lamb_list[-1]
|
||||
beta_end, beta_p_end, X3_end, X4_end = werte[-1]
|
||||
|
||||
d_beta_end_d_beta0 = X3_end
|
||||
delta = beta_end - beta_2
|
||||
@@ -174,16 +202,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
raise RuntimeError("konvergiert nicht.")
|
||||
|
||||
# Z
|
||||
werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N)
|
||||
# werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N, False)
|
||||
lamb_list, werte = rk.rk4(ode_lamb, lamb_1, np.array([beta_1, beta_0, 0.0, 1.0]), dlamb, N, False)
|
||||
|
||||
beta_arr = np.zeros(N + 1)
|
||||
lamb_arr = np.zeros(N + 1)
|
||||
# lamb_arr = np.zeros(N + 1)
|
||||
lamb_arr = np.array(lamb_list)
|
||||
beta_p_arr = np.zeros(N + 1)
|
||||
|
||||
for i, state in enumerate(werte):
|
||||
lamb_arr[i] = state[0]
|
||||
beta_arr[i] = state[1]
|
||||
beta_p_arr[i] = state[2]
|
||||
# lamb_arr[i] = state[0]
|
||||
# beta_arr[i] = state[1]
|
||||
# beta_p_arr[i] = state[2]
|
||||
beta_arr[i] = state[0]
|
||||
beta_p_arr[i] = state[1]
|
||||
|
||||
(_, _, E1, G1,
|
||||
*_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
|
||||
@@ -215,52 +247,80 @@ 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
|
||||
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
|
||||
if all_points:
|
||||
return 0, 0, 0, np.array([]), np.array([])
|
||||
else:
|
||||
return 0, 0, 0
|
||||
|
||||
lamb_0 = 0
|
||||
|
||||
converged = False
|
||||
iterations = 0
|
||||
|
||||
def functions_beta():
|
||||
def g_lamb(beta, lamb, lamb_p, Y3, Y4):
|
||||
return lamb_p
|
||||
# def functions_beta():
|
||||
# def g_lamb(beta, lamb, lamb_p, Y3, Y4):
|
||||
# return lamb_p
|
||||
#
|
||||
# def g_lamb_p(beta, lamb, lamb_p, Y3, Y4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# q_3, q_2, q_1, q_0,
|
||||
# q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
# return q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
#
|
||||
# def g_Y3(beta, lamb, lamb_p, Y3, Y4):
|
||||
# return Y4
|
||||
#
|
||||
# def g_Y4(beta, lamb, lamb_p, Y3, Y4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# q_3, q_2, q_1, q_0,
|
||||
# q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
# return (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
# (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
|
||||
#
|
||||
# return [g_lamb, g_lamb_p, g_Y3, g_Y4]
|
||||
|
||||
def buildODEbeta():
|
||||
def ODE(beta, v):
|
||||
lamb, lamb_p, Y3, Y4 = v
|
||||
|
||||
def g_lamb_p(beta, lamb, lamb_p, Y3, Y4):
|
||||
(BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
return q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
|
||||
def g_Y3(beta, lamb, lamb_p, Y3, Y4):
|
||||
return Y4
|
||||
|
||||
def g_Y4(beta, lamb, lamb_p, Y3, Y4):
|
||||
(BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
return (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
dlamb = lamb_p
|
||||
dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
dY3 = Y4
|
||||
dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
(3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
|
||||
|
||||
return [g_lamb, g_lamb_p, g_Y3, g_Y4]
|
||||
return np.array([dlamb, dlamb_p, dY3, dY4])
|
||||
return ODE
|
||||
|
||||
funcs_beta = functions_beta()
|
||||
# funcs_beta = functions_beta()
|
||||
ode_beta = buildODEbeta()
|
||||
|
||||
for i in range(iter_max):
|
||||
iterations = i + 1
|
||||
|
||||
startwerte = [beta_1, lamb_1, lamb_0, 0.0, 1.0]
|
||||
startwerte = [lamb_1, lamb_0, 0.0, 1.0]
|
||||
|
||||
werte = rk.verfahren(funcs_beta, startwerte, dbeta, N)
|
||||
beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||||
# werte = rk.verfahren(funcs_beta, startwerte, dbeta, N, False)
|
||||
beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
|
||||
|
||||
beta_end = beta_list[-1]
|
||||
# beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||||
lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||||
|
||||
d_lamb_end_d_lambda0 = Y3_end
|
||||
delta = lamb_end - lamb_2
|
||||
@@ -279,16 +339,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
|
||||
|
||||
lamb_0 = lamb_0 - step
|
||||
|
||||
werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N)
|
||||
# werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N, False)
|
||||
beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False)
|
||||
|
||||
beta_arr = np.zeros(N + 1)
|
||||
# beta_arr = np.zeros(N + 1)
|
||||
beta_arr = np.array(beta_list)
|
||||
lamb_arr = np.zeros(N + 1)
|
||||
lambda_p_arr = np.zeros(N + 1)
|
||||
|
||||
for i, state in enumerate(werte):
|
||||
beta_arr[i] = state[0]
|
||||
lamb_arr[i] = state[1]
|
||||
lambda_p_arr[i] = state[2]
|
||||
# beta_arr[i] = state[0]
|
||||
# lamb_arr[i] = state[1]
|
||||
# lambda_p_arr[i] = state[2]
|
||||
lamb_arr[i] = state[0]
|
||||
lambda_p_arr[i] = state[1]
|
||||
|
||||
# Azimute
|
||||
(BETA1, LAMBDA1, E1, G1,
|
||||
@@ -314,26 +378,58 @@ 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
|
||||
if all_points:
|
||||
return alpha_1, alpha_2, s, beta_arr, lamb_arr
|
||||
else:
|
||||
return alpha_1, alpha_2, s
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
# beta1 = np.deg2rad(75)
|
||||
# lamb1 = np.deg2rad(-90)
|
||||
# beta2 = np.deg2rad(75)
|
||||
# lamb2 = np.deg2rad(66)
|
||||
# a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
|
||||
# print(aus.gms("a1", a1, 4))
|
||||
# print(aus.gms("a2", a2, 4))
|
||||
# ell = EllipsoidTriaxial.init_name("Fiction")
|
||||
# # beta1 = np.deg2rad(75)
|
||||
# # lamb1 = np.deg2rad(-90)
|
||||
# # beta2 = np.deg2rad(75)
|
||||
# # lamb2 = np.deg2rad(66)
|
||||
# # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
|
||||
# # print(aus.gms("a1", a1, 4))
|
||||
# # print(aus.gms("a2", a2, 4))
|
||||
# # print(s)
|
||||
# cart1 = ell.para2cart(0, 0)
|
||||
# cart2 = ell.para2cart(0.4, 1.4)
|
||||
# beta1, lamb1 = ell.cart2ell(cart1)
|
||||
# beta2, lamb2 = ell.cart2ell(cart2)
|
||||
#
|
||||
# a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
|
||||
# print(s)
|
||||
cart1 = ell.para2cart(0, 0)
|
||||
cart2 = ell.para2cart(0.4, 0.4)
|
||||
beta1, lamb1 = ell.cart2ell(cart1)
|
||||
beta2, lamb2 = ell.cart2ell(cart2)
|
||||
|
||||
a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=2500)
|
||||
print(s)
|
||||
|
||||
|
||||
# ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
# diffs_panou = []
|
||||
# examples_panou = ne_panou.get_random_examples(4)
|
||||
# for example in examples_panou:
|
||||
# beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example
|
||||
# P0 = ell.ell2cart(beta0, lamb0)
|
||||
# try:
|
||||
# alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
|
||||
# diffs_panou.append(
|
||||
# (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num)))
|
||||
# except:
|
||||
# print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
|
||||
# diffs_panou = np.array(diffs_panou)
|
||||
# print(diffs_panou)
|
||||
#
|
||||
# ell = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
# diffs_karney = []
|
||||
# # examples_karney = ne_karney.get_examples((30500, 40500))
|
||||
# examples_karney = ne_karney.get_random_examples(2)
|
||||
# for example in examples_karney:
|
||||
# beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example
|
||||
#
|
||||
# try:
|
||||
# alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
|
||||
# diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num)))
|
||||
# except:
|
||||
# print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
|
||||
# diffs_karney = np.array(diffs_karney)
|
||||
# print(diffs_karney)
|
||||
|
||||
pass
|
||||
|
||||
@@ -1,84 +0,0 @@
|
||||
def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int, fein: bool = True) -> list:
|
||||
"""
|
||||
Runge-Kutta-Verfahren für ein beliebiges DGLS
|
||||
:param funktionen: Liste mit allen Funktionen
|
||||
:param startwerte: Liste mit allen Startwerten der Variablen
|
||||
:param weite: gesamte Weite über die integriert werden soll
|
||||
:param schritte: Anzahl der Schritte über die gesamte Weite
|
||||
:return: Liste mit Listen für alle Wertepaare
|
||||
"""
|
||||
h = weite / schritte
|
||||
werte = [startwerte]
|
||||
for i in range(schritte):
|
||||
|
||||
zuschlaege_grob = zuschlaege(funktionen, werte[-1], h)
|
||||
werte_grob = [werte[-1][j] if j == 0 else werte[-1][j] + zuschlaege_grob[j - 1]
|
||||
for j in range(len(startwerte))]
|
||||
if fein:
|
||||
zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2)
|
||||
werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1]
|
||||
for j in range(len(startwerte))]
|
||||
|
||||
zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2)
|
||||
werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1]
|
||||
for j in range(len(startwerte))]
|
||||
|
||||
werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j])
|
||||
for j in range(len(startwerte))]
|
||||
|
||||
werte.append(werte_korr)
|
||||
else:
|
||||
werte.append(werte_grob)
|
||||
return werte
|
||||
|
||||
|
||||
def zuschlaege(funktionen: list, startwerte: list, h: float) -> list:
|
||||
"""
|
||||
Berechnung der Zuschläge eines einzelnen Schritts
|
||||
:param funktionen: Liste mit allen Funktionen
|
||||
:param startwerte: Liste mit allen Startwerten der Variablen
|
||||
:param h: Schrittweite
|
||||
:return: Liste mit Zuschlägen für die einzelnen Variablen
|
||||
"""
|
||||
werte = [wert for wert in startwerte]
|
||||
|
||||
k1 = [h * funktion(*werte) for funktion in funktionen]
|
||||
|
||||
werte = [startwerte[i] + (h / 2 if i == 0 else k1[i - 1] / 2)
|
||||
for i in range(len(startwerte))]
|
||||
|
||||
k2 = [h * funktion(*werte) for funktion in funktionen]
|
||||
|
||||
werte = [startwerte[i] + (h / 2 if i == 0 else k2[i - 1] / 2)
|
||||
for i in range(len(startwerte))]
|
||||
|
||||
k3 = [h * funktion(*werte) for funktion in funktionen]
|
||||
|
||||
werte = [startwerte[i] + (h if i == 0 else k3[i - 1])
|
||||
for i in range(len(startwerte))]
|
||||
|
||||
k4 = [h * funktion(*werte) for funktion in funktionen]
|
||||
|
||||
k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))]
|
||||
|
||||
return k_
|
||||
|
||||
def rk_chat(F, v0: list, weite: float, schritte: int):
|
||||
h = weite/schritte
|
||||
v = v0
|
||||
werte = [v]
|
||||
|
||||
for _ in range(schritte):
|
||||
k1 = F(v)
|
||||
k2 = F([v[i] + 0.5 * h * k1[i] for i in range(6)])
|
||||
k3 = F([v[i] + 0.5 * h * k2[i] for i in range(6)])
|
||||
k4 = F([v[i] + h * k3[i] for i in range(6)])
|
||||
|
||||
v = [
|
||||
v[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i])
|
||||
for i in range(6)
|
||||
]
|
||||
|
||||
werte.append(v)
|
||||
|
||||
return werte
|
||||
@@ -1,46 +0,0 @@
|
||||
def f(xi, yi):
|
||||
ys = xi + yi
|
||||
return ys
|
||||
|
||||
|
||||
# Gegeben:
|
||||
x0 = 0
|
||||
y0 = 0
|
||||
h = 0.2
|
||||
xmax = 0.4
|
||||
Ey = 0.0001
|
||||
|
||||
x = [x0]
|
||||
y = [[y0]]
|
||||
|
||||
n = 0
|
||||
while x[-1] < xmax:
|
||||
x.append(x[-1]+h)
|
||||
fn = f(x[n], y[n][-1])
|
||||
if n == 0:
|
||||
y_neu = y[n][-1] + h * fn
|
||||
else:
|
||||
y_neu = y[n-1][-1] + 2*h * fn
|
||||
y.append([y_neu])
|
||||
dy = 1
|
||||
while dy > Ey:
|
||||
y_neu = y[n][-1] + h/2 * (fn + f(x[n+1], y[n+1][-1]))
|
||||
y[-1].append(y_neu)
|
||||
dy = abs(y[-1][-2]-y[-1][-1])
|
||||
n += 1
|
||||
|
||||
print(x)
|
||||
print(y)
|
||||
|
||||
werte = []
|
||||
for i in range(len(x)):
|
||||
werte.append((x[i], y[i][-1]))
|
||||
|
||||
for paar in werte:
|
||||
print(f"({round(paar[0], 5)}, {round(paar[1], 5)})")
|
||||
|
||||
integral = 0
|
||||
for i in range(len(werte)-1):
|
||||
integral += (werte[i+1][1]+werte[i][1])/2 * (werte[i+1][0]-werte[i][0])
|
||||
|
||||
print(f"Integral = {round(integral, 5)}")
|
||||
@@ -1,7 +0,0 @@
|
||||
import num_int_runge_kutta as rk
|
||||
|
||||
f = lambda ti, xi, yi: yi - ti
|
||||
g = lambda ti, xi, yi: xi + yi
|
||||
|
||||
funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 3)
|
||||
print(funktionswerte)
|
||||
@@ -1,6 +0,0 @@
|
||||
import num_int_runge_kutta as rk
|
||||
|
||||
f = lambda xi, yi: xi + yi
|
||||
|
||||
funktionswerte = rk.verfahren([f], [0, 0], 1, 5)
|
||||
print(funktionswerte)
|
||||
@@ -1,10 +0,0 @@
|
||||
import numpy as np
|
||||
import num_int_runge_kutta as rk
|
||||
|
||||
f = lambda ti, ui, phii: -4 * np.sin(phii)
|
||||
g = lambda ti, ui, phii: ui
|
||||
|
||||
funktionswerte = rk.verfahren([f, g], [0, 1, 0], 0.6, 3)
|
||||
|
||||
for wert in funktionswerte:
|
||||
print(f"t = {round(wert[0],1)}s -> phi = {round(wert[2],5)}, phip = {round(wert[1],5)}, v = {round(2.45 * wert[1],5)}")
|
||||
@@ -1,7 +0,0 @@
|
||||
import num_int_runge_kutta as rk
|
||||
|
||||
f = lambda xi, yi, ui: ui
|
||||
g = lambda xi, yi, ui: 4 * ui - 4 * yi
|
||||
|
||||
funktionswerte = rk.verfahren([f, g], [0, 0, 1], 0.6, 2)
|
||||
print(funktionswerte)
|
||||
659
ellipsoide.py
659
ellipsoide.py
@@ -1,8 +1,10 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arctan, arctan2, sqrt
|
||||
from numpy import sin, cos, arctan, arctan2, sqrt, pi, arccos
|
||||
import winkelumrechnungen as wu
|
||||
import ausgaben as aus
|
||||
import jacobian_Ligas
|
||||
import matplotlib.pyplot as plt
|
||||
from typing import Tuple
|
||||
from numpy.typing import NDArray
|
||||
|
||||
|
||||
class EllipsoidBiaxial:
|
||||
@@ -10,8 +12,8 @@ class EllipsoidBiaxial:
|
||||
self.a = a
|
||||
self.b = b
|
||||
self.c = a ** 2 / b
|
||||
self.e = np.sqrt(a ** 2 - b ** 2) / a
|
||||
self.e_ = np.sqrt(a ** 2 - b ** 2) / b
|
||||
self.e = sqrt(a ** 2 - b ** 2) / a
|
||||
self.e_ = sqrt(a ** 2 - b ** 2) / b
|
||||
|
||||
@classmethod
|
||||
def init_name(cls, name: str):
|
||||
@@ -40,28 +42,28 @@ class EllipsoidBiaxial:
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
|
||||
V = lambda self, phi: np.sqrt(1 + self.e_ ** 2 * np.cos(phi) ** 2)
|
||||
V = lambda self, phi: sqrt(1 + self.e_ ** 2 * cos(phi) ** 2)
|
||||
M = lambda self, phi: self.c / self.V(phi) ** 3
|
||||
N = lambda self, phi: self.c / self.V(phi)
|
||||
|
||||
beta2psi = lambda self, beta: np.arctan(self.a / self.b * np.tan(beta))
|
||||
beta2phi = lambda self, beta: np.arctan(self.a ** 2 / self.b ** 2 * np.tan(beta))
|
||||
beta2psi = lambda self, beta: arctan(self.a / self.b * np.tan(beta))
|
||||
beta2phi = lambda self, beta: arctan(self.a ** 2 / self.b ** 2 * np.tan(beta))
|
||||
|
||||
psi2beta = lambda self, psi: np.arctan(self.b / self.a * np.tan(psi))
|
||||
psi2phi = lambda self, psi: np.arctan(self.a / self.b * np.tan(psi))
|
||||
psi2beta = lambda self, psi: arctan(self.b / self.a * np.tan(psi))
|
||||
psi2phi = lambda self, psi: arctan(self.a / self.b * np.tan(psi))
|
||||
|
||||
phi2beta = lambda self, phi: np.arctan(self.b ** 2 / self.a ** 2 * np.tan(phi))
|
||||
phi2psi = lambda self, phi: np.arctan(self.b / self.a * np.tan(phi))
|
||||
phi2beta = lambda self, phi: arctan(self.b ** 2 / self.a ** 2 * np.tan(phi))
|
||||
phi2psi = lambda self, phi: arctan(self.b / self.a * np.tan(phi))
|
||||
|
||||
phi2p = lambda self, phi: self.N(phi) * np.cos(phi)
|
||||
phi2p = lambda self, phi: self.N(phi) * cos(phi)
|
||||
|
||||
def cart2ell(self, Eh, Ephi, x, y, z):
|
||||
p = np.sqrt(x**2+y**2)
|
||||
p = sqrt(x**2+y**2)
|
||||
# print(f"p = {round(p, 5)} m")
|
||||
|
||||
lamb = np.arctan(y/x)
|
||||
lamb = arctan(y/x)
|
||||
|
||||
phi_null = np.arctan(z/p*(1-self.e**2)**-1)
|
||||
phi_null = arctan(z/p*(1-self.e**2)**-1)
|
||||
|
||||
hi = [0]
|
||||
phii = [phi_null]
|
||||
@@ -69,9 +71,9 @@ class EllipsoidBiaxial:
|
||||
i = 0
|
||||
|
||||
while True:
|
||||
N = self.a*(1-self.e**2*np.sin(phii[i])**2)**(-1/2)
|
||||
h = p/np.cos(phii[i])-N
|
||||
phi = np.arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1))
|
||||
N = self.a*(1-self.e**2*sin(phii[i])**2)**(-1/2)
|
||||
h = p/cos(phii[i])-N
|
||||
phi = arctan(z/p*(1-(self.e**2*N)/(N+h))**(-1))
|
||||
hi.append(h)
|
||||
phii.append(phi)
|
||||
dh = abs(hi[i]-h)
|
||||
@@ -86,11 +88,11 @@ class EllipsoidBiaxial:
|
||||
return phi, lamb, h
|
||||
|
||||
def ell2cart(self, phi, lamb, h):
|
||||
W = np.sqrt(1 - self.e**2 * np.sin(phi)**2)
|
||||
W = sqrt(1 - self.e**2 * sin(phi)**2)
|
||||
N = self.a / W
|
||||
x = (N+h) * np.cos(phi) * np.cos(lamb)
|
||||
y = (N+h) * np.cos(phi) * np.sin(lamb)
|
||||
z = (N * (1-self.e**2) + h) * np.sin(lamb)
|
||||
x = (N+h) * cos(phi) * cos(lamb)
|
||||
y = (N+h) * cos(phi) * sin(lamb)
|
||||
z = (N * (1-self.e**2) + h) * sin(lamb)
|
||||
return x, y, z
|
||||
|
||||
class EllipsoidTriaxial:
|
||||
@@ -98,15 +100,15 @@ class EllipsoidTriaxial:
|
||||
self.ax = ax
|
||||
self.ay = ay
|
||||
self.b = b
|
||||
self.ex = np.sqrt((self.ax**2 - self.b**2) / self.ax**2)
|
||||
self.ey = np.sqrt((self.ay**2 - self.b**2) / self.ay**2)
|
||||
self.ee = np.sqrt((self.ax**2 - self.ay**2) / self.ax**2)
|
||||
self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2)
|
||||
self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2)
|
||||
self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2)
|
||||
self.Ex = np.sqrt(self.ax**2 - self.b**2)
|
||||
self.Ey = np.sqrt(self.ay**2 - self.b**2)
|
||||
self.Ee = np.sqrt(self.ax**2 - self.ay**2)
|
||||
self.ex = sqrt((self.ax**2 - self.b**2) / self.ax**2)
|
||||
self.ey = sqrt((self.ay**2 - self.b**2) / self.ay**2)
|
||||
self.ee = sqrt((self.ax**2 - self.ay**2) / self.ax**2)
|
||||
self.ex_ = sqrt((self.ax**2 - self.b**2) / self.b**2)
|
||||
self.ey_ = sqrt((self.ay**2 - self.b**2) / self.b**2)
|
||||
self.ee_ = sqrt((self.ax**2 - self.ay**2) / self.ay**2)
|
||||
self.Ex = sqrt(self.ax**2 - self.b**2)
|
||||
self.Ey = sqrt(self.ay**2 - self.b**2)
|
||||
self.Ee = sqrt(self.ax**2 - self.ay**2)
|
||||
|
||||
@classmethod
|
||||
def init_name(cls, name: str):
|
||||
@@ -156,24 +158,56 @@ class EllipsoidTriaxial:
|
||||
b = 4000000
|
||||
return cls(ax, ay, b)
|
||||
elif name == "KarneyTest2024":
|
||||
ax = np.sqrt(2)
|
||||
ax = sqrt(2)
|
||||
ay = 1
|
||||
b = 1 / np.sqrt(2)
|
||||
b = 1 / sqrt(2)
|
||||
return cls(ax, ay, b)
|
||||
|
||||
def point_on(self, point: np.ndarray) -> bool:
|
||||
def func_H(self, point: NDArray) -> float:
|
||||
"""
|
||||
Test, ob ein Punkt auf dem Ellipsoid liegt.
|
||||
:param point: kartesische 3D-Koordinaten
|
||||
:return: Punkt auf dem Ellispoid?
|
||||
Berechnung H
|
||||
Panou, Korakitis 2019 [43]
|
||||
:param point: Punkt
|
||||
:return: H
|
||||
"""
|
||||
value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2
|
||||
if abs(1-value) < 0.000001:
|
||||
return True
|
||||
else:
|
||||
return False
|
||||
x, y, z = point
|
||||
return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2
|
||||
|
||||
def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
|
||||
def func_n(self, point: NDArray, H: float = None) -> NDArray:
|
||||
"""
|
||||
Berechnung normalen Vektor
|
||||
Panou, Korakitis 2019 [9-12]
|
||||
:param point: Punkt
|
||||
:param H:
|
||||
:return:
|
||||
"""
|
||||
if H is None:
|
||||
H = self.func_H(point)
|
||||
sqrtH = sqrt(H)
|
||||
x, y, z = point
|
||||
return np.array([x / sqrtH,
|
||||
y / ((1 - self.ee ** 2) * sqrtH),
|
||||
z / ((1 - self.ex ** 2) * sqrtH)])
|
||||
|
||||
def func_t12(self, point: NDArray) -> Tuple[float, float]:
|
||||
"""
|
||||
Berechnung Wurzeln
|
||||
Panou, Korakitis 2019 [9-12]
|
||||
:param point: Punkt
|
||||
:return: Wurzeln t1, t2
|
||||
"""
|
||||
x, y, z = point
|
||||
c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
|
||||
c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
|
||||
(self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
|
||||
self.ax ** 2 + self.ay ** 2) * z ** 2)
|
||||
t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
|
||||
if t2 == 0:
|
||||
t2 = 1e-18
|
||||
t1 = c0 / t2
|
||||
return t1, t2
|
||||
|
||||
def ellu2cart(self, beta: float, lamb: float, u: float) -> NDArray:
|
||||
"""
|
||||
Panou 2014 12ff.
|
||||
Elliptische Breite+Länge sind nicht gleich der geodätischen
|
||||
@@ -184,61 +218,13 @@ class EllipsoidTriaxial:
|
||||
:param u: Größe entlang der z-Achse
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
# s1 = u**2 - self.b**2
|
||||
# s2 = -self.ay**2 * np.sin(beta)**2 - self.b**2 * np.cos(beta)**2
|
||||
# s3 = -self.ax**2 * np.sin(lamb)**2 - self.ay**2 * np.cos(lamb)**2
|
||||
# print(s1, s2, s3)
|
||||
|
||||
# xe = np.sqrt(((self.ax**2+s1) * (self.ax**2+s2) * (self.ax**2+s3)) /
|
||||
# ((self.ax**2-self.ay**2) * (self.ax**2-self.b**2)))
|
||||
# ye = np.sqrt(((self.ay**2+s1) * (self.ay**2+s2) * (self.ay**2+s3)) /
|
||||
# ((self.ay**2-self.ax**2) * (self.ay**2-self.b**2)))
|
||||
# ze = np.sqrt(((self.b**2+s1) * (self.b**2+s2) * (self.b**2+s3)) /
|
||||
# ((self.b**2-self.ax**2) * (self.b**2-self.ay**2)))
|
||||
|
||||
x = np.sqrt(u**2 + self.Ex**2) * np.sqrt(np.cos(beta)**2 + self.Ee**2/self.Ex**2 * np.sin(beta)**2) * np.cos(lamb)
|
||||
y = np.sqrt(u**2 + self.Ey**2) * np.cos(beta) * np.sin(lamb)
|
||||
z = u * np.sin(beta) * np.sqrt(1 - self.Ee**2/self.Ex**2 * np.cos(lamb)**2)
|
||||
x = sqrt(u**2 + self.Ex**2) * sqrt(cos(beta)**2 + self.Ee**2/self.Ex**2 * sin(beta)**2) * cos(lamb)
|
||||
y = sqrt(u**2 + self.Ey**2) * cos(beta) * sin(lamb)
|
||||
z = u * sin(beta) * sqrt(1 - self.Ee**2/self.Ex**2 * cos(lamb)**2)
|
||||
|
||||
return np.array([x, y, z])
|
||||
|
||||
def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray:
|
||||
"""
|
||||
Panou, Korakitis 2019 2
|
||||
:param beta: elliptische Breite [rad]
|
||||
:param lamb: elliptische Länge [rad]
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
beta = np.asarray(beta, dtype=float)
|
||||
lamb = np.asarray(lamb, dtype=float)
|
||||
|
||||
beta, lamb = np.broadcast_arrays(beta, lamb)
|
||||
|
||||
B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2
|
||||
L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2
|
||||
|
||||
x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb)
|
||||
y = self.ay * np.cos(beta) * np.sin(lamb)
|
||||
z = self.b / self.Ex * np.sin(beta) * np.sqrt(L)
|
||||
|
||||
xyz = np.stack((x, y, z), axis=-1)
|
||||
|
||||
# Pole
|
||||
mask_south = beta == -np.pi / 2
|
||||
mask_north = beta == np.pi / 2
|
||||
xyz[mask_south] = np.array([0, 0, -self.b])
|
||||
xyz[mask_north] = np.array([0, 0, self.b])
|
||||
|
||||
# Äquator
|
||||
mask_eq = beta == 0
|
||||
xyz[mask_eq & (lamb == -np.pi / 2)] = np.array([0, -self.ay, 0])
|
||||
xyz[mask_eq & (lamb == np.pi / 2)] = np.array([0, self.ay, 0])
|
||||
xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0])
|
||||
xyz[mask_eq & (lamb == np.pi)] = np.array([-self.ax, 0, 0])
|
||||
|
||||
return xyz
|
||||
|
||||
def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]:
|
||||
def cart2ellu(self, point: NDArray) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Panou 2014 15ff.
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
@@ -253,25 +239,137 @@ class EllipsoidTriaxial:
|
||||
|
||||
p = (c2**2 - 3*c1) / 9
|
||||
q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54
|
||||
omega = np.arccos(q / np.sqrt(p**3))
|
||||
omega = arccos(q / sqrt(p**3))
|
||||
|
||||
s1 = 2 * np.sqrt(p) * np.cos(omega/3) - c2/3
|
||||
s2 = 2 * np.sqrt(p) * np.cos(omega/3 - 2*np.pi/3) - c2/3
|
||||
s3 = 2 * np.sqrt(p) * np.cos(omega/3 - 4*np.pi/3) - c2/3
|
||||
s1 = 2 * sqrt(p) * cos(omega/3) - c2/3
|
||||
s2 = 2 * sqrt(p) * cos(omega/3 - 2*pi/3) - c2/3
|
||||
s3 = 2 * sqrt(p) * cos(omega/3 - 4*pi/3) - c2/3
|
||||
# print(s1, s2, s3)
|
||||
|
||||
beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
|
||||
beta = arctan(sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
|
||||
if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7:
|
||||
lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
|
||||
lamb = arctan(sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
|
||||
else:
|
||||
lamb = 0
|
||||
u = np.sqrt(self.b**2 + s1)
|
||||
u = sqrt(self.b**2 + s1)
|
||||
|
||||
return beta, lamb, u
|
||||
|
||||
def cart2ell(self, point: np.ndarray) -> tuple[float, float]:
|
||||
def ell2cart(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray:
|
||||
"""
|
||||
Panou, Korakitis 2019 2f.
|
||||
Panou, Korakitis 2019 2
|
||||
:param beta: elliptische Breite [rad]
|
||||
:param lamb: elliptische Länge [rad]
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
beta = np.asarray(beta, dtype=float)
|
||||
lamb = np.asarray(lamb, dtype=float)
|
||||
|
||||
beta, lamb = np.broadcast_arrays(beta, lamb)
|
||||
|
||||
B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2
|
||||
L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2
|
||||
|
||||
x = self.ax / self.Ex * sqrt(B) * cos(lamb)
|
||||
y = self.ay * cos(beta) * sin(lamb)
|
||||
z = self.b / self.Ex * sin(beta) * sqrt(L)
|
||||
|
||||
xyz = np.stack((x, y, z), axis=-1)
|
||||
|
||||
# Pole
|
||||
mask_south = beta == -pi / 2
|
||||
mask_north = beta == pi / 2
|
||||
xyz[mask_south] = np.array([0, 0, -self.b])
|
||||
xyz[mask_north] = np.array([0, 0, self.b])
|
||||
|
||||
# Äquator
|
||||
mask_eq = beta == 0
|
||||
xyz[mask_eq & (lamb == -pi / 2)] = np.array([0, -self.ay, 0])
|
||||
xyz[mask_eq & (lamb == pi / 2)] = np.array([0, self.ay, 0])
|
||||
xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0])
|
||||
xyz[mask_eq & (lamb == pi)] = np.array([-self.ax, 0, 0])
|
||||
|
||||
return xyz
|
||||
|
||||
def ell2cart_bektas(self, beta: float | NDArray, omega: float | NDArray) -> NDArray:
|
||||
"""
|
||||
Bektas 2015
|
||||
:param beta:
|
||||
:param omega:
|
||||
:return:
|
||||
"""
|
||||
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))
|
||||
y = self.ay * cos(beta) * sin(omega)
|
||||
z = self.b * sin(beta) * sqrt((self.ax**2 * sin(omega)**2 + self.ay**2 * cos(omega)**2 - self.b**2) / (self.ax**2 - self.b**2))
|
||||
|
||||
return np.array([x, y, z])
|
||||
|
||||
def ell2cart_karney(self, beta: float | NDArray, lamb: float | NDArray) -> NDArray:
|
||||
"""
|
||||
Karney 2025 Geographic Lib
|
||||
:param beta:
|
||||
:param lamb:
|
||||
:return:
|
||||
"""
|
||||
k = sqrt(self.ay**2 - self.b**2) / sqrt(self.ax**2 - self.b**2)
|
||||
k_ = sqrt(self.ax**2 - self.ay**2) / sqrt(self.ax**2 - self.b**2)
|
||||
X = self.ax * cos(lamb) * sqrt(k**2*cos(beta)**2+k_**2)
|
||||
Y = self.ay * cos(beta) * sin(lamb)
|
||||
Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2)
|
||||
return np.array([X, Y, Z])
|
||||
|
||||
def cart2ell(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]:
|
||||
"""
|
||||
Panou, Korakitis 2019 3f. (num)
|
||||
:param point:
|
||||
:param eps:
|
||||
:param maxI:
|
||||
:return:
|
||||
"""
|
||||
x, y, z = point
|
||||
beta, lamb = self.cart2ell_panou(point)
|
||||
delta_ell = np.array([np.inf, np.inf]).T
|
||||
|
||||
i = 0
|
||||
while np.sum(delta_ell) > eps and i < maxI:
|
||||
x0, y0, z0 = self.ell2cart(beta, lamb)
|
||||
delta_l = np.array([x-x0, y-y0, z-z0]).T
|
||||
|
||||
B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2
|
||||
L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2
|
||||
|
||||
J = np.array([[(-self.ax * self.Ey ** 2) / (2 * self.Ex) * sin(2 * beta) / sqrt(B) * cos(lamb),
|
||||
-self.ax / self.Ex * sqrt(B) * sin(lamb)],
|
||||
[-self.ay * sin(beta) * sin(lamb),
|
||||
self.ay * cos(beta) * cos(lamb)],
|
||||
[self.b / self.Ex * cos(beta) * sqrt(L),
|
||||
(self.b * self.Ee ** 2) / (2 * self.Ex) * sin(beta) * sin(2 * lamb) / sqrt(L)]])
|
||||
|
||||
N = J.T @ J
|
||||
det = N[0, 0] * N[1, 1] - N[0, 1] * N[1, 0]
|
||||
if abs(det) < eps:
|
||||
det = eps
|
||||
N_inv = 1 / det * np.array([[N[1, 1], -N[0, 1]], [-N[1, 0], N[0, 0]]])
|
||||
delta_ell = N_inv @ J.T @ delta_l
|
||||
beta += delta_ell[0]
|
||||
lamb += delta_ell[1]
|
||||
i += 1
|
||||
|
||||
if i == maxI:
|
||||
raise Exception("Umrechung ist nicht konvergiert")
|
||||
|
||||
point_n = self.ell2cart(beta, lamb)
|
||||
delta_r = np.linalg.norm(point - point_n, axis=-1)
|
||||
|
||||
if delta_r > 1e-4:
|
||||
# raise Exception("Fehler in der Umrechnung cart2ell")
|
||||
print(f"Fehler in der Umrechnung cart2ell, deltaR = {delta_r}m")
|
||||
|
||||
return beta, lamb
|
||||
|
||||
def cart2ell_panou(self, point: NDArray) -> Tuple[float, float]:
|
||||
"""
|
||||
Panou, Korakitis 2019 2f. (analytisch -> Näherung)
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:return: elliptische Breite, elliptische Länge
|
||||
"""
|
||||
@@ -280,35 +378,31 @@ class EllipsoidTriaxial:
|
||||
eps = 1e-9
|
||||
|
||||
if abs(x) < eps and abs(y) < eps: # Punkt in der z-Achse
|
||||
beta = np.pi / 2 if z > 0 else -np.pi / 2
|
||||
beta = pi / 2 if z > 0 else -pi / 2
|
||||
lamb = 0.0
|
||||
return beta, lamb
|
||||
|
||||
elif abs(x) < eps and abs(z) < eps: # Punkt in der y-Achse
|
||||
beta = 0.0
|
||||
lamb = np.pi / 2 if y > 0 else -np.pi / 2
|
||||
lamb = pi / 2 if y > 0 else -pi / 2
|
||||
return beta, lamb
|
||||
|
||||
elif abs(y) < eps and abs(z) < eps: # Punkt in der x-Achse
|
||||
beta = 0.0
|
||||
lamb = 0.0 if x > 0 else np.pi
|
||||
lamb = 0.0 if x > 0 else pi
|
||||
return beta, lamb
|
||||
|
||||
# ---- Allgemeiner Fall -----
|
||||
|
||||
c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
|
||||
c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
|
||||
(self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
|
||||
self.ax ** 2 + self.ay ** 2) * z ** 2)
|
||||
t2 = (-c1 + np.sqrt(c1 ** 2 - 4 * c0)) / 2
|
||||
t1 = c0 / t2
|
||||
t1, t2 = self.func_t12(point)
|
||||
|
||||
num_beta = max(t1 - self.b ** 2, 0)
|
||||
den_beta = max(self.ay ** 2 - t1, 0)
|
||||
num_lamb = max(t2 - self.ay ** 2, 0)
|
||||
den_lamb = max(self.ax ** 2 - t2, 0)
|
||||
beta = np.arctan(np.sqrt(num_beta / den_beta))
|
||||
lamb = np.arctan(np.sqrt(num_lamb / den_lamb))
|
||||
|
||||
beta = arctan(sqrt(num_beta / den_beta))
|
||||
lamb = arctan(sqrt(num_lamb / den_lamb))
|
||||
|
||||
if z < 0:
|
||||
beta = -beta
|
||||
@@ -317,20 +411,63 @@ class EllipsoidTriaxial:
|
||||
lamb = -lamb
|
||||
|
||||
if x < 0:
|
||||
lamb = np.pi - lamb
|
||||
lamb = pi - lamb
|
||||
|
||||
if abs(x) < eps:
|
||||
lamb = -np.pi/2 if y < 0 else np.pi/2
|
||||
lamb = -pi/2 if y < 0 else pi/2
|
||||
|
||||
elif abs(y) < eps:
|
||||
lamb = 0 if x > 0 else np.pi
|
||||
lamb = 0 if x > 0 else pi
|
||||
|
||||
elif abs(z) < eps:
|
||||
beta = 0
|
||||
|
||||
return beta, lamb
|
||||
|
||||
def cart2geod(self, point: np.ndarray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
|
||||
def cart2ell_bektas(self, point: NDArray, eps: float = 1e-12, maxI: int = 100) -> Tuple[float, float]:
|
||||
"""
|
||||
Bektas 2015
|
||||
:param point:
|
||||
:param eps:
|
||||
:param maxI:
|
||||
:return:
|
||||
"""
|
||||
x, y, z = point
|
||||
phi, lamb = self.cart2para(point)
|
||||
p = sqrt((self.ax**2 - self.ay**2) / (self.ax**2 - self.b**2))
|
||||
d_phi = np.inf
|
||||
d_lamb = np.inf
|
||||
|
||||
i = 0
|
||||
while d_phi > eps and d_lamb > eps and i < maxI:
|
||||
lamb_new = arctan2(self.ax * y * sqrt((p**2-1) * sin(phi)**2 + 1), self.ay * x * cos(phi))
|
||||
phi_new = arctan2(self.ay * z * sin(lamb), self.b * y * sqrt(1 - p**2 * cos(lamb)**2))
|
||||
d_phi = abs(phi_new - phi)
|
||||
phi = phi_new
|
||||
d_lamb = abs(lamb_new - lamb)
|
||||
lamb = lamb_new
|
||||
i += 1
|
||||
|
||||
if i == maxI:
|
||||
raise Exception("Umrechung ist nicht konvergiert")
|
||||
|
||||
return phi, lamb
|
||||
|
||||
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
|
||||
:return: kartesische Koordinaten
|
||||
"""
|
||||
v = self.ax / sqrt(1 - self.ex**2*sin(phi)**2-self.ee**2*cos(phi)**2*sin(lamb)**2)
|
||||
xG = (v + h) * cos(phi) * cos(lamb)
|
||||
yG = (v * (1-self.ee**2) + h) * cos(phi) * sin(lamb)
|
||||
zG = (v * (1-self.ex**2) + h) * sin(phi)
|
||||
return np.array([xG, yG, zG])
|
||||
|
||||
def cart2geod(self, point: NDArray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Ligas 2012
|
||||
:param mode: ligas1, ligas2, oder ligas3
|
||||
@@ -344,24 +481,24 @@ class EllipsoidTriaxial:
|
||||
eps = 1e-9
|
||||
|
||||
if abs(xG) < eps and abs(yG) < eps: # Punkt in der z-Achse
|
||||
phi = np.pi / 2 if zG > 0 else -np.pi / 2
|
||||
phi = pi / 2 if zG > 0 else -pi / 2
|
||||
lamb = 0.0
|
||||
h = abs(zG) - self.b
|
||||
return phi, lamb, h
|
||||
|
||||
elif abs(xG) < eps and abs(zG) < eps: # Punkt in der y-Achse
|
||||
phi = 0.0
|
||||
lamb = np.pi / 2 if yG > 0 else -np.pi / 2
|
||||
lamb = pi / 2 if yG > 0 else -pi / 2
|
||||
h = abs(yG) - self.ay
|
||||
return phi, lamb, h
|
||||
|
||||
elif abs(yG) < eps and abs(zG) < eps: # Punkt in der x-Achse
|
||||
phi = 0.0
|
||||
lamb = 0.0 if xG > 0 else np.pi
|
||||
lamb = 0.0 if xG > 0 else pi
|
||||
h = abs(xG) - self.ax
|
||||
return phi, lamb, h
|
||||
|
||||
rG = np.sqrt(xG ** 2 + yG ** 2 + zG ** 2)
|
||||
rG = sqrt(xG ** 2 + yG ** 2 + zG ** 2)
|
||||
pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
|
||||
|
||||
E = 1 / self.ax**2
|
||||
@@ -380,50 +517,110 @@ class EllipsoidTriaxial:
|
||||
invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE)
|
||||
pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1)
|
||||
pEi = pEi.reshape(1, -1).flatten()
|
||||
loa = np.sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
|
||||
loa = sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
|
||||
pE = pEi
|
||||
i += 1
|
||||
|
||||
phi = np.arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / np.sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2))
|
||||
lamb = np.arctan(1/(1-self.ee**2) * pE[1]/pE[0])
|
||||
h = np.sign(zG - pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2)
|
||||
phi = arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2))
|
||||
lamb = arctan(1/(1-self.ee**2) * pE[1]/pE[0])
|
||||
h = np.sign(zG - pE[2]) * np.sign(pE[2]) * sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2)
|
||||
|
||||
if xG < 0 and yG < 0:
|
||||
lamb = -np.pi + lamb
|
||||
lamb = -pi + lamb
|
||||
|
||||
elif xG < 0:
|
||||
lamb = np.pi + lamb
|
||||
lamb = pi + lamb
|
||||
|
||||
if abs(zG) < eps:
|
||||
phi = 0
|
||||
|
||||
return phi, lamb, h
|
||||
|
||||
def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray:
|
||||
def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray:
|
||||
"""
|
||||
Ligas 2012, 250
|
||||
:param phi: geodätische Breite [rad]
|
||||
:param lamb: geodätische Länge [rad]
|
||||
:param h: Höhe über dem Ellipsoid
|
||||
:return: kartesische Koordinaten
|
||||
Panou, Korakitits 2020, 4
|
||||
:param u: Parameter u
|
||||
:param v: Parameter v
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
v = self.ax / np.sqrt(1 - self.ex**2*np.sin(phi)**2-self.ee**2*np.cos(phi)**2*np.sin(lamb)**2)
|
||||
xG = (v + h) * np.cos(phi) * np.cos(lamb)
|
||||
yG = (v * (1-self.ee**2) + h) * np.cos(phi) * np.sin(lamb)
|
||||
zG = (v * (1-self.ex**2) + h) * np.sin(phi)
|
||||
return np.array([xG, yG, zG])
|
||||
x = self.ax * cos(u) * cos(v)
|
||||
y = self.ay * cos(u) * sin(v)
|
||||
z = self.b * sin(u)
|
||||
z = np.broadcast_to(z, np.shape(x))
|
||||
return np.array([x, y, z])
|
||||
|
||||
def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]:
|
||||
def cart2para(self, point: NDArray) -> Tuple[float, float]:
|
||||
"""
|
||||
Panou, Korakitits 2020, 4
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:return: parametrische Koordinaten
|
||||
"""
|
||||
x, y, z = point
|
||||
|
||||
u_check1 = z*sqrt(1 - self.ee**2)
|
||||
u_check2 = sqrt(x**2 * (1-self.ee**2) + y**2) * sqrt(1-self.ex**2)
|
||||
if u_check1 <= u_check2:
|
||||
u = arctan2(u_check1, u_check2)
|
||||
else:
|
||||
u = pi/2 - arctan2(u_check2, u_check1)
|
||||
|
||||
v_check1 = y
|
||||
v_check2 = x*sqrt(1-self.ee**2)
|
||||
v_factor = sqrt(x**2*(1-self.ee**2)+y**2)
|
||||
if v_check1 <= v_check2:
|
||||
v = 2 * arctan2(v_check1, v_check2 + v_factor)
|
||||
else:
|
||||
v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor)
|
||||
|
||||
return u, v
|
||||
|
||||
def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]:
|
||||
cart = self.ell2cart(beta, lamb)
|
||||
return self.cart2para(cart)
|
||||
|
||||
def para2ell(self, u: float, v: float) -> Tuple[float, float]:
|
||||
cart = self.para2cart(u, v)
|
||||
return self.cart2ell(cart)
|
||||
|
||||
def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> Tuple[float, float, float]:
|
||||
cart = self.para2cart(u, v)
|
||||
return self.cart2geod(cart, mode, maxIter, maxLoa)
|
||||
|
||||
def geod2para(self, phi: float, lamb: float, h: float) -> Tuple[float, float]:
|
||||
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]:
|
||||
cart = self.ell2cart(beta, lamb)
|
||||
return self.cart2geod(cart, mode, maxIter, maxLoa)
|
||||
|
||||
def geod2ell(self, phi: float, lamb: float, h: float) -> Tuple[float, float]:
|
||||
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.
|
||||
:param point: kartesische 3D-Koordinaten
|
||||
:return: Punkt auf dem Ellispoid?
|
||||
"""
|
||||
value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2
|
||||
if abs(1-value) < 1e-6:
|
||||
return True
|
||||
else:
|
||||
return False
|
||||
|
||||
def cartonell(self, point: NDArray) -> NDArray:
|
||||
"""
|
||||
Berechnung des Lotpunktes auf einem Ellipsoiden
|
||||
:param point: Punkt in kartesischen Koordinaten, der gelotet werden soll
|
||||
:return: Lotpunkt in kartesischen Koordinaten, geodätische Koordinaten des Punktes
|
||||
:return: Lotpunkt in kartesischen Koordinaten
|
||||
"""
|
||||
phi, lamb, h = self.cart2geod(point, "ligas3")
|
||||
x, y, z = self. geod2cart(phi, lamb, 0)
|
||||
return np.array([x, y, z]), phi, lamb, h
|
||||
p = self. geod2cart(phi, lamb, 0)
|
||||
return p
|
||||
|
||||
def cartellh(self, point: np.ndarray, h: float) -> np.ndarray:
|
||||
def cartellh(self, point: NDArray, h: float) -> NDArray:
|
||||
"""
|
||||
Punkt auf Ellipsoid hoch loten
|
||||
:param point: Punkt auf dem Ellipsoid
|
||||
@@ -434,144 +631,64 @@ class EllipsoidTriaxial:
|
||||
pointH = self. geod2cart(phi, lamb, h)
|
||||
return pointH
|
||||
|
||||
def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray:
|
||||
"""
|
||||
Panou, Korakitits 2020, 4
|
||||
:param u: Parameter u
|
||||
:param v: Parameter v
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
x = self.ax * np.cos(u) * np.cos(v)
|
||||
y = self.ay * np.cos(u) * np.sin(v)
|
||||
z = self.b * np.sin(u)
|
||||
z = np.broadcast_to(z, np.shape(x))
|
||||
return np.array([x, y, z])
|
||||
|
||||
def cart2para(self, point: np.ndarray) -> tuple[float, float]:
|
||||
"""
|
||||
Panou, Korakitits 2020, 4
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:return: parametrische Koordinaten
|
||||
"""
|
||||
x, y, z = point
|
||||
|
||||
u_check1 = z*np.sqrt(1 - self.ee**2)
|
||||
u_check2 = np.sqrt(x**2 * (1-self.ee**2) + y**2) * np.sqrt(1-self.ex**2)
|
||||
if u_check1 <= u_check2:
|
||||
u = np.arctan2(u_check1, u_check2)
|
||||
else:
|
||||
u = np.pi/2 - np.arctan2(u_check2, u_check1)
|
||||
|
||||
v_check1 = y
|
||||
v_check2 = x*np.sqrt(1-self.ee**2)
|
||||
v_factor = np.sqrt(x**2*(1-self.ee**2)+y**2)
|
||||
if v_check1 <= v_check2:
|
||||
v = 2 * np.arctan2(v_check1, v_check2 + v_factor)
|
||||
else:
|
||||
v = np.pi/2 - 2 * np.arctan2(v_check2, v_check1 + v_factor)
|
||||
|
||||
return u, v
|
||||
|
||||
def ell2para(self, beta, lamb) -> tuple[float, float]:
|
||||
cart = self.ell2cart(beta, lamb)
|
||||
return self.cart2para(cart)
|
||||
|
||||
def para2ell(self, u, v) -> tuple[float, float]:
|
||||
cart = self.para2cart(u, v)
|
||||
return self.cart2ell(cart)
|
||||
|
||||
def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
|
||||
cart = self.para2cart(u, v)
|
||||
return self.cart2geod(cart, mode, maxIter, maxLoa)
|
||||
|
||||
def geod2para(self, phi, lamb, h) -> tuple[float, float]:
|
||||
cart = self.geod2cart(phi, lamb, h)
|
||||
return self.cart2para(cart)
|
||||
|
||||
def ell2geod(self, beta, lamb, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
|
||||
cart = self.ell2cart(beta, lamb)
|
||||
return self.cart2geod(cart, mode, maxIter, maxLoa)
|
||||
|
||||
def func_H(self, x, y, z):
|
||||
return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2
|
||||
|
||||
def func_n(self, x, y, z, H=None):
|
||||
if H is None:
|
||||
H = self.func_H(x, y, z)
|
||||
return np.array([x / sqrt(H),
|
||||
y / ((1 - self.ee ** 2) * sqrt(H)),
|
||||
z / ((1 - self.ex ** 2) * sqrt(H))])
|
||||
|
||||
def p_q(self, x, y, z) -> tuple[np.ndarray, np.ndarray]:
|
||||
"""
|
||||
Berechnung sämtlicher Größen
|
||||
:param x: x
|
||||
:param y: y
|
||||
:param z: z
|
||||
:return: Dictionary sämtlicher Größen
|
||||
"""
|
||||
n = self.func_n(x, y, z)
|
||||
|
||||
beta, lamb = self.cart2ell(np.array([x, y, z]))
|
||||
B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2
|
||||
L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2
|
||||
|
||||
c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
|
||||
c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
|
||||
(self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
|
||||
self.ax ** 2 + self.ay ** 2) * z ** 2)
|
||||
t2 = (-c1 + np.sqrt(c1 ** 2 - 4 * c0)) / 2
|
||||
t1 = c0 / t2
|
||||
t2e = self.ax ** 2 * np.sin(lamb) ** 2 + self.ay ** 2 * np.cos(lamb) ** 2
|
||||
t1e = self.ay ** 2 * np.sin(beta) ** 2 + self.b ** 2 * np.cos(beta) ** 2
|
||||
|
||||
F = self.Ey ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(lamb) ** 2
|
||||
p1 = -np.sqrt(L / (F * t2)) * self.ax / self.Ex * np.sqrt(B) * np.sin(lamb)
|
||||
p2 = np.sqrt(L / (F * t2)) * self.ay * np.cos(beta) * np.cos(lamb)
|
||||
p3 = 1 / np.sqrt(F * t2) * (self.b * self.Ee ** 2) / (2 * self.Ex) * np.sin(beta) * np.sin(2 * lamb)
|
||||
# p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * self.ax / (self.Ex * self.Ee) * np.sqrt(B) * np.sqrt(t2 - self.ay ** 2)
|
||||
# p2 = np.sign(x) * np.sqrt(L / (F * t2)) * self.ay / (self.Ey * self.Ee) * np.sqrt((ell.ay ** 2 - t1) * (self.ax ** 2 - t2))
|
||||
# p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * self.b / (self.Ex * self.Ey) * np.sqrt(
|
||||
# (t1 - self.b ** 2) * (t2 - self.ay ** 2) * (self.ax ** 2 - t2))
|
||||
p = np.array([p1, p2, p3])
|
||||
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]])
|
||||
|
||||
return p, q
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("Eitschberger1978")
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
diff_list = []
|
||||
for beta_deg in range(-150, 210, 30):
|
||||
for lamb_deg in range(-150, 210, 30):
|
||||
beta = wu.deg2rad(beta_deg)
|
||||
lamb = wu.deg2rad(lamb_deg)
|
||||
point = ell.ell2cart(beta, lamb)
|
||||
diffs_para = []
|
||||
diffs_ell = []
|
||||
diffs_geod = []
|
||||
points = []
|
||||
for v_deg in range(-180, 181, 5):
|
||||
for u_deg in range(-90, 91, 5):
|
||||
v = wu.deg2rad(v_deg)
|
||||
u = wu.deg2rad(u_deg)
|
||||
point = ell.para2cart(u, v)
|
||||
points.append(point)
|
||||
|
||||
elli = ell.cart2ell(point)
|
||||
cart_elli = ell.ell2cart(elli[0], elli[1])
|
||||
diff_ell = np.sum(np.abs(point-cart_elli))
|
||||
diff_ell = np.linalg.norm(point - cart_elli, axis=-1)
|
||||
|
||||
para = ell.cart2para(point)
|
||||
cart_para = ell.para2cart(para[0], para[1])
|
||||
diff_para = np.sum(np.abs(point-cart_para))
|
||||
|
||||
# geod = ell.cart2geod(point, "ligas1")
|
||||
# cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
||||
# diff_geod1 = np.sum(np.abs(point-cart_geod))
|
||||
#
|
||||
# geod = ell.cart2geod(point, "ligas2")
|
||||
# cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
||||
# diff_geod2 = np.sum(np.abs(point-cart_geod))
|
||||
diff_para = np.linalg.norm(point - cart_para, axis=-1)
|
||||
|
||||
geod = ell.cart2geod(point, "ligas3")
|
||||
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
||||
diff_geod3 = np.sum(np.abs(point-cart_geod))
|
||||
diff_geod3 = np.linalg.norm(point - cart_geod, axis=-1)
|
||||
|
||||
diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod3])
|
||||
diff_list.append([diff_ell])
|
||||
diff_list.append([v_deg, u_deg, diff_ell, diff_para, diff_geod3])
|
||||
diffs_ell.append([diff_ell])
|
||||
diffs_para.append([diff_para])
|
||||
diffs_geod.append([diff_geod3])
|
||||
|
||||
diff_list = np.array(diff_list)
|
||||
pass
|
||||
diffs_ell = np.array(diffs_ell)
|
||||
diffs_para = np.array(diffs_para)
|
||||
diffs_geod = np.array(diffs_geod)
|
||||
|
||||
pass
|
||||
points = np.array(points)
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(projection='3d')
|
||||
|
||||
sc = ax.scatter(
|
||||
points[:, 0],
|
||||
points[:, 1],
|
||||
points[:, 2],
|
||||
c=diffs_ell, # Farbcode = diff
|
||||
cmap='viridis', # Colormap
|
||||
s=10 + 20 * diffs_ell, # optional: Größe abhängig vom diff
|
||||
alpha=0.8
|
||||
)
|
||||
|
||||
# Farbskala
|
||||
cbar = plt.colorbar(sc)
|
||||
cbar.set_label("diff")
|
||||
|
||||
ax.set_xlabel("X")
|
||||
ax.set_ylabel("Y")
|
||||
ax.set_zlabel("Z")
|
||||
|
||||
plt.show()
|
||||
@@ -1,55 +0,0 @@
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
|
||||
def polyapp_tscheby_hayford(s: float) -> float:
|
||||
"""
|
||||
Berechnung der ellipsoidisch geodätischen Breite.
|
||||
Polynomapproximation mittels Tschebyscheff-Polynomen.
|
||||
Auf dem Hayford-Ellipsoid.
|
||||
:param s: Strecke auf einer Ellipse vom Äquator aus
|
||||
:type s: float
|
||||
:return: ellipsoidisch geodätische Breite
|
||||
:rtype: float
|
||||
"""
|
||||
c0 = 1
|
||||
c1 = -0.00837809325
|
||||
c2 = 0.00428127367
|
||||
c3 = -0.00114523986
|
||||
c4 = 0.00023219707
|
||||
c5 = -0.00004421222
|
||||
c6 = 0.00000570244
|
||||
alpha = wu.gms2rad([0, 0, 325643.97199])
|
||||
s90 = 10002288.2990
|
||||
|
||||
xi = s/s90
|
||||
|
||||
phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) +
|
||||
c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6))
|
||||
return phi
|
||||
|
||||
|
||||
def polyapp_tscheby_bessel(s: float) -> float:
|
||||
"""
|
||||
Berechnung der ellipsoidisch geodätischen Breite.
|
||||
Polynomapproximation mittels Tschebyscheff-Polynomen.
|
||||
Auf dem Bessel-Ellipsoid.
|
||||
:param s: Strecke auf einer Ellipse vom Äquator aus
|
||||
:type s: float
|
||||
:return: ellipsoidisch geodätische Breite
|
||||
:rtype: float
|
||||
"""
|
||||
c0 = 1
|
||||
c1 = -0.00831729565
|
||||
c2 = 0.00424914906
|
||||
c3 = -0.00113566119
|
||||
c4 = 0.00022976983
|
||||
c5 = -0.00004363980
|
||||
c6 = 0.00000562025
|
||||
alpha = wu.gms2rad([0, 0, 325632.08677])
|
||||
s90 = 10000855.7644
|
||||
|
||||
xi = s/s90
|
||||
|
||||
phi = alpha * xi * (c0*xi**(2*0) + c1*xi**(2*1) + c2*xi**(2*2) + c3*xi**(2*3) +
|
||||
c4*xi**(2*4) + c5*xi**(2*5) + c6*xi**(2*6))
|
||||
return phi
|
||||
43
runge_kutta.py
Normal file
43
runge_kutta.py
Normal file
@@ -0,0 +1,43 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]:
|
||||
"""
|
||||
Standard Runge-Kutta Verfahren 4. Ordnung
|
||||
:param ode: ODE-System als Funktion
|
||||
:param t0: Startwert der unabhängigen Variable
|
||||
:param v0: Startwerte
|
||||
:param weite: Integrationsweite
|
||||
:param schritte: Schrittzahl
|
||||
:param fein:
|
||||
:return: Variable und Funktionswerte an jedem Stützpunkt
|
||||
"""
|
||||
h = weite/schritte
|
||||
|
||||
t_list = [t0]
|
||||
werte = [v0]
|
||||
|
||||
for _ in range(schritte):
|
||||
t = t_list[-1]
|
||||
v = werte[-1]
|
||||
|
||||
if not fein:
|
||||
v_next = rk4_step(ode, t, v, h)
|
||||
|
||||
else:
|
||||
v_grob = rk4_step(ode, t, v, h)
|
||||
v_half = rk4_step(ode, t, v, 0.5 * h)
|
||||
v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
|
||||
v_next = v_fein + (v_fein - v_grob) / 15.0
|
||||
|
||||
t_list.append(t + h)
|
||||
werte.append(v_next)
|
||||
|
||||
return t_list, werte
|
||||
|
||||
def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray:
|
||||
k1 = ode(t, v)
|
||||
k2 = ode(t + 0.5 * h, v + 0.5 * h * k1)
|
||||
k3 = ode(t + 0.5 * h, v + 0.5 * h * k2)
|
||||
k4 = ode(t + h, v + h * k3)
|
||||
return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
|
||||
78
s_ellipse.py
78
s_ellipse.py
@@ -1,78 +0,0 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def reihenentwicklung(es: float, c: float, phi: float) -> float:
|
||||
"""
|
||||
Berechnung der Strecke auf einer Ellipse.
|
||||
Reihenentwicklung.
|
||||
:param es: zweite numerische Exzentrizität
|
||||
:type es: float
|
||||
:param c: Polkrümmungshalbmesser
|
||||
:type c: float
|
||||
:param phi: ellipsoidisch geodästische Breite in Radiant
|
||||
:type phi: float
|
||||
:return: Strecke auf einer Ellipse vom Äquator aus
|
||||
:rtype: float
|
||||
"""
|
||||
Ass = 1 - 3/4*es**2 + 45/64*es**4 - 175/256*es**6 + 11025/16384*es**8
|
||||
Bss = - 3/4*es**2 + 15/16*es**4 - 525/512*es**6 + 2205/2048*es**8
|
||||
Css = 15/64*es**4 - 105/256*es**6 + 2205/4096*es**8
|
||||
Dss = - 35/512*es**6 + 315/2048*es**8
|
||||
print(f"A'' = {round(Ass, 10):.10f}\nB'' = {round(Bss, 10):.10f}\nC'' = {round(Css, 10):.10f}\nD'' = {round(Dss, 10):.10f}")
|
||||
|
||||
s = c * (Ass*phi + 1/2*Bss * np.sin(2*phi) + 1/4*Css * np.sin(4*phi) + 1/6*Dss * np.sin(6*phi))
|
||||
return s
|
||||
|
||||
|
||||
def polyapp_tscheby_hayford(phi: float) -> float:
|
||||
"""
|
||||
Berechnung der Strecke auf einer Ellipse.
|
||||
Polynomapproximation mittels Tschebyscheff-Polynomen.
|
||||
Auf dem Hayford-Ellipsoid.
|
||||
:param phi: ellipsoidisch geodästische Breite in Radiant
|
||||
:type phi: float
|
||||
:return: Strecke auf einer Ellipse vom Äquator aus
|
||||
:rtype: float
|
||||
"""
|
||||
c1s = 0.00829376218
|
||||
c2s = -0.00398963425
|
||||
c3s = 0.00084200710
|
||||
c4s = -0.0000648906
|
||||
c5s = -0.00001075680
|
||||
c6s = 0.00000396474
|
||||
c7s = -0.00000046347
|
||||
alpha = 9951793.0123
|
||||
|
||||
xi = 2 / np.pi * phi
|
||||
xis = xi ^ 2
|
||||
|
||||
s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3
|
||||
+ c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7)
|
||||
return s
|
||||
|
||||
|
||||
def polyapp_tscheby_bessel(phi: float) -> float:
|
||||
"""
|
||||
Berechnung der Strecke auf einer Ellipse.
|
||||
Polynomapproximation mittels Tschebyscheff-Polynomen.
|
||||
Auf dem Bessel-Ellipsoid.
|
||||
:param phi: ellipsoidisch geodästische Breite in Radiant
|
||||
:type phi: float
|
||||
:return: Strecke auf einer Ellipse vom Äquator aus
|
||||
:rtype: float
|
||||
"""
|
||||
c1s = 0.00823417717
|
||||
c2s = -0.00396170744
|
||||
c3s = 0.00083680249
|
||||
c4s = -0.00006488462
|
||||
c5s = -0.00001053242
|
||||
c6s = 0.00000390854
|
||||
c7s = -0.00000045768
|
||||
alpha = 9950730.8876
|
||||
|
||||
xi = 2 / np.pi * phi
|
||||
xis = xi ** 2
|
||||
|
||||
s = alpha * xi * (1 + c1s * xis + c2s * xis**2 + c3s * xis**3
|
||||
+ c4s * xis**4 + c5s * xis**5 + c6s * xis**6 + c7s * xis**7)
|
||||
return s
|
||||
@@ -1,30 +0,0 @@
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import winkelumrechnungen as wu
|
||||
from dashboard import ellipsoid_figure
|
||||
|
||||
u = np.linspace(0, 2*np.pi, 51)
|
||||
v = np.linspace(0, np.pi, 51)
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
points = []
|
||||
lines = []
|
||||
for u_i, u_value in enumerate(u):
|
||||
for v_i, v_value in enumerate(v):
|
||||
cart = ell.ell2cart(u_value, v_value)
|
||||
if u_i != 0 and v_i != 0:
|
||||
lines.append((points[-1], cart, "red"))
|
||||
points.append(cart)
|
||||
points = []
|
||||
for v_i, v_value in enumerate(v):
|
||||
for u_i, u_value in enumerate(u):
|
||||
cart = ell.ell2cart(u_value, v_value)
|
||||
if u_i != 0 and v_i != 0:
|
||||
lines.append((points[-1], cart, "blue"))
|
||||
points.append(cart)
|
||||
ax = ell.ax
|
||||
ay = ell.ay
|
||||
b = ell.b
|
||||
|
||||
figu = ellipsoid_figure(ax, ay, b, lines=lines)
|
||||
figu.show()
|
||||
31
test.py
31
test.py
@@ -1,29 +1,14 @@
|
||||
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
|
||||
|
||||
J = np.array([
|
||||
[2, 3, 0],
|
||||
[0, 3, 0],
|
||||
[6, 0, 4]
|
||||
])
|
||||
ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978")
|
||||
|
||||
xi = np.array([1, 2, 3])
|
||||
xi_col = xi.reshape(-1, 1)
|
||||
print(xi_col)
|
||||
xi_row = xi_col.reshape(1, -1).flatten()
|
||||
print(xi_row)
|
||||
x, y, z = ell.para2cart(0.5, 0.7)
|
||||
|
||||
# Spaltenvektor-Variante
|
||||
res_col = xi[:, None] - J @ xi[:, None]
|
||||
|
||||
# Zeilenvektor-Variante
|
||||
res_row = xi[None, :] - xi[None, :] @ J
|
||||
|
||||
print("Spaltenvektor:")
|
||||
print(res_col[0,0])
|
||||
print("Zeilenvektor:")
|
||||
print(res_row)
|
||||
t = 5
|
||||
l = 2
|
||||
print(fact(t+1-l) / (fact(t+1-l) * fact(l-1)), comb(l-1, t+1-l))
|
||||
x1 = pq(ell, x, y, z)
|
||||
x2 = ell.p_q(x, y, z)
|
||||
pass
|
||||
@@ -1,84 +0,0 @@
|
||||
import GHA_triaxial.numeric_examples_panou as nep
|
||||
import ellipsoide
|
||||
from GHA_triaxial.panou_2013_2GHA_num import gha2_num
|
||||
from GHA_triaxial.panou import gha1_ana, gha1_num
|
||||
import numpy as np
|
||||
import time
|
||||
|
||||
def test():
|
||||
ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
|
||||
tables = nep.get_tables()
|
||||
|
||||
diffs_gha1_num = []
|
||||
diffs_gha1_ana = []
|
||||
diffs_gha2_num = []
|
||||
times_gha1_num = []
|
||||
times_gha1_ana = []
|
||||
times_gha2_num = []
|
||||
|
||||
for table in tables:
|
||||
diffs_gha1_num.append([])
|
||||
diffs_gha1_ana.append([])
|
||||
diffs_gha2_num.append([])
|
||||
times_gha1_num.append([])
|
||||
times_gha1_ana.append([])
|
||||
times_gha2_num.append([])
|
||||
|
||||
for example in table:
|
||||
beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = example
|
||||
P0 = ell.ell2cart(beta0, lamb0)
|
||||
P1 = ell.ell2cart(beta1, lamb1)
|
||||
|
||||
start = time.perf_counter()
|
||||
try:
|
||||
P1_num = gha1_num(ell, P0, alpha0, s, 10000)
|
||||
end = time.perf_counter()
|
||||
diff_P1_num = np.linalg.norm(P1 - P1_num)
|
||||
except:
|
||||
end = time.perf_counter()
|
||||
diff_P1_num = None
|
||||
time_gha1_num = end - start
|
||||
|
||||
start = time.perf_counter()
|
||||
try:
|
||||
P1_ana = gha1_ana(ell, P0, alpha0, s, 50)
|
||||
end = time.perf_counter()
|
||||
diff_P1_ana = np.linalg.norm(P1 - P1_ana)
|
||||
except:
|
||||
end = time.perf_counter()
|
||||
diff_P1_ana = None
|
||||
time_gha1_ana = end - start
|
||||
|
||||
start = time.perf_counter()
|
||||
try:
|
||||
alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000)
|
||||
end = time.perf_counter()
|
||||
diff_s_num = abs(s - s_num)
|
||||
except:
|
||||
end = time.perf_counter()
|
||||
diff_s_num = None
|
||||
time_gha2_num = None
|
||||
time_gha2_num = end - start
|
||||
|
||||
diffs_gha1_num[-1].append(diff_P1_num)
|
||||
diffs_gha1_ana[-1].append(diff_P1_ana)
|
||||
diffs_gha2_num[-1].append(diff_s_num)
|
||||
times_gha1_num[-1].append(time_gha1_num)
|
||||
times_gha1_ana[-1].append(time_gha1_ana)
|
||||
times_gha2_num[-1].append(time_gha2_num)
|
||||
print(diffs_gha1_num, diffs_gha1_ana, diffs_gha2_num)
|
||||
print(times_gha1_num, times_gha1_ana, times_gha2_num)
|
||||
|
||||
def display():
|
||||
diffs = [[{'gha1_num': np.float64(3.410763124264611e-05), 'gha1_ana': np.float64(3.393273802112796e-05), 'gha2_num': np.float64(3.3931806683540344e-05)}, {'gha1_num': np.float64(0.0008736425000530604), 'gha1_ana': np.float64(0.0008736458415010259), 'gha2_num': None}, {'gha1_num': np.float64(0.0007739730058338136), 'gha1_ana': np.float64(0.0007739621469802854), 'gha2_num': np.float64(1.5832483768463135e-07)}, {'gha1_num': np.float64(0.00010554956741100295), 'gha1_ana': np.float64(8.814246009944831), 'gha2_num': np.float64(4.864111542701721e-05)}, {'gha1_num': np.float64(0.0002135908394614854), 'gha1_ana': np.float64(0.0002138610897967267), 'gha2_num': np.float64(5.0179407158866525)}, {'gha1_num': np.float64(0.00032727226891456654), 'gha1_ana': np.float64(0.00032734569198545905), 'gha2_num': np.float64(9.735533967614174e-05)}, {'gha1_num': np.float64(0.0005195973303787956), 'gha1_ana': np.float64(0.0005197766935509641), 'gha2_num': None}], [{'gha1_num': np.float64(1.780250537652368e-05), 'gha1_ana': np.float64(1.996805145339501e-05), 'gha2_num': np.float64(1.8164515495300293e-05)}, {'gha1_num': np.float64(4.8607540473363564e-05), 'gha1_ana': np.float64(2205539.954949392), 'gha2_num': None}, {'gha1_num': np.float64(0.00017376854985685854), 'gha1_ana': np.float64(328124.1513636429), 'gha2_num': np.float64(0.17443156614899635)}, {'gha1_num': np.float64(5.83429352558999e-05), 'gha1_ana': np.float64(0.01891628037258558), 'gha2_num': np.float64(1.4207654744386673)}, {'gha1_num': np.float64(0.0006421087024666934), 'gha1_ana': np.float64(0.0006420400127297228), 'gha2_num': np.float64(0.12751091085374355)}, {'gha1_num': np.float64(0.0004456207867164434), 'gha1_ana': np.float64(0.0004455649707698245), 'gha2_num': np.float64(0.00922046648338437)}, {'gha1_num': np.float64(0.0002340879908275419), 'gha1_ana': np.float64(0.00023422217242111216), 'gha2_num': np.float64(0.001307751052081585)}], [{'gha1_num': np.float64(976.6580096633622), 'gha1_ana': np.float64(976.6580096562798), 'gha2_num': np.float64(6.96033239364624e-05)}, {'gha1_num': np.float64(2825.2936643258527), 'gha1_ana': np.float64(2794.954866417055), 'gha2_num': np.float64(1.3615936040878296e-05)}, {'gha1_num': np.float64(1248.8942058074501), 'gha1_ana': np.float64(538.5550561841195), 'gha2_num': np.float64(3.722589462995529e-05)}, {'gha1_num': np.float64(2201.1793359793814), 'gha1_ana': np.float64(3735.376499414938), 'gha2_num': np.float64(1.4525838196277618e-05)}, {'gha1_num': np.float64(2262.134819997246), 'gha1_ana': np.float64(25549.567793410763), 'gha2_num': np.float64(9.328126907348633e-06)}, {'gha1_num': np.float64(2673.219788119847), 'gha1_ana': np.float64(21760.866677295206), 'gha2_num': np.float64(8.635222911834717e-06)}, {'gha1_num': np.float64(1708.758419275875), 'gha1_ana': np.float64(3792.1128807063437), 'gha2_num': np.float64(2.4085864424705505e-05)}], [{'gha1_num': np.float64(0.7854659044152204), 'gha1_ana': np.float64(0.785466068424286), 'gha2_num': np.float64(0.785466069355607)}, {'gha1_num': np.float64(237.79878717216718), 'gha1_ana': np.float64(1905080.064324282), 'gha2_num': None}, {'gha1_num': np.float64(55204.601699830164), 'gha1_ana': np.float64(55204.60175211949), 'gha2_num': None}, {'gha1_num': np.float64(12766.348063015519), 'gha1_ana': np.float64(12766.376619517901), 'gha2_num': np.float64(12582.786206113175)}, {'gha1_num': np.float64(29703.049988324146), 'gha1_ana': np.float64(29703.056427749252), 'gha2_num': np.float64(28933.668131249025)}, {'gha1_num': np.float64(43912.03007182513), 'gha1_ana': np.float64(43912.03007528712), 'gha2_num': None}, {'gha1_num': np.float64(28522.29828970693), 'gha1_ana': np.float64(28522.29830145182), 'gha2_num': None}, {'gha1_num': np.float64(17769.115549537233), 'gha1_ana': np.float64(17769.115549483362), 'gha2_num': np.float64(17769.121286311187)}]]
|
||||
arr = []
|
||||
for table in diffs:
|
||||
for example in table:
|
||||
arr.append([example['gha1_num'], example['gha1_ana'], example['gha2_num']])
|
||||
arr = np.array(arr)
|
||||
pass
|
||||
|
||||
if __name__ == "__main__":
|
||||
# test()
|
||||
display()
|
||||
Reference in New Issue
Block a user