diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index f8cf750..0548ec2 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -10,11 +10,17 @@ import ausgaben as aus from utils_angle import cot, arccot, wrap_to_pi -def norm_a(a): - if a < 0.0: - a += np.pi +def norm_a(a: float) -> float: + a = float(a) % (2 * np.pi) 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): dlam = wrap_to_pi(lam2 - lam1) y = np.sin(dlam) * np.cos(beta2) @@ -24,6 +30,7 @@ def sph_azimuth(beta1, lam1, beta2, lam2): a += 2 * np.pi return a + # Panou 2013 def gha2_num( ell: EllipsoidTriaxial, @@ -51,62 +58,62 @@ def gha2_num( # Berechnung Koeffizienten, Gaußschen Fundamentalgrößen 1. Ordnung sowie deren Ableitungen 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 + 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 + 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 ) - BETA_ = (ell.ax ** 2 * ell.Ey ** 2 * np.sin(2 * beta)) / ( - ell.Ex ** 2 - ell.Ey ** 2 * np.sin(beta) ** 2 + 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 + LAMBDA_ = -(ell.b**2 * ell.Ee**2 * np.sin(2 * lamb)) / ( + ell.Ex**2 - ell.Ee**2 * np.cos(lamb) ** 2 ) ** 2 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 + (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 + (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) - G = LAMBDA * (ell.Ey ** 2 * np.cos(beta) ** 2 + ell.Ee ** 2 * np.sin(lamb) ** 2) + E = BETA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2) + G = LAMBDA * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2) 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) + 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) + E_lamb = BETA * ell.Ee**2 * np.sin(2 * lamb) - G_beta = -LAMBDA * ell.Ey ** 2 * np.sin(2 * beta) + 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) + LAMBDA_ * (ell.Ey**2 * np.cos(beta) ** 2 + ell.Ee**2 * np.sin(lamb) ** 2) + + LAMBDA * ell.Ee**2 * np.sin(2 * lamb) ) 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) + 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) + 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_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) + 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 ( @@ -220,10 +227,14 @@ def gha2_num( (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb) 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) if abs(lamb_1 - lamb_0) >= 1e-15: 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) lamb0 = float(lamb_0) @@ -237,7 +248,9 @@ def gha2_num( 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 + 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) alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1) @@ -295,36 +308,38 @@ def gha2_num( beta_p_arr = np.array([st[1] for st in states], dtype=float) (_, _, 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_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_arr[-1])) + alpha_0 = azimut(E_start, G_start, dbeta_du=beta_p_arr[0] * sgn, dlamb_du=1.0 * sgn) + alpha_1 = azimut(E_end, G_end, dbeta_du=beta_p_arr[-1] * sgn, dlamb_du=1.0 * sgn) # Distanz aus Arrays integrand = np.zeros(N + 1, dtype=float) 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) + 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 + 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 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) beta_end, beta_p_end, _, _ = y_end (_, _, 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)) - alpha_2 = norm_a(arccot(np.sqrt(E_end / G_end) * beta_p_end)) + (_, _, E_end, G_end, *_) = BETA_LAMBDA(float(beta_end), lamb1) + 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) N = int(n) @@ -339,15 +354,18 @@ def gha2_num( lamb0 = float(lamb_0) beta1 = float(beta_1) lamb1 = float(lamb_1) + sgn = 1.0 if dbeta >= 0.0 else -1.0 def ode_beta(beta, v): lamb, lamb_p, Y3, Y4 = v (_, _, _, _, 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 + 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 + 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) lamb_p0 = 0.0 @@ -376,36 +394,39 @@ def gha2_num( 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) - (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) - (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1]) + (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0]) + (_, _, 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_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_arr[-1])) + alpha_0 = azimut(E_start, G_start, dbeta_du=1.0 * sgn, dlamb_du=lamb_p_arr[0] * sgn) + 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) for i in range(N + 1): (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i]) - integrand[i] = np.sqrt(Ei + Gi * lamb_p_arr[i]**2) + integrand[i] = np.sqrt(Ei + Gi * lamb_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 + 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 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) lamb_end, lamb_p_end, _, _ = y_end - (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta0, lamb0) - (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta1, lamb1) + (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0) + 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)) - alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_end)) + (_, _, E_end, G_end, *_) = BETA_LAMBDA(beta1, float(lamb_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__": # ell = EllipsoidTriaxial.init_name("BursaSima1980round")