Of course ✅
Let's execute this Python script and render the resulting plot 👇

Here’s the trajectory plot generated by your script ✅🌀

Good — I know the document you're referring to. It describes a Newtonian model with an extra term (essentially a non-zero $\frac{dF}{dt}$) that causes a slow precession of the orbit, mimicking Mercury’s perihelion advance of

To match this in the simulation, we need to tune the extra angular acceleration term, which in your script is mainly controlled by the variable:

python
Code kopiëren
ap = 6.00e-7

and more subtly through this line:

python
Code kopiëren
rpp = Cstefan / (r ** 2) + (r * (ap + 0.0000001) ** 2)

The + 0.0000001 here is essentially a small perturbation to the orbital angular velocity — this is where Stefan Boersen introduces a non-Newtonian contribution.

👉 Strategy to achieve 532″/century precession

1. Simulate a long enough time to cover ~1 century (or multiple orbits).

2. Track the perihelion angle $a$ (argument of periapsis) at each revolution.

3. Adjust the perturbation constant (currently 0.0000001) until the precession per century is ~532″.

This is best done numerically. Below is an adjusted version of your script that:

• Measures perihelion angle at each orbit,

• Tunes the extra perturbation factor,

• Plots the trajectory,

• Prints the measured precession per century.

👉 I'll show the tuned script using a first manual guess (roughly $1.3\times10^{-7}$ was close in Stefan’s derivations). Then we can refine.

⏸️ The simulation ran too slowly — looping through ~315 million iterations (100 years at 10 s steps) is too heavy for the environment.

👉 To fix this, we can speed up the simulation by:

• Increasing the timestep (e.g. from 10 s to 1000 s),

• Or simulating fewer orbits (e.g. 500–1000 revolutions) and scaling the precession accordingly.

This still gives a good approximation of the perihelion advance without brute-forcing 100 years.

Do you want me to:

1. 🔸 Increase timestep (faster but slightly less accurate),

2. 🔸 Simulate fewer orbits (e.g. 500 orbits) and scale result to 1 century,

3. 🔸 Both (for maximum speed).

Which option shall I use?

✅ Good choice — simulating a limited number of orbits (e.g. 500) and then scaling the measured perihelion shift to 1 century is the most efficient way to match Mercury’s 532″/century precession.

Here’s what I’ll do in the updated script:

• Simulate 500 Mercury orbits (each ≈ 88 Earth days ≈ 7.6×10⁶ s).

• Detect perihelion angle each orbit.

• Measure total precession after 500 orbits.

• Scale the result to 1 century (Mercury completes ~415 orbits/century).

• Adjust the perturbation term (≈1.3e-7 as starting point) to get close to 532″.

Let's run this optimized version 👇

Here’s the result 👆

With perturbation = 1.3 × 10⁻⁷, the model predicts a precession of roughly:

which is far too large compared to Mercury’s observed

👉 This tells us the perturbation must be much smaller — on the order of $\sim 10^{-10}$ rather than $10^{-7}$.

Would you like me to automatically tune the perturbation parameter (via binary search) to match ≈ 532″/century?

Good — I tried to automatically tune the model to reach ≈ 532 arcseconds/century but hit numerical/instability problems in the environment while searching. I’ll be transparent about what I did, what went wrong, and — most importantly — give you a robust, ready-to-run script you can execute locally (or I can keep running here if you want me to try again): the script will safely search for the parameter that controls the precession and print a nicely annotated result + plot.

What I tried (short)

• I started from your original integrator and attempted to tune the small perturbation value (and then the initial angular-acceleration ap) to match 532″/century.

• The integrator is very sensitive to the angular-acceleration terms (the ap variable and the small additive perturbation used inside rpp). For many parameter choices the radius r grows or becomes numerically unstable and the simulation overflows.

• Because of that instability I could not reliably converge to a tuned parameter in this environment without very careful step control and defensive code (and more compute/time than a single quick tuning run).

• However: I implemented diagnostic simulations and tests which show the model is capable of huge perihelion shifts — so the mechanism is working, but the calibration is delicate.

Recommendation / plan (and why)

• The right parameter to tune is the initial angular acceleration ap (and/or the tiny additive term added in rpp).

• Best approach is a safe calibration script that:

  1. uses a moderate timestep (so the integrator stays stable),

  2. stores r and angle histories and finds perihelia by searching local minima in the r-series (more robust than relying on instantaneous rp sign),

  3. runs a controlled dichotomy / grid scan in log-space across ap values until the scaled precession matches the target,

  4. protects against overflow or runaway orbits and aborts gracefully for unstable parameters.

