he2851/Masterprojekt
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

ell2

Browse Source
This commit is contained in:
ni3019
2026-02-05 15:56:09 +01:00
parent 39641b5293
commit 09ae06e9b2
1 changed files with 140 additions and 127 deletions
Download Patch File Download Diff File Expand all files Collapse all files

267
GHA_triaxial/gha2_num.py
View File

@@ -17,124 +17,6 @@ def sph_azimuth(beta1, lam1, beta2, lam2):
a += 2 * np.pi
return a
def BETA_LAMBDA(ell, 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(ell, 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(ell, 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,
@@ -155,6 +37,137 @@ 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 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
if lamb_1 != lamb_2:
N = n
@@ -164,7 +177,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
if abs(dlamb) < 1e-15:
beta_0 = 0.0
else:
(_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1)
(_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
beta_0 = np.sqrt(G1 / E1) * cot(alpha0_sph)
ode_lamb = buildODElamb()
@@ -196,7 +209,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
return False, None, None, None
alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
(_, _, E1, G1, *_) = BETA_LAMBDA(ell, beta_1, lamb_1)
(_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph)
guesses = [
@@ -220,7 +233,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
integrand = np.zeros(N + 1)
for i in range(N + 1):
(_, _, Ei, Gi, *_) = BETA_LAMBDA(ell, beta_arr_c[i], lamb_arr_c[i])
(_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr_c[i], lamb_arr_c[i])
integrand[i] = np.sqrt(Ei * beta_p_arr_c[i] ** 2 + Gi)
h = abs(dlamb) / N
@@ -253,9 +266,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
beta_p_arr[i] = state[1]
(_, _, E1, G1,
*_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0])
*_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
(_, _, E2, G2,
*_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1])
*_) = 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])
@@ -263,7 +276,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
integrand = np.zeros(N + 1)
for i in range(N + 1):
(_, _, Ei, Gi,
*_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i])
*_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
integrand[i] = np.sqrt(Ei * beta_p_arr[i] ** 2 + Gi)
h = abs(dlamb) / N
@@ -341,9 +354,9 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
# Azimute
(BETA1, LAMBDA1, E1, G1,
*_) = BETA_LAMBDA(ell, beta_arr[0], lamb_arr[0])
*_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
(BETA2, LAMBDA2, E2, G2,
*_) = BETA_LAMBDA(ell, beta_arr[-1], lamb_arr[-1])
*_) = 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])
@@ -351,7 +364,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1: float, lamb_1: float, beta_2: float
integrand = np.zeros(N + 1)
for i in range(N + 1):
(_, _, Ei, Gi,
*_) = BETA_LAMBDA(ell, beta_arr[i], lamb_arr[i])
*_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
integrand[i] = np.sqrt(Ei + Gi * lambda_p_arr[i] ** 2)
h = abs(dbeta) / N