Umrechnung alpha, Näherungslösung GHA 1

2026-01-11 13:14:43 +01:00
parent 797afdfd6f
commit 4d5b6fcc3e
6 changed files with 205 additions and 71 deletions
--- a/GHA_triaxial/approx_gha1.py
+++ b/GHA_triaxial/approx_gha1.py
@@ -0,0 +1,63 @@
+import numpy as np
+from numpy import sin, cos, arcsin, arccos, arctan2
+from ellipsoide import EllipsoidTriaxial
+import matplotlib.pyplot as plt
+from panou import louville_constant, func_sigma_ell, gha1_ana
+import plotly.graph_objects as go
+import winkelumrechnungen as wu
+
+def gha1(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float, ds: int):
+    l0 = louville_constant(ell, p0, alpha0)
+    points = [p0]
+    alphas = [alpha0]
+    s_curr = 0.0
+    while s_curr < s:
+        ds_step = min(ds, s - s_curr)
+        if ds_step < 1e-8:
+            break
+        p1 = points[-1]
+        alpha1 = alphas[-1]
+        x1, y1, z1 = p1
+        sigma = func_sigma_ell(ell, x1, y1, z1, alpha1)
+        p2 = p1 + ds_step * sigma
+        p2, _, _, _ = ell.cartonell(p2)
+        ds_step = np.linalg.norm(p2 - p1)
+
+        points.append(p2)
+        dalpha = 1e-6
+        l2 = louville_constant(ell, p2, alpha1)
+        dl_dalpha = (louville_constant(ell, p2, alpha1+dalpha) - l2) / dalpha
+        alpha2 = alpha1 + (l0 - l2) / dl_dalpha
+        alphas.append(alpha2)
+        s_curr += ds_step
+    return points[-1], alphas[-1], np.array(points)
+
+def show_points(points, p0, p1):
+    fig = go.Figure()
+
+    fig.add_scatter3d(x=points[:, 0], y=points[:, 1], z=points[:, 2],
+                      mode='lines', line=dict(color="red", width=3), name="Approx")
+    fig.add_scatter3d(x=[p0[0]], y=[p0[1]], z=[p0[2]],
+                      mode='markers', marker=dict(color="green"), name="P0")
+    fig.add_scatter3d(x=[p1[0]], y=[p1[1]], z=[p1[2]],
+                      mode='markers', marker=dict(color="green"), name="P1")
+
+    fig.update_layout(
+        scene=dict(xaxis_title='X [km]',
+                   yaxis_title='Y [km]',
+                   zaxis_title='Z [km]',
+                   aspectmode='data'),
+        title="CHAMP")
+
+    fig.show()
+
+
+if __name__ == '__main__':
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    P0 = ell.para2cart(0, 0)
+    alpha0 = wu.deg2rad(90)
+    s = 1000000
+    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 60, maxPartCircum=32)
+    P1_app, alpha1_app, points = gha1(ell, P0, alpha0, s, 5000)
+    show_points(points, P0, P1_ana)
+    print(np.linalg.norm(P1_app - P1_ana))
--- a/GHA_triaxial/numeric_examples_karney.py
+++ b/GHA_triaxial/numeric_examples_karney.py
@@ -16,7 +16,7 @@ def get_random_examples(num):
    :param num:
    :return:
    """
-    random.seed(42)
+    # random.seed(42)
    with open("Karney_2024_Testset.txt") as datei:
        lines = datei.readlines()
    examples = []
--- a/GHA_triaxial/panou.py
+++ b/GHA_triaxial/panou.py
@@ -68,7 +68,7 @@ def buildODE(ell: EllipsoidTriaxial) -> Callable:
        return np.array([dxds, ddx, dyds, ddy, dzds, ddz])
    return ODE

-def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float]:
+def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, num: int) -> tuple[np.ndarray, float, list]:
    """
    Panou, Korakitits 2019
    :param ell:
@@ -95,15 +95,15 @@ def gha1_num(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,

    p1, q1 = pq_ell(ell, x1, y1, z1)
    sigma = np.array([dx1ds, dy1ds, dz1ds])
-    P = p1 @ sigma
-    Q = q1 @ sigma
+    P = float(p1 @ sigma)
+    Q = float(q1 @ sigma)

    alpha1 = arctan2(P, Q)

    if alpha1 < 0:
        alpha1 += 2 * np.pi

-    return np.array([x1, y1, z1]), alpha1
+    return np.array([x1, y1, z1]), alpha1, werte

 # ---------------------------------------------------------------------------------------------------------------------

@@ -225,8 +225,8 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):

    # 52-53
    sigma = np.array([dx_s, dy_s, dz_s])
