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

Fehlerkorrektur

Browse Source
This commit is contained in:
Tammo.Weber
2026-02-09 11:43:50 +01:00
parent 020d282420
commit 71c13c568a
1 changed files with 86 additions and 65 deletions
Download Patch File Download Diff File Expand all files Collapse all files

71
GHA_triaxial/gha2_num.py
View File

@@ -10,11 +10,17 @@ import ausgaben as aus
from utils_angle import cot, arccot, wrap_to_pi from utils_angle import cot, arccot, wrap_to_pi
def norm_a(a): def norm_a(a: float) -> float:
if a < 0.0: a = float(a) % (2 * np.pi)
a += np.pi
return a return a
def azimut(E: float, G: float, dbeta_du: float, dlamb_du: float) -> float:
north = np.sqrt(E) * dbeta_du
east = np.sqrt(G) * dlamb_du
return norm_a(np.arctan2(east, north))
def sph_azimuth(beta1, lam1, beta2, lam2): def sph_azimuth(beta1, lam1, beta2, lam2):
dlam = wrap_to_pi(lam2 - lam1) dlam = wrap_to_pi(lam2 - lam1)
y = np.sin(dlam) * np.cos(beta2) y = np.sin(dlam) * np.cos(beta2)
@@ -24,6 +30,7 @@ def sph_azimuth(beta1, lam1, beta2, lam2):
a += 2 * np.pi a += 2 * np.pi
return a return a
# Panou 2013 # Panou 2013
def gha2_num( def gha2_num(
ell: EllipsoidTriaxial, ell: EllipsoidTriaxial,
@@ -220,10 +227,14 @@ def gha2_num(
(_, _, E, G, *_) = BETA_LAMBDA(beta, lamb) (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
return np.sqrt(E + G * lamb_p**2) return np.sqrt(E + G * lamb_p**2)
lamb_0 = wrap_to_pi(lamb_0)
lamb_1 = wrap_to_pi(lamb_1)
# Fall 1 (lambda_0 != lambda_1) # Fall 1 (lambda_0 != lambda_1)
if abs(lamb_1 - lamb_0) >= 1e-15: if abs(lamb_1 - lamb_0) >= 1e-15:
N = int(n) N = int(n)
dlamb = float(lamb_1 - lamb_0) dlamb = wrap_to_pi(lamb_1 - lamb_0)
sgn = 1.0 if dlamb >= 0.0 else -1.0
beta0 = float(beta_0) beta0 = float(beta_0)
lamb0 = float(lamb_0) lamb0 = float(lamb_0)
@@ -237,7 +248,9 @@ def gha2_num(
dbeta = beta_p dbeta = beta_p
dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0 dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0
dX3 = X4 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 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], dtype=float) return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1) alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1)
@@ -297,8 +310,8 @@ def gha2_num(
(_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
(_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p_arr[0])) alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p_arr[0] * sgn, dlamb_du=1.0 * sgn)
alpha_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_arr[-1])) alpha_1 = azimut(E_end, G_end, dbeta_du=beta_p_arr[-1] * sgn, dlamb_du=1.0 * sgn)
# Distanz aus Arrays # Distanz aus Arrays
integrand = np.zeros(N + 1, dtype=float) integrand = np.zeros(N + 1, dtype=float)
@@ -308,23 +321,25 @@ def gha2_num(
h = abs(dlamb) / N h = abs(dlamb) / N
if N % 2 == 0: 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 = 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 s = h / 3.0 * S
else: else:
s = np.trapz(integrand, dx=h) s = np.trapz(integrand, dx=h)
return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr
_, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N, integrand_lambda) _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N, integrand_lambda)
beta_end, beta_p_end, _, _ = y_end beta_end, beta_p_end, _, _ = y_end
(_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0) (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
(_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, lamb1) alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p0 * sgn, dlamb_du=1.0 * sgn)
alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p0)) (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), lamb1)
alpha_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_end)) alpha_1 = azimut(E_end, G_end, dbeta_du=float(beta_p_end) * sgn, dlamb_du=1.0 * sgn)
return float(alpha_1), float(alpha_2), float(s) return float(alpha_0), float(alpha_1), float(s)
# Fall 2 (lambda_0 == lambda_1) # Fall 2 (lambda_0 == lambda_1)
N = int(n) N = int(n)
@@ -339,6 +354,7 @@ def gha2_num(
lamb0 = float(lamb_0) lamb0 = float(lamb_0)
beta1 = float(beta_1) beta1 = float(beta_1)
lamb1 = float(lamb_1) lamb1 = float(lamb_1)
sgn = 1.0 if dbeta >= 0.0 else -1.0
def ode_beta(beta, v): def ode_beta(beta, v):
lamb, lamb_p, Y3, Y4 = v lamb, lamb_p, Y3, Y4 = v
@@ -347,7 +363,9 @@ def gha2_num(
dlamb = lamb_p dlamb = lamb_p
dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0 dlamb_p = q_3 * lamb_p**3 + q_2 * lamb_p**2 + q_1 * lamb_p + q_0
dY3 = Y4 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 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], dtype=float) return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
lamb_p0 = 0.0 lamb_p0 = 0.0
@@ -376,11 +394,11 @@ def gha2_num(
lamb_arr = np.array([st[0] for st in states], dtype=float) lamb_arr = np.array([st[0] for st in states], dtype=float)
lamb_p_arr = np.array([st[1] for st in states], dtype=float) lamb_p_arr = np.array([st[1] for st in states], dtype=float)
(BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
(BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p_arr[0])) alpha_0 = azimut(E_start, G_start, dbeta_du=1.0 * sgn, dlamb_du=lamb_p_arr[0] * sgn)
alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_arr[-1])) alpha_1 = azimut(E_end, G_end, dbeta_du=1.0 * sgn, dlamb_du=lamb_p_arr[-1] * sgn)
integrand = np.zeros(N + 1, dtype=float) integrand = np.zeros(N + 1, dtype=float)
for i in range(N + 1): for i in range(N + 1):
@@ -389,23 +407,26 @@ def gha2_num(
h = abs(dbeta) / N h = abs(dbeta) / N
if N % 2 == 0: 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 = 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 s = h / 3.0 * S
else: else:
s = np.trapz(integrand, dx=h) s = np.trapz(integrand, dx=h)
return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr return float(alpha_0), float(alpha_1), float(s), beta_arr, lamb_arr
_, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta) _, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta)
lamb_end, lamb_p_end, _, _ = y_end lamb_end, lamb_p_end, _, _ = y_end
(BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta0, lamb0) (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
(BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta1, lamb1) alpha_0 = azimut(E_start, G_start, dbeta_du=1.0 * sgn, dlamb_du=lamb_p0 * sgn)
alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p0)) (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, float(lamb_end))
alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_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(alpha_1), float(alpha_2), float(s)
if __name__ == "__main__": if __name__ == "__main__":
# ell = EllipsoidTriaxial.init_name("BursaSima1980round") # ell = EllipsoidTriaxial.init_name("BursaSima1980round")