GHA1 numerisch implementiert

# Conflicts:
#	dashboard.py

1,lambda-ES-CMA.py

@@ -0,0 +1,364 @@
+import numpy as np
+import plotly.graph_objects as go
+from ellipsoide import EllipsoidTriaxial
+
+
+def ell2cart(ell: EllipsoidTriaxial, beta, lamb):
+    x = ell.ax * np.cos(beta) * np.cos(lamb)
+    y = ell.ay * np.cos(beta) * np.sin(lamb)
+    z = ell.b  * np.sin(beta)
+    return np.array([x, y, z], dtype=float)
+
+
+def liouville(ell: EllipsoidTriaxial, beta, lamb, alpha):
+    k = (ell.Ee**2) / (ell.Ex**2)
+    c_sq = (np.cos(beta)**2 + k * np.sin(beta)**2) * np.sin(alpha)**2 \
+           + k * np.cos(lamb)**2 * np.cos(alpha)**2
+    return c_sq
+
+
+def normalize_angles(beta, lamb):
+    beta = np.clip(beta, -np.pi / 2.0, np.pi / 2.0)
+    lamb = (lamb + np.pi) % (2 * np.pi) - np.pi
+    return beta, lamb
+
+
+def compute_azimuth(beta0: float, lamb0: float,
+                    beta1: float, lamb1: float) -> float:
+    dlam = lamb1 - lamb0
+    y = np.cos(beta1) * np.sin(dlam)
+    x = np.cos(beta0) * np.sin(beta1) - np.sin(beta0) * np.cos(beta1) * np.cos(dlam)
+    alpha = np.arctan2(y, x)
+    return alpha
+
+
+def sphere_forward_step(beta0: float,
+                        lamb0: float,
+                        alpha0: float,
+                        s: float,
+                        R: float) -> tuple[float, float]:
+
+    delta = s / R
+    sin_beta0 = np.sin(beta0)
+    cos_beta0 = np.cos(beta0)
+
+    sin_beta2 = sin_beta0 * np.cos(delta) + cos_beta0 * np.sin(delta) * np.cos(alpha0)
+    beta2 = np.arcsin(sin_beta2)
+
+    y = np.sin(alpha0) * np.sin(delta) * cos_beta0
+    x = np.cos(delta) - sin_beta0 * sin_beta2
+    dlam = np.arctan2(y, x)
+    lamb2 = lamb0 + dlam
+
+    beta2, lamb2 = normalize_angles(beta2, lamb2)
+    return beta2, lamb2
+
+
+def local_step_objective(candidate: np.ndarray,
+                         beta_start: float,
+                         lamb_start: float,
+                         alpha_prev: float,
+                         c_target: float,
+                         step_length: float,
+                         ell: EllipsoidTriaxial,
+                         beta_pred: float,
+                         lamb_pred: float,
+                         w_L: float = 5.0,
+                         w_d: float = 5.0,
+                         w_p: float = 2.0,
+                         w_a: float = 2.0) -> float:
+
+    beta1, lamb1 = candidate
+    beta1, lamb1 = normalize_angles(beta1, lamb1)
+
+    P0 = ell2cart(ell, beta_start, lamb_start)
+    P1 = ell2cart(ell, beta1, lamb1)
+    d = np.linalg.norm(P1 - P0)
+
+    alpha1 = compute_azimuth(beta_start, lamb_start, beta1, lamb1)
+
+    c1 = liouville(ell, beta1, lamb1, alpha1)
+    J_L = (c1 - c_target) ** 2
+
+    J_d = (d - step_length) ** 2
+
+    d_beta = beta1 - beta_pred
+    d_lamb = lamb1 - lamb_pred
+    d_ang2 = d_beta**2 + (np.cos(beta_pred) * d_lamb)**2
+    J_p = d_ang2
+
+    d_alpha = np.arctan2(np.sin(alpha1 - alpha_prev),
+                         np.cos(alpha1 - alpha_prev))
+    J_a = d_alpha**2
+
+    return w_L * J_L + w_d * J_d + w_p * J_p + w_a * J_a
+
+
+
+
+def ES_CMA_step(beta_start: float,
+                lamb_start: float,
+                alpha_prev: float,
+                c_target: float,
+                step_length: float,
+                ell: EllipsoidTriaxial,
+                beta_pred: float,
+                lamb_pred: float,
+                sigma0: float,
+                stopfitness: float = 1e-18) -> tuple[float, float]:
+
+    N = 2
+    xmean = np.array([beta_pred, lamb_pred], dtype=float)
+    sigma = sigma0
+    stopeval = int(400 * N**2)
+
+    lamb = 30
+    mu = 1
+    weights = np.array([1.0])
+    mueff = 1.0
+
+    cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N)
+    cs = (mueff + 2) / (N + mueff + 5)
+    c1 = 2 / ((N + 1.3)**2 + mueff)
+    cmu = min(1 - c1,
+              2 * (mueff - 2 + 1/mueff) / ((N + 2)**2 + mueff))
+    damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs
+
+    pc = np.zeros(N)
+    ps = np.zeros(N)
+    B = np.eye(N)
+    D = np.eye(N)
+    C = B @ D @ (B @ D).T
+    eigeneval = 0
+    chiN = np.sqrt(N) * (1 - 1 / (4 * N) + 1 / (21 * N**2))
+
+    counteval = 0
+    arx = np.zeros((N, lamb))
+    arz = np.zeros((N, lamb))
+    arfitness = np.zeros(lamb)
+
+    while counteval < stopeval:
+        for k in range(lamb):
+            arz[:, k] = np.random.randn(N)
+            arx[:, k] = xmean + sigma * (B @ D @ arz[:, k])
+
+            arfitness[k] = local_step_objective(
+                arx[:, k],
+                beta_start, lamb_start,
+                alpha_prev,
+                c_target,
+                step_length,
+                ell,
+                beta_pred, lamb_pred
+            )
+            counteval += 1
+
+        idx = np.argsort(arfitness)
+        arfitness = arfitness[idx]
+        arindex = idx
+
+        xold = xmean.copy()
+        xmean = arx[:, arindex[:mu]] @ weights
+        zmean = arz[:, arindex[:mu]] @ weights
+
+        ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * (B @ zmean)
+        norm_ps = np.linalg.norm(ps)
+        hsig = norm_ps / np.sqrt(1 - (1 - cs)**(2 * counteval / lamb)) / chiN < 1.4 + 2 / (N + 1)
+        hsig = 1.0 if hsig else 0.0
+        pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (B @ D @ zmean)
+
+        BDz = B @ D @ arz[:, arindex[:mu]]
+        C = (1 - c1 - cmu) * C \
+            + c1 * (np.outer(pc, pc) + (1 - hsig) * cc * (2 - cc) * C) \
+            + cmu * BDz @ np.diag(weights) @ BDz.T
+
+        sigma = sigma * np.exp((cs / damps) * (norm_ps / chiN - 1))
+
+        if counteval - eigeneval > lamb / ((c1 + cmu) * N * 10):
+            eigeneval = counteval
+            C = (C + C.T) / 2.0
+            eigvals, B = np.linalg.eigh(C)
+            D = np.diag(np.sqrt(np.maximum(eigvals, 1e-20)))
+
+        if arfitness[0] <= stopfitness:
+            break
+
+    xmin = arx[:, arindex[0]]
+    beta1, lamb1 = normalize_angles(xmin[0], xmin[1])
+    return beta1, lamb1
+
+
+def march_geodesic(beta1: float,
+                   lamb1: float,
+                   alpha1: float,
+                   S_total: float,
+                   step_length: float,
+                   ell: EllipsoidTriaxial):
+
+    beta_curr = beta1
+    lamb_curr = lamb1
+    alpha_curr = alpha1
+
+    betas = [beta_curr]
+    lambs = [lamb_curr]
+    alphas = [alpha_curr]
+
+    c_target = liouville(ell, beta_curr, lamb_curr, alpha_curr)
+
+    total_distance = 0.0
+    R_sphere = ell.ax
+    sigma0 = 1e-10
+
+    while total_distance < S_total - 1e-6:
+        remaining = S_total - total_distance
+        this_step = min(step_length, remaining)
+
+        beta_pred, lamb_pred = sphere_forward_step(
+            beta_curr, lamb_curr, alpha_curr,
+            this_step, R_sphere
+        )
+
+        beta_next, lamb_next = ES_CMA_step(
+            beta_curr, lamb_curr, alpha_curr,
+            c_target,
+            this_step,
+            ell,
+            beta_pred,
+            lamb_pred,
+            sigma0=sigma0
+        )
+
+        P_curr = ell2cart(ell, beta_curr, lamb_curr)
+        P_next = ell2cart(ell, beta_next, lamb_next)
+        d_step = np.linalg.norm(P_next - P_curr)
+        total_distance += d_step
+
+        alpha_next = compute_azimuth(beta_curr, lamb_curr, beta_next, lamb_next)
+
+        beta_curr, lamb_curr, alpha_curr = beta_next, lamb_next, alpha_next
+        betas.append(beta_curr)
+        lambs.append(lamb_curr)
+        alphas.append(alpha_curr)
+
+    return np.array(betas), np.array(lambs), np.array(alphas), total_distance
+
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+
+    beta1 = np.deg2rad(10.0)
+    lamb1 = np.deg2rad(10.0)
+    alpha1 = 0.7494562507041596         # Ergebnis 20°, 20°
+
+    STEP_LENGTH = 500.0
+    S_total = 1542703.1458877102
+
+    betas, lambs, alphas, S_real = march_geodesic(
+        beta1, lamb1, alpha1,
+        S_total,
+        STEP_LENGTH,
+        ell
+    )
+
+    print("Anzahl Schritte:", len(betas) - 1)
+    print("Resultierende Gesamtstrecke (Chord, Ellipsoid):", S_real, "m")
+    print("Letzter Punkt (beta, lambda) in Grad:",
+          np.rad2deg(betas[-1]), np.rad2deg(lambs[-1]))
+
+def plot_geodesic_3d(ell: EllipsoidTriaxial,
+                     betas: np.ndarray,
+                     lambs: np.ndarray,
+                     n_beta: int = 60,
+                     n_lamb: int = 120):
+
+    beta_grid = np.linspace(-np.pi/2, np.pi/2, n_beta)
+    lamb_grid = np.linspace(-np.pi, np.pi, n_lamb)
+    B, L = np.meshgrid(beta_grid, lamb_grid, indexing="ij")
+
+    Xs = ell.ax * np.cos(B) * np.cos(L)
+    Ys = ell.ay * np.cos(B) * np.sin(L)
+    Zs = ell.b  * np.sin(B)
+
+    Xp = ell.ax * np.cos(betas) * np.cos(lambs)
+    Yp = ell.ay * np.cos(betas) * np.sin(lambs)
+    Zp = ell.b  * np.sin(betas)
+
+    fig = go.Figure()
+
+    fig.add_trace(
+        go.Surface(
+            x=Xs,
+            y=Ys,
+            z=Zs,
+            opacity=0.6,
+            showscale=False,
+            colorscale="Viridis",
+            name="Ellipsoid"
+        )
+    )
+
+    fig.add_trace(
+        go.Scatter3d(
+            x=Xp,
+            y=Yp,
+            z=Zp,
+            mode="lines+markers",
+            line=dict(width=5, color="red"),
+            marker=dict(size=4, color="black"),
+            name="ES-CMA Schritte"
+        )
+    )
+
+    fig.add_trace(
+        go.Scatter3d(
+            x=[Xp[0]], y=[Yp[0]], z=[Zp[0]],
+            mode="markers",
+            marker=dict(size=6, color="green"),
+            name="Start"
+        )
+    )
+    fig.add_trace(
+        go.Scatter3d(
+            x=[Xp[-1]], y=[Yp[-1]], z=[Zp[-1]],
+            mode="markers",
+            marker=dict(size=6, color="blue"),
+            name="Ende"
+        )
+    )
+
+    fig.update_layout(
+        title="ES-CMA",
+        scene=dict(
+            xaxis_title="X [m]",
+            yaxis_title="Y [m]",
+            zaxis_title="Z [m]",
+            aspectmode="data"
+        )
+    )
+
+    fig.show()
+
+plot_geodesic_3d(ell, betas, lambs)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

GHA/Hallo Welt!.py

@@ -1 +0,0 @@
-print("Hallo Welt!")

GHA/rk.py

@@ -0,0 +1,17 @@
+import Numerische_Integration.num_int_runge_kutta as rk
+from numpy import sin, cos, tan
+import winkelumrechnungen as wu
+from ellipsoide import EllipsoidBiaxial
+
+def gha1(re: EllipsoidBiaxial, x0, y0, z0, A0, s, num):
+    phi0, lamb0, h0 = re.cart2ell(0.001, wu.gms2rad([0, 0, 0.001]), x0, y0, z0)
+
+    f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c
+    f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c)
+    f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c
+
+    funktionswerte = rk.verfahren([f_phi, f_lam, f_A],
+                                  [0, phi0, lamb0, A0],
+                                  s, num)
+    coords = re.ell2cart(funktionswerte[-1][1], funktionswerte[-1][2], h0)
+    return coords

GHA_triaxial/approx_louville.py