-    P = p_s @ sigma
-    Q = q_s @ sigma
+    P = float(p_s @ sigma)
+    Q = float(q_s @ sigma)

    # 51
    alpha1 = arctan2(P, Q)
@@ -237,7 +237,7 @@ def gha1_ana_step(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
    return np.array([x_s, y_s, z_s]), alpha1


-def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4):
+def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float, maxM: int, maxPartCircum: int = 4) -> tuple[np.ndarray, float]:
    if s > np.pi / maxPartCircum * ell.ax:
        s /= 2
        point_step, alpha_step = gha1_ana(ell, point, alpha0, s, maxM, maxPartCircum)
@@ -251,39 +251,103 @@ def gha1_ana(ell: EllipsoidTriaxial, point: np.ndarray, alpha0: float, s: float,

    return point_end, alpha_end

+def alpha_para2ell(ell: EllipsoidTriaxial, u: float, v: float, alpha_para: float) -> tuple[float, float, float]:
+    x, y, z = ell.para2cart(u, v)
+    beta, lamb = ell.para2ell(u, v)
+
+    p_para, q_para = pq_para(ell, x, y, z)
+    sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
+
+    p_ell, q_ell = pq_ell(ell, x, y, z)
+    alpha_ell = arctan2(p_ell @ sigma_para, q_ell @ sigma_para)
+    sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell)
+
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-12:
+        raise Exception("Alpha Umrechnung fehlgeschlagen")
+
+    return beta, lamb, alpha_ell
+
+def alpha_ell2para(ell: EllipsoidTriaxial, beta: float, lamb: float, alpha_ell: float) -> tuple[float, float, float]:
+    x, y, z = ell.ell2cart(beta, lamb)
+    u, v = ell.ell2para(beta, lamb)
+
+    p_ell, q_ell = pq_ell(ell, x, y, z)
+    sigma_ell = p_ell * sin(alpha_ell) + q_ell * cos(alpha_ell)
+
+    p_para, q_para = pq_para(ell, x, y, z)
+    alpha_para = arctan2(p_para @ sigma_ell, q_para @ sigma_ell)
+    sigma_para = p_para * sin(alpha_para) + q_para * cos(alpha_para)
+
+    if np.linalg.norm(sigma_para - sigma_ell) > 1e-12:
+        raise Exception("Alpha Umrechnung fehlgeschlagen")
+        print("Alpha Umrechnung fehlgeschlagen:", np.linalg.norm(sigma_para - sigma_ell))
+
+    return u, v, alpha_para
+
+def func_sigma_ell(ell, x, y, z, alpha):
+    p, q = pq_ell(ell, x, y, z)
+    sigma = p * sin(alpha) + q * cos(alpha)
+    return sigma
+
+def func_sigma_para(ell, x, y, z, alpha):
+    p, q = pq_para(ell, x, y, z)
+    sigma = p * sin(alpha) + q * cos(alpha)
+    return sigma
+
+
+def louville_constant(ell: EllipsoidTriaxial, p0: np.ndarray, alpha: float) -> float:
+    beta, lamb = ell.cart2ell(p0)
+    l = ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 - ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2
+    # x, y, z = p0
+    # t1, t2 = ell.func_t12(x, y, z)
+    # l_cart = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2)
+    # if abs(l - l_cart) > 1e-12:
+    #     # raise Exception("Louville constant fehlgeschlagen")
+    #     print("Diff zwischen constant:", abs(l - l_cart))
+    return l
+
+def louville_l2c(ell, l):
+    return sqrt((l + ell.Ee**2) / ell.Ex**2)
+
+def louville_c2l(ell, c):
+    return ell.Ex**2 * c**2 - ell.Ee**2
+

 if __name__ == "__main__":
-    ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")