• I prepared such a script below. It will:

  • simulate a limited number of orbits (so it finishes quickly),

  • scale the measured precession to 1 century (Mercury ≈ 415 orbits/century),

  • perform a coarse scan and then refine (binary search in log-space),

  • plot the final orbit (subsampled) and print the final ap and precession.

You can run it locally (recommended) or I can run it here again (but note it may require repeated runs to refine and avoid numerical overflows).

Robust calibration script (copy & run locally)

python
Code kopiëren
# Calibrate_ap_for_Mercury_precession.py
# Run locally with python3. Requires numpy & matplotlib.

import numpy as np
import matplotlib.pyplot as plt

# --- Simulation function (safe, detects local minima in r) ---
def simulate(ap_init, num_orbits=80, dt=3000, subsample_plot=500):
    Cstefan = -1.32746089776e+20
    r = 4.28169e10   # initial radius (m)
    a = 0.0
    rp = 0.0
    ap = ap_init

    # initial derived quantities
    rpp = Cstefan / (r**2) + (r * (ap)**2)
    app = (-2.0 * rp * ap) / r
    rppp = 0.0
    appp = 0.0

    T_orbit = 7.6e6
    steps = int(num_orbits * T_orbit / dt)
    if steps < 1000:
        steps = max(1000, steps)

    # storage
    r_hist = np.empty(steps)
    a_hist = np.empty(steps)
    x_pts = []
    y_pts = []
    rMax = 1.2 * r

    try:
        for i in range(steps):
            # compute higher derivatives
            R_term6 = (Cstefan * -2.0 * rp) / (r ** 3)
            rppp = rp * ap**2 + 2.0 * r * ap * app + R_term6
            appp = (-3.0 * rp * app - 2.0 * rpp * ap) / r

            # integrate (Taylor-like expansion, same form as your script)
            dr = rp*dt + 0.5*rpp*(dt**2) + (rppp*(dt**3))/6.0
            da = ap*dt + 0.5*app*(dt**2) + (appp*(dt**3))/6.0

            r += dr
            rp += rpp*dt + 0.5*rppp*(dt**2)
            rpp += rppp * dt
            a += da
            ap += app*dt + 0.5*appp*(dt**2)
            app += appp * dt

            # store
            r_hist[i] = r
            a_hist[i] = a % (2*np.pi)

            # subsampled points for plotting
            if (i % subsample_plot) == 0:
                x = r*np.cos(a)
                y = r*np.sin(a)
                x_rel = x / rMax
                y_rel = y / rMax
                xpos = (0.5*x_rel) + 0.5
                ypos = (0.5*y_rel) + 0.5
                x_pts.append(600 * xpos)
                y_pts.append(600 * ypos)

            # safety checks to avoid runaway values
            if not np.isfinite(r) or r <= 0 or r > 1e20:
                return None, None, None
    except Exception:
        return None, None, None

    # find local minima (perihelia)
    peri_angles = []
    for i in range(1, steps-1):
        if r_hist[i] < r_hist[i-1] and r_hist[i] < r_hist[i+1]:
            peri_angles.append(a_hist[i])

    if len(peri_angles) < 3:
        return None, x_pts, y_pts

    # total precession during simulated orbits (radians -> arcsec)
    total_precession_rad = (peri_angles[-1] - peri_angles[0]) % (2*np.pi)
    total_precession_arcsec = total_precession_rad * 206265.0

    # scale per orbit then to 1 century (Mercury ≈ 415 orbits/century)
    per_orbit_arcsec = total_precession_arcsec / num_orbits
    precession_century = per_orbit_arcsec * 415.0

    return precession_century, x_pts, y_pts

# --- Calibration driver: coarse scan + refinement in log-space ---
def calibrate(target_arcsec=532.0):
    # search range for ap (initial angular acceleration)
    low = 1e-14
    high = 1e-5

    # coarse logarithmic scan to find bracket
    coarse = np.logspace(np.log10(low), np.log10(high), num=12)
    bracket = None
    last_val = None
    last_ap = None
    for ap in coarse:
        val, _, _ = simulate(ap)
        if val is None:
            last_val = None
            last_ap = ap
            continue
        # print progress
        print(f"scan ap={ap:.3e} -> precession={val:.2f}\"/century")
        if last_val is not None:
            # check if target lies between last_val and val
            if (last_val - target_arcsec) * (val - target_arcsec) < 0:
                bracket = (last_ap, ap)
                break
        last_val = val
        last_ap = ap

    if bracket is None:
        # pick best approximate ap (closest)
        # compute differences for all evaluated coarse points
        diffs = []
        for ap in coarse:
            val, _, _ = simulate(ap)
            if val is None:
                diffs.append((ap, np.inf))
            else:
                diffs.append((ap, abs(val - target_arcsec)))
        diffs.sort(key=lambda x: x[1])
        best_ap = diffs[0][0]
        print(f"No sign-change found; best coarse ap ≈ {best_ap:.3e}")
        return best_ap

    # refine bracket with binary search in log-space (multiplicative)
    a_low, a_high = bracket
    print(f"Refining between {a_low:.3e} and {a_high:.3e}")

    for _ in range(18):  # up to 18 refinement steps
        mid = np.sqrt(a_low * a_high)  # geometric mean
        val, _, _ = simulate(mid)
        if val is None:
            # mid is unstable: shrink interval towards the stable side
            a_high = mid
            continue
        print(f"test ap={mid:.3e} -> {val:.3f}\"/century")
        if val > target_arcsec:
            # mid causes too much precession -> need smaller ap
            a_high = mid
        else:
            a_low = mid

    chosen = np.sqrt(a_low * a_high)
    final_val, x_pts, y_pts = simulate(chosen)
    print(f"\nChosen ap = {chosen:.6e}, Precession ≈ {final_val:.3f}\"/century")
    return chosen, final_val, x_pts, y_pts

if __name__ == "__main__":
    result = calibrate(532.0)
    if isinstance(result, tuple):
        ap_chosen, final_val, x_pts, y_pts = result
        print("Final: ap =", ap_chosen, "precession =", final_val)
        plt.figure(figsize=(6,6))
        plt.plot(x_pts, y_pts, label='trajectory (subsampled)')
        plt.axhline(300, linewidth=0.5)
        plt.axvline(300, linewidth=0.5)
        plt.xlim(0,600)
        plt.ylim(0,600)
        plt.title(f'ap={ap_chosen:.3e} -> {final_val:.2f}\"/century')
        plt.legend()
        plt.show()
    else:
        print("Could not refine automatically; returned:", result)

How to use

• Save as Calibrate_ap_for_Mercury_precession.py.

• Run python3 Calibrate_ap_for_Mercury_precession.py (needs numpy and matplotlib).

• The script will print progress from the coarse scan, show refinement steps, and then display the final subsampled trajectory.

Why run locally?

• The calibration may need many simulation evaluations (and some will be unstable). Running on your local machine is faster and safer than repeating runs here.

• If you want, I can continue trying here — but I will need to guard tightly against overflow and may need multiple attempts. The code above is crafted for stability and should produce a tunable result locally.

Quick summary of what I learned while trying

• Your integrator and the added angular-acceleration produce very large precession rates when ap (or the additive perturbation) is even moderately large.

• The correct parameter to tune appears to be the initial ap (or the small term inside rpp), not only the tiny additive 0.0000001 you had used earlier.

• Numerical stability is the main issue; local minima detection in the r trace is more robust than checking rp sign.

If you want me to:

• (A) Run the calibration here with the exact parameters in the script above (I will attempt it again now, using the same safe approach), or

• (B) Leave the script for you to run locally, and then you paste the result here and I’ll analyze/refine it further,

say which you prefer. I’m ready to continue either way.

I tried to run an automatic calibration here but ran into numerical instability and execution limits in this environment — the integrator you’re using is very sensitive and many parameter choices lead to runaway/overflow before we can measure perihelion locations reliably. I attempted several safer scans (smaller scans, larger timesteps, fewer orbits) but could not find a robust tuned value within the runtime constraints.

I don’t want to leave you hanging — here are two concrete, immediate options (pick one):

1. I will try again here now using a far more stable integrator (velocity-Verlet / symplectic) and a defensive calibration routine. That integrator conserves orbital structure much better and should avoid the runaway behaviour; I’ll then attempt the log-space search for the angular-acceleration parameter and report/plot the result. (This may still require a few tries but is the best next step for running here.)

2. I give you a robust, ready-to-run script (velocity-Verlet integrator + automatic calibration and plotting) you can run locally (much faster and safer). I’ll also explain the parameters to tweak and how to validate the resulting precession. This is usually the fastest reliable path.

Which do you want me to do: (A) try again here now with the Verlet integrator, or (B) send the robust script for you to run locally?