.

GHA_triaxial/gha1_approx.py

@@ -21,6 +21,7 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
     alphas = [alpha0]
     s_curr = 0.0
     last_sigma = None
+    last_p = None
 
     while s_curr < s:
         ds_step = min(ds, s - s_curr)
@@ -30,9 +31,13 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
         p1 = points[-1]
         alpha1 = alphas[-1]
 
-        sigma = func_sigma_ell(ell, p1, alpha1)
-        if last_sigma is not None and np.linalg.norm(sigma - last_sigma) > 0.2:
-            raise Exception("GHA1_approx: Plötzlicher Richtungswechsel")
+        p, q = pq_ell(ell, p1)
+        if last_p is not None and np.dot(p, last_p) < 0:
+            p = -p
+            q = -q
+        sigma = p * sin(alpha1) + q * cos(alpha1)
+        if last_sigma is not None and np.dot(sigma, last_sigma) < 0:
+            sigma = -sigma
         p2 = p1 + ds_step * sigma
         p2 = ell.point_onto_ellipsoid(p2)
 
@@ -82,10 +87,10 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray):
 
 if __name__ == '__main__':
     ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    P0 = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(1))
+    P0 = ell.ell2cart(wu.deg2rad(89), wu.deg2rad(1))
     alpha0 = wu.deg2rad(2)
-    s = 100000
-    P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
+    s = 200000
+    P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=100, all_points=True)
     P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=20, maxPartCircum=2)
-    show_points(points, P0, P1_ana)
     print(np.linalg.norm(P1_app - P1_ana))
+    show_points(points, P0, P1_ana)

GHA_triaxial/utils.py

@@ -128,21 +128,27 @@ def pq_ell(ell: EllipsoidTriaxial, point: NDArray) -> Tuple[NDArray, NDArray]:
     n = ell.func_n(point)
 
     beta, lamb = ell.cart2ell(point)
-    B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2
-    L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2
+    if abs(cos(beta)) < 1e-12 and abs(np.sin(lamb)) < 1e-12:
+        if beta > 0:
+            p = np.array([0, -1, 0])
+        else:
+            p = np.array([0, 1, 0])
+    else:
+        B = ell.Ex ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(beta) ** 2
+        L = ell.Ex ** 2 - ell.Ee ** 2 * cos(lamb) ** 2
 
-    c1 = x ** 2 + y ** 2 + z ** 2 - (ell.ax ** 2 + ell.ay ** 2 + ell.b ** 2)
-    c0 = (ell.ax ** 2 * ell.ay ** 2 + ell.ax ** 2 * ell.b ** 2 + ell.ay ** 2 * ell.b ** 2 -
-          (ell.ay ** 2 + ell.b ** 2) * x ** 2 - (ell.ax ** 2 + ell.b ** 2) * y ** 2 - (
-                  ell.ax ** 2 + ell.ay ** 2) * z ** 2)
-    t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
+        c1 = x ** 2 + y ** 2 + z ** 2 - (ell.ax ** 2 + ell.ay ** 2 + ell.b ** 2)
+        c0 = (ell.ax ** 2 * ell.ay ** 2 + ell.ax ** 2 * ell.b ** 2 + ell.ay ** 2 * ell.b ** 2 -
+              (ell.ay ** 2 + ell.b ** 2) * x ** 2 - (ell.ax ** 2 + ell.b ** 2) * y ** 2 - (
+                      ell.ax ** 2 + ell.ay ** 2) * z ** 2)
+        t2 = (-c1 + sqrt(c1 ** 2 - 4 * c0)) / 2
 
-    F = ell.Ey ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(lamb) ** 2
-    p1 = -sqrt(L / (F * t2)) * ell.ax / ell.Ex * sqrt(B) * sin(lamb)
-    p2 = sqrt(L / (F * t2)) * ell.ay * cos(beta) * cos(lamb)
-    p3 = 1 / sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * sin(beta) * sin(2 * lamb)
-    p = np.array([p1, p2, p3])
-    p = p / np.linalg.norm(p)
+        F = ell.Ey ** 2 * cos(beta) ** 2 + ell.Ee ** 2 * sin(lamb) ** 2
+        p1 = -sqrt(L / (F * t2)) * ell.ax / ell.Ex * sqrt(B) * sin(lamb)
+        p2 = sqrt(L / (F * t2)) * ell.ay * cos(beta) * cos(lamb)
+        p3 = 1 / sqrt(F * t2) * (ell.b * ell.Ee ** 2) / (2 * ell.Ex) * sin(beta) * sin(2 * lamb)
+        p = np.array([p1, p2, p3])
+        p = p / np.linalg.norm(p)
     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]])