@@ -0,0 +1,108 @@
+import numpy as np
+from numpy import sin, cos, arcsin, arccos, arctan2
+from ellipsoide import EllipsoidTriaxial
+import matplotlib.pyplot as plt
+
+dbeta_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ey**2 * sin(alpha)**2 * cos(beta) * sin(beta)
+dlamb_dc = lambda ell, beta, lamb, alpha: -2 * ell.Ee**2 * cos(alpha)**2 * sin(lamb) * cos(lamb)
+dalpha_dc = lambda ell, beta, lamb, alpha:  (2 * sin(alpha) * cos(alpha) *
+                                             (ell.Ey**2 * cos(beta)**2 + ell.Ee**2 * sin(lamb)**2))
+
+lamb2_sphere = lambda r, phi1, lamb1, a12, s: lamb1 + arctan2(sin(s/r) * sin(a12),
+                                                              cos(phi1) * cos(s/r) - sin(s/r) * cos(a12))
+phi2_sphere = lambda r, phi1, lamb1, a12, s: arcsin(sin(phi1) * cos(s/r) + cos(phi1) * sin(s/r) * cos(a12))
+
+a12_sphere = lambda phi1, lamb1, phi2, lamb2: arctan2(cos(phi2) * sin(lamb2 - lamb1),
+                                                      cos(phi1) * cos(phi2) -
+                                                      sin(phi1) * cos(phi2) * cos(lamb2 - lamb1))
+s_sphere = lambda r, phi1, lamb1, phi2, lamb2: r * arccos(sin(phi1) * sin(phi2) +
+                                                          cos(phi1) * cos(phi2) * cos(lamb2 - lamb1))
+
+louville = lambda beta, lamb, alpha: (ell.Ey**2 * cos(beta)**2 * sin(alpha)**2 -
+                                          ell.Ee**2 * sin(lamb)**2 * cos(alpha)**2)
+
+def points_approx_gha2(r: float, phi1: np.ndarray, lamb1: np.ndarray, phi2: np.ndarray, lamb2: np.ndarray, num: int = None, step_size: float = 10000):
+    s_approx = s_sphere(r, phi1, lamb1, phi2, lamb2)
+    if num is not None:
+        step_size = s_approx / (num+1)
+    a_approx = a12_sphere(phi1, lamb1, phi2, lamb2)
+
+    points = [np.array([phi1, lamb1, a_approx])]
+    current_s = step_size
+    while current_s < s_approx:
+        phi_n = phi2_sphere(r, phi1, lamb1, a_approx, current_s)
+        lamb_n = lamb2_sphere(r, phi1, lamb1, a_approx, current_s)
+        points.append(np.array([phi_n, lamb_n, a_approx]))
+        current_s += step_size
+    points.append(np.array([phi2, lamb2, a_approx]))
+    return points
+
+
+def num_update(ell: EllipsoidTriaxial, points, diffs):
+    for i, (beta, lamb, alpha) in enumerate(points):
+        dalpha = dalpha_dc(ell, beta, lamb, alpha)
+        if i == 0 or i == len(points) - 1:
+            grad = np.array([0, 0, dalpha])
+        else:
+            dbeta = dbeta_dc(ell, beta, lamb, alpha)
+            dlamb = dlamb_dc(ell, beta, lamb, alpha)
+            grad = np.array([dbeta, dlamb, dalpha])
+
+        delta = -diffs[i] * grad / np.dot(grad, grad)
+        points[i] += delta
+    return points
+
+
+def gha2(ell: EllipsoidTriaxial, p1: np.ndarray, p2: np.ndarray, maxI: int):
+    beta1, lamb1 = ell.cart2ell(p1)
+    beta2, lamb2 = ell.cart2ell(p2)
+    points = points_approx_gha2(ell.ax, beta1, lamb1, beta2, lamb2, 5)
+
+    for j in range(maxI):
+        constants = [louville(point[0], point[1], point[2]) for point in points]
+        mean_constant = np.mean(constants)
+        diffs = constants - mean_constant
+        if np.mean(np.abs(diffs)) > 10:
+            points = num_update(ell, points, diffs)
+        else:
+            break
+        for k in range(maxI):
+
+
+            last_diff_alpha = points[-2][-1] - points[-3][-1]
+            alpha_extrap = points[-2][-1] + last_diff_alpha
+            if abs(alpha_extrap - points[-1][-1]) > 0.0005:
+                pass
+            else:
+                break
+            pass
+    pass
+    return points
+
+def show_points(ell: EllipsoidTriaxial, points):
+    points_cart = []
+    for point in points:
+        points_cart.append(ell.ell2cart(point[0], point[1]))
+    points_cart = np.array(points_cart)
+
+    fig = plt.figure()
+    ax = fig.add_subplot(111, projection='3d')
+
+    ax.plot(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2])
+    ax.scatter(points_cart[:, 0], points_cart[:, 1], points_cart[:, 2])
+
+    ax.set_xlabel('X')
+    ax.set_ylabel('Y')
+    ax.set_zlabel('Z')
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("Eitschberger1978")
+    p1 = np.array([4189000, 812000, 4735000])
+    p2 = np.array([4090000, 868000, 4808000])
+    p1, phi1, lamb1, h1 = ell.cartonell(p1)
+    p2, phi2, lamb2, h2 = ell.cartonell(p2)
+    points = gha2(ell, p1, p2, 10)
+    show_points(ell, points)

GHA_triaxial/numeric_examples_panou.py

@@ -0,0 +1,111 @@
+import winkelumrechnungen as wu
+
+table1 = [
+    (wu.deg2rad(0),  wu.deg2rad(0),   wu.deg2rad(0),   wu.deg2rad(90),   1.00000000000,
+     wu.gms2rad([90,0,0.0000]),     wu.gms2rad([90,0,0.0000]),     10018754.9569),
+
+    (wu.deg2rad(1),  wu.deg2rad(0),   wu.deg2rad(-80), wu.deg2rad(5),    0.05883743460,
+     wu.gms2rad([179,7,12.2719]),   wu.gms2rad([174,40,13.8487]),  8947130.7221),
+
+    (wu.deg2rad(5),  wu.deg2rad(0),   wu.deg2rad(-60), wu.deg2rad(40),   0.34128138370,
+     wu.gms2rad([160,13,24.5001]),  wu.gms2rad([137,26,47.0036]),  8004762.4330),
+
+    (wu.deg2rad(30), wu.deg2rad(0),   wu.deg2rad(-30), wu.deg2rad(175),  0.86632464962,
+     wu.gms2rad([91,7,30.9337]),    wu.gms2rad([91,7,30.8672]),    19547128.7971),
+
+    (wu.deg2rad(60), wu.deg2rad(0),   wu.deg2rad(60),  wu.deg2rad(175),  0.06207487624,
+     wu.gms2rad([2,52,26.2393]),    wu.gms2rad([177,4,13.6373]),   6705715.1610),
+
+    (wu.deg2rad(75), wu.deg2rad(0),   wu.deg2rad(80),  wu.deg2rad(120),  0.11708984898,
+     wu.gms2rad([23,20,34.7823]),   wu.gms2rad([140,55,32.6385]),  2482501.2608),
+
+    (wu.deg2rad(80), wu.deg2rad(0),   wu.deg2rad(60),  wu.deg2rad(90),   0.17478427424,
+     wu.gms2rad([72,26,50.4024]),   wu.gms2rad([159,38,30.3547]),  3519745.1283)
+]
+
+table2 = [
+    (wu.deg2rad(0),  wu.deg2rad(-90), wu.deg2rad(0),  wu.deg2rad(89.5), 1.00000000000,
+     wu.gms2rad([90,0,0.0000]),      wu.gms2rad([90,0,0.0000]),     19981849.8629),
+
+    (wu.deg2rad(1),  wu.deg2rad(-90), wu.deg2rad(1),  wu.deg2rad(89.5), 0.18979826428,
+     wu.gms2rad([10,56,33.6952]),    wu.gms2rad([169,3,26.4359]),   19776667.0342),
+
+    (wu.deg2rad(5),  wu.deg2rad(-90), wu.deg2rad(5),  wu.deg2rad(89),   0.09398403161,
+     wu.gms2rad([5,24,48.3899]),     wu.gms2rad([174,35,12.6880]),  18889165.0873),
+
+    (wu.deg2rad(30), wu.deg2rad(-90), wu.deg2rad(30), wu.deg2rad(86),   0.06004022935,
+     wu.gms2rad([3,58,23.8038]),     wu.gms2rad([176,2,7.2825]),    13331814.6078),
+
+    (wu.deg2rad(60), wu.deg2rad(-90), wu.deg2rad(60), wu.deg2rad(78),   0.06076096484,
+     wu.gms2rad([6,56,46.4585]),     wu.gms2rad([173,11,5.9592]),   6637321.6350),
+
+    (wu.deg2rad(75), wu.deg2rad(-90), wu.deg2rad(75), wu.deg2rad(66),   0.05805851008,
+     wu.gms2rad([12,40,34.9009]),    wu.gms2rad([168,20,26.7339]),  3267941.2812),
+
+    (wu.deg2rad(80), wu.deg2rad(-90), wu.deg2rad(80), wu.deg2rad(55),   0.05817384452,
+     wu.gms2rad([18,35,40.7848]),    wu.gms2rad([164,25,34.0017]),  2132316.9048)
+]
+
+table3 = [
+    (wu.deg2rad(0),   wu.deg2rad(0.5),  wu.deg2rad(80),  wu.deg2rad(0.5),  0.05680316848,
+     wu.gms2rad([0,-0,16.0757]),      wu.gms2rad([0,1,32.5762]),   8831874.3717),
+
+    (wu.deg2rad(-1),  wu.deg2rad(5),    wu.deg2rad(75),  wu.deg2rad(5),    0.05659149555,
+     wu.gms2rad([0,-1,47.2105]),      wu.gms2rad([0,6,54.0958]),   8405370.4947),
+
+    (wu.deg2rad(-5),  wu.deg2rad(30),   wu.deg2rad(60),  wu.deg2rad(30),   0.04921108945,
+     wu.gms2rad([0,-4,22.3516]),      wu.gms2rad([0,8,42.0756]),   7204083.8568),
+
+    (wu.deg2rad(-30), wu.deg2rad(45),   wu.deg2rad(30),  wu.deg2rad(45),   0.04017812574,
+     wu.gms2rad([0,-3,41.2461]),      wu.gms2rad([0,3,41.2461]),   6652788.1287),
+
+    (wu.deg2rad(-60), wu.deg2rad(60),   wu.deg2rad(5),   wu.deg2rad(60),   0.02843082609,
+     wu.gms2rad([0,-8,40.4575]),      wu.gms2rad([0,4,22.1675]),   7213412.4477),
+
+    (wu.deg2rad(-75), wu.deg2rad(85),   wu.deg2rad(1),   wu.deg2rad(85),   0.00497802414,
+     wu.gms2rad([0,-6,44.6115]),      wu.gms2rad([0,1,47.0474]),   8442938.5899),
+
+    (wu.deg2rad(-80), wu.deg2rad(89.5), wu.deg2rad(0),   wu.deg2rad(89.5), 0.00050178253,
+     wu.gms2rad([0,-1,27.9705]),      wu.gms2rad([0,0,16.0490]),   8888783.7815)
+]
+
+table4 = [
+    (wu.deg2rad(0),   wu.deg2rad(0),   wu.deg2rad(0),   wu.deg2rad(90),    1.00000000000,
+     wu.gms2rad([90,0,0.0000]),      wu.gms2rad([90,0,0.0000]),    10018754.1714),
+
+    (wu.deg2rad(1),   wu.deg2rad(0),   wu.deg2rad(0),   wu.deg2rad(179.5), 0.30320665822,
+     wu.gms2rad([17,39,11.0942]),    wu.gms2rad([162,20,58.9032]), 19884417.8083),
+
+    (wu.deg2rad(5),   wu.deg2rad(0),   wu.deg2rad(-80), wu.deg2rad(170),   0.03104258442,
+     wu.gms2rad([178,12,51.5083]),   wu.gms2rad([10,17,52.6423]),  11652530.7514),
+
+    (wu.deg2rad(30),  wu.deg2rad(0),   wu.deg2rad(-75), wu.deg2rad(120),   0.24135347134,
+     wu.gms2rad([163,49,4.4615]),    wu.gms2rad([68,49,50.9617]),  14057886.8752),
+
+    (wu.deg2rad(60),  wu.deg2rad(0),   wu.deg2rad(-60), wu.deg2rad(40),    0.19408499032,
+     wu.gms2rad([157,9,33.5589]),    wu.gms2rad([157,9,33.5589]),  13767414.8267),
+
+    (wu.deg2rad(75),  wu.deg2rad(0),   wu.deg2rad(-30), wu.deg2rad(0.5),   0.00202789418,
+     wu.gms2rad([179,33,3.8613]),    wu.gms2rad([179,51,57.0077]), 11661713.4496),
+
+    (wu.deg2rad(80),  wu.deg2rad(0),   wu.deg2rad(-5),  wu.deg2rad(120),   0.15201222384,
+     wu.gms2rad([61,5,33.9600]),     wu.gms2rad([171,13,22.0148]), 11105138.2902),
+
+    (wu.deg2rad(0),   wu.deg2rad(0),   wu.deg2rad(60),  wu.deg2rad(0),     0.00000000000,
+     wu.gms2rad([0,0,0.0000]),       wu.gms2rad([0,0,0.0000]),     6663348.2060)
+]
+
+tables = [table1, table2, table3, table4]
+
+def get_example(table, example):
+    table -= 1
+    example -= 1
+    return tables[table][example]
+
+def get_tables():
+    return tables
+
+
+if __name__ == "__main__":
+    test = get_example(1, 4)
+    pass

GHA_triaxial/panou.py

