Merge remote-tracking branch 'origin/main'
This commit is contained in:
Tammo.Weber
@@ -12,6 +12,33 @@ from numpy.typing import NDArray
|
||||
from GHA_triaxial.utils import alpha_ell2para, pq_ell
|
||||
|
||||
|
||||
def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
"""
|
||||
Aufbau des DGL-Systems
|
||||
:param ell: Ellipsoid
|
||||
:return: DGL-System
|
||||
"""
|
||||
|
||||
def ODE(s: float, v: NDArray) -> NDArray:
|
||||
"""
|
||||
DGL-System
|
||||
:param s: unabhängige Variable
|
||||
:param v: abhängige Variablen
|
||||
:return: Ableitungen der abhängigen Variablen
|
||||
"""
|
||||
x, dxds, y, dyds, z, dzds = v
|
||||
|
||||
H = ell.func_H(np.array([x, y, z]))
|
||||
h = dxds ** 2 + 1 / (1 - ell.ee ** 2) * dyds ** 2 + 1 / (1 - ell.ex ** 2) * dzds ** 2
|
||||
|
||||
ddx = -(h / H) * x
|
||||
ddy = -(h / H) * y / (1 - ell.ee ** 2)
|
||||
ddz = -(h / H) * z / (1 - ell.ex ** 2)
|
||||
|
||||
return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
|
||||
|
||||
return ODE
|
||||
|
||||
def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, num: int, all_points: bool = False) -> Tuple[NDArray, float] | Tuple[NDArray, float, List]:
|
||||
"""
|
||||
Panou, Korakitits 2019
|
||||
@@ -34,33 +61,6 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu
|
||||
|
||||
v_init = np.array([x0, dxds0, y0, dyds0, z0, dzds0])
|
||||
|
||||
def buildODE(ell: EllipsoidTriaxial) -> Callable:
|
||||
"""
|
||||
Aufbau des DGL-Systems
|
||||
:param ell: Ellipsoid
|
||||
:return: DGL-System
|
||||
"""
|
||||
|
||||
def ODE(s: float, v: NDArray) -> NDArray:
|
||||
"""
|
||||
DGL-System
|
||||
:param s: unabhängige Variable
|
||||
:param v: abhängige Variablen
|
||||
:return: Ableitungen der abhängigen Variablen
|
||||
"""
|
||||
x, dxds, y, dyds, z, dzds = v
|
||||
|
||||
H = ell.func_H(np.array([x, y, z]))
|
||||
h = dxds ** 2 + 1 / (1 - ell.ee ** 2) * dyds ** 2 + 1 / (1 - ell.ex ** 2) * dzds ** 2
|
||||
|
||||
ddx = -(h / H) * x
|
||||
ddy = -(h / H) * y / (1 - ell.ee ** 2)
|
||||
ddz = -(h / H) * z / (1 - ell.ex ** 2)
|
||||
|
||||
return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
|
||||
|
||||
return ODE
|
||||
|
||||
ode = buildODE(ell)
|
||||
|
||||
_, werte = rk.rk4(ode, 0, v_init, s, num)
|
||||
|
||||
@@ -1,12 +1,141 @@
|
||||
import numpy as np
|
||||
from ellipsoide import EllipsoidTriaxial
|
||||
import runge_kutta as rk
|
||||
import GHA_triaxial.numeric_examples_karney as ne_karney
|
||||
import GHA_triaxial.numeric_examples_panou as ne_panou
|
||||
import winkelumrechnungen as wu
|
||||
from typing import Tuple
|
||||
from numpy.typing import NDArray
|
||||
|
||||
from utils_angle import arccot, cot, wrap_to_pi
|
||||
|
||||
|
||||
def sph_azimuth(beta1, lam1, beta2, lam2):
|
||||
# sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß)
|
||||
dlam = wrap_to_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) # (-pi, pi]
|
||||
if a < 0:
|
||||
a += 2 * np.pi
|
||||
return a
|
||||
|
||||
def BETA_LAMBDA(beta, lamb):
|
||||
|
||||
BETA = (ell.ay**2 * np.sin(beta)**2 + ell.b**2 * np.cos(beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)
|
||||
LAMBDA = (ell.ax**2 * np.sin(lamb)**2 + ell.ay**2 * np.cos(lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)
|
||||
|
||||
# Erste Ableitungen von ΒETA und LAMBDA
|
||||
BETA_ = (ell.ax**2 * ell.Ey**2 * np.sin(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2
|
||||
LAMBDA_ = - (ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2
|
||||
|
||||
# Zweite Ableitungen von ΒETA und LAMBDA
|
||||
BETA__ = ((2 * ell.ax**2 * ell.Ey**4 * np.sin(2*beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**3) + ((2 * ell.ax**2 * ell.Ey**2 * np.cos(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2)
|
||||
LAMBDA__ = (((2 * ell.b**2 * ell.Ee**4 * np.sin(2*lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**3) -
|
||||
((2 * ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2))
|
||||
|
||||
E = BETA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
|
||||
F = 0
|
||||
G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
|
||||
|
||||
# Erste Ableitungen von E und G
|
||||
E_beta = BETA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - BETA * ell.Ey**2 * np.sin(2*beta)
|
||||
E_lamb = BETA * ell.Ee**2 * np.sin(2*lamb)
|
||||
|
||||
G_beta = - LAMBDA * ell.Ey**2 * np.sin(2*beta)
|
||||
G_lamb = LAMBDA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + LAMBDA * ell.Ee**2 * np.sin(2*lamb)
|
||||
|
||||
# Zweite Ableitungen von E und G
|
||||
E_beta_beta = BETA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - 2 * BETA_ * ell.Ey**2 * np.sin(2*beta) - 2 * BETA * ell.Ey**2 * np.cos(2*beta)
|
||||
E_beta_lamb = BETA_ * ell.Ee**2 * np.sin(2*lamb)
|
||||
E_lamb_lamb = 2 * BETA * ell.Ee**2 * np.cos(2*lamb)
|
||||
|
||||
G_beta_beta = - 2 * LAMBDA * ell.Ey**2 * np.cos(2*beta)
|
||||
G_beta_lamb = - LAMBDA_ * ell.Ey**2 * np.sin(2*beta)
|
||||
G_lamb_lamb = LAMBDA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + 2 * LAMBDA_ * ell.Ee**2 * np.sin(2*lamb) + 2 * LAMBDA * ell.Ee**2 * np.cos(2*lamb)
|
||||
|
||||
return (BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb)
|
||||
|
||||
def p_coef(beta, lamb):
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
|
||||
|
||||
p_3 = - 0.5 * (E_lamb / G)
|
||||
p_2 = (G_beta / G) - 0.5 * (E_beta / E)
|
||||
p_1 = 0.5 * (G_lamb / G) - (E_lamb / E)
|
||||
p_0 = 0.5 * (G_beta / E)
|
||||
|
||||
p_33 = - 0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G ** 2))
|
||||
p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G ** 2)) - 0.5 * ((E * E_beta_beta - E_beta * E_beta) / (E ** 2))
|
||||
p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2))
|
||||
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)
|
||||
|
||||
def buildODElamb():
|
||||
def ODE(lamb, v):
|
||||
beta, beta_p, X3, X4 = v
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
p_3, p_2, p_1, p_0,
|
||||
p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
|
||||
dbeta = beta_p
|
||||
dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
dX3 = X4
|
||||
dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
(3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
|
||||
return np.array([dbeta, dbeta_p, dX3, dX4])
|
||||
|
||||
return ODE
|
||||
|
||||
def q_coef(beta, lamb):
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
|
||||
|
||||
q_3 = - 0.5 * (G_beta / E)
|
||||
q_2 = (E_lamb / E) - 0.5 * (G_lamb / G)
|
||||
q_1 = 0.5 * (E_beta / E) - (G_beta / G)
|
||||
q_0 = 0.5 * (E_lamb / G)
|
||||
|
||||
q_33 = - 0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E ** 2))
|
||||
q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E ** 2)) - 0.5 * ((G * G_lamb_lamb - G_lamb * G_lamb) / (G ** 2))
|
||||
q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2))
|
||||
q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G ** 2))
|
||||
|
||||
return (BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00)
|
||||
|
||||
def buildODEbeta():
|
||||
def ODE(beta, v):
|
||||
lamb, lamb_p, Y3, Y4 = v
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
|
||||
dlamb = lamb_p
|
||||
dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
dY3 = Y4
|
||||
dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
(3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
|
||||
|
||||
return np.array([dlamb, dlamb_p, dY3, dY4])
|
||||
|
||||
return ODE
|
||||
|
||||
# Panou 2013
|
||||
def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float, lamb_2: float,
|
||||
n: int = 16000, epsilon: float = 10**-12, iter_max: int = 30, all_points: bool = False
|
||||
@@ -27,152 +156,8 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
|
||||
|
||||
# h_x, h_y, h_e entsprechen E_x, E_y, E_e
|
||||
|
||||
def arccot(x):
|
||||
return np.arctan2(1.0, x)
|
||||
|
||||
def cot(a):
|
||||
return np.cos(a) / np.sin(a)
|
||||
|
||||
def wrap_to_pi(x):
|
||||
return (x + np.pi) % (2 * np.pi) - np.pi
|
||||
|
||||
def sph_azimuth(beta1, lam1, beta2, lam2):
|
||||
# sphärischer Anfangsazimut (von Norden/meridian, im Bogenmaß)
|
||||
dlam = wrap_to_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) # (-pi, pi]
|
||||
if a < 0:
|
||||
a += 2 * np.pi
|
||||
return a
|
||||
|
||||
def BETA_LAMBDA(beta, lamb):
|
||||
|
||||
BETA = (ell.ay**2 * np.sin(beta)**2 + ell.b**2 * np.cos(beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)
|
||||
LAMBDA = (ell.ax**2 * np.sin(lamb)**2 + ell.ay**2 * np.cos(lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)
|
||||
|
||||
# Erste Ableitungen von ΒETA und LAMBDA
|
||||
BETA_ = (ell.ax**2 * ell.Ey**2 * np.sin(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2
|
||||
LAMBDA_ = - (ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2
|
||||
|
||||
# Zweite Ableitungen von ΒETA und LAMBDA
|
||||
BETA__ = ((2 * ell.ax**2 * ell.Ey**4 * np.sin(2*beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**3) + ((2 * ell.ax**2 * ell.Ey**2 * np.cos(2*beta)) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)**2)
|
||||
LAMBDA__ = (((2 * ell.b**2 * ell.Ee**4 * np.sin(2*lamb)**2) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**3) -
|
||||
((2 * ell.b**2 * ell.Ee**2 * np.sin(2*lamb)) / (ell.Ex**2 - ell.Ee**2 * np.cos(lamb)**2)**2))
|
||||
|
||||
E = BETA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
|
||||
F = 0
|
||||
G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2)
|
||||
|
||||
# Erste Ableitungen von E und G
|
||||
E_beta = BETA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - BETA * ell.Ey**2 * np.sin(2*beta)
|
||||
E_lamb = BETA * ell.Ee**2 * np.sin(2*lamb)
|
||||
|
||||
G_beta = - LAMBDA * ell.Ey**2 * np.sin(2*beta)
|
||||
G_lamb = LAMBDA_ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + LAMBDA * ell.Ee**2 * np.sin(2*lamb)
|
||||
|
||||
# Zweite Ableitungen von E und G
|
||||
E_beta_beta = BETA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) - 2 * BETA_ * ell.Ey**2 * np.sin(2*beta) - 2 * BETA * ell.Ey**2 * np.cos(2*beta)
|
||||
E_beta_lamb = BETA_ * ell.Ee**2 * np.sin(2*lamb)
|
||||
E_lamb_lamb = 2 * BETA * ell.Ee**2 * np.cos(2*lamb)
|
||||
|
||||
G_beta_beta = - 2 * LAMBDA * ell.Ey**2 * np.cos(2*beta)
|
||||
G_beta_lamb = - LAMBDA_ * ell.Ey**2 * np.sin(2*beta)
|
||||
G_lamb_lamb = LAMBDA__ * (ell.Ey**2 * np.cos(beta)**2 + ell.Ee**2 * np.sin(lamb)**2) + 2 * LAMBDA_ * ell.Ee**2 * np.sin(2*lamb) + 2 * LAMBDA * ell.Ee**2 * np.cos(2*lamb)
|
||||
|
||||
return (BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb)
|
||||
|
||||
def p_coef(beta, lamb):
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
|
||||
|
||||
p_3 = - 0.5 * (E_lamb / G)
|
||||
p_2 = (G_beta / G) - 0.5 * (E_beta / E)
|
||||
p_1 = 0.5 * (G_lamb / G) - (E_lamb / E)
|
||||
p_0 = 0.5 * (G_beta / E)
|
||||
|
||||
p_33 = - 0.5 * ((E_beta_lamb * G - E_lamb * G_beta) / (G ** 2))
|
||||
p_22 = ((G * G_beta_beta - G_beta * G_beta) / (G ** 2)) - 0.5 * ((E * E_beta_beta - E_beta * E_beta) / (E ** 2))
|
||||
p_11 = 0.5 * ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2)) - ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2))
|
||||
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)
|
||||
|
||||
def q_coef(beta, lamb):
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
BETA_, LAMBDA_, BETA__, LAMBDA__,
|
||||
E_beta, E_lamb, G_beta, G_lamb,
|
||||
E_beta_beta, E_beta_lamb, E_lamb_lamb,
|
||||
G_beta_beta, G_beta_lamb, G_lamb_lamb) = BETA_LAMBDA(beta, lamb)
|
||||
|
||||
q_3 = - 0.5 * (G_beta / E)
|
||||
q_2 = (E_lamb / E) - 0.5 * (G_lamb / G)
|
||||
q_1 = 0.5 * (E_beta / E) - (G_beta / G)
|
||||
q_0 = 0.5 * (E_lamb / G)
|
||||
|
||||
q_33 = - 0.5 * ((E * G_beta_lamb - E_lamb * G_lamb) / (E ** 2))
|
||||
q_22 = ((E * E_lamb_lamb - E_lamb * E_lamb) / (E ** 2)) - 0.5 * ((G * G_lamb_lamb - G_lamb * G_lamb) / (G ** 2))
|
||||
q_11 = 0.5 * ((E * E_beta_lamb - E_beta * E_lamb) / (E ** 2)) - ((G * G_beta_lamb - G_beta * G_lamb) / (G ** 2))
|
||||
q_00 = 0.5 * ((E_lamb_lamb * G - E_lamb * G_lamb) / (G ** 2))
|
||||
|
||||
return (BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00)
|
||||
|
||||
if lamb_1 != lamb_2:
|
||||
# def functions():
|
||||
# def f_beta(lamb, beta, beta_p, X3, X4):
|
||||
# return beta_p
|
||||
#
|
||||
# def f_beta_p(lamb, beta, beta_p, X3, X4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# p_3, p_2, p_1, p_0,
|
||||
# p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
# return p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
#
|
||||
# def f_X3(lamb, beta, beta_p, X3, X4):
|
||||
# return X4
|
||||
#
|
||||
# def f_X4(lamb, beta, beta_p, X3, X4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# p_3, p_2, p_1, p_0,
|
||||
# p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
# return (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
# (3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
|
||||
#
|
||||
# return [f_beta, f_beta_p, f_X3, f_X4]
|
||||
|
||||
def buildODElamb():
|
||||
def ODE(lamb, v):
|
||||
beta, beta_p, X3, X4 = v
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
p_3, p_2, p_1, p_0,
|
||||
p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
|
||||
|
||||
dbeta = beta_p
|
||||
dbeta_p = p_3 * beta_p ** 3 + p_2 * beta_p ** 2 + p_1 * beta_p + p_0
|
||||
dX3 = X4
|
||||
dX4 = (p_33 * beta_p ** 3 + p_22 * beta_p ** 2 + p_11 * beta_p + p_00) * X3 + \
|
||||
(3 * p_3 * beta_p ** 2 + 2 * p_2 * beta_p + p_1) * X4
|
||||
|
||||
return np.array([dbeta, dbeta_p, dX3, dX4])
|
||||
|
||||
return ODE
|
||||
|
||||
N = n
|
||||
|
||||
dlamb = lamb_2 - lamb_1
|
||||
alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
|
||||
|
||||
@@ -182,10 +167,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
|
||||
(_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
|
||||
beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph)
|
||||
|
||||
converged = False
|
||||
iterations = 0
|
||||
|
||||
# funcs = functions()
|
||||
ode_lamb = buildODElamb()
|
||||
|
||||
def solve_newton(beta_p0_init: float):
|
||||
@@ -307,11 +288,10 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
|
||||
return alpha_1, alpha_2, s
|
||||
|
||||
else: # lamb_1 == lamb_2
|
||||
|
||||
N = n
|
||||
dbeta = beta_2 - beta_1
|
||||
|
||||
if abs(dbeta) < 10**-15:
|
||||
if abs(dbeta) < 1e-15:
|
||||
if all_points:
|
||||
return 0, 0, 0, np.array([]), np.array([])
|
||||
else:
|
||||
@@ -319,68 +299,20 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
|
||||
|
||||
lamb_0 = 0
|
||||
|
||||
converged = False
|
||||
iterations = 0
|
||||
|
||||
# def functions_beta():
|
||||
# def g_lamb(beta, lamb, lamb_p, Y3, Y4):
|
||||
# return lamb_p
|
||||
#
|
||||
# def g_lamb_p(beta, lamb, lamb_p, Y3, Y4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# q_3, q_2, q_1, q_0,
|
||||
# q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
# return q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
#
|
||||
# def g_Y3(beta, lamb, lamb_p, Y3, Y4):
|
||||
# return Y4
|
||||
#
|
||||
# def g_Y4(beta, lamb, lamb_p, Y3, Y4):
|
||||
# (BETA, LAMBDA, E, G,
|
||||
# q_3, q_2, q_1, q_0,
|
||||
# q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
# return (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
# (3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
|
||||
#
|
||||
# return [g_lamb, g_lamb_p, g_Y3, g_Y4]
|
||||
|
||||
def buildODEbeta():
|
||||
def ODE(beta, v):
|
||||
lamb, lamb_p, Y3, Y4 = v
|
||||
|
||||
(BETA, LAMBDA, E, G,
|
||||
q_3, q_2, q_1, q_0,
|
||||
q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
|
||||
|
||||
dlamb = lamb_p
|
||||
dlamb_p = q_3 * lamb_p ** 3 + q_2 * lamb_p ** 2 + q_1 * lamb_p + q_0
|
||||
dY3 = Y4
|
||||
dY4 = (q_33 * lamb_p ** 3 + q_22 * lamb_p ** 2 + q_11 * lamb_p + q_00) * Y3 + \
|
||||
(3 * q_3 * lamb_p ** 2 + 2 * q_2 * lamb_p + q_1) * Y4
|
||||
|
||||
return np.array([dlamb, dlamb_p, dY3, dY4])
|
||||
return ODE
|
||||
|
||||
# funcs_beta = functions_beta()
|
||||
ode_beta = buildODEbeta()
|
||||
|
||||
for i in range(iter_max):
|
||||
iterations = i + 1
|
||||
|
||||
startwerte = [lamb_1, lamb_0, 0.0, 1.0]
|
||||
|
||||
# werte = rk.verfahren(funcs_beta, startwerte, dbeta, N, False)
|
||||
beta_list, werte = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
|
||||
|
||||
beta_end = beta_list[-1]
|
||||
# beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||||
lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||||
|
||||
d_lamb_end_d_lambda0 = Y3_end
|
||||
delta = lamb_end - lamb_2
|
||||
|
||||
if abs(delta) < epsilon:
|
||||
converged = True
|
||||
break
|
||||
|
||||
if abs(d_lamb_end_d_lambda0) < 1e-20:
|
||||
@@ -393,7 +325,6 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
|
||||
|
||||
lamb_0 = lamb_0 - step
|
||||
|
||||
# werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N, False)
|
||||
beta_list, werte = rk.rk4(ode_beta, beta_1, np.array([lamb_1, lamb_0, 0.0, 1.0]), dbeta, N, False)
|
||||
|
||||
# beta_arr = np.zeros(N + 1)
|
||||
|
||||
@@ -335,14 +335,14 @@ class EllipsoidTriaxial:
|
||||
Z = self.b * sin(beta) * sqrt(k**2 + k_**2 * sin(lamb)**2)
|
||||
return np.array([X, Y, Z])
|
||||
|
||||
def cart2ell_yFake(self, point: NDArray) -> Tuple[float, float]:
|
||||
def cart2ell_yFake(self, point: NDArray, delta_y: float = 1e-4) -> Tuple[float, float]:
|
||||
"""
|
||||
Bei Fehlschlagen von cart2ell
|
||||
:param point: Punkt in kartesischen Koordinaten
|
||||
:param delta_y: Startwert für Suche nach kleinstmöglichem delta_y
|
||||
:return: ellipsoidische Breite und Länge
|
||||
"""
|
||||
x, y, z = point
|
||||
delta_y = 1e-4
|
||||
best_delta = np.inf
|
||||
while True:
|
||||
try:
|
||||
|
||||
13
utils_angle.py
Normal file
13
utils_angle.py
Normal file
@@ -0,0 +1,13 @@
|
||||
import numpy as np
|
||||
|
||||
|
||||
def arccot(x):
|
||||
return np.arctan2(1.0, x)
|
||||
|
||||
|
||||
def cot(a):
|
||||
return np.cos(a) / np.sin(a)
|
||||
|
||||
|
||||
def wrap_to_pi(x):
|
||||
return (x + np.pi) % (2 * np.pi) - np.pi
|
||||
Reference in New Issue
Block a user