Exceptions einheitlich

GHA_triaxial/gha1_ES.py

@@ -245,7 +245,7 @@ def gha1_ES(ell: EllipsoidTriaxial, beta0: float, omega0: float, alpha0: float,
         P_all.append(P)
         alpha_end.append(alpha)
         if step > nsteps_est + 50:
-            raise RuntimeError("Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
+            raise RuntimeError("GHA1_ES: Zu viele Schritte – vermutlich Konvergenzproblem / falsche Azimut-Konvention.")
     Pk = P_all[-1]
     alpha1 = float(alpha_end[-1])
 

GHA_triaxial/gha1_ana.py

@@ -132,7 +132,7 @@ def gha1_ana(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, ma
 
     _, _, h = ell.cart2geod(point_end, "ligas3")
     if h > 1e-5:
-        raise Exception("Analytische Methode ist explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
+        raise Exception("GHA1_ana: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
 
     return point_end, alpha_end
 

GHA_triaxial/gha1_approx.py

@@ -20,6 +20,7 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
     points = [p0]
     alphas = [alpha0]
     s_curr = 0.0
+    last_sigma = None
 
     while s_curr < s:
         ds_step = min(ds, s - s_curr)
@@ -30,6 +31,8 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
         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")
         p2 = p1 + ds_step * sigma
         p2 = ell.point_onto_ellipsoid(p2)
 
@@ -43,7 +46,7 @@ def gha1_approx(ell: EllipsoidTriaxial, p0: np.ndarray, alpha0: float, s: float,
 
         ds_step = np.linalg.norm(p2 - p1)
         s_curr += ds_step
-        if s_curr > 10000000:
+        last_sigma = sigma
         pass
 
     if all_points:
@@ -79,10 +82,10 @@ def show_points(points: NDArray, p0: NDArray, p1: NDArray):
 
 if __name__ == '__main__':
     ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    P0 = ell.para2cart(0.2, 0.3)
+    P0 = ell.ell2cart(wu.deg2rad(10), wu.deg2rad(1))
-    alpha0 = wu.deg2rad(35)
+    alpha0 = wu.deg2rad(2)
-    s = 13000000
+    s = 100000
     P1_app, alpha1_app, points, alphas = gha1_approx(ell, P0, alpha0, s, ds=10000, all_points=True)
-    P1_ana, alpha1_ana = gha1_ana(ell, P0, alpha0, s, maxM=60, maxPartCircum=16)
+    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))

GHA_triaxial/gha1_num.py

@@ -78,6 +78,10 @@ def gha1_num(ell: EllipsoidTriaxial, point: NDArray, alpha0: float, s: float, nu
     if alpha1 < 0:
         alpha1 += 2 * np.pi
 
+    _, _, h = ell.cart2geod(point1, "ligas3")
+    if h > 1e-5:
+        raise Exception("GHA1_num: explodiert, Punkt liegt nicht mehr auf dem Ellipsoid")
+
     if all_points:
         return point1, alpha1, werte
     else:

GHA_triaxial/gha2_ES.py

@@ -109,7 +109,7 @@ def gha2_ES(ell: EllipsoidTriaxial, P0: NDArray, Pk: NDArray, maxSegLen: float =
                 startIter += 1
                 maxIter = 10000
                 if startIter > maxIter:
-                    raise RuntimeError("Abbruch: maximale Iterationen überschritten.")
+                    raise RuntimeError("GHA2_ES: maximale Iterationen überschritten")
             new_points.append(B)
 
         points = new_points

GHA_triaxial/gha2_num.py

@@ -282,7 +282,7 @@ def gha2_num(
                 if (best is None) or (s_quick < best[0]):
                     best = (s_quick, sol)
             if best is None:
-                raise RuntimeError("Keine Startwert-Variante konvergiert (lambda-Fall).")
+                raise RuntimeError("GHA2_num: Keine Startwert-Variante konvergiert (lambda-Fall)")
             beta_p0_sol = best[1]
 
         beta_p0 = float(beta_p0_sol)
@@ -362,7 +362,7 @@ def gha2_num(
             break
 
         if abs(Y3_end) < 1e-20:
-            raise RuntimeError("Abbruch (Ableitung ~ 0) im beta-Fall.")
+            raise RuntimeError("GHA2_num: Ableitung ~ 0 im beta-Fall")
 
         step = delta / Y3_end
         step = np.clip(step, -1.0, 1.0)

Tests/algorithms_test.ipynb

@@ -114,6 +114,7 @@
     "# test = \"Karney\"\n",
     "# test = \"Panou\"\n",
     "test = \"Random\"\n",
+    "\n",
     "if test == \"Karney\":\n",
     "    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n",
     "    examples = get_examples_karney(4, seed=42)\n",
@@ -128,9 +129,9 @@
     "            examples.append(example)\n",
     "elif test == \"Random\":\n",
     "    ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n",
-    "    examples = [] # beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s\n",
+    "    examples = []\n",
     "    random.seed(42)\n",
-    "    for _ in range(30):\n",
+    "    for _ in range(50):\n",
     "        beta0 = wu.deg2rad(random.randint(-90, 90))\n",
     "        lamb0 = wu.deg2rad(random.randint(-179, 180))\n",
     "        alpha0_ell = wu.deg2rad(random.randint(0, 359))\n",
@@ -151,16 +152,36 @@
    },
    "cell_type": "code",
    "source": [
-    "if test != \"Random\":\n",
     "results = {}\n",
     "for i, example in enumerate(examples):\n",
     "    print(f\"----- Beispiel {i+1}/{len(examples)}\")\n",
     "    example_results = {}\n",
     "\n",
+    "    if test != \"Random\":\n",
     "        beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n",
     "        P0 = ell.ell2cart(beta0, lamb0)\n",
     "        P1 = ell.ell2cart(beta1, lamb1)\n",
     "        _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
+    "    else:\n",
+    "        beta0, lamb0, alpha0_ell, s = example\n",
+    "        P0 = ell.ell2cart(beta0, lamb0)\n",
+    "        _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
+    "\n",
+    "        # try:\n",
+    "            # P1, alpha1_para = gha1_ana(ell, P0, alpha0_para, s, maxM=80, maxPartCircum=8)\n",
+    "            # beta1, lamb1 = ell.cart2ell(P1)\n",
+    "            # u1, v1 = ell.cart2para(P1)\n",
+    "            # _, _, alpha1_ell = alpha_para2ell(ell, u1, v1, alpha1_para)\n",
+    "        # except Exception as e:\n",
+    "        #     print(\"Referenz-Berechnung analytisch fehlgeschlagen: \", e)\n",
+    "        #     continue\n",
+    "\n",
+    "        try:\n",
+    "            P1, alpha1_ell = gha1_num(ell, P0, alpha0_ell, s, num=10000)\n",
+    "            beta1, lamb1 = ell.cart2ell(P1)\n",
+    "        except Exception as e:\n",
+    "            print(\"Referenz-Berechnung numerisch fehlgeschlagen: \", e)\n",
+    "            continue\n",
     "\n",
     "    for steps in steps_gha1_num:\n",
     "        start = time.perf_counter()\n",
@@ -303,131 +324,6 @@
    ],
    "execution_count": 15
   },
-  {
-   "metadata": {},
-   "cell_type": "code",
-   "source": [
-    "if test == \"Random\":\n",
-    "    results = {}\n",
-    "    for i, example in enumerate(examples):\n",
-    "        print(f\"----- Beispiel {i+1}/{len(examples)}\")\n",
-    "        example_results = {}\n",
-    "\n",
-    "        beta0, lamb0, alpha0_ell, s = example\n",
-    "        P0 = ell.ell2cart(beta0, lamb0)\n",
-    "        _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
-    "\n",
-    "        # P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=60, maxPartCircum=4)\n",
-    "        # beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
-    "        # _, _, alpha1_ana = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
-    "\n",
-    "        P1, alpha1_ell = gha1_num(ell, P0, alpha0_ell, s, num=10000)\n",
-    "        beta1, lamb1 = ell.cart2ell(P1)\n",
-    "\n",
-    "        for maxM in maxM_gha1_ana:\n",
-    "            for parts in parts_gha1_ana:\n",
-    "                start = time.perf_counter()\n",
-    "                try:\n",
-    "                    P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n",
-    "                    end = time.perf_counter()\n",
-    "                    beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
-    "                    _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
-    "                    d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n",
-    "                    d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n",
-    "                    d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n",
-    "                    d_time = end - start\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "                except Exception as e:\n",
-    "                    print(e)\n",
-    "                    example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha1_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n",
-    "                d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n",
-    "                d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha1_ES:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                P1_ES, alpha1_ES = gha1_ES(ell, beta0, lamb0, alpha0_ell, s, maxSegLen=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                beta1_ES, lamb1_ES = ell.cart2ell(P1_ES)\n",
-    "                d_beta1 = abs(beta1_ES - beta1)\n",
-    "                d_lamb1 = abs(lamb1_ES - lamb1)\n",
-    "                d_alpha1 = abs(alpha1_ES - alpha1_ell)\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA1_ES_{dsPart}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA1_ES_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for steps in steps_gha2_num:\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                with warnings.catch_warnings():\n",
-    "                    warnings.simplefilter(\"ignore\", RuntimeWarning)\n",
-    "                    alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n",
-    "                d_s = abs(s_num - s) / 1000\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha2_ES:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                with suppress_print():\n",
-    "                    alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n",
-    "                d_s = abs(s_ES - s) / 1000\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_ES_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_ES_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        for dsPart in dsPart_gha2_approx:\n",
-    "            ds = ell.ax/dsPart\n",
-    "            start = time.perf_counter()\n",
-    "            try:\n",
-    "                alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n",
-    "                end = time.perf_counter()\n",
-    "                d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n",
-    "                d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
-    "                d_s = abs(s_approx - s) / 1000\n",
-    "                d_time = end - start\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
-    "            except Exception as e:\n",
-    "                print(e)\n",
-    "                example_results[f\"GHA2_approx_{dsPart}\"] = (nan, nan, nan, nan)\n",
-    "\n",
-    "        results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results\n",
-    "\n",
-    "\n"
-   ],
-   "id": "6e20a077d522d5",
-   "outputs": [],
-   "execution_count": null
-  },
   {
    "metadata": {},
    "cell_type": "code",

ellipsoide.py

@@ -412,12 +412,12 @@ class EllipsoidTriaxial:
                 i += 1
 
             if i == maxI:
-                raise Exception("Umrechnung ist nicht konvergiert")
+                raise Exception("Umrechnung cart2ell: nicht konvergiert")
 
             point_n = self.ell2cart(beta, lamb)
             delta_r = np.linalg.norm(point - point_n, axis=-1)
             if delta_r > 1e-6:
-                raise Exception("Fehler in der Umrechnung cart2ell")
+                raise Exception("Umrechnung cart2ell: Punktdifferenz")
 
             return beta, lamb
 
@@ -511,7 +511,7 @@ class EllipsoidTriaxial:
             i += 1
 
         if i == maxI:
-            raise Exception("Umrechung ist nicht konvergiert")
+            raise Exception("Umrechnung cart2ell: nicht konvergiert")
 
         return phi, lamb