@@ -0,0 +1,191 @@
+import numpy as np
+from numpy import sin, cos, sqrt, arctan2
+import ellipsoide
+import Numerische_Integration.num_int_runge_kutta as rk
+import winkelumrechnungen as wu
+import ausgaben as aus
+import GHA.rk as ghark
+from scipy.special import factorial as fact
+from math import comb
+import GHA_triaxial.numeric_examples_panou as nep
+
+# Panou, Korakitits 2019
+
+
+def gha1_num_old(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, num):
+    phi, lamb, h = ell.cart2geod(point, "ligas3")
+    x, y, z = ell.geod2cart(phi, lamb, 0)
+    p, q = ell.p_q(x, y, z)
+
+    dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0)
+    dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0)
+    dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0)
+
+    h = lambda dxds, dyds, dzds: dxds**2 + 1/(1-ell.ee**2)*dyds**2 + 1/(1-ell.ex**2)*dzds**2
+
+    f1 = lambda x, dxds, y, dyds, z, dzds: dxds
+    f2 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * x
+    f3 = lambda x, dxds, y, dyds, z, dzds: dyds
+    f4 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * y/(1-ell.ee**2)
+    f5 = lambda x, dxds, y, dyds, z, dzds: dzds
+    f6 = lambda x, dxds, y, dyds, z, dzds: -h(dxds, dyds, dzds) / ell.func_H(x, y, z) * z/(1-ell.ex**2)
+
+    funktionswerte = rk.verfahren([f1, f2, f3, f4, f5, f6], [x, dxds0, y, dyds0, z, dzds0], s, num, fein=False)
+    P2 = funktionswerte[-1]
+    P2 = (P2[0], P2[2], P2[4])
+    return P2
+
+def buildODE(ell):
+    def ODE(v):
+        x, dxds, y, dyds, z, dzds = v
+
+        H = ell.func_H(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 [dxds, ddx, dyds, ddy, dzds, ddz]
+    return ODE
+
+def gha1_num(ell, point, alpha0, s, num):
+    phi, lam, _ = ell.cart2geod(point, "ligas3")
+    x0, y0, z0 = ell.geod2cart(phi, lam, 0)
+
+    p, q = ell.p_q(x0, y0, z0)
+    dxds0 = p[0] * sin(alpha0) + q[0] * cos(alpha0)
+    dyds0 = p[1] * sin(alpha0) + q[1] * cos(alpha0)
+    dzds0 = p[2] * sin(alpha0) + q[2] * cos(alpha0)
+
+    v_init = [x0, dxds0, y0, dyds0, z0, dzds0]
+
+    F = buildODE(ell)
+
+    werte = rk.rk_chat(F, v_init, s, num)
+    x1, _, y1, _, z1, _ = werte[-1]
+
+    return x1, y1, z1
+
+
+def checkLiouville(ell: ellipsoide.EllipsoidTriaxial, points):
+    constantValues = []
+    for point in points:
+        x = point[1]
+        dxds = point[2]
+        y = point[3]
+        dyds = point[4]
+        z = point[5]
+        dzds = point[6]
+
+        values = ell.p_q(x, y, z)
+        p = values["p"]
+        q = values["q"]
+        t1 = values["t1"]
+        t2 = values["t2"]
+
+        P = p[0]*dxds + p[1]*dyds + p[2]*dzds
+        Q = q[0]*dxds + q[1]*dyds + q[2]*dzds
+        alpha = arctan2(P, Q)
+
+        c = ell.ay**2 - (t1 * sin(alpha)**2 + t2 * cos(alpha)**2)
+        constantValues.append(c)
+    pass
+
+
+def gha1_ana(ell: ellipsoide.EllipsoidTriaxial, point, alpha0, s, maxM):
+    """
+    Panou, Korakitits 2020, 5ff.
+    :param ell:
+    :param point:
+    :param alpha0:
+    :param s:
+    :param maxM:
+    :return:
+    """
+    x, y, z = point
+    x_m = [x]
+    y_m = [y]
+    z_m = [z]
+
+    # erste Ableitungen (7-8)
+    H = x ** 2 + y ** 2 / (1 - ell.ee ** 2) ** 2 + z ** 2 / (1 - ell.ex ** 2) ** 2
+    sqrtH = sqrt(H)
+    n = np.array([x / sqrtH,
+                  y / ((1-ell.ee**2) * sqrtH),
+                  z / ((1-ell.ex**2) * sqrtH)])
+    u, v = ell.cart2para(np.array([x, y, z]))
+    G = sqrt(1 - ell.ex**2 * cos(u)**2 - ell.ee**2 * sin(u)**2 * sin(v)**2)
+    q = np.array([-1/G * sin(u) * cos(v),
+                  -1/G * sqrt(1-ell.ee**2) * sin(u) * sin(v),
+                  1/G * sqrt(1-ell.ex**2) * cos(u)])
+    p = np.array([q[1]*n[2] - q[2]*n[1],
+                  q[2]*n[0] - q[0]*n[2],
+                  q[0]*n[1] - q[1]*n[0]])
+    x_m.append(p[0] * sin(alpha0) + q[0] * cos(alpha0))
+    y_m.append(p[1] * sin(alpha0) + q[1] * cos(alpha0))
+    z_m.append(p[2] * sin(alpha0) + q[2] * cos(alpha0))
+
+    # H Ableitungen (7)
+    H_ = lambda p: np.sum([comb(p, i) * (x_m[p - i] * x_m[i] +
+                                         1 / (1-ell.ee**2) ** 2 * y_m[p-i] * y_m[i] +
+                                         1 / (1-ell.ex**2) ** 2 * z_m[p-i] * z_m[i]) for i in range(0, p+1)])
+
+    # h Ableitungen (7)
+    h_ = lambda q: np.sum([comb(q, j) * (x_m[q-j+1] * x_m[j+1] +
+                                         1 / (1 - ell.ee ** 2) * y_m[q-j+1] * y_m[j+1] +
+                                         1 / (1 - ell.ex ** 2) * z_m[q-j+1] * z_m[j+1]) for j in range(0, q+1)])
+
+    # h/H Ableitungen (6)
+    hH_ = lambda t: 1/H_(0) * (h_(t) - fact(t) *
+                               np.sum([H_(t+1-l) / (fact(t+1-l) * fact(l-1)) * hH_t[l-1] for l in range(1, t+1)]))
+
+    # xm, ym, zm Ableitungen (6)
+    x_ = lambda m: -np.sum([comb(m-2, k) * hH_t[m-2-k] * x_m[k] for k in range(0, m-2+1)])
+    y_ = lambda m: -1 / (1-ell.ee**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * y_m[k] for k in range(0, m-2+1)])
+    z_ = lambda m: -1 / (1-ell.ex**2) * np.sum([comb(m-2, k) * hH_t[m-2-k] * z_m[k] for k in range(0, m-2+1)])
+
+    hH_t = []
+    a_m = []
+    b_m = []
+    c_m = []
+    for m in range(0, maxM+1):
+        if m >= 2:
+            hH_t.append(hH_(m-2))
+            x_m.append(x_(m))
+            y_m.append(y_(m))
+            z_m.append(z_(m))
+        a_m.append(x_m[m] / fact(m))
+        b_m.append(y_m[m] / fact(m))
+        c_m.append(z_m[m] / fact(m))
+
+    # am, bm, cm (6)
+    x_s = 0
+    for a in reversed(a_m):
+        x_s = x_s * s + a
+    y_s = 0
+    for b in reversed(b_m):
+        y_s = y_s * s + b
+    z_s = 0
+    for c in reversed(c_m):
+        z_s = z_s * s + c
+
+    return x_s, y_s, z_s
+    pass
+
+
+if __name__ == "__main__":
+    # ell = ellipsoide.EllipsoidTriaxial.init_name("Eitschberger1978")
+    ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
+    # ellbi = ellipsoide.EllipsoidTriaxial.init_name("Bessel-biaxial")
+    re = ellipsoide.EllipsoidBiaxial.init_name("Bessel")
+
+    # Panou 2013, 7, Table 1, beta0=60°
+    beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = nep.get_example(table=1, example=5)
+    P0 = ell.ell2cart(beta0, lamb0)
+    P1 = ell.ell2cart(beta1, lamb1)
+
+    # P1_num = gha1_num(ell, P0, alpha0, s, 1000)
+    P1_num = gha1_num(ell, P0, alpha0, s, 10000)
+    P1_ana = gha1_ana(ell, P0, alpha0, s, 30)
+    pass

GHA_triaxial/panou_2013_2GHA_num.py

@@ -0,0 +1,339 @@
+import numpy as np
+from ellipsoide import EllipsoidTriaxial
+import Numerische_Integration.num_int_runge_kutta as rk
+import ausgaben as aus
+
+# Panou 2013
+def gha2_num(ell: EllipsoidTriaxial, beta_1, lamb_1, beta_2, lamb_2, n=16000, epsilon=10**-12, iter_max=30):
+    """
+
+    :param ell: triaxiales Ellipsoid
+    :param beta_1: reduzierte ellipsoidische Breite Punkt 1
+    :param lamb_1: elllipsoidische Länge Punkt 1
+    :param beta_2: reduzierte ellipsoidische Breite Punkt 2
+    :param lamb_2: elllipsoidische Länge Punkt 2
+    :param n: Anzahl Schritte
+    :param epsilon:
+    :param iter_max: Maximale Anzhal Iterationen
+    :return:
+    """
+
+    # h_x, h_y, h_e entsprechen E_x, E_y, E_e
+
+    def arccot(x):
+        return np.arctan2(1.0, x)
+
+
+    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]
+
+        N = n
+
+        dlamb = lamb_2 - lamb_1
+
+        if abs(dlamb) < 1e-15:
+            beta_0 = 0.0
+        else:
+            beta_0 = (beta_2 - beta_1) / (lamb_2 - lamb_1)
+
+        converged = False
+        iterations = 0
+
+        funcs = functions()
+
+        for i in range(iter_max):
+            iterations = i + 1
+
+            startwerte = [lamb_1, beta_1, beta_0, 0.0, 1.0]
+
+            werte = rk.verfahren(funcs, startwerte, dlamb, N)
+            lamb_end, beta_end, beta_p_end, X3_end, X4_end = werte[-1]
+
+            d_beta_end_d_beta0 = X3_end
+            delta = beta_end - beta_2
+
+            if abs(delta) < epsilon:
+                converged = True
+                break
+
+            if abs(d_beta_end_d_beta0) < 1e-20:
+                raise RuntimeError("Abbruch.")
+
+            max_step = 0.5
+            step = delta / d_beta_end_d_beta0
+            if abs(step) > max_step:
+                step = np.sign(step) * max_step
+            beta_0 = beta_0 - step
+
+        if not converged:
+            raise RuntimeError("konvergiert nicht.")
+
+        # Z
+        werte = rk.verfahren(funcs, [lamb_1, beta_1, beta_0, 0.0, 1.0], dlamb, N)
+
+        beta_arr = np.zeros(N + 1)
+        lamb_arr = np.zeros(N + 1)
+        beta_p_arr = np.zeros(N + 1)
+
+        for i, state in enumerate(werte):
+            lamb_arr[i] = state[0]
+            beta_arr[i] = state[1]
+            beta_p_arr[i] = state[2]
+
+        (_, _, E1, G1,
+         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+        (_, _, E2, G2,
+         *_) = 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])
+
+        integrand = np.zeros(N + 1)
+        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)
+
+        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
+        else:
+            s = np.trapz(integrand, dx=h)
+
+        beta0 = beta_arr[0]
+        lamb0 = lamb_arr[0]
+        c = np.sqrt(
+            (np.cos(beta0) ** 2 + (ell.Ee**2 / ell.Ex**2) * np.sin(beta0) ** 2) * np.sin(alpha_1) ** 2
+            + (ell.Ee**2 / ell.Ex**2) * np.cos(lamb0) ** 2 * np.cos(alpha_1) ** 2
+        )
+
+        return alpha_1, alpha_2, s
+
+    if lamb_1 == lamb_2:
+
+        N = n
+        dbeta = beta_2 - beta_1
+
+        if abs(dbeta) < 10**-15:
+            return 0, 0, 0
+
+        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]
+
+        funcs_beta = functions_beta()
+
+        for i in range(iter_max):
+            iterations = i + 1
+
+            startwerte = [beta_1, lamb_1, lamb_0, 0.0, 1.0]
+
+            werte = rk.verfahren(funcs_beta, startwerte, dbeta, N)
+            beta_end, 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:
+                raise RuntimeError("Abbruch (Ableitung ~ 0).")
+
+            max_step = 1.0
+            step = delta / d_lamb_end_d_lambda0
+            if abs(step) > max_step:
+                step = np.sign(step) * max_step
+
+            lamb_0 = lamb_0 - step
+
+        werte = rk.verfahren(funcs_beta, [beta_1, lamb_1, lamb_0, 0.0, 1.0], dbeta, N)
+
+        beta_arr = np.zeros(N + 1)
+        lamb_arr = np.zeros(N + 1)
+        lambda_p_arr = np.zeros(N + 1)
+
+        for i, state in enumerate(werte):
+            beta_arr[i] = state[0]
+            lamb_arr[i] = state[1]
+            lambda_p_arr[i] = state[2]
+
+        # Azimute
+        (BETA1, LAMBDA1, E1, G1,
+         *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+        (BETA2, LAMBDA2, E2, G2,
+         *_) = 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])
+
+        integrand = np.zeros(N + 1)
+        for i in range(N + 1):
+            (_, _, Ei, Gi,
+             *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
+            integrand[i] = np.sqrt(Ei + Gi * lambda_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
+        else:
+            s = np.trapz(integrand, dx=h)
+
+        return alpha_1, alpha_2, s
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    # beta1 = np.deg2rad(75)
+    # lamb1 = np.deg2rad(-90)
+    # beta2 = np.deg2rad(75)
+    # lamb2 = np.deg2rad(66)
+    # a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2)
+    # print(aus.gms("a1", a1, 4))
+    # print(aus.gms("a2", a2, 4))
+    # print(s)
+    cart1 = ell.para2cart(0, 0)
+    cart2 = ell.para2cart(0.4, 0.4)
+    beta1, lamb1 = ell.cart2ell(cart1)
+    beta2, lamb2 = ell.cart2ell(cart2)
+
+    a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=2500)
+    print(s)
+
+
+

Hansen_ES_CMA.py

@@ -0,0 +1,119 @@
+import numpy as np
+
+
+def felli(x):
+    N = x.shape[0]
+    if N < 2:
+        raise ValueError("dimension must be greater than one")
+    exponents = np.arange(N) / (N - 1)
+    return float(np.sum((1e6 ** exponents) * (x ** 2)))
+
+
+def escma():
+    #Initialization
+
+    # User defined input parameters
+    N = 10                       
+    xmean = np.random.rand(N)
+    sigma = 0.5
+    stopfitness = 1e-10
+    stopeval = int(1e3 * N**2)
+
+    # Strategy parameter setting: Selection
+    lambda_ = 4 + int(np.floor(3 * np.log(N)))
+    mu = lambda_ / 2.0
+
+    # muXone recombination weights
+    weights = np.log(mu + 0.5) - np.log(np.arange(1, int(mu) + 1))
+    mu = int(np.floor(mu))
+    weights = weights / np.sum(weights)
+    mueff = np.sum(weights)**2 / np.sum(weights**2)
+
+    # Strategy parameter setting: Adaptation
+    cc = (4 + mueff / N) / (N + 4 + 2 * mueff / N)
+    cs = (mueff + 2) / (N + mueff + 5)
+    c1 = 2 / ((N + 1.3)**2 + mueff)
+    cmu = min(1 - c1,
+              2 * (mueff - 2 + 1 / mueff) / ((N + 2)**2 + 2 * mueff))
+    damps = 1 + 2 * max(0, np.sqrt((mueff - 1) / (N + 1)) - 1) + cs
+
+    # Initialize dynamic (internal) strategy parameters and constants
+    pc = np.zeros(N)
+    ps = np.zeros(N)
+    B = np.eye(N)
+    D = np.eye(N)
+    C = B @ D @ (B @ D).T
+    eigeneval = 0
+    chiN = np.sqrt(N) * (1 - 1/(4*N) + 1/(21 * N**2))
+
+    #Generation Loop
+    counteval = 0
+    arx = np.zeros((N, lambda_))
+    arz = np.zeros((N, lambda_))
+    arfitness = np.zeros(lambda_)
+
+    while counteval < stopeval:
+
+        # Generate and evaluate lambda offspring
+        for k in range(lambda_):
+            arz[:, k] = np.random.randn(N)
+            arx[:, k] = xmean + sigma * (B @ D @ arz[:, k])
+            arfitness[k] = felli(arx[:, k])
+            counteval += 1
+
+        # Sort by fitness and compute weighted mean into xmean
+        idx = np.argsort(arfitness)
+        arfitness = arfitness[idx]
+        arindex = idx
+
+        xold = xmean.copy()
+        xmean = arx[:, arindex[:mu]] @ weights
+        zmean = arz[:, arindex[:mu]] @ weights
+
+        # Cumulation: Update evolution paths
+        ps = (1 - cs) * ps + np.sqrt(cs * (2 - cs) * mueff) * (B @ zmean)
+        norm_ps = np.linalg.norm(ps)
+        hsig = norm_ps / np.sqrt(1 - (1 - cs)**(2 * counteval / lambda_)) / chiN < \
+               (1.4 + 2 / (N + 1))
+        hsig = 1.0 if hsig else 0.0
+
+        pc = (1 - cc) * pc + hsig * np.sqrt(cc * (2 - cc) * mueff) * (B @ D @ zmean)
+
+        # Adapt covariance matrix C
+        BDz = B @ D @ arz[:, arindex[:mu]]
+        C = (1 - c1 - cmu) * C \
+            + c1 * (np.outer(pc, pc) + (1 - hsig) * cc * (2 - cc) * C) \
+            + cmu * BDz @ np.diag(weights) @ BDz.T
+
+        # Adapt step-size sigma
+        sigma = sigma * np.exp((cs / damps) * (norm_ps / chiN - 1))
+
+        # Update B and D from C (Eigenzerlegung, O(N^2))
+        if counteval - eigeneval > lambda_ / ((c1 + cmu) * N * 10):
+            eigeneval = counteval
+            # enforce symmetry
+            C = (C + C.T) / 2.0
+            eigvals, B = np.linalg.eigh(C)
+            D = np.diag(np.sqrt(eigvals))
+
+        # Break, if fitness is good enough
+        if arfitness[0] <= stopfitness:
+            break
+
+        # Escape flat fitness, or better terminate?
+        if arfitness[0] == arfitness[int(np.ceil(0.7 * lambda_)) - 1]:
+            sigma = sigma * np.exp(0.2 + cs / damps)
+            print("warning: flat fitness, consider reformulating the objective")
+
+        print(f"{counteval}: {arfitness[0]}")
+
+    #Final Message
+    print(f"{counteval}: {arfitness[0]}")
+    xmin = arx[:, arindex[0]]
+    return xmin
+
+
+if __name__ == "__main__":
+    xmin = escma()
+    print("Bestes gefundenes x:", xmin)
+    print("f(xmin) =", felli(xmin))

Numerische_Integration/num_int_runge_kutta.py

@@ -1,4 +1,4 @@
-def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) -> list:
+def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int, fein: bool = True) -> list:
     """
     Runge-Kutta-Verfahren für ein beliebiges DGLS
     :param funktionen: Liste mit allen Funktionen
@@ -14,19 +14,21 @@ def verfahren(funktionen: list, startwerte: list, weite: float, schritte: int) -
         zuschlaege_grob = zuschlaege(funktionen, werte[-1], h)
         werte_grob = [werte[-1][j] if j == 0 else werte[-1][j] + zuschlaege_grob[j - 1]
                       for j in range(len(startwerte))]
+        if fein:
+            zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2)
+            werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1]
+                            for j in range(len(startwerte))]
 
-        zuschlaege_fein_1 = zuschlaege(funktionen, werte[-1], h / 2)
+            zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2)
-        werte_fein_1 = [werte[-1][j] + h/2 if j == 0 else werte[-1][j]+zuschlaege_fein_1[j-1]
+            werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1]
-                        for j in range(len(startwerte))]
+                            for j in range(len(startwerte))]
 
-        zuschlaege_fein_2 = zuschlaege(funktionen, werte_fein_1, h / 2)
+            werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j])
-        werte_fein_2 = [werte_fein_1[j] + h/2 if j == 0 else werte_fein_1[j]+zuschlaege_fein_2[j-1]
+                          for j in range(len(startwerte))]
-                        for j in range(len(startwerte))]
 
-        werte_korr = [werte_fein_2[j] if j == 0 else werte_fein_2[j] + 1/15 * (werte_fein_2[j] - werte_grob[j])
+            werte.append(werte_korr)
-                      for j in range(len(startwerte))]
+        else:
-
+            werte.append(werte_grob)
-        werte.append(werte_korr)
     return werte
 
 
@@ -60,3 +62,23 @@ def zuschlaege(funktionen: list, startwerte: list, h: float) -> list:
     k_ = [(k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6 for i in range(len(k1))]
 
     return k_
+
+def rk_chat(F, v0: list, weite: float, schritte: int):
+    h = weite/schritte
+    v = v0
+    werte = [v]
+
+    for _ in range(schritte):
+        k1 = F(v)
+        k2 = F([v[i] + 0.5 * h * k1[i] for i in range(6)])
+        k3 = F([v[i] + 0.5 * h * k2[i] for i in range(6)])
+        k4 = F([v[i] + h * k3[i] for i in range(6)])
+
+        v = [
+            v[i] + (h / 6) * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i])
+            for i in range(6)
+        ]
+
+        werte.append(v)
+
+    return werte

Richtungsreduktion.py

@@ -1,14 +0,0 @@
-from numpy import sqrt, cos, sin
-import ellipsoide
-import ausgaben as aus
-import winkelumrechnungen as wu
-
-
-re = ellipsoide.Ellipsoid.init_name("Bessel")
-eta = lambda phi: sqrt(re.e_ ** 2 * cos(phi) ** 2)
-
-A = wu.gms2rad([327,0,0])
-phi = wu.gms2rad([35,0,0])
-h = 3500
-dA = eta(phi)**2 * h * sin(A)*cos(A) / re.N(phi)
-print(aus.gms("dA", dA, 10))

ausgaben.py

@@ -10,9 +10,7 @@ def xyz(x: float, y: float, z: float, stellen: int) -> str:
     :param stellen: Anzahl Nachkommastellen
     :return: String zur Ausgabe der Koordinaten
     """
