From 77c7a6f9ab8e3ac81f0a0704019a8f68d65221f4 Mon Sep 17 00:00:00 2001 From: Hendrik Date: Thu, 5 Feb 2026 11:12:17 +0100 Subject: [PATCH] Umstrukturierung --- GHA_triaxial/gha1_num.py | 54 +++---- GHA_triaxial/gha2_num.py | 335 ++++++++++++++++----------------------- utils_angle.py | 13 ++ 3 files changed, 173 insertions(+), 229 deletions(-) create mode 100644 utils_angle.py diff --git a/GHA_triaxial/gha1_num.py b/GHA_triaxial/gha1_num.py index c527216..de50d74 100644 --- a/GHA_triaxial/gha1_num.py +++ b/GHA_triaxial/gha1_num.py @@ -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) diff --git a/GHA_triaxial/gha2_num.py b/GHA_triaxial/gha2_num.py index 4d4d229..734dc22 100644 --- a/GHA_triaxial/gha2_num.py +++ b/GHA_triaxial/gha2_num.py @@ -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) diff --git a/utils_angle.py b/utils_angle.py new file mode 100644 index 0000000..e560739 --- /dev/null +++ b/utils_angle.py @@ -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