Compare commits
9 Commits
db05f7b6db
...
main
| Author | SHA1 | Date | |
|---|---|---|---|
| d90ff5df69 | |||
| 59ad560f36 | |||
| 5a293a823a | |||
| 36b62059fc | |||
| 57a086f6cb | |||
| b8d07307aa | |||
| 798cace25d | |||
| 4f7b9aaef0 | |||
| 1fbfb555a4 |
@@ -1,7 +1,8 @@
|
||||
import numpy as np
|
||||
from numpy.typing import NDArray
|
||||
|
||||
|
||||
def felli(x):
|
||||
def felli(x: NDArray) -> float:
|
||||
N = x.shape[0]
|
||||
if N < 2:
|
||||
raise ValueError("dimension must be greater than one")
|
||||
@@ -11,8 +12,7 @@ def felli(x):
|
||||
|
||||
def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000,
|
||||
func_args=(), func_kwargs=None, seed=0,
|
||||
bestEver = np.inf, noImproveGen = 0, absTolImprove = 1e-12, maxNoImproveGen = 100, sigmaImprove = 1e-12):
|
||||
|
||||
bestEver=np.inf, noImproveGen=0, absTolImprove=1e-12, maxNoImproveGen=100, sigmaImprove=1e-12):
|
||||
if func_kwargs is None:
|
||||
func_kwargs = {}
|
||||
|
||||
@@ -27,7 +27,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
N = xmean.shape[0]
|
||||
|
||||
if stopeval is None:
|
||||
stopeval = int(1e3 * N**2)
|
||||
stopeval = int(1e3 * N ** 2)
|
||||
|
||||
# Strategy parameter setting: Selection
|
||||
lambda_ = 4 + int(np.floor(3 * np.log(N)))
|
||||
@@ -37,14 +37,14 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
weights = np.log(mu + 0.5) - np.log(np.arange(1, int(mu) + 1))
|
||||
mu = int(np.floor(mu))
|
||||
weights = weights / np.sum(weights)
|
||||
mueff = np.sum(weights)**2 / np.sum(weights**2)
|
||||
mueff = np.sum(weights) ** 2 / np.sum(weights ** 2)
|
||||
|
||||
# Strategy parameter setting: Adaptation
|
||||
cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N)
|
||||
cs = (mueff + 2) / (N + mueff + 5)
|
||||
c1 = 2 / ((N + 1.3)**2 + mueff)
|
||||
c1 = 2 / ((N + 1.3) ** 2 + mueff)
|
||||
cmu = min(1 - c1,
|
||||
2 * (mueff - 2 + 1 / mueff) / ((N + 2)**2 + 2 * mueff))
|
||||
2 * (mueff - 2 + 1 / mueff) / ((N + 2) ** 2 + 2 * mueff))
|
||||
damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs
|
||||
|
||||
# Initialize dynamic (internal) strategy parameters and constants
|
||||
@@ -54,7 +54,7 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
D = np.eye(N)
|
||||
C = B @ D @ (B @ D).T
|
||||
eigeneval = 0
|
||||
chiN = np.sqrt(N) * (1 - 1/(4*N) + 1/(21 * N**2))
|
||||
chiN = np.sqrt(N) * (1 - 1 / (4 * N) + 1 / (21 * N ** 2))
|
||||
|
||||
# Generation Loop
|
||||
counteval = 0
|
||||
@@ -91,10 +91,9 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
bestEver = fbest
|
||||
noImproveGen = 0
|
||||
else:
|
||||
noImproveGen = noImproveGen + 1
|
||||
noImproveGen += 1
|
||||
|
||||
|
||||
if gen == 1 or gen%50==0:
|
||||
if gen == 1 or gen % 50 == 0:
|
||||
# print(f' [CMA-ES] Gen {gen}, best = {round(fbest, 6)}, sigma = {sigma:.3g}')
|
||||
pass
|
||||
|
||||
@@ -106,13 +105,10 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
# print(f' [CMA-ES] Abbruch: sigma zu klein {sigma:.3g}')
|
||||
break
|
||||
|
||||
|
||||
|
||||
# Cumulation: Update evolution paths
|
||||
ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * (B @ zmean)
|
||||
norm_ps = np.linalg.norm(ps)
|
||||
hsig = norm_ps / np.sqrt(1 - (1 - cs)**(2 * counteval / lambda_)) / chiN < \
|
||||
(1.4 + 2 / (N + 1))
|
||||
hsig = norm_ps / np.sqrt(1 - (1 - cs) ** (2 * counteval / lambda_)) / chiN < (1.4 + 2 / (N + 1))
|
||||
hsig = 1.0 if hsig else 0.0
|
||||
|
||||
pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (B @ D @ zmean)
|
||||
@@ -142,13 +138,13 @@ def escma(func, *, N=10, xmean=None, sigma=0.5, stopfitness=1e-14, stopeval=2000
|
||||
if arfitness[0] == arfitness[int(np.ceil(0.7 * lambda_)) - 1]:
|
||||
sigma = sigma * np.exp(0.2 + cs / damps)
|
||||
# print(' [CMA-ES] stopfitness erreicht.')
|
||||
#print("warning: flat fitness, consider reformulating the objective")
|
||||
# print("warning: flat fitness, consider reformulating the objective")
|
||||
break
|
||||
|
||||
#print(f"{counteval}: {arfitness[0]}")
|
||||
# print(f"{counteval}: {arfitness[0]}")
|
||||
|
||||
#Final Message
|
||||
#print(f"{counteval}: {arfitness[0]}")
|
||||
# Final Message
|
||||
# print(f"{counteval}: {arfitness[0]}")
|
||||
xmin = arx[:, arindex[0]]
|
||||
bestValue = arfitness[0]
|
||||
# print(f' [CMA-ES] Ende: Gen = {gen}, best = {round(bestValue, 6)}')
|
||||
@@ -1,29 +1,17 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from codeop import PyCF_ALLOW_INCOMPLETE_INPUT
|
||||
from typing import List, Optional, Tuple
|
||||
from typing import List, Tuple
|
||||
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
from ES.Hansen_ES_CMA import escma
|
||||
from GHA_triaxial.gha1_ana import gha1_ana
|
||||
from GHA_triaxial.gha1_approx import gha1_approx
|
||||
from Hansen_ES_CMA import escma
|
||||
from utils_angle import wrap_to_pi
|
||||
from numpy.typing import NDArray
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
|
||||
def ellipsoid_formparameter(ell: EllipsoidTriaxial):
|
||||
"""
|
||||
Berechnet die Formparameter des dreiachsigen Ellipsoiden nach Karney (2025), Gl. (2)
|
||||
:param ell: Ellipsoid
|
||||
:return: e, k und k'
|
||||
"""
|
||||
nenner = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0))
|
||||
k = np.sqrt(max(ell.ay * ell.ay - ell.b * ell.b, 0.0)) / nenner
|
||||
k_ = np.sqrt(max(ell.ax * ell.ax - ell.ay * ell.ay, 0.0)) / nenner
|
||||
e = np.sqrt(max(ell.ax * ell.ax - ell.b * ell.b, 0.0)) / ell.ay
|
||||
|
||||
return e, k, k_
|
||||
from GHA_triaxial.utils import jacobi_konstante
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from utils_angle import wrap_mpi_pi
|
||||
|
||||
|
||||
def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
|
||||
@@ -40,7 +28,7 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
|
||||
R (XYZ) = Punkt in XYZ
|
||||
"""
|
||||
# Berechnungshilfen
|
||||
omega = wrap_to_pi(omega)
|
||||
omega = wrap_mpi_pi(omega)
|
||||
cb = np.cos(beta)
|
||||
sb = np.sin(beta)
|
||||
co = np.cos(omega)
|
||||
@@ -61,47 +49,32 @@ def ENU_beta_omega(beta: float, omega: float, ell: EllipsoidTriaxial) \
|
||||
# --- Ableitungen - Karney Gl. (5a,b,c)---
|
||||
# E = ∂R/∂ω
|
||||
dX_dw = -ell.ax * so * Sx / D
|
||||
dY_dw = ell.ay * cb * co
|
||||
dY_dw = ell.ay * cb * co
|
||||
dZ_dw = ell.b * sb * (so * co * (ell.ax*ell.ax - ell.ay*ell.ay) / Sz) / D
|
||||
E = np.array([dX_dw, dY_dw, dZ_dw], dtype=float)
|
||||
|
||||
# N = ∂R/∂β
|
||||
dX_db = ell.ax * co * (sb * cb * (ell.b*ell.b - ell.ay*ell.ay) / Sx) / D
|
||||
dY_db = -ell.ay * sb * so
|
||||
dZ_db = ell.b * cb * Sz / D
|
||||
dZ_db = ell.b * cb * Sz / D
|
||||
N = np.array([dX_db, dY_db, dZ_db], dtype=float)
|
||||
|
||||
# U = Grad(x^2/a^2 + y^2/b^2 + z^2/c^2 - 1)
|
||||
U = np.array([X/(ell.ax*ell.ax), Y/(ell.ay*ell.ay), Z/(ell.b*ell.b)], dtype=float)
|
||||
|
||||
En = np.linalg.norm(E)
|
||||
Nn = np.linalg.norm(N)
|
||||
Un = np.linalg.norm(U)
|
||||
En = float(np.linalg.norm(E))
|
||||
Nn = float(np.linalg.norm(N))
|
||||
Un = float(np.linalg.norm(U))
|
||||
|
||||
N_hat = N / Nn
|
||||
E_hat = E / En
|
||||
U_hat = U / Un
|
||||
E_hat = E_hat - float(np.dot(E_hat, N_hat)) * N_hat
|
||||
E_hat -= float(np.dot(E_hat, N_hat)) * N_hat
|
||||
E_hat = E_hat / max(np.linalg.norm(E_hat), 1e-18)
|
||||
|
||||
return E_hat, N_hat, U_hat, En, Nn, R
|
||||
|
||||
|
||||
def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTriaxial) -> float:
|
||||
"""
|
||||
Jacobi-Konstante nach Karney (2025), Gl. (14)
|
||||
:param beta: Beta Koordinate
|
||||
:param omega: Omega Koordinate
|
||||
:param alpha: Azimut alpha
|
||||
:param ell: Ellipsoid
|
||||
:return: Jacobi-Konstante
|
||||
"""
|
||||
e, k, k_ = ellipsoid_formparameter(ell)
|
||||
gamma_jacobi = float((k ** 2) * (np.cos(beta) ** 2) * (np.sin(alpha) ** 2) - (k_ ** 2) * (np.sin(omega) ** 2) * (np.cos(alpha) ** 2))
|
||||
|
||||
return gamma_jacobi
|
||||
|
||||
|
||||
def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_hat_curr: NDArray, U_hat_curr: NDArray) -> float:
|
||||
"""
|
||||
Berechnet das Azimut in der lokalen Tangentialebene am aktuellen Punkt P_curr, gemessen
|
||||
@@ -117,11 +90,10 @@ def azimuth_at_ESpoint(P_prev: NDArray, P_curr: NDArray, E_hat_curr: NDArray, N_
|
||||
vT = v - float(np.dot(v, U_hat_curr)) * U_hat_curr
|
||||
vTn = max(np.linalg.norm(vT), 1e-18)
|
||||
vT_hat = vT / vTn
|
||||
#vT_hat = vT / np.linalg.norm(vT)
|
||||
sE = float(np.dot(vT_hat, E_hat_curr))
|
||||
sN = float(np.dot(vT_hat, N_hat_curr))
|
||||
|
||||
return wrap_to_pi(float(np.arctan2(sE, sN)))
|
||||
return wrap_mpi_pi(float(np.arctan2(sE, sN)))
|
||||
|
||||
|
||||
def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float, gamma0: float,
|
||||
@@ -154,11 +126,10 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
|
||||
d_beta = float(np.clip(d_beta, -0.2, 0.2)) # rad
|
||||
d_omega = float(np.clip(d_omega, -0.2, 0.2)) # rad
|
||||
|
||||
|
||||
#d_beta = ds * float(np.cos(alpha_i)) / Nn_i
|
||||
#d_omega = ds * float(np.sin(alpha_i)) / En_i
|
||||
# d_beta = ds * float(np.cos(alpha_i)) / Nn_i
|
||||
# d_omega = ds * float(np.sin(alpha_i)) / En_i
|
||||
beta_pred = beta_i + d_beta
|
||||
omega_pred = wrap_to_pi(omega_i + d_omega)
|
||||
omega_pred = wrap_mpi_pi(omega_i + d_omega)
|
||||
|
||||
xmean = np.array([beta_pred, omega_pred], dtype=float)
|
||||
|
||||
@@ -175,10 +146,10 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
|
||||
:return: Fitnesswert (f)
|
||||
"""
|
||||
beta = x[0]
|
||||
omega = wrap_to_pi(x[1])
|
||||
omega = wrap_mpi_pi(x[1])
|
||||
|
||||
P = ell.ell2cart_karney(beta, omega) # in kartesischer Koordinaten
|
||||
d = float(np.linalg.norm(P - P_i)) # Distanz zwischen
|
||||
P = ell.ell2cart_karney(beta, omega) # in kartesischer Koordinaten
|
||||
d = float(np.linalg.norm(P - P_i)) # Distanz zwischen
|
||||
|
||||
# maxSegLen einhalten
|
||||
J_len = ((d - ds) / ds) ** 2
|
||||
@@ -197,11 +168,10 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
|
||||
|
||||
return f
|
||||
|
||||
|
||||
xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES
|
||||
xb = escma(fitness, N=2, xmean=xmean, sigma=sigma0) # Aufruf CMA-ES
|
||||
|
||||
beta_best = xb[0]
|
||||
omega_best = wrap_to_pi(xb[1])
|
||||
omega_best = wrap_mpi_pi(xb[1])
|
||||
P_best = ell.ell2cart_karney(beta_best, omega_best)
|
||||
E_j, N_j, U_j, _, _, _ = ENU_beta_omega(beta_best, omega_best, ell)
|
||||
alpha_end = azimuth_at_ESpoint(P_i, P_best, E_j, N_j, U_j)
|
||||
@@ -209,7 +179,7 @@ def optimize_next_point(beta_i: float, omega_i: float, alpha_i: float, ds: float
|
||||
return beta_best, omega_best, P_best, alpha_end
|
||||
|
||||
|
||||
def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000, all_points: boolean = False)\
|
||||
def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float, s_total: float, maxSegLen: float = 1000, all_points: bool = False)\
|
||||
-> Tuple[NDArray, float, NDArray] | Tuple[NDArray, float]:
|
||||
"""
|
||||
Aufruf der 1. GHA mittels CMA-ES
|
||||
@@ -223,10 +193,10 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
|
||||
:return: Zielpunkt Pk, Azimut am Zielpunkt und Punktliste
|
||||
"""
|
||||
beta = float(beta0)
|
||||
omega = wrap_to_pi(float(omega0))
|
||||
alpha = wrap_to_pi(float(alpha0))
|
||||
omega = wrap_mpi_pi(float(omega0))
|
||||
alpha = wrap_mpi_pi(float(alpha0))
|
||||
|
||||
gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0
|
||||
gamma0 = jacobi_konstante(beta, omega, alpha, ell) # Referenz-γ0
|
||||
|
||||
P_all: List[NDArray] = [ell.ell2cart_karney(beta, omega)]
|
||||
alpha_end: List[float] = [alpha]
|
||||
@@ -243,7 +213,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
|
||||
ell=ell, maxSegLen=maxSegLen)
|
||||
s_acc += ds
|
||||
P_all.append(P)
|
||||
alpha_end.append(alpha)
|
||||
alpha_end.append(wrap_mpi_pi(alpha))
|
||||
if step > nsteps_est + 50:
|
||||
raise RuntimeError("GHA1_ES: Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
|
||||
Pk = P_all[-1]
|
||||
@@ -258,8 +228,8 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
s = 180000
|
||||
#alpha0 = 3
|
||||
alpha0 = wu.gms2rad([5 ,0 ,0])
|
||||
# alpha0 = 3
|
||||
alpha0 = wu.gms2rad([5, 0, 0])
|
||||
beta = 0
|
||||
omega = 0
|
||||
P0 = ell.ell2cart(beta, omega)
|
||||
@@ -272,6 +242,6 @@ if __name__ == "__main__":
|
||||
print(res)
|
||||
print(alpha)
|
||||
print(points)
|
||||
#print("alpha1 (am Endpunkt):", res.alpha1)
|
||||
# print("alpha1 (am Endpunkt):", res.alpha1)
|
||||
print(res - point1)
|
||||
print(point1app - point1, "approx")
|
||||
@@ -1,14 +1,13 @@
|
||||
import numpy as np
|
||||
from Hansen_ES_CMA import escma
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from numpy.typing import NDArray
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from numpy.typing import NDArray
|
||||
|
||||
from ES.Hansen_ES_CMA import escma
|
||||
from GHA_triaxial.gha2_num import gha2_num
|
||||
from GHA_triaxial.utils import sigma2alpha, pq_ell
|
||||
|
||||
|
||||
|
||||
from GHA_triaxial.utils import sigma2alpha
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
|
||||
|
||||
def Sehne(P1: NDArray, P2: NDArray) -> float:
|
||||
@@ -19,14 +18,11 @@ def Sehne(P1: NDArray, P2: NDArray) -> float:
|
||||
:return: Bogenlänge s
|
||||
"""
|
||||
R12 = P2-P1
|
||||
s = np.linalg.norm(R12)
|
||||
s = float(np.linalg.norm(R12))
|
||||
|
||||
return s
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float = None, all_points: bool = False) -> Tuple[float, float, float, NDArray] | Tuple[float, float, float]:
|
||||
"""
|
||||
Berechnen der 2. GHA mithilfe der CMA-ES.
|
||||
@@ -70,7 +66,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
|
||||
|
||||
R0 = (ell.ax + ell.ay + ell.b) / 3
|
||||
if maxSegLen is None:
|
||||
maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius
|
||||
maxSegLen = R0 * 1 / (637.4*2) # 10km Segment bei mittleren Erdradius
|
||||
|
||||
sigma_uv_nom = 1e-3 * (maxSegLen / R0) # ~1e-5
|
||||
points: list[NDArray] = [P0, Pk]
|
||||
@@ -102,7 +98,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
|
||||
|
||||
sigmaStep = sigma_uv_nom * (Sehne(A, B) / maxSegLen)
|
||||
|
||||
u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) # Aufruf CMA-ES
|
||||
u, v = escma(midpoint_fitness, N=2, xmean=xmean, sigma=sigmaStep) # Aufruf CMA-ES
|
||||
|
||||
P_next = ell.para2cart(u, v)
|
||||
new_points.append(P_next)
|
||||
@@ -3,12 +3,13 @@ from math import comb
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arctan2
|
||||
from numpy._typing import NDArray
|
||||
import winkelumrechnungen as wu
|
||||
from numpy import arctan2, cos, sin
|
||||
from numpy.typing import NDArray
|
||||
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import winkelumrechnungen as wu
|
||||
from GHA_triaxial.utils import pq_para
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from utils_angle import wrap_0_2pi
|
||||
|
||||
|
||||
def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int) -> Tuple[NDArray, float]:
|
||||
@@ -110,7 +111,7 @@ def gha1_ana_step(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: floa
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
|
||||
return p1, alpha1
|
||||
return p1, wrap_0_2pi(alpha1)
|
||||
|
||||
|
||||
def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> Tuple[NDArray, float]:
|
||||
@@ -134,7 +135,7 @@ def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, ma
|
||||
if h > 1e-5:
|
||||
raise Exception("GHA1_ana: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
|
||||
|
||||
return point_end, alpha_end
|
||||
return point_end, wrap_0_2pi(alpha_end)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
@@ -142,4 +143,3 @@ if __name__ == "__main__":
|
||||
p0 = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(20))
|
||||
p1, alpha1 = gha1_ana(ell, p0, wu.deg2rad(36), 200000, 70)
|
||||
print(p1, wu.rad2gms(alpha1))
|
||||
|
||||
|
||||
@@ -1,13 +1,18 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos
|
||||
from numpy.typing import NDArray
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from GHA_triaxial.gha1_ana import gha1_ana
|
||||
from GHA_triaxial.utils import func_sigma_ell, louville_constant, pq_ell
|
||||
import plotly.graph_objects as go
|
||||
import winkelumrechnungen as wu
|
||||
from typing import Tuple
|
||||
|
||||
def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray]:
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from numpy import cos, sin
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
from GHA_triaxial.utils import louville_constant, pq_ell
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from utils_angle import wrap_0_2pi
|
||||
|
||||
|
||||
def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: float, all_points: bool = False) \
|
||||
-> Tuple[NDArray, float] | Tuple[NDArray, float, NDArray, NDArray]:
|
||||
"""
|
||||
Berechung einer Näherungslösung der ersten Hauptaufgabe
|
||||
:param ell: Ellipsoid
|
||||
@@ -37,19 +42,24 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
|
||||
if last_p is not None and np.dot(p, last_p) < 0:
|
||||
p = -p
|
||||
q = -q
|
||||
last_p = p
|
||||
sigma = p * sin(alpha1) + q * cos(alpha1)
|
||||
if last_sigma is not None and np.dot(sigma, last_sigma) < 0:
|
||||
sigma = -sigma
|
||||
alpha1 += np.pi
|
||||
alpha1 = wrap_0_2pi(alpha1)
|
||||
p2 = p1 + ds_step * sigma
|
||||
p2 = ell.point_onto_ellipsoid(p2)
|
||||
|
||||
dalpha = 1e-6
|
||||
dalpha = 1e-9
|
||||
l2 = louville_constant(ell, p2, alpha1)
|
||||
dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
|
||||
alpha2 = alpha1 + (l0 - l2) / dl_dalpha
|
||||
|
||||
if abs(dl_dalpha) < 1e-20:
|
||||
alpha2 = alpha1 + 0
|
||||
else:
|
||||
alpha2 = alpha1 + (l0 - l2) / dl_dalpha
|
||||
points.append(p2)
|
||||
alphas.append(alpha2)
|
||||
alphas.append(wrap_0_2pi(alpha2))
|
||||
|
||||
ds_step = np.linalg.norm(p2 - p1)
|
||||
s_curr += ds_step
|
||||
@@ -88,11 +98,11 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray):
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
P0 = ell.ell2cart(wu.deg2rad(89), wu.deg2rad(1))
|
||||
alpha0 = wu.deg2rad(2)
|
||||
s = 200000
|
||||
P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=100, all_points=True)
|
||||
P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=20, maxPartCircum=2)
|
||||
print(np.linalg.norm(P1_app - P1_ana))
|
||||
show_points(points, P0, P1_ana)
|
||||
ell = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
P0 = ell.ell2cart(wu.deg2rad(15), wu.deg2rad(15))
|
||||
alpha0 = wu.deg2rad(270)
|
||||
s = 1
|
||||
P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=0.1, all_points=True)
|
||||
# P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=40, maxPartCircum=32)
|
||||
# print(np.linalg.norm(P1_app - P1_ana))
|
||||
# show_points(points, P0, P0)
|
||||
|
||||
@@ -1,15 +1,16 @@
|
||||
from typing import Callable, List, Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arctan2
|
||||
import ellipsoide
|
||||
import runge_kutta as rk
|
||||
import winkelumrechnungen as wu
|
||||
import GHA_triaxial.numeric_examples_karney as ne_karney
|
||||
from GHA_triaxial.gha1_ana import gha1_ana
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from typing import Callable, Tuple, List
|
||||
from numpy import arctan2, cos, sin
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import GHA_triaxial.numeric_examples_karney as ne_karney
|
||||
import runge_kutta as rk
|
||||
import winkelumrechnungen as wu
|
||||
from GHA_triaxial.gha1_ana import gha1_ana
|
||||
from GHA_triaxial.utils import alpha_ell2para, pq_ell
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from utils_angle import wrap_0_2pi
|
||||
|
||||
|
||||
def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
@@ -75,8 +76,7 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu
|
||||
|
||||
alpha1 = arctan2(P, Q)
|
||||
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
alpha1 = wrap_0_2pi(alpha1)
|
||||
|
||||
_, _, h = ell.cart2geod(point1, "ligas3")
|
||||
if h > 1e-5:
|
||||
@@ -108,7 +108,7 @@ if __name__ == "__main__":
|
||||
# diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
|
||||
# print(diffs_panou)
|
||||
|
||||
ell: EllipsoidTriaxial = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
diffs_karney = []
|
||||
# examples_karney = ne_karney.get_examples((30499, 30500, 40500))
|
||||
examples_karney = ne_karney.get_random_examples(20)
|
||||
|
||||
@@ -1,12 +1,13 @@
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from GHA_triaxial.gha2_num import gha2_num
|
||||
import plotly.graph_objects as go
|
||||
import winkelumrechnungen as wu
|
||||
from numpy.typing import NDArray
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
from GHA_triaxial.gha2_num import gha2_num
|
||||
from GHA_triaxial.utils import sigma2alpha
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
|
||||
|
||||
def gha2_approx(ell: EllipsoidTriaxial, p0: NDArray, p1: NDArray, ds: float, all_points: bool = False) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray]:
|
||||
|
||||
@@ -1,13 +1,15 @@
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from runge_kutta import rk4, rk4_step, rk4_end, rk4_integral
|
||||
from numpy.typing import NDArray
|
||||
|
||||
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
|
||||
import ausgaben as aus
|
||||
from utils_angle import cot, arccot, wrap_to_pi
|
||||
import winkelumrechnungen as wu
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from runge_kutta import rk4, rk4_end, rk4_integral
|
||||
from utils_angle import cot, wrap_0_2pi, wrap_mpi_pi
|
||||
|
||||
|
||||
def norm_a(a: float) -> float:
|
||||
@@ -22,7 +24,7 @@ def azimut(E: float, G: float, dbeta_du: float, dlamb_du: float) -> float:
|
||||
|
||||
|
||||
def sph_azimuth(beta1, lam1, beta2, lam2):
|
||||
dlam = wrap_to_pi(lam2 - lam1)
|
||||
dlam = wrap_mpi_pi(lam2 - lam1)
|
||||
y = np.sin(dlam) * np.cos(beta2)
|
||||
x = np.cos(beta1) * np.sin(beta2) - np.sin(beta1) * np.cos(beta2) * np.cos(dlam)
|
||||
a = np.arctan2(y, x)
|
||||
@@ -176,7 +178,7 @@ def gha2_num(
|
||||
)
|
||||
p_00 = 0.5 * ((E * G_beta_beta - E_beta * G_beta) / (E**2))
|
||||
|
||||
return (BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00)
|
||||
return BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00
|
||||
|
||||
# Berechnung der ODE Koeffizienten für Fall 2 (lambda_0 == lambda_1)
|
||||
def q_coef(beta, lamb):
|
||||
@@ -347,8 +349,8 @@ def gha2_num(
|
||||
|
||||
return best[0], best[1], sgn, dbeta, ode_beta
|
||||
|
||||
lamb0 = float(wrap_to_pi(lamb_0))
|
||||
lamb1 = float(wrap_to_pi(lamb_1))
|
||||
lamb0 = float(wrap_mpi_pi(lamb_0))
|
||||
lamb1 = float(wrap_mpi_pi(lamb_1))
|
||||
beta0 = float(beta_0)
|
||||
beta1 = float(beta_1)
|
||||
|
||||
@@ -491,7 +493,7 @@ def gha2_num(
|
||||
else:
|
||||
s = np.trapz(integrand, dx=h)
|
||||
|
||||
return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr
|
||||
return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr
|
||||
|
||||
_, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N_full, integrand_lambda)
|
||||
beta_end, beta_p_end, _, _ = y_end
|
||||
@@ -502,7 +504,7 @@ def gha2_num(
|
||||
(_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), float(lamb0 + dlamb))
|
||||
alpha_1 = azimut(E_end, G_end, dbeta_du=float(beta_p_end) * sgn, dlamb_du=1.0 * sgn)
|
||||
|
||||
return float(alpha_0), float(alpha_1), float(s)
|
||||
return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s)
|
||||
|
||||
# Fall 2 (lambda_0 == lambda_1)
|
||||
N = int(n)
|
||||
@@ -574,7 +576,7 @@ def gha2_num(
|
||||
else:
|
||||
s = np.trapz(integrand, dx=h)
|
||||
|
||||
return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr
|
||||
return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(s), beta_arr, lamb_arr
|
||||
|
||||
_, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta)
|
||||
lamb_end, lamb_p_end, _, _ = y_end
|
||||
@@ -585,57 +587,57 @@ def gha2_num(
|
||||
(_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, float(lamb_end))
|
||||
alpha_1 = azimut(E_end, G_end, dbeta_du=1.0 * sgn, dlamb_du=float(lamb_p_end) * sgn)
|
||||
|
||||
return float(alpha_0), float(alpha_1), float(s)
|
||||
return float(wrap_0_2pi(alpha_0)), float(wrap_0_2pi(alpha_1)), float(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)
|
||||
# a0, a1, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
|
||||
# print(aus.gms("a0", a0, 4))
|
||||
# print(aus.gms("a1", a1, 4))
|
||||
# print("s: ", s)
|
||||
# # 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)
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
beta1 = np.deg2rad(75)
|
||||
lamb1 = np.deg2rad(-90)
|
||||
beta2 = np.deg2rad(75)
|
||||
lamb2 = np.deg2rad(66)
|
||||
a0, a1, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=100)
|
||||
print(aus.gms("a0", a0, 4))
|
||||
print(aus.gms("a1", a1, 4))
|
||||
print("s: ", s)
|
||||
# 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)
|
||||
|
||||
# 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)
|
||||
a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
|
||||
print(s)
|
||||
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||||
diffs_panou = []
|
||||
examples_panou = ne_panou.get_random_examples(4)
|
||||
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,11 +1,12 @@
|
||||
import random
|
||||
from typing import List
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
from typing import List, Tuple
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from GHA_triaxial.gha1_ES import jacobi_konstante
|
||||
from GHA_triaxial.utils import jacobi_konstante
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
|
||||
ell = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
file_path = r"Karney_2024_Testset.txt"
|
||||
|
||||
def line2example(line: str) -> List:
|
||||
"""
|
||||
@@ -31,7 +32,7 @@ def get_random_examples(num: int, seed: int = None) -> List:
|
||||
"""
|
||||
if seed is not None:
|
||||
random.seed(seed)
|
||||
with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
|
||||
with open(file_path) as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
for i in range(num):
|
||||
@@ -46,7 +47,7 @@ def get_examples(l_i: List) -> List:
|
||||
:param l_i: Liste von Indizes
|
||||
:return: Liste mit Beispielen
|
||||
"""
|
||||
with open("Karney_2024_Testset.txt") as datei:
|
||||
with open(file_path) as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
for i in l_i:
|
||||
@@ -54,53 +55,21 @@ def get_examples(l_i: List) -> List:
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
# beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s
|
||||
|
||||
def get_random_examples_simple_short(num: int, seed: int = None) -> List:
|
||||
if seed is not None:
|
||||
random.seed(seed)
|
||||
with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
while len(examples) < num:
|
||||
example = line2example(lines[random.randint(0, len(lines) - 1)])
|
||||
beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
|
||||
if s < 1 and abs(abs(beta0) - np.pi/2) > 1e-5 and lamb0 != 0 and abs(abs(lamb0) - np.pi) > 1e-5:
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
def get_random_examples_umbilics_start(num: int, seed: int = None) -> List:
|
||||
if seed is not None:
|
||||
random.seed(seed)
|
||||
with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
while len(examples) < num:
|
||||
example = line2example(lines[random.randint(0, len(lines) - 1)])
|
||||
beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
|
||||
if abs(abs(beta0) - np.pi/2) < 1e-5 and (lamb0 == 0 or abs(abs(lamb0) - np.pi) < 1e-5):
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
def get_random_examples_umbilics_end(num: int, seed: int = None) -> List:
|
||||
if seed is not None:
|
||||
random.seed(seed)
|
||||
with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
while len(examples) < num:
|
||||
example = line2example(lines[random.randint(0, len(lines) - 1)])
|
||||
beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
|
||||
if abs(abs(beta1) - np.pi/2) < 1e-5 and (lamb1 == 0 or abs(abs(lamb1) - np.pi) < 1e-5):
|
||||
examples.append(example)
|
||||
return examples
|
||||
|
||||
def get_random_examples_gamma(group: str, num: int, seed: int = None, length: str = None) -> List:
|
||||
"""
|
||||
Zufällige Beispiele aus Karney in Gruppen nach Einteilung anhand der Jacobi-Konstanten
|
||||
:param group: Gruppe
|
||||
:param num: Anzahl
|
||||
:param seed: Random-Seed
|
||||
:param length: long oder short, sond egal
|
||||
:return: Liste mit Beispielen
|
||||
"""
|
||||
eps = 1e-20
|
||||
long_short = 2
|
||||
if seed is not None:
|
||||
random.seed(seed)
|
||||
with open(r"C:\Users\moell\OneDrive\Desktop\Vorlesungen\Master-Projekt\Python_Masterprojekt\GHA_triaxial\Karney_2024_Testset.txt") as datei:
|
||||
with open(file_path) as datei:
|
||||
lines = datei.readlines()
|
||||
examples = []
|
||||
i = 0
|
||||
@@ -113,7 +82,7 @@ def get_random_examples_gamma(group: str, num: int, seed: int = None, length: st
|
||||
beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
|
||||
gamma = jacobi_konstante(beta0, lamb0, alpha0_ell, ell)
|
||||
|
||||
if group not in ["a", "b", "c", "d", "e"]:
|
||||
if group not in ["a", "b", "c", "d", "e", "de"]:
|
||||
break
|
||||
elif group == "a" and not 1 >= gamma >= 0.01:
|
||||
continue
|
||||
@@ -125,6 +94,8 @@ def get_random_examples_gamma(group: str, num: int, seed: int = None, length: st
|
||||
continue
|
||||
elif group == "e" and not -1e-17 >= gamma >= -1:
|
||||
continue
|
||||
elif group == "de" and not -eps > gamma > -1:
|
||||
continue
|
||||
|
||||
if length == "short":
|
||||
if example[6] < long_short:
|
||||
|
||||
@@ -1,11 +1,13 @@
|
||||
from __future__ import annotations
|
||||
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import arctan2, sin, cos, sqrt
|
||||
from numpy._typing import NDArray
|
||||
from numpy import arctan2, cos, sin, sqrt
|
||||
from numpy.typing import NDArray
|
||||
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
from utils_angle import wrap_0_2pi
|
||||
|
||||
|
||||
def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float:
|
||||
@@ -21,7 +23,7 @@ def sigma2alpha(ell: EllipsoidTriaxial, sigma: NDArray, point: NDArray) -> float
|
||||
Q = float(q @ sigma)
|
||||
|
||||
alpha = arctan2(P, Q)
|
||||
return alpha
|
||||
return wrap_0_2pi(alpha)
|
||||
|
||||
|
||||
def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> Tuple[float, float, float]:
|
||||
@@ -43,10 +45,10 @@ def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float
|
||||
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-9:
|
||||
if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
|
||||
raise Exception("alpha_para2ell: Differenz in den Richtungsableitungen")
|
||||
|
||||
return beta, lamb, alpha_ell
|
||||
return beta, lamb, wrap_0_2pi(alpha_ell)
|
||||
|
||||
|
||||
def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> Tuple[float, float, float]:
|
||||
@@ -68,10 +70,10 @@ def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell:
|
||||
alpha_para = arctan2(p_para @ sigma_ell, q_para @ sigma_ell)
|
||||
sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
|
||||
|
||||
if np.linalg.norm(sigma_para - sigma_ell) > 1e-9:
|
||||
if np.linalg.norm(sigma_para - sigma_ell) > 1e-7:
|
||||
raise Exception("alpha_ell2para: Differenz in den Richtungsableitungen")
|
||||
|
||||
return u, v, alpha_para
|
||||
return u, v, wrap_0_2pi(alpha_para)
|
||||
|
||||
|
||||
def func_sigma_ell(ell: EllipsoidTriaxial, point: NDArray, alpha_ell: float) -> NDArray:
|
||||
@@ -124,11 +126,10 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
:param point: Punkt
|
||||
:return: p und q
|
||||
"""
|
||||
x, y, z = point
|
||||
n = ell.func_n(point)
|
||||
|
||||
beta, lamb = ell.cart2ell(point)
|
||||
if abs(cos(beta)) < 1e-12 and abs(np.sin(lamb)) < 1e-12:
|
||||
if abs(cos(beta)) < 1e-15 and abs(np.sin(lamb)) < 1e-15:
|
||||
if beta > 0:
|
||||
p = np.array([0, -1, 0])
|
||||
else:
|
||||
@@ -137,11 +138,7 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2
|
||||
L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2
|
||||
|
||||
c1 = x ** 2 + y ** 2 + z ** 2 - (ell.ax ** 2 + ell.ay ** 2 + ell.b ** 2)
|
||||
c0 = (ell.ax ** 2 * ell.ay ** 2 + ell.ax ** 2 * ell.b ** 2 + ell.ay ** 2 * ell.b ** 2 -
|
||||
(ell.ay ** 2 + ell.b ** 2) * x ** 2 - (ell.ax ** 2 + ell.b ** 2) * y ** 2 - (
|
||||
ell.ax ** 2 + ell.ay ** 2) * z ** 2)
|
||||
t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
|
||||
_, t2 = ell.func_t12(point)
|
||||
|
||||
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)
|
||||
@@ -181,3 +178,24 @@ def pq_para(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
q = q / np.linalg.norm(q)
|
||||
|
||||
return p, q
|
||||
|
||||
|
||||
def jacobi_konstante(beta: float, omega: float, alpha: float, ell: EllipsoidTriaxial) -> float:
|
||||
"""
|
||||
Jacobi-Konstante nach Karney (2025), Gl. (14)
|
||||
:param beta: Beta Koordinate
|
||||
:param omega: Omega Koordinate
|
||||
:param alpha: Azimut alpha
|
||||
:param ell: Ellipsoid
|
||||
:return: Jacobi-Konstante
|
||||
"""
|
||||
gamma_jacobi = float((ell.k ** 2) * (np.cos(beta) ** 2) * (np.sin(alpha) ** 2) - (ell.k_ ** 2) * (np.sin(omega) ** 2) * (np.cos(alpha) ** 2))
|
||||
return gamma_jacobi
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
alpha_para = 0
|
||||
u, v = ell.ell2para(np.pi/2, 0)
|
||||
alpha_ell = alpha_para2ell(ell, u, v, alpha_para)
|
||||
pass
|
||||
|
||||
Binary file not shown.
@@ -1,34 +0,0 @@
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
from GHA_biaxial.bessel import gha1 as gha1_bessel
|
||||
from GHA_biaxial.gauss import gha1 as gha1_gauss
|
||||
from GHA_biaxial.rk import gha1 as gha1_rk
|
||||
from GHA_biaxial.gauss import gha2 as gha2_gauss
|
||||
|
||||
re = EllipsoidBiaxial.init_name("Bessel")
|
||||
|
||||
# phi0 = 0.6
|
||||
# lamb0 = 1.2
|
||||
# alpha0 = 0.45
|
||||
# s = 123456
|
||||
#
|
||||
# values_bessel = gha1_bessel(re, phi0, lamb0, alpha0, s)
|
||||
# alpha1_bessel = values_bessel[-1]
|
||||
# p1_bessel = re.bi_ell2cart(values_bessel[0], values_bessel[1], 0)
|
||||
#
|
||||
# values_gauss1 = gha1_gauss(re, phi0, lamb0, alpha0, s)
|
||||
# alpha1_gauss1 = values_gauss1[-1]
|
||||
# p1_gauss = re.bi_ell2cart(values_gauss1[0], values_gauss1[1], 0)
|
||||
#
|
||||
# values_rk = gha1_rk(re, phi0, lamb0 , alpha0, s, 10000)
|
||||
# alpha1_rk = values_rk[-1]
|
||||
# p1_rk = re.bi_ell2cart(values_rk[0], values_rk[1], 0)
|
||||
#
|
||||
# alpha0_gauss, alpha1_gauss2, s_gauss = gha2_gauss(re, phi0, lamb0, values_gauss1[0], values_gauss1[1])
|
||||
|
||||
phi0 = 0.6
|
||||
lamb0 = 1.2
|
||||
|
||||
cart = re.bi_ell2cart(phi0, lamb0, 0)
|
||||
ell = re.bi_cart2ell(cart)
|
||||
pass
|
||||
@@ -10,7 +10,7 @@ def xyz(x: float, y: float, z: float, stellen: int) -> str:
|
||||
:param stellen: Anzahl Nachkommastellen
|
||||
:return: String zur Ausgabe der Koordinaten
|
||||
"""
|
||||
return f"""x = {(round(x,stellen))} m y = {(round(y,stellen))} m z = {(round(z,stellen))} m"""
|
||||
return f"""x = {(round(x, stellen))} m y = {(round(y, stellen))} m z = {(round(z, stellen))} m"""
|
||||
|
||||
|
||||
def gms(name: str, rad: float, stellen: int) -> str:
|
||||
@@ -21,5 +21,5 @@ def gms(name: str, rad: float, stellen: int) -> str:
|
||||
:param stellen: Anzahl Nachkommastellen
|
||||
:return: String zur Ausgabe des Winkels
|
||||
"""
|
||||
gms = wu.rad2gms(rad)
|
||||
return f"{name} = {int(gms[0])}° {int(gms[1])}' {round(gms[2],stellen):.{stellen}f}''"
|
||||
values = wu.rad2gms(rad)
|
||||
return f"{name} = {int(values[0])}° {int(values[1])}' {round(values[2], stellen):.{stellen}f}''"
|
||||
|
||||
88
dashboard.py
88
dashboard.py
@@ -1,34 +1,31 @@
|
||||
from dash import Dash, dash, html, dcc, Input, Output, State, no_update, ctx
|
||||
import plotly.graph_objects as go
|
||||
import numpy as np
|
||||
import dash_bootstrap_components as dbc
|
||||
|
||||
import builtins
|
||||
from dash.exceptions import PreventUpdate
|
||||
import traceback
|
||||
|
||||
import webbrowser
|
||||
from threading import Timer
|
||||
import dash_bootstrap_components as dbc
|
||||
import numpy as np
|
||||
import plotly.graph_objects as go
|
||||
from dash import Dash, Input, Output, State, dcc, html, no_update
|
||||
from dash.exceptions import PreventUpdate
|
||||
from numpy import pi
|
||||
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import winkelumrechnungen as wu
|
||||
import ausgaben as aus
|
||||
from GHA_triaxial.utils import alpha_ell2para, alpha_para2ell
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
from ES.gha1_ES import gha1_ES
|
||||
from ES.gha2_ES import gha2_ES
|
||||
from GHA_triaxial.gha1_ana import gha1_ana
|
||||
from GHA_triaxial.gha1_num import gha1_num
|
||||
from GHA_triaxial.gha1_ES import gha1_ES
|
||||
from GHA_triaxial.gha1_approx import gha1_approx
|
||||
|
||||
from GHA_triaxial.gha2_num import gha2_num
|
||||
from GHA_triaxial.gha2_ES import gha2_ES
|
||||
from GHA_triaxial.gha1_num import gha1_num
|
||||
from GHA_triaxial.gha2_approx import gha2_approx
|
||||
from GHA_triaxial.gha2_num import gha2_num
|
||||
from GHA_triaxial.utils import alpha_ell2para, alpha_para2ell
|
||||
from ellipsoid_triaxial import EllipsoidTriaxial
|
||||
|
||||
|
||||
# Prints von importierten Funktionen unterdücken
|
||||
def _no_print(*args, **kwargs):
|
||||
pass
|
||||
|
||||
|
||||
builtins.print = _no_print
|
||||
|
||||
|
||||
@@ -39,7 +36,7 @@ app.title = "Geodätische Hauptaufgaben"
|
||||
|
||||
|
||||
# Erzeugen der Eingabefelder
|
||||
def inputfeld(left_text, input_id, right_text="", width=200, min=None, max=None):
|
||||
def inputfeld(left_text, input_id, right_text="", width=200, mini=None, maxi=None):
|
||||
return html.Div(
|
||||
children=[
|
||||
html.Span(f"{left_text} =", style={"minWidth": 36, "textAlign": "right", "marginRight": 5}),
|
||||
@@ -142,7 +139,7 @@ def method_failed(method_label: str, exc: Exception):
|
||||
return html.Div([
|
||||
html.Strong(f"{method_label}: "),
|
||||
html.Span("konnte nicht berechnet werden. ", style={"color": "red"}),
|
||||
#html.Span(f"({type(exc).__name__}: {exc})", style={"color": "#b02a37"}),
|
||||
# html.Span(f"({type(exc).__name__}: {exc})", style={"color": "#b02a37"}),
|
||||
|
||||
html.Details([
|
||||
html.Summary("Details"),
|
||||
@@ -177,7 +174,7 @@ def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
|
||||
scene=dict(
|
||||
xaxis=dict(
|
||||
range=[-rx, rx],
|
||||
#title="X [m]",
|
||||
# title="X [m]",
|
||||
title="",
|
||||
showgrid=False,
|
||||
zeroline=False,
|
||||
@@ -186,7 +183,7 @@ def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
|
||||
),
|
||||
yaxis=dict(
|
||||
range=[-ry, ry],
|
||||
#title="Y [m]",
|
||||
# title="Y [m]",
|
||||
title="",
|
||||
showgrid=False,
|
||||
zeroline=False,
|
||||
@@ -195,7 +192,7 @@ def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
|
||||
),
|
||||
zaxis=dict(
|
||||
range=[-rz, rz],
|
||||
#title="Z [m]",
|
||||
# title="Z [m]",
|
||||
title="",
|
||||
showgrid=False,
|
||||
zeroline=False,
|
||||
@@ -209,8 +206,8 @@ def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
|
||||
)
|
||||
|
||||
# Ellipsoid
|
||||
u = np.linspace(-np.pi/2, np.pi/2, 80)
|
||||
v = np.linspace(-np.pi, np.pi, 160)
|
||||
u = np.linspace(-pi/2, pi/2, 80)
|
||||
v = np.linspace(-pi, pi, 160)
|
||||
U, V = np.meshgrid(u, v)
|
||||
X, Y, Z = ell.para2cart(U, V)
|
||||
fig.add_trace(go.Surface(
|
||||
@@ -260,7 +257,7 @@ def figure_constant_lines(fig, ell: EllipsoidTriaxial, coordsystem: str = "para"
|
||||
all_beta[-1] -= 1e-8
|
||||
constants_lamb = wu.deg2rad(np.arange(-180, 180, 15))
|
||||
for lamb in constants_lamb:
|
||||
if lamb != 0 and abs(lamb) != np.pi:
|
||||
if lamb != 0 and abs(lamb) != pi:
|
||||
xyz = ell.ell2cart(all_beta, lamb)
|
||||
fig.add_trace(go.Scatter3d(
|
||||
x=xyz[:, 0], y=xyz[:, 1], z=xyz[:, 2], mode="lines",
|
||||
@@ -335,8 +332,10 @@ def figure_lines(fig, line, name, color):
|
||||
))
|
||||
return fig
|
||||
|
||||
|
||||
# HTML der beiden Tabs
|
||||
# Tab 1
|
||||
|
||||
pane_gha1 = html.Div(
|
||||
[
|
||||
html.Div(
|
||||
@@ -468,7 +467,7 @@ app.layout = html.Div(
|
||||
style={"fontFamily": "Arial", "padding": "10px", "width": "95%", "margin": "0 auto"},
|
||||
children=[
|
||||
html.H2("Geodätische Hauptaufgaben für dreiachsige Ellipsoide"),
|
||||
#html.H2("für dreiachsige Ellipsoide"),
|
||||
# html.H2("für dreiachsige Ellipsoide"),
|
||||
|
||||
html.Div(
|
||||
style={
|
||||
@@ -488,15 +487,13 @@ app.layout = html.Div(
|
||||
dcc.Dropdown(
|
||||
id="dropdown-ellipsoid",
|
||||
options=[
|
||||
{"label": "BursaFialova1993", "value": "BursaFialova1993"},
|
||||
{"label": "BursaSima1980", "value": "BursaSima1980"},
|
||||
{"label": "BursaSima1980round", "value": "BursaSima1980round"},
|
||||
{"label": "Eitschberger1978", "value": "Eitschberger1978"},
|
||||
{"label": "Bursa1972", "value": "Bursa1972"},
|
||||
{"label": "Bursa1970", "value": "Bursa1970"},
|
||||
{"label": "BesselBiaxial", "value": "BesselBiaxial"},
|
||||
{"label": "KarneyTest2024", "value": "KarneyTest2024"},
|
||||
{"label": "Fiction", "value": "Fiction"},
|
||||
{"label": "Burša und Šíma (1980)", "value": "BursaSima1980round"},
|
||||
{"label": "Karney (2024)", "value": "KarneyTest2024"},
|
||||
{"label": "Fictional", "value": "Fiction"},
|
||||
{"label": "Burša und Fialová (1993)", "value": "BursaFialova1993"},
|
||||
{"label": "Eitschberger (1978)", "value": "Eitschberger1978"},
|
||||
{"label": "Burša (1972)", "value": "Bursa1972"},
|
||||
{"label": "Burša (1970)", "value": "Bursa1970"}
|
||||
|
||||
],
|
||||
value="",
|
||||
@@ -507,9 +504,9 @@ app.layout = html.Div(
|
||||
|
||||
html.Div(
|
||||
[
|
||||
inputfeld("aₓ", "input-ax", "m", min=0, width="clamp(80px, 7vw, 200px)"),
|
||||
inputfeld("aᵧ", "input-ay", "m", min=0, width="clamp(80px, 7vw, 200px)"),
|
||||
inputfeld("b", "input-b", "m", min=0, width="clamp(80px, 7vw, 200px)"),
|
||||
inputfeld("aₓ", "input-ax", "m", mini=0, width="clamp(80px, 7vw, 200px)"),
|
||||
inputfeld("aᵧ", "input-ay", "m", mini=0, width="clamp(80px, 7vw, 200px)"),
|
||||
inputfeld("b", "input-b", "m", mini=0, width="clamp(80px, 7vw, 200px)"),
|
||||
],
|
||||
style={
|
||||
"display": "grid",
|
||||
@@ -520,7 +517,7 @@ app.layout = html.Div(
|
||||
},
|
||||
),
|
||||
|
||||
#html.Br(),
|
||||
# html.Br(),
|
||||
|
||||
dcc.Tabs(
|
||||
id="tabs-GHA",
|
||||
@@ -572,7 +569,7 @@ app.layout = html.Div(
|
||||
dcc.Store(id="calc-token-gha1", data=0),
|
||||
dcc.Store(id="calc-token-gha2", data=0),
|
||||
|
||||
#html.P("© 2026", style={"fontSize": "10px", "color": "gray", "textAlign": "center", "marginTop": "16px"}),
|
||||
# html.P("© 2026", style={"fontSize": "10px", "color": "gray", "textAlign": "center", "marginTop": "16px"}),
|
||||
],
|
||||
|
||||
|
||||
@@ -662,10 +659,8 @@ def toggle_ds(v):
|
||||
return "on" not in (v or [])
|
||||
|
||||
|
||||
|
||||
# Abfrage ob Berechnungsverfahren gewählt
|
||||
from dash.exceptions import PreventUpdate
|
||||
from dash import no_update, html
|
||||
|
||||
|
||||
@app.callback(
|
||||
Output("calc-token-gha1", "data"),
|
||||
@@ -961,7 +956,7 @@ def compute_gha1_ana(n1, cb_ana, max_M, maxPartCircum, beta0, lamb0, s, a0, ax,
|
||||
def compute_gha1_num(n1, cb_num, n_in, beta0, lamb0, s, a0, ax, ay, b):
|
||||
if not n1:
|
||||
return no_update, no_update
|
||||
if "on" not in (cb_num or []):
|
||||
if "on" not in (cb_num or []):
|
||||
return "", None
|
||||
|
||||
n_in = int(n_in) if n_in else 2000
|
||||
@@ -974,7 +969,6 @@ def compute_gha1_num(n1, cb_num, n_in, beta0, lamb0, s, a0, ax, ay, b):
|
||||
alpha_rad = wu.deg2rad(float(a0))
|
||||
s_val = float(s)
|
||||
|
||||
|
||||
P0 = ell.ell2cart(beta_rad, lamb_rad)
|
||||
|
||||
P1_num, alpha1, werte = gha1_num(ell, P0, alpha_rad, s_val, n_in, all_points=True)
|
||||
@@ -1402,7 +1396,7 @@ def clear_all_stores_on_ellipsoid_change(ax, ay, b):
|
||||
if None in (ax, ay, b):
|
||||
return (no_update,)*7
|
||||
|
||||
return (None, None, None, None, None, None, None)
|
||||
return None, None, None, None, None, None, None
|
||||
|
||||
# Funktionen zur Erzeugung der Überschriften
|
||||
@app.callback(
|
||||
@@ -1483,6 +1477,6 @@ if __name__ == "__main__":
|
||||
# Automatisiertes Öffnen der Seite im Browser
|
||||
HOST = "127.0.0.1"
|
||||
PORT = 8050
|
||||
#Timer(1.0, webbrowser.open_new_tab(f"http://{HOST}:{PORT}/")).start
|
||||
# Timer(1.0, webbrowser.open_new_tab(f"http://{HOST}:{PORT}/")).start
|
||||
|
||||
app.run(host=HOST, port=PORT, debug=False)
|
||||
|
||||
@@ -1,116 +1,20 @@
|
||||
import numpy as np
|
||||
from numpy import sin, cos, arctan, arctan2, sqrt, pi, arccos
|
||||
import winkelumrechnungen as wu
|
||||
import jacobian_Ligas
|
||||
import matplotlib.pyplot as plt
|
||||
from typing import Tuple
|
||||
from numpy.typing import NDArray
|
||||
import math
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import arccos, arctan, arctan2, cos, pi, sin, sqrt
|
||||
from numpy.typing import NDArray
|
||||
|
||||
class EllipsoidBiaxial:
|
||||
def __init__(self, a: float, b: float):
|
||||
self.a = a
|
||||
self.b = b
|
||||
self.c = a ** 2 / b
|
||||
self.e = sqrt(a ** 2 - b ** 2) / a
|
||||
self.e_ = sqrt(a ** 2 - b ** 2) / b
|
||||
import jacobian_Ligas
|
||||
from utils_angle import wrap_mhalfpi_halfpi, wrap_mpi_pi
|
||||
|
||||
@classmethod
|
||||
def init_name(cls, name: str):
|
||||
if name == "Bessel":
|
||||
a = 6377397.15508
|
||||
b = 6356078.96290
|
||||
return cls(a, b)
|
||||
elif name == "Hayford":
|
||||
a = 6378388
|
||||
f = 1/297
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
elif name == "Krassowski":
|
||||
a = 6378245
|
||||
f = 298.3
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
elif name == "WGS84":
|
||||
a = 6378137
|
||||
f = 298.257223563
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
|
||||
@classmethod
|
||||
def init_af(cls, a: float, f: float):
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
|
||||
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.arctan2(self.a * np.sin(beta), self.b * np.cos(beta))
|
||||
beta2phi = lambda self, beta: np.arctan2(self.a ** 2 * np.sin(beta), self.b ** 2 * np.cos(beta))
|
||||
|
||||
psi2beta = lambda self, psi: np.arctan2(self.b * np.sin(psi), self.a * np.cos(psi))
|
||||
psi2phi = lambda self, psi: np.arctan2(self.a * np.sin(psi), self.b * np.cos(psi))
|
||||
|
||||
phi2beta = lambda self, phi: np.arctan2(self.b**2 * np.sin(phi), self.a**2 * np.cos(phi))
|
||||
phi2psi = lambda self, phi: np.arctan2(self.b * np.sin(phi), self.a * np.cos(phi))
|
||||
|
||||
phi2p = lambda self, phi: self.N(phi) * cos(phi)
|
||||
|
||||
def bi_cart2ell(self, point: NDArray, Eh: float = 0.001, Ephi: float = wu.gms2rad([0, 0, 0.001])) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Umrechnung von kartesischen in ellipsoidische Koordinaten auf einem Rotationsellipsoid
|
||||
# TODO: Quelle
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:param Eh: Grenzwert für die Höhe
|
||||
:param Ephi: Grenzwert für die Breite
|
||||
:return: ellipsoidische Breite, Länge, geodätische Höhe
|
||||
"""
|
||||
x, y, z = point
|
||||
|
||||
lamb = arctan2(y, x)
|
||||
|
||||
p = sqrt(x**2+y**2)
|
||||
|
||||
phi_null = arctan2(z, p*(1 - self.e**2))
|
||||
|
||||
hi = [0]
|
||||
phii = [phi_null]
|
||||
|
||||
i = 0
|
||||
|
||||
while True:
|
||||
N = self.a / sqrt(1 - self.e**2 * sin(phii[i])**2)
|
||||
h = p / cos(phii[i]) - N
|
||||
phi = arctan2(z, p * (1-(self.e**2*N) / (N+h)))
|
||||
hi.append(h)
|
||||
phii.append(phi)
|
||||
dh = abs(hi[i]-h)
|
||||
dphi = abs(phii[i]-phi)
|
||||
i = i+1
|
||||
if dh < Eh:
|
||||
if dphi < Ephi:
|
||||
break
|
||||
return phi, lamb, h
|
||||
|
||||
def bi_ell2cart(self, phi: float, lamb: float, h: float) -> NDArray:
|
||||
"""
|
||||
Umrechnung von ellipsoidischen in kartesische Koordinaten auf einem Rotationsellipsoid
|
||||
# TODO: Quelle
|
||||
:param phi: ellipsoidische Breite
|
||||
:param lamb: ellipsoidische Länge
|
||||
:param h: geodätische Höhe
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
W = sqrt(1 - self.e**2 * sin(phi)**2)
|
||||
N = self.a / W
|
||||
x = (N+h) * cos(phi) * cos(lamb)
|
||||
y = (N+h) * cos(phi) * sin(lamb)
|
||||
z = (N * (1-self.e**2) + h) * sin(phi)
|
||||
return np.array([x, y, z])
|
||||
|
||||
class EllipsoidTriaxial:
|
||||
"""
|
||||
Klasse für dreiachsige Ellipsoide
|
||||
Parameter: Formparameter
|
||||
Funktionen: Koordinatenumrechnungen
|
||||
"""
|
||||
def __init__(self, ax: float, ay: float, b: float):
|
||||
self.ax = ax
|
||||
self.ay = ay
|
||||
@@ -124,14 +28,19 @@ class EllipsoidTriaxial:
|
||||
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)
|
||||
nenner = sqrt(max(self.ax * self.ax - self.b * self.b, 0.0))
|
||||
self.k = sqrt(max(self.ay * self.ay - self.b * self.b, 0.0)) / nenner
|
||||
self.k_ = sqrt(max(self.ax * self.ax - self.ay * self.ay, 0.0)) / nenner
|
||||
self.e = sqrt(max(self.ax * self.ax - self.b * self.b, 0.0)) / self.ay
|
||||
|
||||
@classmethod
|
||||
def init_name(cls, name: str):
|
||||
def init_name(cls, name: str) -> EllipsoidTriaxial:
|
||||
"""
|
||||
Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, BursaSima1980round, Eitschberger1978, Bursa1972,
|
||||
Bursa1970, BesselBiaxial, Fiction, KarneyTest2024
|
||||
Mögliche Ellipsoide: BursaSima1980round, KarneyTest2024, Fiction, BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972,
|
||||
Bursa1970
|
||||
Panou et al (2020)
|
||||
:param name: Name des dreiachsigen Ellipsoids
|
||||
:return: dreiachsiger Ellipsoid
|
||||
"""
|
||||
if name == "BursaFialova1993":
|
||||
ax = 6378171.36
|
||||
@@ -164,11 +73,6 @@ class EllipsoidTriaxial:
|
||||
ay = 6378105
|
||||
b = 6356754
|
||||
return cls(ax, ay, b)
|
||||
elif name == "BesselBiaxial":
|
||||
ax = 6377397.15509
|
||||
ay = 6377397.15508
|
||||
b = 6356078.96290
|
||||
return cls(ax, ay, b)
|
||||
elif name == "Fiction":
|
||||
ax = 6000000
|
||||
ay = 4000000
|
||||
@@ -179,6 +83,8 @@ class EllipsoidTriaxial:
|
||||
ay = 1
|
||||
b = 1 / sqrt(2)
|
||||
return cls(ax, ay, b)
|
||||
else:
|
||||
raise Exception(f"EllipsoidTriaxial.init_name: Name {name} unbekannt")
|
||||
|
||||
def func_H(self, point: NDArray) -> float:
|
||||
"""
|
||||
@@ -218,8 +124,10 @@ class EllipsoidTriaxial:
|
||||
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)
|
||||
if c1 ** 2 - 4 * c0 < 0:
|
||||
t2 = np.nan
|
||||
if c1 ** 2 - 4 * c0 < -1e-9:
|
||||
raise Exception("t1, t2: Negativer Wurzelterm")
|
||||
elif c1 ** 2 - 4 * c0 < 0:
|
||||
t2 = 0
|
||||
else:
|
||||
t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
|
||||
if t2 == 0:
|
||||
@@ -261,7 +169,6 @@ class EllipsoidTriaxial:
|
||||
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 = arctan(sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
|
||||
if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7:
|
||||
@@ -284,6 +191,11 @@ class EllipsoidTriaxial:
|
||||
|
||||
beta, lamb = np.broadcast_arrays(beta, lamb)
|
||||
|
||||
beta = np.where(
|
||||
np.isclose(np.abs(beta), pi / 2, atol=1e-15),
|
||||
beta * 8999999999999999 / 9000000000000000,
|
||||
beta
|
||||
)
|
||||
B = self.Ex ** 2 * cos(beta) ** 2 + self.Ee ** 2 * sin(beta) ** 2
|
||||
L = self.Ex ** 2 - self.Ee ** 2 * cos(lamb) ** 2
|
||||
|
||||
@@ -419,7 +331,7 @@ class EllipsoidTriaxial:
|
||||
if delta_r > 1e-6:
|
||||
raise Exception("Umrechnung cart2ell: Punktdifferenz")
|
||||
|
||||
return beta, lamb
|
||||
return wrap_mhalfpi_halfpi(beta), wrap_mpi_pi(lamb)
|
||||
|
||||
except Exception as e:
|
||||
# Wenn die Berechnung fehlschlägt auf Grund von sehr kleinem y, solange anpassen, bis Umrechnung ohne Fehler
|
||||
@@ -577,6 +489,8 @@ class EllipsoidTriaxial:
|
||||
invJ, fxE = jacobian_Ligas.case2(E, F, G, np.array([xG, yG, zG]), pE)
|
||||
elif mode == "ligas3":
|
||||
invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE)
|
||||
else:
|
||||
raise Exception(f"cart2geod: Modus {mode} nicht bekannt")
|
||||
pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1)
|
||||
pEi = pEi.reshape(1, -1).flatten()
|
||||
loa = sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
|
||||
@@ -598,15 +512,15 @@ class EllipsoidTriaxial:
|
||||
phi, lamb, h = self.cart2geod(point, f"ligas{new_mode}", maxIter, maxLoa)
|
||||
else:
|
||||
if xG < 0 and yG < 0:
|
||||
lamb = -pi + lamb
|
||||
lamb += -pi
|
||||
|
||||
elif xG < 0:
|
||||
lamb = pi + lamb
|
||||
lamb += pi
|
||||
|
||||
if abs(zG) < eps:
|
||||
phi = 0
|
||||
|
||||
return phi, lamb, h
|
||||
wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb)
|
||||
return wrap_mhalfpi_halfpi(phi), wrap_mpi_pi(lamb), h
|
||||
|
||||
def para2cart(self, u: float | NDArray, v: float | NDArray) -> NDArray:
|
||||
"""
|
||||
@@ -643,8 +557,8 @@ class EllipsoidTriaxial:
|
||||
v = 2 * arctan2(v_check1, v_check2 + v_factor)
|
||||
else:
|
||||
v = pi/2 - 2 * arctan2(v_check2, v_check1 + v_factor)
|
||||
|
||||
return u, v
|
||||
wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v)
|
||||
return wrap_mhalfpi_halfpi(u), wrap_mpi_pi(v)
|
||||
|
||||
def ell2para(self, beta: float, lamb: float) -> Tuple[float, float]:
|
||||
"""
|
||||
@@ -749,63 +663,71 @@ class EllipsoidTriaxial:
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
ell = EllipsoidTriaxial.init_name("BursaSima1980")
|
||||
diff_list = []
|
||||
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)
|
||||
ell = EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
# cart = ell.ell2cart(pi/2, 0)
|
||||
# print(cart)
|
||||
# cart = ell.ell2cart(pi/2*8999999999999999/9000000000000000, 0)
|
||||
# print(cart)
|
||||
elli = ell.cart2ell(np.array([0, 0.0, 1/sqrt(2)]))
|
||||
print(elli)
|
||||
|
||||
elli = ell.cart2ell(point)
|
||||
cart_elli = ell.ell2cart(elli[0], elli[1])
|
||||
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.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.linalg.norm(point - cart_geod, axis=-1)
|
||||
|
||||
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)
|
||||
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()
|
||||
# ell = EllipsoidTriaxial.init_name("BursaSima1980")
|
||||
# diff_list = []
|
||||
# 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.linalg.norm(point - cart_elli, axis=-1)
|
||||
#
|
||||
# para = ell.cart2para(point)
|
||||
# cart_para = ell.para2cart(para[0], para[1])
|
||||
# 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.linalg.norm(point - cart_geod, axis=-1)
|
||||
#
|
||||
# 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)
|
||||
# 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,6 +1,8 @@
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy.typing import NDArray
|
||||
from typing import Tuple
|
||||
|
||||
|
||||
def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
"""
|
||||
@@ -34,7 +36,7 @@ def case1(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArr
|
||||
|
||||
return invJ, fxE
|
||||
|
||||
def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, NDArray]:
|
||||
def case2(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
"""
|
||||
Aufstellen des Gleichungssystem für den zweiten Fall
|
||||
:param E: Konstante E
|
||||
@@ -68,7 +70,7 @@ def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple
|
||||
|
||||
return invJ, fxE
|
||||
|
||||
def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray) -> Tuple[NDArray, NDArray]:
|
||||
def case3(E: float, F: float, G: float, pG: NDArray, pE: NDArray) -> Tuple[NDArray, NDArray]:
|
||||
"""
|
||||
Aufstellen des Gleichungssystem für den dritten Fall
|
||||
:param E: Konstante E
|
||||
|
||||
@@ -1,10 +1,20 @@
|
||||
from numpy import *
|
||||
import scipy as sp
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
from typing import Tuple
|
||||
|
||||
import scipy as sp
|
||||
from ellipsoid_biaxial import EllipsoidBiaxial
|
||||
from numpy import *
|
||||
|
||||
|
||||
def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0:float, s: float) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Berechnung der 1.GHA auf einem Rotationsellipsoid nach Bessel
|
||||
:param re:
|
||||
:param phi0:
|
||||
:param lamb0:
|
||||
:param alpha0:
|
||||
:param s:
|
||||
:return:
|
||||
"""
|
||||
psi0 = re.phi2psi(phi0)
|
||||
clairant = arcsin(cos(psi0) * sin(alpha0))
|
||||
sigma0 = arcsin(sin(psi0) / cos(clairant))
|
||||
@@ -1,8 +1,8 @@
|
||||
from numpy import sin, cos, pi, sqrt, tan, arcsin, arccos, arctan
|
||||
import ausgaben as aus
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
from typing import Tuple
|
||||
|
||||
from ellipsoid_biaxial import EllipsoidBiaxial
|
||||
from numpy import arctan, cos, sin, sqrt, tan
|
||||
|
||||
|
||||
def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, eps: float = 1e-12) -> Tuple[float, float, float]:
|
||||
"""
|
||||
@@ -1,13 +1,26 @@
|
||||
import runge_kutta as rk
|
||||
from numpy import sin, cos, tan
|
||||
import winkelumrechnungen as wu
|
||||
from ellipsoide import EllipsoidBiaxial
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from ellipsoid_biaxial import EllipsoidBiaxial
|
||||
from numpy import cos, sin, tan
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import runge_kutta as rk
|
||||
|
||||
|
||||
def gha1(re: EllipsoidBiaxial, phi0: float, lamb0: float, alpha0: float, s: float, num: int) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Berechnung der 1. GHA auf einem Rotationsellipsoid mittels RK4
|
||||
:param re:
|
||||
:param phi0:
|
||||
:param lamb0:
|
||||
:param alpha0:
|
||||
:param s:
|
||||
:param num:
|
||||
:return:
|
||||
"""
|
||||
def buildODE():
|
||||
def ODE(s, v):
|
||||
def ODE(s: float, v: NDArray):
|
||||
phi, lam, A = v
|
||||
V = re.V(phi)
|
||||
dphi = cos(A) * V ** 3 / re.c
|
||||
0
nicht abgeben/Tests/__init__.py
Normal file
0
nicht abgeben/Tests/__init__.py
Normal file
File diff suppressed because it is too large
Load Diff
@@ -17,10 +17,11 @@
|
||||
"source": [
|
||||
"%reload_ext autoreload\n",
|
||||
"%autoreload 2\n",
|
||||
"import numpy as np\n",
|
||||
"\n",
|
||||
"import winkelumrechnungen as wu\n",
|
||||
"from ellipsoide import EllipsoidTriaxial\n",
|
||||
"from GHA_triaxial.utils import alpha_para2ell, alpha_ell2para\n",
|
||||
"import numpy as np"
|
||||
"from GHA_triaxial.utils import alpha_ell2para, alpha_para2ell\n",
|
||||
"from ellipsoid_triaxial import EllipsoidTriaxial"
|
||||
],
|
||||
"id": "46aa84a937fea491",
|
||||
"outputs": [],
|
||||
@@ -20,13 +20,14 @@
|
||||
"source": [
|
||||
"%reload_ext autoreload\n",
|
||||
"%autoreload 2\n",
|
||||
"import pickle\n",
|
||||
"import numpy as np\n",
|
||||
"import winkelumrechnungen as wu\n",
|
||||
"from itertools import product\n",
|
||||
"\n",
|
||||
"import numpy as np\n",
|
||||
"import pandas as pd\n",
|
||||
"from ellipsoide import EllipsoidTriaxial\n",
|
||||
"import plotly.graph_objects as go"
|
||||
"import plotly.graph_objects as go\n",
|
||||
"\n",
|
||||
"import winkelumrechnungen as wu\n",
|
||||
"from ellipsoid_triaxial import EllipsoidTriaxial"
|
||||
],
|
||||
"outputs": [],
|
||||
"execution_count": null
|
||||
BIN
nicht abgeben/Tests/gha_resultsKarney.pkl
Normal file
BIN
nicht abgeben/Tests/gha_resultsKarney.pkl
Normal file
Binary file not shown.
BIN
nicht abgeben/Tests/gha_resultsRandom_num.pkl
Normal file
BIN
nicht abgeben/Tests/gha_resultsRandom_num.pkl
Normal file
Binary file not shown.
0
nicht abgeben/__init__.py
Normal file
0
nicht abgeben/__init__.py
Normal file
126
nicht abgeben/ellipsoid_biaxial.py
Normal file
126
nicht abgeben/ellipsoid_biaxial.py
Normal file
@@ -0,0 +1,126 @@
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import arctan2, cos, sin, sqrt
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
|
||||
class EllipsoidBiaxial:
|
||||
"""
|
||||
Klasse für Rotationsellipdoide
|
||||
"""
|
||||
def __init__(self, a: float, b: float):
|
||||
self.a = a
|
||||
self.b = b
|
||||
self.c = a ** 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) -> EllipsoidBiaxial:
|
||||
"""
|
||||
Erstellen eines Rotationsellipdoids nach Namen
|
||||
:param name: Name des Rotationsellipsoids
|
||||
:return: Rotationsellipsoid
|
||||
"""
|
||||
if name == "Bessel":
|
||||
a = 6377397.15508
|
||||
b = 6356078.96290
|
||||
return cls(a, b)
|
||||
elif name == "Hayford":
|
||||
a = 6378388
|
||||
f = 1/297
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
elif name == "Krassowski":
|
||||
a = 6378245
|
||||
f = 298.3
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
elif name == "WGS84":
|
||||
a = 6378137
|
||||
f = 298.257223563
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
else:
|
||||
raise Exception(f"EllipsoidBiaxial.init_name: Name {name} unbekannt")
|
||||
|
||||
@classmethod
|
||||
def init_af(cls, a: float, f: float) -> EllipsoidBiaxial:
|
||||
"""
|
||||
Erstellen eines Rotationsellipdoids aus der großen Halbachse und der Abplattung
|
||||
:param a: große Halbachse
|
||||
:param f: großen Halbachse
|
||||
:return: Rotationsellipsoid
|
||||
"""
|
||||
b = a - a * f
|
||||
return cls(a, b)
|
||||
|
||||
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: arctan2(self.a * sin(beta), self.b * cos(beta))
|
||||
beta2phi = lambda self, beta: arctan2(self.a ** 2 * sin(beta), self.b ** 2 * cos(beta))
|
||||
|
||||
psi2beta = lambda self, psi: arctan2(self.b * sin(psi), self.a * cos(psi))
|
||||
psi2phi = lambda self, psi: arctan2(self.a * sin(psi), self.b * cos(psi))
|
||||
|
||||
phi2beta = lambda self, phi: arctan2(self.b**2 * sin(phi), self.a**2 * cos(phi))
|
||||
phi2psi = lambda self, phi: arctan2(self.b * sin(phi), self.a * cos(phi))
|
||||
|
||||
phi2p = lambda self, phi: self.N(phi) * cos(phi)
|
||||
|
||||
def bi_cart2ell(self, point: NDArray, Eh: float = 0.001, Ephi: float = wu.gms2rad([0, 0, 0.001])) -> Tuple[float, float, float]:
|
||||
"""
|
||||
Umrechnung von kartesischen in ellipsoidische Koordinaten auf einem Rotationsellipsoid
|
||||
# TODO: Quelle
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:param Eh: Grenzwert für die Höhe
|
||||
:param Ephi: Grenzwert für die Breite
|
||||
:return: ellipsoidische Breite, Länge, geodätische Höhe
|
||||
"""
|
||||
x, y, z = point
|
||||
|
||||
lamb = arctan2(y, x)
|
||||
|
||||
p = sqrt(x**2+y**2)
|
||||
|
||||
phi_null = arctan2(z, p*(1 - self.e**2))
|
||||
|
||||
hi = [0]
|
||||
phii = [phi_null]
|
||||
|
||||
i = 0
|
||||
|
||||
while True:
|
||||
N = self.a / sqrt(1 - self.e**2 * sin(phii[i])**2)
|
||||
h = p / cos(phii[i]) - N
|
||||
phi = arctan2(z, p * (1-(self.e**2*N) / (N+h)))
|
||||
hi.append(h)
|
||||
phii.append(phi)
|
||||
dh = abs(hi[i]-h)
|
||||
dphi = abs(phii[i]-phi)
|
||||
i += 1
|
||||
if dh < Eh:
|
||||
if dphi < Ephi:
|
||||
break
|
||||
return phi, lamb, h
|
||||
|
||||
def bi_ell2cart(self, phi: float, lamb: float, h: float) -> NDArray:
|
||||
"""
|
||||
Umrechnung von ellipsoidischen in kartesische Koordinaten auf einem Rotationsellipsoid
|
||||
# TODO: Quelle
|
||||
:param phi: ellipsoidische Breite
|
||||
:param lamb: ellipsoidische Länge
|
||||
:param h: geodätische Höhe
|
||||
:return: Punkt in kartesischen Koordinaten
|
||||
"""
|
||||
W = sqrt(1 - self.e**2 * sin(phi)**2)
|
||||
N = self.a / W
|
||||
x = (N+h) * cos(phi) * cos(lamb)
|
||||
y = (N+h) * cos(phi) * sin(lamb)
|
||||
z = (N * (1-self.e**2) + h) * sin(phi)
|
||||
return np.array([x, y, z])
|
||||
@@ -1,7 +1,9 @@
|
||||
import numpy as np
|
||||
from numpy import sqrt, arctan2, sin, cos, arcsin, arccos
|
||||
from numpy.typing import NDArray
|
||||
from typing import Tuple
|
||||
|
||||
import numpy as np
|
||||
from numpy import arccos, arcsin, arctan2, cos, pi, sin, sqrt
|
||||
from numpy.typing import NDArray
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
|
||||
@@ -77,7 +79,7 @@ def gha2(R: float, phi0: float, lamb0: float, phi1: float, lamb1: float) -> Tupl
|
||||
alpha1 = arctan2(-cos(phi0) * sin(lamb1 - lamb0),
|
||||
cos(phi1) * sin(phi0) - sin(phi1) * cos(phi0) * cos(lamb1 - lamb0))
|
||||
if alpha1 < 0:
|
||||
alpha1 += 2 * np.pi
|
||||
alpha1 += 2 * pi
|
||||
|
||||
return alpha0, alpha1, s
|
||||
|
||||
@@ -9,10 +9,11 @@
|
||||
},
|
||||
"cell_type": "code",
|
||||
"source": [
|
||||
"import plotly.graph_objects as go\n",
|
||||
"import numpy as np\n",
|
||||
"from ellipsoide import EllipsoidTriaxial\n",
|
||||
"import winkelumrechnungen as wu"
|
||||
"import plotly.graph_objects as go\n",
|
||||
"\n",
|
||||
"import winkelumrechnungen as wu\n",
|
||||
"from ellipsoid_triaxial import EllipsoidTriaxial"
|
||||
],
|
||||
"id": "731173e4745cfe7c",
|
||||
"outputs": [],
|
||||
8
nicht abgeben/test.py
Normal file
8
nicht abgeben/test.py
Normal file
@@ -0,0 +1,8 @@
|
||||
import numpy as np
|
||||
|
||||
import ellipsoid_triaxial
|
||||
|
||||
ell = ellipsoid_triaxial.EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
|
||||
cart = ell.para2cart(0, np.pi/2)
|
||||
print(cart)
|
||||
7
requirements.txt
Normal file
7
requirements.txt
Normal file
@@ -0,0 +1,7 @@
|
||||
numpy~=2.3.4
|
||||
plotly~=6.4.0
|
||||
pandas~=2.3.3
|
||||
scipy~=1.16.3
|
||||
dash-bootstrap-components~=2.0.4
|
||||
dash~=4.0.0
|
||||
matplotlib~=3.10.7
|
||||
@@ -1,7 +1,10 @@
|
||||
from typing import Callable
|
||||
|
||||
import numpy as np
|
||||
from numpy.typing import NDArray
|
||||
|
||||
|
||||
def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]:
|
||||
def rk4(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, fein: bool = False) -> tuple[list, list]:
|
||||
"""
|
||||
Standard Runge-Kutta Verfahren 4. Ordnung
|
||||
:param ode: ODE-System als Funktion
|
||||
@@ -9,7 +12,7 @@ def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool
|
||||
:param v0: Startwerte
|
||||
:param weite: Integrationsweite
|
||||
:param schritte: Schrittzahl
|
||||
:param fein:
|
||||
:param fein: Fein-Rechnung?
|
||||
:return: Variable und Funktionswerte an jedem Stützpunkt
|
||||
"""
|
||||
h = weite/schritte
|
||||
@@ -35,14 +38,32 @@ def rk4(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool
|
||||
|
||||
return t_list, werte
|
||||
|
||||
def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray:
|
||||
def rk4_step(ode: Callable, t: float, v: NDArray, h: float) -> NDArray:
|
||||
"""
|
||||
Ein Schritt des Runge-Kutta Verfahrens 4. Ordnung
|
||||
:param ode: ODE-System als Funktion
|
||||
:param t: unabhängige Variable
|
||||
:param v: abhängige Variablen
|
||||
:param h: Schrittweite
|
||||
:return: abhängige Variablen nach einem Schritt
|
||||
"""
|
||||
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)
|
||||
|
||||
def rk4_end(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False):
|
||||
def rk4_end(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, fein: bool = False):
|
||||
"""
|
||||
Standard Runge-Kutta Verfahren 4. Ordnung, nur Ausgabe der letzten Variablenwerte
|
||||
:param ode: ODE-System als Funktion
|
||||
:param t0: Startwert der unabhängigen Variable
|
||||
:param v0: Startwerte
|
||||
:param weite: Integrationsweite
|
||||
:param schritte: Schrittzahl
|
||||
:param fein: Fein-Rechnung?
|
||||
:return: Variable und Funktionswerte am letzten Stützpunkt
|
||||
"""
|
||||
h = weite / schritte
|
||||
t = float(t0)
|
||||
v = np.array(v0, dtype=float, copy=True)
|
||||
@@ -62,8 +83,19 @@ def rk4_end(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: b
|
||||
return t, v
|
||||
|
||||
# RK4 mit Simpson bzw. Trapez
|
||||
def rk4_integral( ode, t0: float, v0: np.ndarray, weite: float, schritte: int, integrand_at, fein: bool = False, simpson: bool = True, ):
|
||||
|
||||
def rk4_integral(ode: Callable, t0: float, v0: NDArray, weite: float, schritte: int, integrand_at: Callable, fein: bool = False, simpson: bool = True):
|
||||
"""
|
||||
Runge-Kutta Verfahren 4. Ordnung mit Simpson bzw. Trapez
|
||||
:param ode: ODE-System als Funktion
|
||||
:param t0: Startwert der unabhängigen Variable
|
||||
:param v0: Startwerte
|
||||
:param weite: Integrationsweite
|
||||
:param integrand_at: Funktion
|
||||
:param schritte: Schrittzahl
|
||||
:param fein: Fein-Rechnung?
|
||||
:param simpson: Simpson? Wenn nein, dann Trapez
|
||||
:return: Variable und Funktionswerte am letzten Stützpunkt
|
||||
"""
|
||||
h = weite / schritte
|
||||
habs = abs(h)
|
||||
|
||||
|
||||
7
test.py
7
test.py
@@ -1,7 +0,0 @@
|
||||
import numpy as np
|
||||
import ellipsoide
|
||||
|
||||
ell = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
|
||||
|
||||
cart = ell.para2cart(0, np.pi/2)
|
||||
print(cart)
|
||||
@@ -1,13 +1,54 @@
|
||||
import numpy as np
|
||||
|
||||
import winkelumrechnungen as wu
|
||||
|
||||
def arccot(x):
|
||||
|
||||
def arccot(x: float) -> float:
|
||||
"""
|
||||
Berechnung von arccot eines Winkels
|
||||
:param x: Winkel
|
||||
:return: arccot(Winkel)
|
||||
"""
|
||||
return np.arctan2(1.0, x)
|
||||
|
||||
|
||||
def cot(a):
|
||||
return np.cos(a) / np.sin(a)
|
||||
def cot(x: float) -> float:
|
||||
"""
|
||||
Berechnung von cot eines Winkels
|
||||
:param x: Winkel
|
||||
:return: cot(Winkel)
|
||||
"""
|
||||
return np.cos(x) / np.sin(x)
|
||||
|
||||
|
||||
def wrap_to_pi(x):
|
||||
def wrap_mpi_pi(x: float) -> float:
|
||||
"""
|
||||
Wrap eines Winkels in den Wertebereich [-π, π)
|
||||
:param x: Winkel
|
||||
:return: Winkel in [-π, π)
|
||||
"""
|
||||
return (x + np.pi) % (2 * np.pi) - np.pi
|
||||
|
||||
|
||||
def wrap_mhalfpi_halfpi(x: float) -> float:
|
||||
"""
|
||||
Wrap eines Winkels in den Wertebereich [-π/2, π/2)
|
||||
:param x: Winkel
|
||||
:return: Winkel in [-π/2, π/2)
|
||||
"""
|
||||
return (x + np.pi / 2) % np.pi - np.pi / 2
|
||||
|
||||
|
||||
def wrap_0_2pi(x: float) -> float:
|
||||
"""
|
||||
Wrap eines Winkels in den Wertebereich [0, 2π)
|
||||
:param x: Winkel
|
||||
:return: Winkel in [0, 2π)
|
||||
"""
|
||||
return x % (2 * np.pi)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
print(wu.rad2deg(wrap_mhalfpi_halfpi(wu.deg2rad(181))))
|
||||
print(wu.rad2deg(wrap_0_2pi(wu.deg2rad(181))))
|
||||
print(wu.rad2deg(wrap_mpi_pi(wu.deg2rad(181))))
|
||||
|
||||
Reference in New Issue
Block a user