-    diffs_panou = []
-    examples_panou = ne_panou.get_random_examples(5)
-    for example in examples_panou:
-        beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example
-        P0 = ell.ell2cart(beta0, lamb0)
-
-        P1_num, alpha1_num = gha1_num(ell, P0, alpha0, s, 100)
-        beta1_num, lamb1_num = ell.cart2ell(P1_num)
-        P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 60)
-        beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
-        diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana)))
-    diffs_panou = np.array(diffs_panou)
-    mask_360 = (diffs_panou > 359) & (diffs_panou < 361)
-    diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
-    print(diffs_panou)
+    # ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
+    # diffs_panou = []
+    # examples_panou = ne_panou.get_random_examples(5)
+    # for example in examples_panou:
+    #     beta0, lamb0, beta1, lamb1, _, alpha0_ell, alpha1_ell, s = example
+    #     P0 = ell.ell2cart(beta0, lamb0)
+    #
+    #     P1_num, alpha1_num, _ = gha1_num(ell, P0, alpha0_ell, s, 100)
+    #     beta1_num, lamb1_num = ell.cart2ell(P1_num)
+    #
+    #     _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
+    #     P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 60)
+    #     beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
+    #     diffs_panou.append((abs(beta1-beta1_num), abs(lamb1-lamb1_num), abs(beta1-beta1_ana), abs(lamb1-lamb1_ana)))
+    # diffs_panou = np.array(diffs_panou)
+    # mask_360 = (diffs_panou > 359) & (diffs_panou < 361)
+    # diffs_panou[mask_360] = np.abs(diffs_panou[mask_360] - 360)
+    # print(diffs_panou)

    ell = ellipsoide.EllipsoidTriaxial.init_name("KarneyTest2024")
-
    diffs_karney = []
-    examples_karney = ne_karney.get_examples((30499, 30500, 40500))
-    # examples_karney = ne_karney.get_random_examples(5)
+    # examples_karney = ne_karney.get_examples((30499, 30500, 40500))
+    examples_karney = ne_karney.get_random_examples(20)
    for example in examples_karney:
-        beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example
+        beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example
        P0 = ell.ell2cart(beta0, lamb0)

-        P1_num, alpha1_num = gha1_num(ell, P0, alpha0, s, 100)
+        P1_num, alpha1_num, _ = gha1_num(ell, P0, alpha0_ell, s, 5000)
        beta1_num, lamb1_num = ell.cart2ell(P1_num)
+
        try:
-            P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, 40)
+            _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)
+            P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0_para, s, 30, maxPartCircum=16)
            beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)
        except:
            beta1_ana, lamb1_ana = np.inf, np.inf
@@ -293,4 +357,5 @@ if __name__ == "__main__":
    mask_360 = (diffs_karney > 359) & (diffs_karney < 361)
    diffs_karney[mask_360] = np.abs(diffs_karney[mask_360] - 360)
    print(diffs_karney)
+
    pass
--- a/GHA_triaxial/panou_2013_2GHA_num.py
+++ b/GHA_triaxial/panou_2013_2GHA_num.py
@@ -244,7 +244,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
            + (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2
        )

-        return alpha_1, alpha_2, s
+        return alpha_1, alpha_2, s, beta_arr, lamb_arr

    if lamb_1 == lamb_2:

@@ -252,7 +252,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
        dbeta = beta_2 - beta_1

        if abs(dbeta) < 10**-15:
-            return 0, 0, 0
+            return 0, 0, 0, np.array([]), np.array([])

        lamb_0 = 0

@@ -369,7 +369,7 @@ def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, ep
        else:
            s = np.trapz(integrand, dx=h)

-        return alpha_1, alpha_2, s
+        return alpha_1, alpha_2, s, beta_arr, lamb_arr


 if __name__ == "__main__":
@@ -391,35 +391,37 @@ if __name__ == "__main__":
    # print(s)