-    return f"""x = {(round(x,stellen))} m
+    return f"""x = {(round(x,stellen))} m  y = {(round(y,stellen))} m  z = {(round(z,stellen))} m"""
-y = {(round(y,stellen))} m
-z = {(round(z,stellen))} m"""
 
 
 def gms(name: str, rad: float, stellen: int) -> str:

dashboard.py

@@ -0,0 +1,516 @@
+from dash import Dash, html, dcc, Input, Output, State, no_update
+import dash
+import plotly.graph_objects as go
+import numpy as np
+
+from GHA_triaxial.panou import gha1_ana
+from GHA_triaxial.panou import gha1_num
+from GHA_triaxial.panou_2013_2GHA_num import gha2_num
+from ellipsoide import EllipsoidTriaxial
+import winkelumrechnungen as wu
+import ausgaben as aus
+
+
+app = Dash(__name__, suppress_callback_exceptions=True)
+app.title = "Geodätische Hauptaufgaben"
+
+
+def abplattung(a, b):
+    return (a - b) / a
+
+def ellipsoid_figure(ell: EllipsoidTriaxial, title="Dreiachsiges Ellipsoid"):
+    fig = go.Figure()
+
+    # Darstellung
+    rx, ry, rz = 1.05*ell.ax, 1.05*ell.ay, 1.05*ell.b
+    fig.update_layout(
+        title=title,
+        scene=dict(
+            xaxis=dict(range=[-rx, rx], title="X [m]"),
+            yaxis=dict(range=[-ry, ry], title="Y [m]"),
+            zaxis=dict(range=[-rz, rz], title="Z [m]"),
+            aspectmode="data"
+        ),
+        margin=dict(l=0, r=0, t=40, b=0),
+    )
+
+    # Ellipsoid
+    u = np.linspace(-np.pi/2, np.pi/2, 80)
+    v = np.linspace(-np.pi,     np.pi,   160)
+    U, V = np.meshgrid(u, v)
+    X, Y, Z = ell.para2cart(U, V)
+    fig.add_trace(go.Surface(
+        x=X, y=Y, z=Z, showscale=False, opacity=0.7,
+        surfacecolor=np.zeros_like(X),
+        colorscale=[[0, "rgb(200,220,255)"], [1, "rgb(200,220,255)"]],
+        name="Ellipsoid"
+    ))
+
+    return fig
+
+def figure_constant_lines(fig, ell: EllipsoidTriaxial, coordsystem: str = "para"):
+    if coordsystem == "para":
+        constants_u = wu.deg2rad(np.arange(0, 360, 15))
+        all_v = np.linspace(-np.pi / 2, np.pi / 2, 361)
+        for u in constants_u:
+            xm, ym, zm = ell.para2cart(u, all_v)
+            fig.add_trace(go.Scatter3d(
+                x=xm, y=ym, z=zm, mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+        all_u = np.linspace(0, 2 * np.pi, 361)
+        constants_v = wu.deg2rad(np.arange(-75, 90, 15))
+        for v in constants_v:
+            x, y, z = ell.para2cart(all_u, v)
+            fig.add_trace(go.Scatter3d(
+                x=x, y=y, z=z, mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+    elif coordsystem == "ell":
+        constants_beta = wu.deg2rad(np.arange(-75, 90, 15))
+        all_lamb = np.linspace(0, 2 * np.pi, 361)
+        for beta in constants_beta:
+            xyz = ell.ell2cart(beta, all_lamb)
+            fig.add_trace(go.Scatter3d(
+                x=xyz[:, 0], y=xyz[:, 1], z=xyz[:, 2], mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+        all_beta = np.linspace(-np.pi / 2, np.pi / 2, 361)
+        constants_lamb = wu.deg2rad(np.arange(0, 360, 15))
+        for lamb in constants_lamb:
+            xyz = ell.ell2cart(all_beta, lamb)
+            fig.add_trace(go.Scatter3d(
+                x=xyz[:, 0], y=xyz[:, 1], z=xyz[:, 2], mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+    elif coordsystem == "geod":
+        constants_phi = wu.deg2rad(np.arange(-75, 90, 15))
+        all_lamb = np.linspace(0, 2 * np.pi, 361)
+        for phi in constants_phi:
+            x, y, z = ell.geod2cart(phi, all_lamb, 0)
+            fig.add_trace(go.Scatter3d(
+                x=x, y=y, z=z, mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+        all_phi = np.linspace(-np.pi / 2, np.pi / 2, 361)
+        constants_lamb = wu.deg2rad(np.arange(0, 360, 15))
+        for lamb in constants_lamb:
+            x, y, z = ell.geod2cart(all_phi, lamb, 0)
+            fig.add_trace(go.Scatter3d(
+                x=x, y=y, z=z, mode="lines",
+                line=dict(width=1, color="black"),
+                showlegend=False
+            ))
+
+    return fig
+
+def figure_points(fig, points):
+    """
+
+    :param fig: plotly.graph_objects.Figure
+    :param points: Punktliste [(name, (x,y,z), color)]
+    :return: plotly.graph_objects.Figure
+    """
+    for name, (px, py, pz), color in points:
+        fig.add_trace(go.Scatter3d(
+            x=[px], y=[py], z=[pz],
+            mode="markers+text",
+            marker=dict(size=6, color=color),
+            text=[name], textposition="top center",
+            name=name, showlegend=False
+        ))
+    return fig
+
+def figure_lines(fig, lines):
+    """
+
+    :param fig: plotly.graph_objects.Figure
+    :param lines: Linienliste [((x1,y1,z1), (x2,y2,z2), color)]
+    :return: plotly.graph_objects.Figure
+    """
+    for (p1, p2, color) in lines:
+        xline = [p1[0], p2[0]]
+        yline = [p1[1], p2[1]]
+        zline = [p1[2], p2[2]]
+        fig.add_trace(go.Scatter3d(
+            x=xline, y=yline, z=zline,
+            mode="lines",
+            line=dict(width=4, color=color),
+            showlegend=False
+        ))
+    return fig
+
+
+
+app.layout = html.Div(
+    style={"fontFamily": "Arial", "padding": "5px", "width": "70%", "margin-left": "auto"},
+    children=[
+        html.H1("Geodätische Hauptaufgaben"),
+        html.H2("für dreiachsige Ellipsoide"),
+
+        html.Label("Ellipsoid wählen:"),
+        dcc.Dropdown(
+            id="dropdown-ellipsoid",
+            options=[
+                {"label": "BursaFialova1993", "value": "BursaFialova1993"},
+                {"label": "BursaSima1980", "value": "BursaSima1980"},
+                {"label": "BursaSima1980round", "value": "BursaSima1980round"},
+                {"label": "Eitschberger1978", "value": "Eitschberger1978"},
+                {"label": "Bursa1972", "value": "Bursa1972"},
+                {"label": "Bursa1970", "value": "Bursa1970"},
+                {"label": "BesselBiaxial", "value": "BesselBiaxial"},
+                {"label": "Fiction", "value": "Fiction"},
+                #{"label": "Ei", "value": "Ei"},
+            ],
+            value="",
+            style={"width": "300px", "marginBottom": "20px"},
+        ),
+
+        html.Label("Halbachsen:"),
+        dcc.Input(
+            id="input-ax",
+            type="number",
+            min=0,
+            placeholder="ax...[m]",
+            style={"marginBottom": "10px", "display": "block", "width": "300px"},
+        ),
+        dcc.Input(
+            id="input-ay",
+            type="number",
+            min=0,
+            placeholder="ay...[m]",
+            style={"marginBottom": "10px", "display": "block", "width": "300px"},
+        ),
+        dcc.Input(
+            id="input-b",
+            type="number",
+            min=0,
+            placeholder="b...[m]",
+            style={"marginBottom": "20px", "display": "block", "width": "300px"},
+        ),
+
+        html.Button(
+            "Ellipsoid berechnen",
+            id="calc-ell",
+            n_clicks=0,
+            style={"marginRight": "10px", "marginBottom": "20px"},
+        ),
+
+        html.Div(id="output-area", style={"marginBottom": "20px"}),
+
+        dcc.Tabs(
+            id="tabs-GHA",
+            value="tab-GHA1",
+            style={"marginRight": "10px", "marginBottom": "20px", "width": "50%"},
+            children=[
+                dcc.Tab(label="Erste Hauptaufgabe", value="tab-GHA1"),
+                dcc.Tab(label="Zweite Hauptaufgabe", value="tab-GHA2"),
+            ],
+        ),
+
+        html.Div(
+            id="tabs-GHA-out",
+            style={"marginRight": "10px", "marginBottom": "20px", "width": "50%"},
+        ),
+
+        html.Div(id="output-gha1", style={"marginBottom": "20px"}),
+
+        html.Div(id="output-gha2", style={"marginBottom": "20px"}),
+
+        dcc.Graph(
+            id="ellipsoid-plot",
+            style={"height": "500px", "width": "700px"},
+        ),
+
+        html.P(
+            "© 2025",
+            style={
+                "fontSize": "12px",
+                "color": "gray",
+                "textAlign": "center",
+            },
+        ),
+    ],
+)
+
+
+@app.callback(
+    Output("input-ax", "value"),
+    Output("input-ay", "value"),
+    Output("input-b", "value"),
+    Input("dropdown-ellipsoid", "value"),
+)
+def fill_inputs_from_dropdown(selected_ell):
+    if not selected_ell:
+        return None, None, None
+
+
+    ell = EllipsoidTriaxial.init_name(selected_ell)
+    ax = ell.ax
+    ay = ell.ay
+    b = ell.b
+    return ax, ay, b
+
+
+@app.callback(
+    Output("output-area", "children"),
+    Input("calc-ell", "n_clicks"),
+    State("input-ax", "value"),
+    State("input-ay", "value"),
+    State("input-b", "value"),
+)
+def update_output(n_clicks, ax, ay, b):
+    if not n_clicks:
+        return ""
+    if n_clicks and ax is None or ay is None or b is None:
+        return html.Span("Bitte Ellipsoid auswählen!", style={"color": "red"})
+    if ay >= ax or b >= ay or ax <= 0 or ay <= 0 or b <= 0:
+        return html.Span("Eingabe inkorrekt.", style={"color": "red"})
+    ell = EllipsoidTriaxial(ax, ay, b)
+    return f"ex = {round(ell.ex, 6)}, ", f"ey = {round(ell.ey, 6)}, ", f"ee = {round(ell.ee, 6)}"
+
+@app.callback(
+    Output("tabs-GHA-out", "children"),
+    Input("tabs-GHA", "value"),
+)
+def render_content(tab):
+    show1 = {"display": "block"} if tab == "tab-GHA1" else {"display": "none"}
+    show2 = {"display": "block"} if tab == "tab-GHA2" else {"display": "none"}
+
+    pane_gha1 = html.Div(
+        [
+            dcc.Input(id="input-GHA1-beta1", type="number", min=-90, max=90, placeholder="β1...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA1-lamb1", type="number", min=-180, max=180, placeholder="λ1...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA1-s", type="number", min=0, placeholder="s...[m]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA1-a", type="number", min=0, max=360, placeholder="α...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+
+            dcc.Checklist(
+                id="method-checklist-1",
+                options=[
+                    {"label": "Analytisch", "value": "analytisch"},
+                    {"label": "Numerisch", "value": "numerisch"},
+                    {"label": "Stochastisch (ES)", "value": "stochastisch"},
+                ],
+                value=[],
+                style={"marginBottom": "20px"},
+            ),
+            html.Div(
+                [
+                    html.Button(
+                        "Berechnen",
+                        id="button-calc-gha1",
+                        n_clicks=0,
+                        style={"marginRight": "10px"},
+                    ),
+                ],
+                style={"marginBottom": "20px"},
+            ),
+        ],
+        id="pane-gha1",
+        style=show1,
+    )
+
+    pane_gha2 = html.Div(
+        [
+            dcc.Input(id="input-GHA2-beta1", type="number", min=-90, max=90, placeholder="β1...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA2-lamb1", type="number", min=-180, max=180, placeholder="λ1...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA2-beta2", type="number", min=-90, max=90, placeholder="β2...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+            dcc.Input(id="input-GHA2-lamb2", type="number", min=-180, max=180, placeholder="λ2...[°]",
+                      style={"marginBottom": "20px", "display": "block", "width": "300px"}),
+
+            dcc.Checklist(
+                id="method-checklist-2",
+                options=[
+                    {"label": "Numerisch", "value": "numerisch"},
+                    {"label": "Stochastisch (ES)", "value": "stochastisch"},
+                ],
+                value=[],
+                style={"marginBottom": "20px"},
+            ),
+            html.Div(
+                [
+                    html.Button(
+                        "Berechnen",
+                        id="button-calc-gha2",
+                        n_clicks=0,
+                        style={"marginRight": "10px"},
+                    ),
+                ],
+                style={"marginBottom": "20px"},
+            ),
+        ],
+        id="pane-gha2",
+        style=show2,
+    )
+
+    return html.Div([pane_gha1, pane_gha2])
+
+
+@app.callback(
+    Output("output-gha1", "children"),
+    Output("output-gha2", "children"),
+    Output("ellipsoid-plot", "figure"),
+    Input("button-calc-gha1", "n_clicks"),
+    Input("button-calc-gha2", "n_clicks"),
+    State("input-GHA1-beta1", "value"),
+    State("input-GHA1-lamb1", "value"),
+    State("input-GHA1-s", "value"),
+    State("input-GHA1-a", "value"),
+    State("input-GHA2-beta1", "value"),
+    State("input-GHA2-lamb1", "value"),
+    State("input-GHA2-beta2", "value"),
+    State("input-GHA2-lamb2", "value"),
+    State("input-ax", "value"),
+    State("input-ay", "value"),
+    State("input-b", "value"),
+    State("method-checklist-1", "value"),
+    State("method-checklist-2", "value"),
+
+    prevent_initial_call=True,
+)
+def calc_and_plot(n1, n2,
+                  beta11, lamb11, s, a_deg,
+                  beta21, lamb21, beta22, lamb22,
+                  ax, ay, b, method1, method2):
+
+    if not (n1 or n2):
+        return no_update, no_update, no_update
+
+    if not ax or not ay or not b:
+        return html.Span("Bitte Ellipsoid auswählen!", style={"color": "red"}), "", go.Figure()
+
+    ell = EllipsoidTriaxial(ax, ay, b)
+
+    if dash.ctx.triggered_id == "button-calc-gha1":
+        if None in (beta11, lamb11, s, a_deg):
+
+            return html.Span("Bitte β₁, λ₁, s und α eingeben.", style={"color": "red"}), "",  go.Figure()
+
+        beta_rad = wu.deg2rad(float(beta11))
+        lamb_rad = wu.deg2rad(float(lamb11))
+        alpha_rad = wu.deg2rad(float(a_deg))
+        s_val = float(s)
+
+        p1 = tuple(map(float, ell.ell2cart(beta_rad, lamb_rad)))
+        out1 = []
+
+        if "analytisch" in method1:
+            # ana
+            x2, y2, z2 = gha1_ana(ell, p1, alpha_rad, s_val, 70)
+            p2_ana = (float(x2), float(y2), float(z2))
+            beta2, lamb2 = ell.cart2ell([x2, y2, z2])
+
+            #out1 += f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)},"
+            out1.append(
+                html.Div([
+                    html.Strong("Analytisch: "),
+                    html.Br(),
+                    html.Span(f"kartesisch: x₂={x2:.4f} m, y₂={y2:.4f} m, z₂={z2:.4f} m"),
+                    html.Br(),
+                    html.Span(f"ellipsoidisch: {aus.gms('β₂', beta2, 4)}, {aus.gms('λ₂', lamb2, 4)}")
+                ])
+            )
+
+        if "numerisch" in method1:
+            # num
+            p2_num = gha1_num(ell, p1, alpha_rad, s_val, 10000)
+            beta2_num, lamb2_num = ell.cart2ell(p2_num)
+
+            out1.append(
+                html.Div([
+                    html.Strong("Numerisch: "),
+                    html.Br(),
+                    html.Span(f"kartesisch: x₂={p2_num[0]:.4f} m, y₂={p2_num[1]:.4f} m, z₂={p2_num[2]:.4f} m"),
+                    html.Br(),
+                    html.Span(f"ellipsoidisch: {aus.gms('β₂', beta2_num, 4)}, {aus.gms('λ₂', lamb2_num, 4)}")
+                ])
+            )
+
+        if "stochastisch" in method1:
+            # stoch
+            p2_stoch = "noch nicht implementiert.."
+
+            out1.append(
+                html.Div([
+                    html.Strong("Stochastisch (ES): "),
+                    html.Span(f"{p2_stoch}")
+                ])
+            )
+
+        if not method1:
+            return html.Span("Bitte Berechnungsverfahren auswählen!", style={"color": "red"}), "", go.Figure()
+
+        fig = ellipsoid_figure(ell, title="Erste Hauptaufgabe - analystisch")
+        #fig = figure_constant_lines(fig, ell, "geod")
+        fig = figure_constant_lines(fig, ell, "ell")
+        #fig = figure_constant_lines(fig, ell, "para")
+        fig = figure_points(fig, [("P1", p1, "black"), ("P2", p2_ana, "red")])
+
+        #out1 = f"kartesisch: x₂={p2[0]:.5f} m, y₂={p2[1]:.5f} m, z₂={p2[2]:.5f} m; ellipsoidisch: {aus.gms("β₂", beta2, 5)}, {aus.gms("λ₂", lamb2, 5)}, {p2_num}"
+        return out1, "", fig
+
+    if dash.ctx.triggered_id == "button-calc-gha2":
+        if None in (beta21, lamb21, beta22, lamb22):
+            return html.Span("Bitte β₁, λ₁, β₂, λ₂ eingeben.", style={"color": "red"}), "", go.Figure()
+
+        p1 = tuple(ell.ell2cart(np.deg2rad(float(beta21)), np.deg2rad(float(lamb21))))
+        p2 = tuple(ell.ell2cart(np.deg2rad(float(beta22)), np.deg2rad(float(lamb22))))
+
+        out2 = []
+
+        if "numerisch" in method2:
+            alpha_1, alpha_2, s12 = gha2_num(
+                ell,
+                np.deg2rad(float(beta21)), np.deg2rad(float(lamb21)),
+                np.deg2rad(float(beta22)), np.deg2rad(float(lamb22))
+            )
+
+            out2.append(
+                html.Div([
+                    html.Strong("Numerisch: "),
+                    html.Span(f"{aus.gms('α₁₂', alpha_1, 4)}, {aus.gms('α₂₁', alpha_2, 4)}, s = {s12:.4f} m"),
+                ])
+            )
+
+        if "stochastisch" in method2:
+            # stoch
+            a_stoch = "noch nicht implementiert.."
+
+            out2.append(
+                html.Div([
+                    html.Strong("Stochastisch (ES): "),
+                    html.Span(f"{a_stoch}")
+                ])
+            )
+
+        if not method2:
+            return html.Span("Bitte Berechnungsverfahren auswählen!", style={"color": "red"}), "", go.Figure()
+
+        fig = ellipsoid_figure(ell, title="Zweite Hauptaufgabe")
+        fig = figure_constant_lines(fig, ell, "ell")
+        fig = figure_points(fig, [("P1", p1, "black"), ("P2", p2, "red")])
+        
+        return "", out2, fig
+
+    return no_update, no_update, no_update
+
+
+if __name__ == "__main__":
+    app.run(debug=False)

ellipsoide.py

@@ -1,9 +1,11 @@
 import numpy as np
+from numpy import sin, cos, arctan, arctan2, sqrt
 import winkelumrechnungen as wu
 import ausgaben as aus
+import jacobian_Ligas
 
 
-class Ellipsoid:
+class EllipsoidBiaxial:
     def __init__(self, a: float, b: float):
         self.a = a
         self.b = b
@@ -53,9 +55,9 @@ class Ellipsoid:
 
     phi2p = lambda self, phi: self.N(phi) * np.cos(phi)
 
-    def ellipsoidische_Koords (self, Eh, Ephi, x, y, z):
+    def cart2ell(self, Eh, Ephi, x, y, z):
         p = np.sqrt(x**2+y**2)
-        print(f"p = {round(p, 5)} m")
+        # print(f"p = {round(p, 5)} m")
 
         lamb = np.arctan(y/x)
 
@@ -79,5 +81,497 @@ class Ellipsoid:
                 if dphi < Ephi:
                     break
         for i in range(len(phii)):
-            print(f"P3[{i}]: {aus.gms('phi', phii[i], 5)}\th = {round(hi[i], 5)} m")
+            # print(f"P3[{i}]: {aus.gms('phi', phii[i], 5)}\th = {round(hi[i], 5)} m")
+            pass
         return phi, lamb, h
+
+    def ell2cart(self, phi, lamb, h):
+        W = np.sqrt(1 - self.e**2 * np.sin(phi)**2)
+        N = self.a / W
+        x = (N+h) * np.cos(phi) * np.cos(lamb)
+        y = (N+h) * np.cos(phi) * np.sin(lamb)
+        z = (N * (1-self.e**2) + h) * np.sin(lamb)
+        return x, y, z
+
+class EllipsoidTriaxial:
+    def __init__(self, ax: float, ay: float, b: float):
+        self.ax = ax
+        self.ay = ay
+        self.b = b
+        self.ex = np.sqrt((self.ax**2 - self.b**2) / self.ax**2)
+        self.ey = np.sqrt((self.ay**2 - self.b**2) / self.ay**2)
+        self.ee = np.sqrt((self.ax**2 - self.ay**2) / self.ax**2)
+        self.ex_ = np.sqrt((self.ax**2 - self.b**2) / self.b**2)
+        self.ey_ = np.sqrt((self.ay**2 - self.b**2) / self.b**2)
+        self.ee_ = np.sqrt((self.ax**2 - self.ay**2) / self.ay**2)
+        self.Ex = np.sqrt(self.ax**2 - self.b**2)
+        self.Ey = np.sqrt(self.ay**2 - self.b**2)
+        self.Ee = np.sqrt(self.ax**2 - self.ay**2)
+
+    @classmethod
+    def init_name(cls, name: str):
+        """
+        Mögliche Ellipsoide: BursaFialova1993, BursaSima1980, Eitschberger1978, Bursa1972, Bursa1970, BesselBiaxial
+        :param name: Name des dreiachsigen Ellipsoids
+        """
+        if name == "BursaFialova1993":
+            ax = 6378171.36
+            ay = 6378101.61
+            b = 6356751.84
+            return cls(ax, ay, b)
+        elif name == "BursaSima1980":
+            ax = 6378172
+            ay = 6378102.7
+            b = 6356752.6
+            return cls(ax, ay, b)
+        elif name == "BursaSima1980round":
+            # Panou 2013
+            ax = 6378172
+            ay = 6378103
+            b = 6356753
+            return cls(ax, ay, b)
+        elif name == "Eitschberger1978":
+            ax = 6378173.435
+            ay = 6378103.9
+            b = 6356754.4
+            return cls(ax, ay, b)
+        elif name == "Bursa1972":
+            ax = 6378173
+            ay = 6378104
+            b = 6356754
+            return cls(ax, ay, b)
+        elif name == "Bursa1970":
+            ax = 6378173
+            ay = 6378105
+            b = 6356754
+            return cls(ax, ay, b)
+        elif name == "BesselBiaxial":
+            ax = 6377397.15509
+            ay = 6377397.15508
+            b = 6356078.96290
+            return cls(ax, ay, b)
+        elif name == "Fiction":
+            ax = 6000000
+            ay = 5000000
+            b = 4000000
+            return cls(ax, ay, b)
+        elif name == "KarneyTest2024":
+            ax = np.sqrt(2)
+            ay = 1
+            b = 1 / np.sqrt(2)
+            return cls(ax, ay, b)
+
+    def point_on(self, point: np.ndarray) -> bool:
+        """
+        Test, ob ein Punkt auf dem Ellipsoid liegt.
+        :param point: kartesische 3D-Koordinaten
+        :return: Punkt auf dem Ellispoid?
+        """
+        value = point[0]**2/self.ax**2 + point[1]**2/self.ay**2 + point[2]**2/self.b**2
+        if abs(1-value) < 0.000001:
+            return True
+        else:
+            return False
+
+    def ellu2cart(self, beta: float, lamb: float, u: float) -> np.ndarray:
+        """
+        Panou 2014 12ff.
+        Elliptische Breite+Länge sind nicht gleich der geodätischen
+        Verhältnisse des Ellipsoids bekannt, Größe verändern bis Punkt erreicht,
+           dann ist u die Größe entlang der z-Achse
+        :param beta: ellipsoidische Breite [rad]
+        :param lamb: ellipsoidische Länge [rad]
+        :param u: Größe entlang der z-Achse
+        :return: Punkt in kartesischen Koordinaten
+        """
+        # s1 = u**2 - self.b**2
+        # s2 = -self.ay**2 * np.sin(beta)**2 - self.b**2 * np.cos(beta)**2
+        # s3 = -self.ax**2 * np.sin(lamb)**2 - self.ay**2 * np.cos(lamb)**2
+        # print(s1, s2, s3)
+
+        # xe = np.sqrt(((self.ax**2+s1) * (self.ax**2+s2) * (self.ax**2+s3)) /
+        #              ((self.ax**2-self.ay**2) * (self.ax**2-self.b**2)))
+        # ye = np.sqrt(((self.ay**2+s1) * (self.ay**2+s2) * (self.ay**2+s3)) /
+        #              ((self.ay**2-self.ax**2) * (self.ay**2-self.b**2)))
+        # ze = np.sqrt(((self.b**2+s1) * (self.b**2+s2) * (self.b**2+s3)) /
+        #              ((self.b**2-self.ax**2) * (self.b**2-self.ay**2)))
+
+        x = np.sqrt(u**2 + self.Ex**2) * np.sqrt(np.cos(beta)**2 + self.Ee**2/self.Ex**2 * np.sin(beta)**2) * np.cos(lamb)
+        y = np.sqrt(u**2 + self.Ey**2) * np.cos(beta) * np.sin(lamb)
+        z = u * np.sin(beta) * np.sqrt(1 - self.Ee**2/self.Ex**2 * np.cos(lamb)**2)
+
+        return np.array([x, y, z])
+
+    def ell2cart(self, beta: float | np.ndarray, lamb: float | np.ndarray) -> np.ndarray:
+        """
+        Panou, Korakitis 2019 2
+        :param beta: elliptische Breite [rad]
+        :param lamb: elliptische Länge [rad]
+        :return: Punkt in kartesischen Koordinaten
+        """
+        beta = np.asarray(beta, dtype=float)
+        lamb = np.asarray(lamb, dtype=float)
+
+        beta, lamb = np.broadcast_arrays(beta, lamb)
+
+        B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2
+        L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2
+
+        x = self.ax / self.Ex * np.sqrt(B) * np.cos(lamb)
+        y = self.ay * np.cos(beta) * np.sin(lamb)
+        z = self.b / self.Ex * np.sin(beta) * np.sqrt(L)
+
+        xyz = np.stack((x, y, z), axis=-1)
+
+        # Pole
+        mask_south = beta == -np.pi / 2
+        mask_north = beta == np.pi / 2
+        xyz[mask_south] = np.array([0, 0, -self.b])
+        xyz[mask_north] = np.array([0, 0, self.b])
+
+        # Äquator
+        mask_eq = beta == 0
+        xyz[mask_eq & (lamb == -np.pi / 2)] = np.array([0, -self.ay, 0])
+        xyz[mask_eq & (lamb == np.pi / 2)] = np.array([0, self.ay, 0])
+        xyz[mask_eq & (lamb == 0)] = np.array([self.ax, 0, 0])
+        xyz[mask_eq & (lamb == np.pi)] = np.array([-self.ax, 0, 0])
+
+        return xyz
+
+    def cart2ellu(self, point: np.ndarray) -> tuple[float, float, float]:
+        """
+        Panou 2014 15ff.
+        :param point: Punkt in kartesischen Koordinaten
+        :return: elliptische Breite, elliptische Länge, Größe entlang der z-Achse
+        """
+        x, y, z = point
+        c2 = self.ax**2 + self.ay**2 + self.b**2 - x**2 - y**2 - z**2
+        c1 = (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)
+        c0 = (self.ax**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)
+
+        p = (c2**2 - 3*c1) / 9
+        q = (9*c1*c2 - 27*c0 - 2*c2**3) / 54
+        omega = np.arccos(q / np.sqrt(p**3))
+
+        s1 = 2 * np.sqrt(p) * np.cos(omega/3) - c2/3
+        s2 = 2 * np.sqrt(p) * np.cos(omega/3 - 2*np.pi/3) - c2/3
+        s3 = 2 * np.sqrt(p) * np.cos(omega/3 - 4*np.pi/3) - c2/3
+        # print(s1, s2, s3)
+
+        beta = np.arctan(np.sqrt((-self.b**2 - s2) / (self.ay**2 + s2)))
+        if abs((-self.ay**2 - s3) / (self.ax**2 + s3)) > 1e-7:
+            lamb = np.arctan(np.sqrt((-self.ay**2 - s3) / (self.ax**2 + s3)))
+        else:
+            lamb = 0
+        u = np.sqrt(self.b**2 + s1)
+
+        return beta, lamb, u
+
+    def cart2ell(self, point: np.ndarray) -> tuple[float, float]:
+        """
+        Panou, Korakitis 2019 2f.
+        :param point: Punkt in kartesischen Koordinaten
+        :return: elliptische Breite, elliptische Länge
+        """
+        x, y, z = point
+
+        eps = 1e-9
+
+        if abs(x) < eps and abs(y) < eps:  # Punkt in der z-Achse
+            beta = np.pi / 2 if z > 0 else -np.pi / 2
+            lamb = 0.0
+            return beta, lamb
+
+        elif abs(x) < eps and abs(z) < eps:  # Punkt in der y-Achse
+            beta = 0.0
+            lamb = np.pi / 2 if y > 0 else -np.pi / 2
+            return beta, lamb
+
+        elif abs(y) < eps and abs(z) < eps:  # Punkt in der x-Achse
+            beta = 0.0
+            lamb = 0.0 if x > 0 else np.pi
+            return beta, lamb
+
+        # ---- 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 + np.sqrt(c1 ** 2 - 4 * c0)) / 2
+        t1 = c0 / t2
+
+        num_beta = max(t1 - self.b ** 2, 0)
+        den_beta = max(self.ay ** 2 - t1, 0)
+        num_lamb = max(t2 - self.ay ** 2, 0)
+        den_lamb = max(self.ax ** 2 - t2, 0)
+        beta = np.arctan(np.sqrt(num_beta / den_beta))
+        lamb = np.arctan(np.sqrt(num_lamb / den_lamb))
+
+        if z < 0:
+            beta = -beta
+
+        if y < 0:
+            lamb = -lamb
+
+        if x < 0:
+            lamb = np.pi - lamb
+
+        if abs(x) < eps:
+            lamb = -np.pi/2 if y < 0 else np.pi/2
+
+        elif abs(y) < eps:
+            lamb = 0 if x > 0 else np.pi
+
+        elif abs(z) < eps:
+            beta = 0
+
+        return beta, lamb
+
+    def cart2geod(self, point: np.ndarray, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
+        """
+        Ligas 2012
+        :param mode: ligas1, ligas2, oder ligas3
+        :param point: Punkt in kartesischen Koordinaten
+        :param maxIter: maximale Anzahl Iterationen
+        :param maxLoa: Level of Accuracy, das erreicht werden soll
+        :return: phi, lambda, h
+        """
+        xG, yG, zG = point
+
+        eps = 1e-9
+
+        if abs(xG) < eps and abs(yG) < eps:  # Punkt in der z-Achse
+            phi = np.pi / 2 if zG > 0 else -np.pi / 2
+            lamb = 0.0
+            h = abs(zG) - self.b
+            return phi, lamb, h
+
+        elif abs(xG) < eps and abs(zG) < eps:  # Punkt in der y-Achse
+            phi = 0.0
+            lamb = np.pi / 2 if yG > 0 else -np.pi / 2
+            h = abs(yG) - self.ay
+            return phi, lamb, h
+
+        elif abs(yG) < eps and abs(zG) < eps:  # Punkt in der x-Achse
+            phi = 0.0
+            lamb = 0.0 if xG > 0 else np.pi
+            h = abs(xG) - self.ax
+            return phi, lamb, h
+
+        rG = np.sqrt(xG ** 2 + yG ** 2 + zG ** 2)
+        pE = np.array([self.ax * xG / rG, self.ax * yG / rG, self.ax * zG / rG], dtype=np.float64)
+
+        E = 1 / self.ax**2
+        F = 1 / self.ay**2
+        G = 1 / self.b**2
+
+        i = 0
+        loa = np.inf
+
+        while i < maxIter and loa > maxLoa:
+            if mode == "ligas1":
+                invJ, fxE = jacobian_Ligas.case1(E, F, G, np.array([xG, yG, zG]), pE)
+            elif mode == "ligas2":
+                invJ, fxE = jacobian_Ligas.case2(E, F, G, np.array([xG, yG, zG]), pE)
+            elif mode == "ligas3":
+                invJ, fxE = jacobian_Ligas.case3(E, F, G, np.array([xG, yG, zG]), pE)
+            pEi = pE.reshape(-1, 1) - invJ @ fxE.reshape(-1, 1)
+            pEi = pEi.reshape(1, -1).flatten()
+            loa = np.sqrt((pEi[0]-pE[0])**2 + (pEi[1]-pE[1])**2 + (pEi[2]-pE[2])**2)
+            pE = pEi
+            i += 1
+
+        phi = np.arctan((1-self.ee**2) / (1-self.ex**2) * pE[2] / np.sqrt((1-self.ee**2)**2 * pE[0]**2 + pE[1]**2))
+        lamb = np.arctan(1/(1-self.ee**2) * pE[1]/pE[0])
+        h = np.sign(zG - pE[2]) * np.sign(pE[2]) * np.sqrt((pE[0] - xG) ** 2 + (pE[1] - yG) ** 2 + (pE[2] - zG) ** 2)
+
+        if xG < 0 and yG < 0:
+            lamb = -np.pi + lamb
+
+        elif xG < 0:
+            lamb = np.pi + lamb
+
+        if abs(zG) < eps:
+            phi = 0
+
+        return phi, lamb, h
+
+    def geod2cart(self, phi: float | np.ndarray, lamb: float | np.ndarray, h: float) -> np.ndarray:
+        """
+        Ligas 2012, 250
+        :param phi: geodätische Breite [rad]
+        :param lamb: geodätische Länge [rad]
+        :param h: Höhe über dem Ellipsoid
+        :return: kartesische Koordinaten
+        """
+        v = self.ax / np.sqrt(1 - self.ex**2*np.sin(phi)**2-self.ee**2*np.cos(phi)**2*np.sin(lamb)**2)
+        xG = (v + h) * np.cos(phi) * np.cos(lamb)
+        yG = (v * (1-self.ee**2) + h) * np.cos(phi) * np.sin(lamb)
+        zG = (v * (1-self.ex**2) + h) * np.sin(phi)
+        return np.array([xG, yG, zG])
+
+    def cartonell(self, point: np.ndarray) -> tuple[np.ndarray, float, float, float]:
+        """
+        Berechnung des Lotpunktes auf einem Ellipsoiden
+        :param point: Punkt in kartesischen Koordinaten, der gelotet werden soll
+        :return: Lotpunkt in kartesischen Koordinaten, geodätische Koordinaten des Punktes
+        """
+        phi, lamb, h = self.cart2geod(point, "ligas3")
+        x, y, z = self. geod2cart(phi, lamb, 0)
+        return np.array([x, y, z]), phi, lamb, h
+
+    def cartellh(self, point: np.ndarray, h: float) -> np.ndarray:
+        """
+        Punkt auf Ellipsoid hoch loten
+        :param point: Punkt auf dem Ellipsoid
+        :param h: Höhe über dem Ellipsoid
+        :return: hochgeloteter Punkt
+        """
+        phi, lamb, _ = self.cart2geod(point, "ligas3")
+        pointH = self. geod2cart(phi, lamb, h)
+        return pointH
+
+    def para2cart(self, u: float | np.ndarray, v: float | np.ndarray) -> np.ndarray:
+        """
+        Panou, Korakitits 2020, 4
+        :param u: Parameter u
+        :param v: Parameter v
+        :return: Punkt in kartesischen Koordinaten
+        """
+        x = self.ax * np.cos(u) * np.cos(v)
+        y = self.ay * np.cos(u) * np.sin(v)
+        z = self.b * np.sin(u)
+        z = np.broadcast_to(z, np.shape(x))
+        return np.array([x, y, z])
+
+    def cart2para(self, point: np.ndarray) -> tuple[float, float]:
+        """
+        Panou, Korakitits 2020, 4
+        :param point: Punkt in kartesischen Koordinaten
+        :return: parametrische Koordinaten
+        """
+        x, y, z = point
+
+        u_check1 = z*np.sqrt(1 - self.ee**2)
+        u_check2 = np.sqrt(x**2 * (1-self.ee**2) + y**2) * np.sqrt(1-self.ex**2)
+        if u_check1 <= u_check2:
+            u = np.arctan2(u_check1, u_check2)
+        else:
+            u = np.pi/2 - np.arctan2(u_check2, u_check1)
+
+        v_check1 = y
+        v_check2 = x*np.sqrt(1-self.ee**2)
+        v_factor = np.sqrt(x**2*(1-self.ee**2)+y**2)
+        if v_check1 <= v_check2:
+            v = 2 * np.arctan2(v_check1, v_check2 + v_factor)
+        else:
+            v = np.pi/2 - 2 * np.arctan2(v_check2, v_check1 + v_factor)
+
+        return u, v
+
+    def ell2para(self, beta, lamb) -> tuple[float, float]:
+        cart = self.ell2cart(beta, lamb)
+        return self.cart2para(cart)
+
+    def para2ell(self, u, v) -> tuple[float, float]:
+        cart = self.para2cart(u, v)
+        return self.cart2ell(cart)
+
+    def para2geod(self, u: float, v: float, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
+        cart = self.para2cart(u, v)
+        return self.cart2geod(cart, mode, maxIter, maxLoa)
+
+    def geod2para(self, phi, lamb, h) -> tuple[float, float]:
+        cart = self.geod2cart(phi, lamb, h)
+        return self.cart2para(cart)
+
+    def ell2geod(self, beta, lamb, mode: str = "ligas3", maxIter: int = 30, maxLoa: float = 0.005) -> tuple[float, float, float]:
+        cart = self.ell2cart(beta, lamb)
+        return self.cart2geod(cart, mode, maxIter, maxLoa)
+
+    def func_H(self, x, y, z):
+        return x ** 2 + y ** 2 / (1 - self.ee ** 2) ** 2 + z ** 2 / (1 - self.ex ** 2) ** 2
+
+    def func_n(self, x, y, z, H=None):
+        if H is None:
+            H = self.func_H(x, y, z)
+        return np.array([x / sqrt(H),
+                         y / ((1 - self.ee ** 2) * sqrt(H)),
+                         z / ((1 - self.ex ** 2) * sqrt(H))])
+
+    def p_q(self, x, y, z) -> tuple[np.ndarray, np.ndarray]:
+        """
+        Berechnung sämtlicher Größen
+        :param x: x
+        :param y: y
+        :param z: z
+        :return: Dictionary sämtlicher Größen
+        """
+        n = self.func_n(x, y, z)
+
+        beta, lamb = self.cart2ell(np.array([x, y, z]))
+        B = self.Ex ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(beta) ** 2
+        L = self.Ex ** 2 - self.Ee ** 2 * np.cos(lamb) ** 2
+
+        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 + np.sqrt(c1 ** 2 - 4 * c0)) / 2
+        t1 = c0 / t2
+        t2e = self.ax ** 2 * np.sin(lamb) ** 2 + self.ay ** 2 * np.cos(lamb) ** 2
+        t1e = self.ay ** 2 * np.sin(beta) ** 2 + self.b ** 2 * np.cos(beta) ** 2
+
+        F = self.Ey ** 2 * np.cos(beta) ** 2 + self.Ee ** 2 * np.sin(lamb) ** 2
+        p1 = -np.sqrt(L / (F * t2)) * self.ax / self.Ex * np.sqrt(B) * np.sin(lamb)
+        p2 = np.sqrt(L / (F * t2)) * self.ay * np.cos(beta) * np.cos(lamb)
+        p3 = 1 / np.sqrt(F * t2) * (self.b * self.Ee ** 2) / (2 * self.Ex) * np.sin(beta) * np.sin(2 * lamb)
+        # p1 = -np.sign(y) * np.sqrt(L / (F * t2)) * self.ax / (self.Ex * self.Ee) * np.sqrt(B) * np.sqrt(t2 - self.ay ** 2)
+        # p2 = np.sign(x) * np.sqrt(L / (F * t2)) * self.ay / (self.Ey * self.Ee) * np.sqrt((ell.ay ** 2 - t1) * (self.ax ** 2 - t2))
+        # p3 = np.sign(x) * np.sign(y) * np.sign(z) * 1 / np.sqrt(F * t2) * self.b / (self.Ex * self.Ey) * np.sqrt(
+        #     (t1 - self.b ** 2) * (t2 - self.ay ** 2) * (self.ax ** 2 - t2))
+        p = np.array([p1, p2, p3])
+        q = np.array([n[1] * p[2] - n[2] * p[1],
+                      n[2] * p[0] - n[0] * p[2],
+                      n[0] * p[1] - n[1] * p[0]])
+
+        return p, q
+
+
+if __name__ == "__main__":
+    ell = EllipsoidTriaxial.init_name("Eitschberger1978")
+    diff_list = []
+    for beta_deg in range(-150, 210, 30):
+        for lamb_deg in range(-150, 210, 30):
+            beta = wu.deg2rad(beta_deg)
+            lamb = wu.deg2rad(lamb_deg)
+            point = ell.ell2cart(beta, lamb)
+
+            elli = ell.cart2ell(point)
+            cart_elli = ell.ell2cart(elli[0], elli[1])
+            diff_ell = np.sum(np.abs(point-cart_elli))
+
+            para = ell.cart2para(point)
+            cart_para = ell.para2cart(para[0], para[1])
+            diff_para = np.sum(np.abs(point-cart_para))
+
+            # geod = ell.cart2geod(point, "ligas1")
+            # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
+            # diff_geod1 = np.sum(np.abs(point-cart_geod))
+            #
+            # geod = ell.cart2geod(point, "ligas2")
+            # cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
+            # diff_geod2 = np.sum(np.abs(point-cart_geod))
+
+            geod = ell.cart2geod(point, "ligas3")
+            cart_geod = ell.geod2cart(geod[0], geod[1], geod[2])
+            diff_geod3 = np.sum(np.abs(point-cart_geod))
+
+            diff_list.append([beta_deg, lamb_deg, diff_ell, diff_para, diff_geod3])
+            diff_list.append([diff_ell])
+
+    diff_list = np.array(diff_list)
+    pass

hausarbeit.py

@@ -1,78 +0,0 @@
-from numpy import cos, sin, tan
-import winkelumrechnungen as wu
-import s_ellipse as s_ell
-import ausgaben as aus
-from ellipsoide import Ellipsoid
-import Numerische_Integration.num_int_runge_kutta as rk
-import GHA.gauss as gauss
-import GHA.bessel as bessel
-
-matrikelnummer = "6044051"
-print(f"Matrikelnummer: {matrikelnummer}")
-
-m4 = matrikelnummer[-4]
-m3 = matrikelnummer[-3]
-m2 = matrikelnummer[-2]
-m1 = matrikelnummer[-1]
-print(f"m1={m1}\tm2={m2}\tm3={m3}\tm4={m4}")
-
-re = Ellipsoid.init_name("Bessel")
-nks = 3
-
-print(f"\na = {re.a} m\nb = {re.b} m\nc = {re.c} m\ne' = {re.e_}")
-
-print("\n\nAufgabe 1 (via Reihenentwicklung)")
-phi1 = re.beta2phi(wu.gms2rad([int("1"+m1), int("1"+m2), int("1"+m3)]))
-print(f"{aus.gms('phi_P1', phi1, 5)}")
-s = s_ell.reihenentwicklung(re.e_, re.c, phi1)
-print(f'\ns_P1 = {round(s,nks)} m')
-# s_poly = s_ell.polyapp_tscheby_bessel(phi1)
-# print(f'Meridianbogenlänge: {round(s_poly,nks)} m (Polynomapproximation)')
-
-print("\n\nAufgabe 2 (via Gauß´schen Mittelbreitenformeln)")
-
-lambda1 = wu.gms2rad([int("1"+m3), int("1"+m2), int("1"+m4)])
-A12 = wu.gms2rad([90, 0, 0])
-s12 = float("1"+m4+m3+m2+"."+m1+"0")
-
-# print("\nvia Runge-Kutta-Verfahren")
-# re = Ellipsoid.init_name("Bessel")
-# f_phi = lambda s, phi, lam, A: cos(A) * re.V(phi) ** 3 / re.c
-# f_lam = lambda s, phi, lam, A: sin(A) * re.V(phi) / (cos(phi) * re.c)
-# f_A = lambda s, phi, lam, A: tan(phi) * sin(A) * re.V(phi) / re.c
-#
-# funktionswerte = rk.verfahren([f_phi, f_lam, f_A],
-#                               [0, phi1, lambda1, A12],
-#                               s12, 1)
-#
-# for s, phi, lam, A in funktionswerte:
-#     print(f"{s} m, {aus.gms('phi', phi, nks)}, {aus.gms('lambda', lam, nks)}, {aus.gms('A', A, nks)}")
-
-# print("via Gauß´schen Mittelbreitenformeln")
-phi_p2, lambda_p2, A_p2 = gauss.gha1(Ellipsoid.init_name("Bessel"),
-                                     phi_p1=phi1,
-                                     lambda_p1=lambda1,
-                                     A_p1=A12,
-                                     s=s12, eps=wu.gms2rad([1*10**-8, 0, 00]))
-
-print(f"\nP2: {aus.gms('phi', phi_p2[-1], nks)}\t\t{aus.gms('lambda', lambda_p2[-1], nks)}\t{aus.gms('A', A_p2[-1], nks)}")
-
-# print("\nvia Verfahren nach Bessel")
-# phi_p2, lambda_p2, A_p2 = bessel.gha1(Ellipsoid.init_name("Bessel"),
-#                                       phi_p1=phi1,
-#                                       lambda_p1=lambda1,
-#                                       A_p1=A12,
-#                                       s=s12)
-# print(f"P2: {aus.gms('phi', phi_p2, nks)}, {aus.gms('lambda', lambda_p2, nks)}, {aus.gms('A', A_p2, nks)}")
-
-print("\n\nAufgabe 3")
-p = re.phi2p(phi_p2[-1])
-print(f"p = {round(p,2)} m")
-
-
-print("\n\nAufgabe 4")
-x = float("4308"+m3+"94.556")
-y = float("1214"+m2+"88.242")
-z = float("4529"+m4+"03.878")
-phi_p3, lambda_p3, h_p3 = re.ellipsoidische_Koords(0.001, wu.gms2rad([0, 0, 0.001]), x, y, z)
-print(f"\nP3: {aus.gms('phi', phi_p3, nks)}, {aus.gms('lambda', lambda_p3, 5)}, h = {round(h_p3,nks)} m")

jacobian_Ligas.py

@@ -0,0 +1,74 @@
+import numpy as np
+
+def case1(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+    j11 = 2 * E * pE[0]
+    j12 = 2 * F * pE[1]
+    j13 = 2 * G * pE[2]
+    j21 = F * pE[1] - (pE[1] - pG[1]) * E
+    j22 = (pE[0] - pG[0]) * F - E * pE[0]
+    j23 = 0
+    j31 = G * pE[2] - (pE[2] - pG[2]) * E
+    j32 = 0
+    j33 = (pE[0] - pG[0]) * G - E * pE[0]
+
+    detJ = j11 * j22 * j33 - j21 * j12 * j33 - j31 * j13 * j22
+
+    invJ = 1/detJ * np.array([[j22*j33, -j12*j33, -j13*j22],
+                              [-j21*j33, j11*j33-j13*j31, j13*j21],
+                              [-j22*j31, j12*j31, j11*j22-j12*j21]])
+
+    fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
+                    (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
+                    (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0]])
+
+    return invJ, fxE
+
+def case2(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+    j11 = 2 * E * pE[0]
+    j12 = 2 * F * pE[1]
+    j13 = 2 * G * pE[2]
+    j21 = F * pE[1] - (pE[1] - pG[1]) * E
+    j22 = (pE[0] - pG[0]) * F - E * pE[0]
+    j23 = 0
+    j31 = 0
+    j32 = G * pE[2] - (pE[2] - pG[2]) * F
+    j33 = (pE[1] - pG[1]) * G - F * pE[1]
+
+    detJ = j11 * j22 * j33 - j21 * j12 * j33 + j21 * j13 * j32
+    if detJ == 0:
+        invJ, fxE = case3(E, F, G, pG, pE)
+    else:
+        invJ = 1/detJ * np.array([[j22*j33, -(j12*j33-j13*j32), -j13*j22],
+                                  [-j21*j33, j11*j33, j13*j21],
+                                  [j21*j32, -j11*j32, j11*j22-j12*j21]])
+
+        fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
+                        (pE[0]-pG[0]) * F*pE[1] - (pE[1]-pG[1]) * E*pE[0],
+                        (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
+
+    return invJ, fxE
+
+def case3(E: float, F: float, G: float, pG: np.ndarray, pE: np.ndarray):
+    j11 = 2 * E * pE[0]
+    j12 = 2 * F * pE[1]
+    j13 = 2 * G * pE[2]
+    j21 = G * pE[2] - (pE[2] - pG[2]) * E
+    j22 = 0
+    j23 = (pE[0] - pG[0]) * G - E * pE[0]
+    j31 = 0
+    j32 = G * pE[2] - (pE[2] - pG[2]) * F
+    j33 = (pE[1] - pG[1]) * G - F * pE[1]
+
+    detJ = -j11 * j23 * j32 - j21 * j12 * j33 + j21 * j13 * j32
+    if detJ == 0:
+        invJ, fxE = case2(E, F, G, pG, pE)
+    else:
+        invJ = 1/detJ * np.array([[-j23*j32, -(j12*j33-j13*j32), j12*j23],
+                                  [-j21*j33, j11*j33, -(j11*j23-j13*j21)],
+                                  [j21*j32, -j11*j32, -j12*j21]])
+
+        fxE = np.array([E*pE[0]**2 + F*pE[1]**2 + G*pE[2]**2 - 1,
+                        (pE[0]-pG[0]) * G*pE[2] - (pE[2]-pG[2]) * E*pE[0],
+                        (pE[1]-pG[1]) * G*pE[2] - (pE[2]-pG[2]) * F*pE[1]])
+
+    return invJ, fxE

kugel.py

@@ -0,0 +1,78 @@
+import numpy as np
+
+
+def cart2sph(x, y, z):
+    r = np.sqrt(x**2 + y**2 + z**2)
+    phi = np.atan2(z, np.sqrt(x**2 + y**2))
+    lamb = np.atan2(y, x)
+
+    return r, np.rad2deg(phi), np.rad2deg(lamb)
+
+
+def sph2cart(r, phi, lamb):
+    phi_rad = np.deg2rad(phi)
+    lamb_rad = np.deg2rad(lamb)
+
+    x = r * np.cos(phi_rad) * np.cos(lamb_rad)
+    y = r * np.cos(phi_rad) * np.sin(lamb_rad)
+    z = r * np.sin(phi_rad)
+
+    return x, y, z
+
+
+def kugel_erste_gha(R, phi_1, lamb_1, s_d, a_12):
+
+    s = s_d / R
+    lamb_1_rad = np.deg2rad(lamb_1)
+    phi_1_rad = np.deg2rad(phi_1)
+    a_12_rad = np.deg2rad(a_12)
+
+    sin_satz = np.sin(s) * np.sin(a_12_rad)
+    sin_kos = np.cos(phi_1_rad) * np.cos(s) - np.sin(phi_1_rad) * np.sin(s) * np.cos(a_12_rad)
+    delta_lam = np.atan2(sin_satz, sin_kos)
+    lamb_2_rad = lamb_1_rad + delta_lam
+
+    phi_2_rad = np.asin(np.sin(phi_1_rad) * np.cos(s) + np.cos(phi_1_rad) * np.sin(s) * np.cos(a_12_rad))
+
+    sin_satz_2 = -np.cos(phi_1_rad) * np.sin(a_12_rad)
+    sin_kos_2 = np.sin(s) * np.sin(phi_1_rad) - np.cos(s) * np.cos(phi_1_rad) * np.cos(a_12_rad)
+    a_21_rad = np.atan2(sin_satz_2, sin_kos_2)
+    if a_21_rad < 0:
+        a_21_rad += 2 * np.pi
+
+    return np.rad2deg(phi_2_rad), np.rad2deg(lamb_2_rad), np.rad2deg(a_21_rad)
+
+
+def kugel_zweite_gha(R, phi_1, lamb_1, phi_2, lamb_2):
+    phi_1_rad = np.deg2rad(phi_1)
+    lamb_1_rad = np.deg2rad(lamb_1)
+    phi_2_rad = np.deg2rad(phi_2)
+    lamb_2_rad = np.deg2rad(lamb_2)
+
+    sin_satz_1 = np.cos(phi_2_rad) * np.sin(lamb_2_rad - lamb_1_rad)
+    sin_kos_1 = np.cos(phi_1_rad) * np.sin(phi_2_rad) - np.sin(phi_1_rad) * np.cos(phi_2_rad) * np.cos(lamb_2_rad - lamb_1_rad)
+    a_12_rad = np.atan2(sin_satz_1, sin_kos_1)
+
+    kos = np.sin(phi_1_rad) * np.sin(phi_2_rad) + np.cos(phi_1_rad) * np.cos(phi_2_rad) * np.cos(lamb_2_rad - lamb_1_rad)
+    s = np.acos(kos)
+    s_d = R * s
+
+    sin_satz_2 = -np.cos(phi_1_rad) * np.sin(lamb_2_rad - lamb_1_rad)
+    sin_kos_2 = np.cos(phi_2_rad) * np.sin(phi_1_rad) - np.sin(phi_2_rad) * np.cos(phi_1_rad) * np.cos(lamb_2_rad - lamb_1_rad)
+    a_21_rad = np.atan2(sin_satz_2, sin_kos_2)
+    if a_21_rad < 0:
+        a_21_rad += 2 * np.pi
+
+    return s_d, np.rad2deg(a_12_rad), np.rad2deg(a_21_rad)
+
+
+if __name__ == "__main__":
+    R = 6378815.904                                             # Bern
+    phi_1 = 10
+    lamb_1 = 40
+    a_12 = 100
+    s = 10000
+    phi_2, lamb_2, _ = kugel_erste_gha(R, phi_1, lamb_1, s, a_12)
+
+    print(kugel_erste_gha(R, phi_1, lamb_1, s, a_12))
+    print(kugel_zweite_gha(R, phi_1, lamb_1, phi_2, lamb_2))

show_constant_lines.py

@@ -0,0 +1,30 @@
+import numpy as np
+import plotly.graph_objects as go
+from ellipsoide import EllipsoidTriaxial
+import winkelumrechnungen as wu
+from dashboard import ellipsoid_figure
+
+u = np.linspace(0, 2*np.pi, 51)
+v = np.linspace(0, np.pi, 51)
+ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+points = []
+lines = []
+for u_i, u_value in enumerate(u):
+    for v_i, v_value in enumerate(v):
+        cart = ell.ell2cart(u_value, v_value)
+        if u_i != 0 and v_i != 0:
+            lines.append((points[-1], cart, "red"))
+        points.append(cart)
+points = []
+for v_i, v_value in enumerate(v):
+    for u_i, u_value in enumerate(u):
+        cart = ell.ell2cart(u_value, v_value)
+        if u_i != 0 and v_i != 0:
+            lines.append((points[-1], cart, "blue"))
+        points.append(cart)
+ax = ell.ax
+ay = ell.ay
+b = ell.b
+
+figu = ellipsoid_figure(ax, ay, b, lines=lines)
+figu.show()

test.py

@@ -0,0 +1,29 @@
+import numpy as np
+from scipy.special import factorial as fact
+from math import comb
+
+J = np.array([
+    [2, 3, 0],
+    [0, 3, 0],
+    [6, 0, 4]
+])
+
+xi = np.array([1, 2, 3])
+xi_col = xi.reshape(-1, 1)
+print(xi_col)
+xi_row = xi_col.reshape(1, -1).flatten()
+print(xi_row)
+
+# Spaltenvektor-Variante
+res_col = xi[:, None] - J @ xi[:, None]
+
+# Zeilenvektor-Variante
+res_row = xi[None, :] - xi[None, :] @ J
+
+print("Spaltenvektor:")
+print(res_col[0,0])
+print("Zeilenvektor:")
+print(res_row)
+t = 5
+l = 2
+print(fact(t+1-l) / (fact(t+1-l) * fact(l-1)), comb(l-1, t+1-l))

test_algorithms.py

@@ -0,0 +1,84 @@
+import GHA_triaxial.numeric_examples_panou as nep
+import ellipsoide
+from GHA_triaxial.panou_2013_2GHA_num import gha2_num
+from GHA_triaxial.panou import gha1_ana, gha1_num
+import numpy as np
+import time
+
+def test():
+    ell = ellipsoide.EllipsoidTriaxial.init_name("BursaSima1980round")
+
+    tables = nep.get_tables()
+
+    diffs_gha1_num = []
+    diffs_gha1_ana = []
+    diffs_gha2_num = []
+    times_gha1_num = []
+    times_gha1_ana = []
+    times_gha2_num = []
+
+    for table in tables:
+        diffs_gha1_num.append([])
+        diffs_gha1_ana.append([])
+        diffs_gha2_num.append([])
+        times_gha1_num.append([])
+        times_gha1_ana.append([])
+        times_gha2_num.append([])
+
+        for example in table:
+            beta0, lamb0, beta1, lamb1, c, alpha0, alpha1, s = example
+            P0 = ell.ell2cart(beta0, lamb0)
+            P1 = ell.ell2cart(beta1, lamb1)
+
+            start = time.perf_counter()
+            try:
+                P1_num = gha1_num(ell, P0, alpha0, s, 10000)
+                end = time.perf_counter()
+                diff_P1_num = np.linalg.norm(P1 - P1_num)
+            except:
+                end = time.perf_counter()
+                diff_P1_num = None
+            time_gha1_num = end - start
+
+            start = time.perf_counter()
+            try:
+                P1_ana = gha1_ana(ell, P0, alpha0, s, 50)
+                end = time.perf_counter()
+                diff_P1_ana = np.linalg.norm(P1 - P1_ana)
+            except:
+                end = time.perf_counter()
+                diff_P1_ana = None
+            time_gha1_ana = end - start
+
+            start = time.perf_counter()
+            try:
+                alpha0_num, alpha1_num, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=1000)
+                end = time.perf_counter()
+                diff_s_num = abs(s - s_num)
+            except:
+                end = time.perf_counter()
+                diff_s_num = None
+                time_gha2_num = None
+            time_gha2_num = end - start
+
+            diffs_gha1_num[-1].append(diff_P1_num)
+            diffs_gha1_ana[-1].append(diff_P1_ana)
+            diffs_gha2_num[-1].append(diff_s_num)
+            times_gha1_num[-1].append(time_gha1_num)
+            times_gha1_ana[-1].append(time_gha1_ana)
+            times_gha2_num[-1].append(time_gha2_num)
+    print(diffs_gha1_num, diffs_gha1_ana, diffs_gha2_num)
+    print(times_gha1_num, times_gha1_ana, times_gha2_num)
+
+def display():
+    diffs = [[{'gha1_num': np.float64(3.410763124264611e-05), 'gha1_ana': np.float64(3.393273802112796e-05), 'gha2_num': np.float64(3.3931806683540344e-05)}, {'gha1_num': np.float64(0.0008736425000530604), 'gha1_ana': np.float64(0.0008736458415010259), 'gha2_num': None}, {'gha1_num': np.float64(0.0007739730058338136), 'gha1_ana': np.float64(0.0007739621469802854), 'gha2_num': np.float64(1.5832483768463135e-07)}, {'gha1_num': np.float64(0.00010554956741100295), 'gha1_ana': np.float64(8.814246009944831), 'gha2_num': np.float64(4.864111542701721e-05)}, {'gha1_num': np.float64(0.0002135908394614854), 'gha1_ana': np.float64(0.0002138610897967267), 'gha2_num': np.float64(5.0179407158866525)}, {'gha1_num': np.float64(0.00032727226891456654), 'gha1_ana': np.float64(0.00032734569198545905), 'gha2_num': np.float64(9.735533967614174e-05)}, {'gha1_num': np.float64(0.0005195973303787956), 'gha1_ana': np.float64(0.0005197766935509641), 'gha2_num': None}], [{'gha1_num': np.float64(1.780250537652368e-05), 'gha1_ana': np.float64(1.996805145339501e-05), 'gha2_num': np.float64(1.8164515495300293e-05)}, {'gha1_num': np.float64(4.8607540473363564e-05), 'gha1_ana': np.float64(2205539.954949392), 'gha2_num': None}, {'gha1_num': np.float64(0.00017376854985685854), 'gha1_ana': np.float64(328124.1513636429), 'gha2_num': np.float64(0.17443156614899635)}, {'gha1_num': np.float64(5.83429352558999e-05), 'gha1_ana': np.float64(0.01891628037258558), 'gha2_num': np.float64(1.4207654744386673)}, {'gha1_num': np.float64(0.0006421087024666934), 'gha1_ana': np.float64(0.0006420400127297228), 'gha2_num': np.float64(0.12751091085374355)}, {'gha1_num': np.float64(0.0004456207867164434), 'gha1_ana': np.float64(0.0004455649707698245), 'gha2_num': np.float64(0.00922046648338437)}, {'gha1_num': np.float64(0.0002340879908275419), 'gha1_ana': np.float64(0.00023422217242111216), 'gha2_num': np.float64(0.001307751052081585)}], [{'gha1_num': np.float64(976.6580096633622), 'gha1_ana': np.float64(976.6580096562798), 'gha2_num': np.float64(6.96033239364624e-05)}, {'gha1_num': np.float64(2825.2936643258527), 'gha1_ana': np.float64(2794.954866417055), 'gha2_num': np.float64(1.3615936040878296e-05)}, {'gha1_num': np.float64(1248.8942058074501), 'gha1_ana': np.float64(538.5550561841195), 'gha2_num': np.float64(3.722589462995529e-05)}, {'gha1_num': np.float64(2201.1793359793814), 'gha1_ana': np.float64(3735.376499414938), 'gha2_num': np.float64(1.4525838196277618e-05)}, {'gha1_num': np.float64(2262.134819997246), 'gha1_ana': np.float64(25549.567793410763), 'gha2_num': np.float64(9.328126907348633e-06)}, {'gha1_num': np.float64(2673.219788119847), 'gha1_ana': np.float64(21760.866677295206), 'gha2_num': np.float64(8.635222911834717e-06)}, {'gha1_num': np.float64(1708.758419275875), 'gha1_ana': np.float64(3792.1128807063437), 'gha2_num': np.float64(2.4085864424705505e-05)}], [{'gha1_num': np.float64(0.7854659044152204), 'gha1_ana': np.float64(0.785466068424286), 'gha2_num': np.float64(0.785466069355607)}, {'gha1_num': np.float64(237.79878717216718), 'gha1_ana': np.float64(1905080.064324282), 'gha2_num': None}, {'gha1_num': np.float64(55204.601699830164), 'gha1_ana': np.float64(55204.60175211949), 'gha2_num': None}, {'gha1_num': np.float64(12766.348063015519), 'gha1_ana': np.float64(12766.376619517901), 'gha2_num': np.float64(12582.786206113175)}, {'gha1_num': np.float64(29703.049988324146), 'gha1_ana': np.float64(29703.056427749252), 'gha2_num': np.float64(28933.668131249025)}, {'gha1_num': np.float64(43912.03007182513), 'gha1_ana': np.float64(43912.03007528712), 'gha2_num': None}, {'gha1_num': np.float64(28522.29828970693), 'gha1_ana': np.float64(28522.29830145182), 'gha2_num': None}, {'gha1_num': np.float64(17769.115549537233), 'gha1_ana': np.float64(17769.115549483362), 'gha2_num': np.float64(17769.121286311187)}]]
+    arr = []
+    for table in diffs:
+        for example in table:
+            arr.append([example['gha1_num'], example['gha1_ana'], example['gha2_num']])
+    arr = np.array(arr)
+    pass
+
+if __name__ == "__main__":
+    # test()
+    display()

winkelumrechnungen.py

@@ -1,4 +1,5 @@
 from numpy import *
+import numpy as np
 
 
 def deg2gms(deg: float) -> list:
@@ -10,13 +11,13 @@ def deg2gms(deg: float) -> list:
     :rtype: list
     """
     gra = deg // 1
-    min = gra % 1
+    minu = gra % 1
     gra = gra // 1
-    min *= 60
+    minu *= 60
-    sek = min % 1
+    sek = minu % 1
-    min = min // 1
+    minu = minu // 1
     sek *= 60
-    return [gra, min, sek]
+    return [gra, minu, sek]
 
 
 def deg2gra(deg: float) -> float:
@@ -30,13 +31,13 @@ def deg2gra(deg: float) -> float:
     return deg * 10/9
 
 
-def deg2rad(deg: float) -> float:
+def deg2rad(deg: float | np.ndarray) -> float | np.ndarray:
     """
     Umrechnung von Grad in Radiant
     :param deg: Winkel in Grad
-    :type deg: float
+    :type deg: float or np.ndarray
     :return: Winkel in Radiant
-    :rtype: float
+    :rtype: float or np.ndarray
     """
     return deg * pi / 180
 
@@ -51,13 +52,13 @@ def gra2gms(gra: float) -> list:
     """
     deg = gra2deg(gra)
     gra = deg // 1
-    min = gra % 1
+    minu = gra % 1
     gra = gra // 1
-    min *= 60
+    minu *= 60
-    sek = min % 1
+    sek = minu % 1
-    min = min // 1
+    minu = minu // 1
     sek *= 60
-    return [gra, min, sek]
+    return [gra, minu, sek]
 
 
 def gra2rad(gra: float) -> float:
@@ -113,13 +114,13 @@ def rad2gms(rad: float) -> list:
     :rtype: list
     """
     deg = rad2deg(rad)
-    min = deg % 1
+    minu = deg % 1
     gra = deg // 1
-    min *= 60
+    minu *= 60
-    sek = min % 1
+    sek = minu % 1
-    min = min // 1
+    minu = minu // 1
     sek *= 60
-    return [gra, min, sek]
+    return [gra, minu, sek]
 
 
 def gms2rad(gms: list) -> float: