Merge remote-tracking branch 'origin/main'

GHA_triaxial/gha2_num.py

@@ -1,24 +1,19 @@
 import numpy as np
 from ellipsoide import EllipsoidTriaxial
-import runge_kutta as rk
+from runge_kutta import rk4, rk4_step, rk4_end, rk4_integral
 import GHA_triaxial.numeric_examples_karney as ne_karney
 import GHA_triaxial.numeric_examples_panou as ne_panou
 import winkelumrechnungen as wu
 from typing import Tuple
 from numpy.typing import NDArray
 import ausgaben as aus
-from utils_angle import cot, wrap_to_pi
+from utils_angle import cot, arccot, wrap_to_pi
 
 
-def arccot(x):
+def norm_a(a):
-    x = np.asarray(x)
+    if a < 0.0:
-    a = np.arctan2(1.0, x)
+        a += np.pi
-    return np.where(x < 0.0, a - np.pi, a)
+    return a
-
-def normalize_alpha_0_pi(alpha):
-    if alpha < 0.0:
-        alpha += np.pi
-    return alpha
 
 def sph_azimuth(beta1, lam1, beta2, lam2):
     dlam = wrap_to_pi(lam2 - lam1)
@@ -32,17 +27,29 @@ def sph_azimuth(beta1, lam1, beta2, lam2):
 # Panou 2013
 def gha2_num(
     ell: EllipsoidTriaxial,
+    beta_0: float,
+    lamb_0: float,
     beta_1: float,
     lamb_1: float,
-    beta_2: float,
-    lamb_2: float,
     n: int = 16000,
     epsilon: float = 10**-12,
     iter_max: int = 30,
     all_points: bool = False,
 ) -> Tuple[float, float, float] | Tuple[float, float, float, NDArray, NDArray]:
+    """
+    :param ell: Ellipsoid
+    :param beta_0: Beta Punkt 0
+    :param lamb_0: Lambda Punkt 0
+    :param beta_1: Beta Punkt 1
+    :param lamb_1: Lambda Punkt 1
+    :param n: Anzahl Schritte
+    :param epsilon: Genauigkeit
+    :param iter_max: Maximale Iterationen
+    :param all_points: Ausgabe aller Punkte
+    :return: Azimut Startpunkt, Azumit Zielpunkt, Strecke
+    """
 
-
+    # Berechnung Koeffizienten, Gaußschen Fundamentalgrößen 1. Ordnung sowie deren Ableitungen
     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
@@ -123,6 +130,7 @@ def gha2_num(
             G_lamb_lamb,
         )
 
+    # Berechnung der ODE Koeffizienten für Fall 1 (lambda_0 != lambda_1)
     def p_coef(beta, lamb):
         (
             BETA,
@@ -161,6 +169,7 @@ def gha2_num(
 
         return (BETA, LAMBDA, E, G, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00)
 
+    # Berechnung der ODE Koeffizienten für Fall 2 (lambda_0 == lambda_1)
     def q_coef(beta, lamb):
         (
             BETA,
@@ -197,69 +206,7 @@ def gha2_num(
         )
         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)
+        return BETA, LAMBDA, E, G, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00
-
-    def rk4_last(f, t0, y0, dt, N):
-        h = dt / N
-        t = t0
-        y = np.array(y0, dtype=float, copy=True)
-        for _ in range(N):
-            k1 = f(t, y)
-            k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
-            k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
-            k4 = f(t + h, y + h * k3)
-            y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
-            t = t + h
-        return t, y
-
-    def rk4_last_with_integral(f, t0, y0, dt, N, integrand_at):
-
-        h = dt / N
-        habs = abs(h)
-        t = t0
-        y = np.array(y0, dtype=float, copy=True)
-
-        if N % 2 == 0:
-            # Simpson streaming
-            f0 = integrand_at(t, y)
-            odd_sum = 0.0
-            even_sum = 0.0
-
-            for i in range(1, N + 1):
-                k1 = f(t, y)
-                k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
-                k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
-                k4 = f(t + h, y + h * k3)
-                y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
-                t = t + h
-
-                fi = integrand_at(t, y)
-                if i == N:
-                    fN = fi
-                elif i % 2 == 1:
-                    odd_sum += fi
-                else:
-                    even_sum += fi
-
-            S = f0 + fN + 4.0 * odd_sum + 2.0 * even_sum
-            s = (habs / 3.0) * S
-            return t, y, s
-
-        # Trapez streaming
-        f_prev = integrand_at(t, y)
-        acc = 0.0
-        for _ in range(N):
-            k1 = f(t, y)
-            k2 = f(t + 0.5 * h, y + 0.5 * h * k1)
-            k3 = f(t + 0.5 * h, y + 0.5 * h * k2)
-            k4 = f(t + h, y + h * k3)
-            y = y + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
-            t = t + h
-            f_cur = integrand_at(t, y)
-            acc += 0.5 * (f_prev + f_cur)
-            f_prev = f_cur
-        s = habs * acc
-        return t, y, s
 
     def integrand_lambda(lamb, y):
         beta = y[0]
@@ -273,42 +220,41 @@ def gha2_num(
         (_, _, E, G, *_) = BETA_LAMBDA(beta, lamb)
         return np.sqrt(E + G * lamb_p**2)
 
-    if lamb_1 != lamb_2:
+    # Fall 1 (lambda_0 != lambda_1)
-        N = n
+    if abs(lamb_1 - lamb_0) >= 1e-15:
-        dlamb = lamb_2 - lamb_1
+        N = int(n)
+        dlamb = float(lamb_1 - lamb_0)
 
-        def buildODElamb():
+        beta0 = float(beta_0)
-            def ODE(lamb, v):
+        lamb0 = float(lamb_0)
+        beta1 = float(beta_1)
+        lamb1 = float(lamb_1)
+
+        def ode_lamb(lamb, v):
             beta, beta_p, X3, X4 = v
             (_, _, _, _, p_3, p_2, p_1, p_0, p_33, p_22, p_11, p_00) = p_coef(beta, lamb)
 
             dbeta = beta_p
             dbeta_p = p_3 * beta_p**3 + p_2 * beta_p**2 + p_1 * beta_p + p_0
             dX3 = X4
-                dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (
+            dX4 = (p_33 * beta_p**3 + p_22 * beta_p**2 + p_11 * beta_p + p_00) * X3 + (3*p_3*beta_p**2 + 2*p_2*beta_p + p_1) * X4
-                    3 * p_3 * beta_p**2 + 2 * p_2 * beta_p + p_1
-                ) * X4
             return np.array([dbeta, dbeta_p, dX3, dX4], dtype=float)
 
-            return ODE
+        alpha0_sph = sph_azimuth(beta0, lamb0, beta1, lamb1)
-
+        (_, _, E0, G0, *_) = BETA_LAMBDA(beta0, lamb0)
-        ode_lamb = buildODElamb()
+        beta_p0_sph = np.sqrt(G0 / E0) * cot(alpha0_sph)
-
-        alpha0_sph = sph_azimuth(beta_1, lamb_1, beta_2, lamb_2)
-        (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
-        beta_p0_sph = np.sqrt(G1 / E1) * cot(alpha0_sph) if abs(dlamb) >= 1e-15 else 0.0
-
 
         N_newton = min(N, 4000)
 
         def solve_newton(beta_p0_init: float):
             beta_p0 = float(beta_p0_init)
             for _ in range(iter_max):
-                startwerte = np.array([beta_1, beta_p0, 0.0, 1.0], dtype=float)
+                v0 = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
-                _, y_end = rk4_last(ode_lamb, lamb_1, startwerte, dlamb, N_newton)
+                _, y_end = rk4_end(ode_lamb, lamb0, v0, dlamb, N_newton)
-                beta_end, _, X3_end, _ = y_end
+
+                beta_end, _, X3_end, _ = y_end
+                delta = beta_end - beta1
 
-                delta = beta_end - beta_2
                 if abs(delta) < epsilon:
                     return True, beta_p0
 
@@ -316,59 +262,45 @@ def gha2_num(
                     return False, None
 
                 step = delta / X3_end
-                max_step = 0.5
+                step = np.clip(step, -0.5, 0.5)
-                if abs(step) > max_step:
+                beta_p0 -= step
-                    step = np.sign(step) * max_step
-                beta_p0 = beta_p0 - step
 
             return False, None
 
         ok, beta_p0_sol = solve_newton(beta_p0_sph)
 
         if not ok:
-
             candidates = [-beta_p0_sph, 0.5 * beta_p0_sph, 2.0 * beta_p0_sph]
             N_quick = min(N, 2000)
             best = None
-
             for g in candidates:
-                ok_g, beta_p0_sol_g = solve_newton(g)
+                ok_g, sol = solve_newton(g)
                 if not ok_g:
                     continue
-
+                v0_g = np.array([beta0, sol, 0.0, 1.0], dtype=float)
-                startwerte_g = np.array([beta_1, beta_p0_sol_g, 0.0, 1.0], dtype=float)
+                _, _, s_quick = rk4_integral(ode_lamb, lamb0, v0_g, dlamb, N_quick, integrand_lambda)
-                _, _, s_quick = rk4_last_with_integral(
-                    ode_lamb, lamb_1, startwerte_g, dlamb, N_quick, integrand_lambda
-                )
-
                 if (best is None) or (s_quick < best[0]):
-                    best = (s_quick, beta_p0_sol_g)
+                    best = (s_quick, sol)
-
             if best is None:
-                raise RuntimeError("Keine Startwert-Variante konvergiert.")
+                raise RuntimeError("Keine Startwert-Variante konvergiert (lambda-Fall).")
-
             beta_p0_sol = best[1]
 
-        beta_0 = beta_p0_sol
+        beta_p0 = float(beta_p0_sol)
-        startwerte_final = np.array([beta_1, beta_0, 0.0, 1.0], dtype=float)
+        v0_final = np.array([beta0, beta_p0, 0.0, 1.0], dtype=float)
 
         if all_points:
-            lamb_list, states = rk.rk4(ode_lamb, lamb_1, startwerte_final, dlamb, N, False)
+            lamb_list, states = rk4(ode_lamb, lamb0, v0_final, dlamb, N, False)
-
             lamb_arr = np.array(lamb_list, dtype=float)
             beta_arr = np.array([st[0] for st in states], dtype=float)
             beta_p_arr = np.array([st[1] for st in states], dtype=float)
 
-            (_, _, E1, G1, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+            (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
-            (_, _, E2, G2, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
+            (_, _, E_end,   G_end,   *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
 
-            alpha_1 = arccot(np.sqrt(E1 / G1) * beta_p_arr[0])
+            alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p_arr[0]))
-            alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_arr[-1])
+            alpha_2 = norm_a(arccot(np.sqrt(E_end   / G_end)   * beta_p_arr[-1]))
 
-            alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+            # Distanz aus Arrays
-            alpha_2 = normalize_alpha_0_pi(float(alpha_2))
-
-            # Distanz s aus Arrays (Simpson/Trapz)
             integrand = np.zeros(N + 1, dtype=float)
             for i in range(N + 1):
                 (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
@@ -381,91 +313,75 @@ def gha2_num(
             else:
                 s = np.trapz(integrand, dx=h)
 
-            return alpha_1, alpha_2, s, beta_arr, lamb_arr
+            return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr
 
-        _, y_end, s = rk4_last_with_integral(ode_lamb, lamb_1, startwerte_final, dlamb, N, integrand_lambda)
+        _, y_end, s = rk4_integral(ode_lamb, lamb0, v0_final, dlamb, N, integrand_lambda)
         beta_end, beta_p_end, _, _ = y_end
 
-        (_, _, E1, G1, *_) = BETA_LAMBDA(beta_1, lamb_1)
+        (_, _, E_start, G_start, *_) = BETA_LAMBDA(beta0, lamb0)
-        (_, _, E2, G2, *_) = BETA_LAMBDA(beta_2, lamb_2)
+        (_, _, E_end,   G_end,   *_) = BETA_LAMBDA(beta1, lamb1)
 
-        alpha_1 = arccot(np.sqrt(E1 / G1) * beta_0)
+        alpha_1 = norm_a(arccot(np.sqrt(E_start / G_start) * beta_p0))
-        alpha_2 = arccot(np.sqrt(E2 / G2) * beta_p_end)
+        alpha_2 = norm_a(arccot(np.sqrt(E_end   / G_end)   * beta_p_end))
 
-        alpha_1 = normalize_alpha_0_pi(float(alpha_1))
+        return float(alpha_1), float(alpha_2), float(s)
-        alpha_2 = normalize_alpha_0_pi(float(alpha_2))
 
-        return alpha_1, alpha_2, s
+    # Fall 2 (lambda_0 == lambda_1)
-
+    N = int(n)
-    N = n
+    dbeta = float(beta_1 - beta_0)
-    dbeta = beta_2 - beta_1
 
     if abs(dbeta) < 1e-15:
         if all_points:
             return 0.0, 0.0, 0.0, np.array([]), np.array([])
         return 0.0, 0.0, 0.0
 
-    # ODE-System (lambda, lambda', Y3, Y4) in Abhängigkeit von beta
+    beta0 = float(beta_0)
-    def buildODEbeta():
+    lamb0 = float(lamb_0)
-        def ODE(beta, v):
+    beta1 = float(beta_1)
+    lamb1 = float(lamb_1)
+
+    def ode_beta(beta, v):
         lamb, lamb_p, Y3, Y4 = v
         (_, _, _, _, q_3, q_2, q_1, q_0, q_33, q_22, q_11, q_00) = q_coef(beta, lamb)
 
         dlamb = lamb_p
         dlamb_p = q_3*lamb_p**3 + q_2*lamb_p**2 + q_1*lamb_p + q_0
         dY3 = Y4
-            dY4 = (q_33 * lamb_p**3 + q_22 * lamb_p**2 + q_11 * lamb_p + q_00) * Y3 + (
+        dY4 = (q_33*lamb_p**3 + q_22*lamb_p**2 + q_11*lamb_p + q_00)*Y3 + (3*q_3*lamb_p**2 + 2*q_2*lamb_p + q_1)*Y4
-                3 * q_3 * lamb_p**2 + 2 * q_2 * lamb_p + q_1
-            ) * Y4
-
         return np.array([dlamb, dlamb_p, dY3, dY4], dtype=float)
 
-        return ODE
+    lamb_p0 = 0.0
-
-    ode_beta = buildODEbeta()
-
-    # Newton auf lambda'_0
-    lamb_0 = 0.0
     for _ in range(iter_max):
-        startwerte = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float)
+        v0 = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
-        beta_list, states = rk.rk4(ode_beta, beta_1, startwerte, dbeta, N, False)
+        _, y_end = rk4_end(ode_beta, beta0, v0, dbeta, N)
 
-        lamb_end, lamb_p_end, Y3_end, _ = states[-1]
+        lamb_end, _, Y3_end, _ = y_end
-        delta = lamb_end - lamb_2
+        delta = lamb_end - lamb1
 
         if abs(delta) < epsilon:
             break
 
         if abs(Y3_end) < 1e-20:
-            raise RuntimeError("Abbruch (Ableitung ~ 0).")
+            raise RuntimeError("Abbruch (Ableitung ~ 0) im beta-Fall.")
 
         step = delta / Y3_end
-        max_step = 1.0
+        step = np.clip(step, -1.0, 1.0)
-        if abs(step) > max_step:
+        lamb_p0 -= step
-            step = np.sign(step) * max_step
-        lamb_0 = lamb_0 - step
 
-    startwerte_final = np.array([lamb_1, lamb_0, 0.0, 1.0], dtype=float)
+    v0_final = np.array([lamb0, lamb_p0, 0.0, 1.0], dtype=float)
 
     if all_points:
-        beta_list, states = rk.rk4(ode_beta, beta_1, startwerte_final, dbeta, N, False)
+        beta_list, states = rk4(ode_beta, beta0, v0_final, dbeta, N, False)
-
         beta_arr = np.array(beta_list, dtype=float)
         lamb_arr = np.array([st[0] for st in states], dtype=float)
         lamb_p_arr = np.array([st[1] for st in states], dtype=float)
 
-        # Azimute
+        (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta_arr[0],  lamb_arr[0])
-        (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_arr[0], lamb_arr[0])
+        (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
-        (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_arr[-1], lamb_arr[-1])
 
-        alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_p_arr[0])
+        alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p_arr[0]))
-        alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_arr[-1])
+        alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_arr[-1]))
 
-        # optionaler Quadrantenfix (robust)
-        alpha_1 = normalize_alpha_0_pi(float(alpha_1))
-        alpha_2 = normalize_alpha_0_pi(float(alpha_2))
-
-        # Distanz
         integrand = np.zeros(N + 1, dtype=float)
         for i in range(N + 1):
             (_, _, Ei, Gi, *_) = BETA_LAMBDA(beta_arr[i], lamb_arr[i])
@@ -478,32 +394,29 @@ def gha2_num(
         else:
             s = np.trapz(integrand, dx=h)
 
-        return alpha_1, alpha_2, s, beta_arr, lamb_arr
+        return float(alpha_1), float(alpha_2), float(s), beta_arr, lamb_arr
 
-    # all_points == False: streaming Integral
+    _, y_end, s = rk4_integral(ode_beta, beta0, v0_final, dbeta, N, integrand_beta)
-    _, y_end, s = rk4_last_with_integral(ode_beta, beta_1, startwerte_final, dbeta, N, integrand_beta)
     lamb_end, lamb_p_end, _, _ = y_end
 
-    (BETA1, LAMBDA1, _, _, *_) = BETA_LAMBDA(beta_1, lamb_1)
+    (BETA_s, LAMBDA_s, _, _, *_) = BETA_LAMBDA(beta0, lamb0)
-    (BETA2, LAMBDA2, _, _, *_) = BETA_LAMBDA(beta_2, lamb_2)
+    (BETA_e, LAMBDA_e, _, _, *_) = BETA_LAMBDA(beta1, lamb1)
 
-    alpha_1 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA1 / BETA1) * lamb_0)
+    alpha_1 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_s / BETA_s) * lamb_p0))
-    alpha_2 = (np.pi / 2.0) - arccot(np.sqrt(LAMBDA2 / BETA2) * lamb_p_end)
+    alpha_2 = norm_a((np.pi/2.0) - arccot(np.sqrt(LAMBDA_e / BETA_e) * lamb_p_end))
-
-    alpha_1 = normalize_alpha_0_pi(float(alpha_1))
-    alpha_2 = normalize_alpha_0_pi(float(alpha_2))
-
-    return alpha_1, alpha_2, s
 
+    return float(alpha_1), float(alpha_2), float(s)
 
 if __name__ == "__main__":
-    ell = EllipsoidTriaxial.init_name("BursaSima1980round")
+    # ell = EllipsoidTriaxial.init_name("BursaSima1980round")
-    beta1 = np.deg2rad(75)
+    # beta1 = np.deg2rad(75)
-    lamb1 = np.deg2rad(-90)
+    # lamb1 = np.deg2rad(-90)
-    beta2 = np.deg2rad(75)
+    # beta2 = np.deg2rad(75)
-    lamb2 = np.deg2rad(66)
+    # lamb2 = np.deg2rad(66)
-    a1, a2, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
+    # a0, a1, s = gha2_num(ell, beta1, lamb1, beta2, lamb2, n=5000)
-    print(aus.gms("a1", a1, 4))
+    # print(aus.gms("a0", a0, 4))
+    # print(aus.gms("a1", a1, 4))
+    # print("s: ", s)
     # # print(aus.gms("a2", a2, 4))
     # # print(s)
     # cart1 = ell.para2cart(0, 0)

runge_kutta.py

@@ -41,3 +41,83 @@ def rk4_step(ode, t: float, v: np.ndarray, h: float) -> np.ndarray:
     k3 = ode(t + 0.5 * h, v + 0.5 * h * k2)
     k4 = ode(t + h, v + h * k3)
     return v + (h / 6.0) * (k1 + 2 * k2 + 2 * k3 + k4)
+
+def rk4_end(ode, t0: float, v0: np.ndarray, weite: float, schritte: int, fein: bool = False):
+    h = weite / schritte
+    t = float(t0)
+    v = np.array(v0, dtype=float, copy=True)
+
+    for _ in range(schritte):
+        if not fein:
+            v_next = rk4_step(ode, t, v, h)
+        else:
+            v_grob = rk4_step(ode, t, v, h)
+            v_half = rk4_step(ode, t, v, 0.5 * h)
+            v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+            v_next = v_fein + (v_fein - v_grob) / 15.0
+
+        t += h
+        v = v_next
+
+    return t, v
+
+# RK4 mit Simpson bzw. Trapez
+def rk4_integral( ode, t0: float, v0: np.ndarray, weite: float, schritte: int, integrand_at, fein: bool = False, simpson: bool = True, ):
+
+    h = weite / schritte
+    habs = abs(h)
+
+    t = float(t0)
+    v = np.array(v0, dtype=float, copy=True)
+
+    if simpson and (schritte % 2 == 0):
+        f0 = float(integrand_at(t, v))
+        odd_sum = 0.0
+        even_sum = 0.0
+        fN = None
+
+        for i in range(1, schritte + 1):
+            if not fein:
+                v_next = rk4_step(ode, t, v, h)
+            else:
+                v_grob = rk4_step(ode, t, v, h)
+                v_half = rk4_step(ode, t, v, 0.5 * h)
+                v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+                v_next = v_fein + (v_fein - v_grob) / 15.0
+
+            t += h
+            v = v_next
+
+            fi = float(integrand_at(t, v))
+            if i == schritte:
+                fN = fi
+            elif i % 2 == 1:
+                odd_sum += fi
+            else:
+                even_sum += fi
+
+        S = f0 + fN + 4.0 * odd_sum + 2.0 * even_sum
+        s = (habs / 3.0) * S
+        return t, v, s
+
+    f_prev = float(integrand_at(t, v))
+    acc = 0.0
+
+    for _ in range(schritte):
+        if not fein:
+            v_next = rk4_step(ode, t, v, h)
+        else:
+            v_grob = rk4_step(ode, t, v, h)
+            v_half = rk4_step(ode, t, v, 0.5 * h)
+            v_fein = rk4_step(ode, t + 0.5 * h, v_half, 0.5 * h)
+            v_next = v_fein + (v_fein - v_grob) / 15.0
+
+        t += h
+        v = v_next
+
+        f_cur = float(integrand_at(t, v))
+        acc += 0.5 * (f_prev + f_cur)
+        f_prev = f_cur
+
+    s = habs * acc
+    return t, v, s