333 lines
12 KiB
Python
333 lines
12 KiB
Python
import numpy as np
|
||
from ellipsoide import EllipsoidTriaxial
|
||
import Numerische_Integration.num_int_runge_kutta as rk
|
||
import ausgaben as aus
|
||
|
||
# Panou 2013
|
||
def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30):
|
||
"""
|
||
|
||
:param ell: triaxiales Ellipsoid
|
||
:param beta_1: reduzierte ellipsoidische Breite Punkt 1
|
||
:param lamb_1: elllipsoidische Länge Punkt 1
|
||
:param beta_2: reduzierte ellipsoidische Breite Punkt 2
|
||
:param lamb_2: elllipsoidische Länge Punkt 2
|
||
:param n: Anzahl Schritte
|
||
:param epsilon:
|
||
:param iter_max: Maximale Anzhal Iterationen
|
||
:return:
|
||
"""
|
||
|
||
# h_x, h_y, h_e entsprechen E_x, E_y, E_e
|
||
|
||
def arccot(x):
|
||
return np.arctan2(1.0, x)
|
||
|
||
|
||
def BETA_LAMBDA(beta, lamb):
|
||
|
||
BETA = (ell.ay**2 * np.sin(beta)**2 + ell.b**2 * np.cos(beta)**2) / (ell.Ex**2 - ell.Ey**2 * np.sin(beta)**2)
|
||
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]
|
||
|
||
N = n
|
||
|
||
dlamb = lamb_2 - lamb_1
|
||
|
||
if abs(dlamb) < 1e-15:
|
||
beta_0 = 0.0
|
||
else:
|
||
beta_0 = (beta_2 - beta_1) / (lamb_2 - lamb_1)
|
||
|
||
converged = False
|
||
iterations = 0
|
||
|
||
funcs = functions()
|
||
|
||
for i in range(iter_max):
|
||
iterations = i + 1
|
||
|
||
startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0]
|
||
|
||
werte = rk.verfahren(funcs, startwerte, dlamb, N)
|
||
lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1]
|
||
|
||
d_beta_end_d_beta0 = X3_end
|
||
delta = beta_end - beta_2
|
||
|
||
if abs(delta) < epsilon:
|
||
converged = True
|
||
break
|
||
|
||
if abs(d_beta_end_d_beta0) < 1e-20:
|
||
raise RuntimeError("Abbruch.")
|
||
|
||
max_step = 0.5
|
||
step = delta / d_beta_end_d_beta0
|
||
if abs(step) > max_step:
|
||
step = np.sign(step) * max_step
|
||
beta_0 = beta_0 - step
|
||
|
||
if not converged:
|
||
raise RuntimeError("konvergiert nicht.")
|
||
|
||
# Z
|
||
werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N)
|
||
|
||
beta_arr = np.zeros(N + 1)
|
||
lamb_arr = np.zeros(N + 1)
|
||
beta_p_arr = np.zeros(N + 1)
|
||
|
||
for i, state in enumerate(werte):
|
||
lamb_arr[i] = state[0]
|
||
beta_arr[i] = state[1]
|
||
beta_p_arr[i] = state[2]
|
||
|
||
(_, _, E1, G1,
|
||
*_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
|
||
(_, _, E2, G2,
|
||
*_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
|
||
|
||
alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0])
|
||
alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1])
|
||
|
||
integrand = np.zeros(N + 1)
|
||
for i in range(N + 1):
|
||
(_, _, Ei, Gi,
|
||
*_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
|
||
integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi)
|
||
|
||
h = abs(dlamb) / N
|
||
if N % 2 == 0:
|
||
S = integrand[0] + integrand[-1] \
|
||
+ 4.0 * np.sum(integrand[1:-1:2]) \
|
||
+ 2.0 * np.sum(integrand[2:-1:2])
|
||
s = h / 3.0 * S
|
||
else:
|
||
s = np.trapz(integrand, dx=h)
|
||
|
||
beta0 = beta_arr[0]
|
||
lamb0 = lamb_arr[0]
|
||
c = np.sqrt(
|
||
(np.cos(beta0) ** 2 + (ell.Ee**2 / ell.Ex**2) * np.sin(beta0) ** 2) * np.sin(alpha_1) ** 2
|
||
+ (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2
|
||
)
|
||
|
||
return alpha_1, alpha_2, s
|
||
|
||
if lamb_1 == lamb_2:
|
||
|
||
N = n
|
||
dbeta = beta_2 - beta_1
|
||
|
||
if abs(dbeta) < 10**-15:
|
||
return 0, 0, 0
|
||
|
||
lamb_0 = 0
|
||
|
||
converged = False
|
||
iterations = 0
|
||
|
||
def functions_beta():
|
||
def g_lamb(beta, lamb, lamb_p, Y3, Y4):
|
||
return lamb_p
|
||
|
||
def 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]
|
||
|
||
funcs_beta = functions_beta()
|
||
|
||
for i in range(iter_max):
|
||
iterations = i + 1
|
||
|
||
startwerte = [beta_1, lamb_1, lamb_0, 0.0, 1.0]
|
||
|
||
werte = rk.verfahren(funcs_beta, startwerte, dbeta, N)
|
||
beta_end, lamb_end, lamb_p_end, Y3_end, Y4_end = werte[-1]
|
||
|
||
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:
|
||
raise RuntimeError("Abbruch (Ableitung ~ 0).")
|
||
|
||
max_step = 1.0
|
||
step = delta / d_lamb_end_d_lambda0
|
||
if abs(step) > max_step:
|
||
step = np.sign(step) * max_step
|
||
|
||
lamb_0 = lamb_0 - step
|
||
|
||
werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N)
|
||
|
||
beta_arr = np.zeros(N + 1)
|
||
lamb_arr = np.zeros(N + 1)
|
||
lambda_p_arr = np.zeros(N + 1)
|
||
|
||
for i, state in enumerate(werte):
|
||
beta_arr[i] = state[0]
|
||
lamb_arr[i] = state[1]
|
||
lambda_p_arr[i] = state[2]
|
||
|
||
# Azimute
|
||
(BETA1, LAMBDA1, E1, G1,
|
||
*_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
|
||
(BETA2, LAMBDA2, E2, G2,
|
||
*_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
|
||
|
||
alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lambda_p_arr[0])
|
||
alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lambda_p_arr[-1])
|
||
|
||
integrand = np.zeros(N + 1)
|
||
for i in range(N + 1):
|
||
(_, _, Ei, Gi,
|
||
*_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
|
||
integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2)
|
||
|
||
h = abs(dbeta) / N
|
||
if N % 2 == 0:
|
||
S = integrand[0] + integrand[-1] \
|
||
+ 4.0 * np.sum(integrand[1:-1:2]) \
|
||
+ 2.0 * np.sum(integrand[2:-1:2])
|
||
s = h / 3.0 * S
|
||
else:
|
||
s = np.trapz(integrand, dx=h)
|
||
|
||
return alpha_1, alpha_2, s
|
||
|
||
|
||
if __name__ == "__main__":
|
||
ell = EllipsoidTriaxial.init_name("BursaSima1980round")
|
||
beta1 = np.deg2rad(75)
|
||
lamb1 = np.deg2rad(-90)
|
||
beta2 = np.deg2rad(75)
|
||
lamb2 = np.deg2rad(66)
|
||
a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
|
||
print(aus.gms("a1", a1, 4))
|
||
print(aus.gms("a2", a2, 4))
|
||
print(s)
|
||
|
||
|
||
|