549 lines
20 KiB
Python
549 lines
20 KiB
Python
import numpy as np
|
|
import winkelumrechnungen as wu
|
|
import ausgaben as aus
|
|
import jacobian_Ligas
|
|
|
|
|
|
class EllipsoidBiaxial:
|
|
def __init__(self, a: float, b: float):
|
|
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
|
|
|
|
@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: np.sqrt(1 + self.e_ ** 2 * np.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))
|
|
|
|
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))
|
|
|
|
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))
|
|
|
|
phi2p = lambda self, phi: self.N(phi) * np.cos(phi)
|
|
|
|
def cart2ell(self, Eh, Ephi, x, y, z):
|
|
p = np.sqrt(x**2+y**2)
|
|
# print(f"p = {round(p, 5)} m")
|
|
|
|
lamb = np.arctan(y/x)
|
|
|
|
phi_null = np.arctan(z/p*(1-self.e**2)**-1)
|
|
|
|
hi = [0]
|
|
phii = [phi_null]
|
|
|
|
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))
|
|
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
|
|
for i in range(len(phii)):
|
|
# print(f"P3[{i}]: {aus.gms('phi', phii[i], 5)}\th = {round(hi[i], 5)} m")
|
|
pass
|
|
return phi, lamb, h
|
|
|
|
def ell2cart(self, phi, lamb, h):
|
|
W = np.sqrt(1 - self.e**2 * np.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)
|
|
return x, y, z
|
|
|
|
class EllipsoidTriaxial:
|
|
def __init__(self, ax: float, ay: float, b: float):
|
|
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)
|
|
|
|
@classmethod
|
|
def init_name(cls, name: str):
|
|
"""
|
|
Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972, Bursa1970, BesselBiaxial
|
|
:param name: Name des dreiachsigen Ellipsoids
|
|
"""
|
|
if name == "BursaFialova1993":
|
|
ax = 6378171.36
|
|
ay = 6378101.61
|
|
b = 6356751.84
|
|
return cls(ax, ay, b)
|
|
elif name == "BursaSima1980":
|
|
ax = 6378172
|
|
ay = 6378102.7
|
|
b = 6356752.6
|
|
return cls(ax, ay, b)
|
|
elif name == "BursaSima1980round":
|
|
# Panou 2013
|
|
ax = 6378172
|
|
ay = 6378103
|
|
b = 6356753
|
|
return cls(ax, ay, b)
|
|
elif name == "Eitschberger1978":
|
|
ax = 6378173.435
|
|
ay = 6378103.9
|
|
b = 6356754.4
|
|
return cls(ax, ay, b)
|
|
elif name == "Bursa1972":
|
|
ax = 6378173
|
|
ay = 6378104
|
|
b = 6356754
|
|
return cls(ax, ay, b)
|
|
elif name == "Bursa1970":
|
|
ax = 6378173
|
|
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 = 5500000
|
|
ay = 4500000
|
|
b = 4000000
|
|
return cls(ax, ay, b)
|
|
|
|
def point_on(self, point: np.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) < 0.000001:
|
|
return True
|
|
else:
|
|
return False
|
|
|
|
def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
|
|
"""
|
|
Panou 2014 12ff.
|
|
Ellipsoidische Breite+Länge sind nicht gleich der geodätischen
|
|
Verhältnisse des Ellipsoids bekannt, Größe verändern bis Punkt erreicht,
|
|
dann ist u die Größe entlang der z-Achse
|
|
:param beta: ellipsoidische Breite [rad]
|
|
:param lamb: ellipsoidische Länge [rad]
|
|
: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)
|
|
|
|
return np.array([x, y, z])
|
|
|
|
def ell2cart(self, beta: float, lamb: float) -> np.ndarray:
|
|
"""
|
|
Panou, Korakitis 2019 2
|
|
:param beta: ellipsoidische Breite [rad]
|
|
:param lamb: ellipsoidische Länge [rad]
|
|
:return: Punkt in kartesischen Koordinaten
|
|
"""
|
|
if beta == -np.pi/2:
|
|
return np.array([0, 0, -self.b])
|
|
elif beta == np.pi/2:
|
|
return np.array([0, 0, self.b])
|
|
elif beta == 0 and lamb == -np.pi/2:
|
|
return np.array([0, -self.ay, 0])
|
|
elif beta == 0 and lamb == np.pi/2:
|
|
return np.array([0, self.ay, 0])
|
|
elif beta == 0 and lamb == 0:
|
|
return np.array([self.ax, 0, 0])
|
|
elif beta == 0 and lamb == np.pi:
|
|
return np.array([-self.ax, 0, 0])
|
|
else:
|
|
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)
|
|
return np.array([x, y, z])
|
|
|
|
def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]:
|
|
"""
|
|
Panou 2014 15ff.
|
|
:param point: Punkt in kartesischen Koordinaten
|
|
:return: ellipsoidische Breite, ellipsoidische Länge, Größe entlang der z-Achse
|
|
"""
|
|
x, y, z = point
|
|
c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2
|
|
c1 = (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)
|
|
c0 = (self.ax**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)
|
|
|
|
p = (c2**2 - 3*c1) / 9
|
|
q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54
|
|
omega = np.arccos(q / np.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
|
|
# print(s1, s2, s3)
|
|
|
|
beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
|
|
lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
|
|
u = np.sqrt(self.b**2 + s1)
|
|
|
|
return beta, lamb, u
|
|
|
|
def cart2ell(self, point: np.ndarray) -> tuple[float, float]:
|
|
"""
|
|
Panou, Korakitis 2019 2f.
|
|
:param point: Punkt in kartesischen Koordinaten
|
|
:return: ellipsoidische Breite, ellipsoidische Länge
|
|
"""
|
|
x, y, z = point
|
|
|
|
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
|
|
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
|
|
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
|
|
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
|
|
|
|
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))
|
|
|
|
if z < 0:
|
|
beta = -beta
|
|
|
|
if y < 0:
|
|
lamb = -lamb
|
|
|
|
if x < 0:
|
|
lamb = np.pi - lamb
|
|
|
|
if abs(x) < eps:
|
|
lamb = -np.pi/2 if y < 0 else np.pi/2
|
|
|
|
elif abs(y) < eps:
|
|
lamb = 0 if x > 0 else np.pi
|
|
|
|
elif abs(z) < eps:
|
|
beta = 0
|
|
|
|
return beta, lamb
|
|
|
|
def cart2geod(self, mode: str, point: np.ndarray, maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
|
|
"""
|
|
Ligas 2012
|
|
:param mode: ligas1, ligas2, oder ligas3
|
|
:param point: Punkt in kartesischen Koordinaten
|
|
:param maxIter: maximale Anzahl Iterationen
|
|
:param maxLoa: Level of Accuracy, das erreicht werden soll
|
|
:return: phi, lambda, h
|
|
"""
|
|
xG, yG, zG = point
|
|
|
|
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
|
|
lamb = 0.0
|
|
h = abs(zG) - ell.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
|
|
h = abs(yG) - ell.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
|
|
h = abs(xG) - ell.ax
|
|
return phi, lamb, h
|
|
|
|
# elif abs(zG) < eps: # Punkt in der xy-Ebene
|
|
# phi = 0
|
|
# lamb = np.arctan2(yG / ell.ay**2, xG / ell.ax**2)
|
|
# rG = np.sqrt(xG ** 2 + yG ** 2)
|
|
# pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
|
|
# rE = np.sqrt(pE[0] ** 2 + pE[1] ** 2)
|
|
# h = rG - rE
|
|
# return phi, lamb, h
|
|
#
|
|
# elif abs(yG) < eps: # Punkt in der xz-Ebene
|
|
# phi = np.arctan2(zG / ell.b**2, xG / ell.ax**2)
|
|
# lamb = 0 if xG > 0 else np.pi
|
|
# rG = np.sqrt(xG ** 2 + zG ** 2)
|
|
# pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
|
|
# rE = np.sqrt(pE[0] ** 2 + pE[2] ** 2)
|
|
# h = rG - rE
|
|
# return phi, lamb, h
|
|
|
|
rG = np.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
|
|
F = 1 / self.ay**2
|
|
G = 1 / self.b**2
|
|
|
|
i = 0
|
|
loa = np.inf
|
|
|
|
while i < maxIter and loa > maxLoa:
|
|
if mode == "ligas1":
|
|
invJ, fxE = jacobian_Ligas.case1(E, F, G, np.array([xG, yG, zG]), pE)
|
|
elif mode == "ligas2":
|
|
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)
|
|
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)
|
|
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)
|
|
|
|
if xG < 0 and yG < 0:
|
|
lamb = -np.pi + lamb
|
|
|
|
elif xG < 0:
|
|
lamb = np.pi + lamb
|
|
|
|
if abs(zG) < eps:
|
|
phi = 0
|
|
|
|
return phi, lamb, h
|
|
|
|
def geod2cart(self, phi: float, lamb: float, h: float) -> np.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 / 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])
|
|
|
|
def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]:
|
|
"""
|
|
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
|
|
"""
|
|
phi, lamb, h = self.cart2geod("ligas3", point)
|
|
x, y, z = self. geod2cart(phi, lamb, 0)
|
|
return np.array([x, y, z]), phi, lamb, h
|
|
|
|
def cartellh(self, point: np.ndarray, h: float) -> np.ndarray:
|
|
"""
|
|
Punkt auf Ellipsoid hoch loten
|
|
:param point: Punkt auf dem Ellipsoid
|
|
:param h: Höhe über dem Ellipsoid
|
|
:return: hochgeloteter Punkt
|
|
"""
|
|
phi, lamb, _ = self.cart2geod("ligas3", point)
|
|
pointH = self. geod2cart(phi, lamb, h)
|
|
return pointH
|
|
|
|
def para2cart(self, u: float, v: float) -> 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)
|
|
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 p_q(self, x, y, z) -> dict:
|
|
"""
|
|
Berechnung sämtlicher Größen
|
|
:param x: x
|
|
:param y: y
|
|
:param z: z
|
|
:return: Dictionary sämtlicher Größen
|
|
"""
|
|
H = x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2
|
|
|
|
n = np.array([x / np.sqrt(H), y / ((1 - self.ee ** 2) * np.sqrt(H)), z / ((1 - self.ex ** 2) * np.sqrt(H))])
|
|
|
|
beta, lamb, u = self.cart2ellu(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[1] * p[1] - n[1] * p[0]])
|
|
|
|
return {"H": H, "n": n, "beta": beta, "lamb": lamb, "u": u, "B": B, "L": L, "c1": c1, "c0": c0, "t1": t1,
|
|
"t2": t2,
|
|
"F": F, "p": p, "q": q}
|
|
|
|
|
|
if __name__ == "__main__":
|
|
ell = EllipsoidTriaxial.init_name("Eitschberger1978")
|
|
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)
|
|
|
|
elli = ell.cart2ell(point)
|
|
cart_elli = ell.ell2cart(elli[0], elli[1])
|
|
diff_ell = np.sum(np.abs(point-cart_elli))
|
|
|
|
para = ell.cart2para(point)
|
|
cart_para = ell.para2cart(para[0], para[1])
|
|
diff_para = np.sum(np.abs(point-cart_para))
|
|
|
|
geod = ell.cart2geod("ligas1", point)
|
|
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
|
diff_geod1 = np.sum(np.abs(point-cart_geod))
|
|
|
|
geod = ell.cart2geod("ligas2", point)
|
|
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
|
diff_geod2 = np.sum(np.abs(point-cart_geod))
|
|
|
|
geod = ell.cart2geod("ligas3", point)
|
|
cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
|
|
diff_geod3 = np.sum(np.abs(point-cart_geod))
|
|
|
|
diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod1, diff_geod2, diff_geod3])
|
|
|
|
diff_list = np.array(diff_list)
|
|
pass |