-    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    diffs_panou = []
-    examples_panou = ne_panou.get_random_examples(4)
-    for example in examples_panou:
-        beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example
-        P0 = ell.ell2cart(beta0, lamb0)
-        try:
-            alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
-            diffs_panou.append(
-                (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num)))
-        except:
-            print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
-    diffs_panou = np.array(diffs_panou)
-    print(diffs_panou)
+    # ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    # diffs_panou = []
+    # examples_panou = ne_panou.get_random_examples(4)
+    # for example in examples_panou:
+    #     beta0, lamb0, beta1, lamb1, _, alpha0, alpha1, s = example
+    #     P0 = ell.ell2cart(beta0, lamb0)
+    #     try:
+    #         alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
+    #         diffs_panou.append(
+    #             (wu.rad2deg(abs(alpha0 - alpha0_num)), wu.rad2deg(abs(alpha1 - alpha1_num)), abs(s - s_num)))
+    #     except:
+    #         print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
+    # diffs_panou = np.array(diffs_panou)
+    # print(diffs_panou)
+    #
+    # ell = EllipsoidTriaxial.init_name("KarneyTest2024")
+    # diffs_karney = []
+    # # examples_karney = ne_karney.get_examples((30500, 40500))
+    # examples_karney = ne_karney.get_random_examples(2)
+    # for example in examples_karney:
+    #     beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example
+    #
+    #     try:
+    #         alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
+    #         diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num)))
+    #     except:
+    #         print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
+    # diffs_karney = np.array(diffs_karney)
+    # print(diffs_karney)

-    ell = EllipsoidTriaxial.init_name("KarneyTest2024")
-    diffs_karney = []
-    # examples_karney = ne_karney.get_examples((30500, 40500))
-    examples_karney = ne_karney.get_random_examples(2)
-    for example in examples_karney:
-        beta0, lamb0, alpha0, beta1, lamb1, alpha1, s = example

-        try:
-            alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=4000, iter_max=10)
-            diffs_karney.append((wu.rad2deg(abs(alpha0-alpha0_num)), wu.rad2deg(abs(alpha1-alpha1_num)), abs(s-s_num)))
-        except:
-            print(f"Fehler für {beta0}, {lamb0}, {beta1}, {lamb1}")
-    diffs_karney = np.array(diffs_karney)
-    print(diffs_karney)
    pass


--- a/dashboard.py
+++ b/dashboard.py
@@ -413,7 +413,7 @@ def calc_and_plot(n1, n2,

        if "analytisch" in method1:
            # ana
-            x2, y2, z2 = gha1_ana(ell, p1, alpha_rad, s_val, 70)
+            (x2, y2, z2), alpha2 = gha1_ana(ell, p1, alpha_rad, s_val, 70)
            p2_ana = (float(x2), float(y2), float(z2))
            beta2, lamb2 = ell.cart2ell([x2, y2, z2])

@@ -430,7 +430,7 @@ def calc_and_plot(n1, n2,

        if "numerisch" in method1:
            # num
-            x1, y1, z1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000)
+            (x1, y1, z1), alpha1, werte = gha1_num(ell, p1, alpha_rad, s_val, 10000)
            p2_num = x1, y1, z1
            beta2_num, lamb2_num = ell.cart2ell(p2_num)

--- a/ellipsoide.py
+++ b/ellipsoide.py
@@ -173,6 +173,17 @@ class EllipsoidTriaxial:
                         y / ((1 - self.ee ** 2) * sqrtH),
                         z / ((1 - self.ex ** 2) * sqrtH)])

+    def func_t12(self, x, y, z):
+        c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
+        c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
+              (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
+                      self.ax ** 2 + self.ay ** 2) * z ** 2)
+        t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
+        if t2 == 0:
+            t2 = 1e-18
+        t1 = c0 / t2
+        return t1, t2
+
    def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
        """
        Panou 2014 12ff.
@@ -360,14 +371,7 @@ class EllipsoidTriaxial:

        # ---- Allgemeiner Fall -----

-        c1 = x ** 2 + y ** 2 + z ** 2 - (self.ax ** 2 + self.ay ** 2 + self.b ** 2)
-        c0 = (self.ax ** 2 * self.ay ** 2 + self.ax ** 2 * self.b ** 2 + self.ay ** 2 * self.b ** 2 -
-              (self.ay ** 2 + self.b ** 2) * x ** 2 - (self.ax ** 2 + self.b ** 2) * y ** 2 - (
-                      self.ax ** 2 + self.ay ** 2) * z ** 2)
-        t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
-        if t2 == 0:
-            t2 = 1e-14
-        t1 = c0 / t2
+        t1, t2 = self.func_t12(x, y, z)

        num_beta = max(t1 - self.b ** 2, 0)
        den_beta = max(self.ay ** 2 - t1, 0)