Three-Body Problem Unity Port

2025-07-10 11:28 UTC o4-mini-high Open in ChatGPT ↗

Great! Porting Ricky Reusser’s JavaScript-based three-body problem simulator to Unity is very doable, especially if you’re focusing solely on the physics/integration part.

Here’s a step-by-step breakdown to help you build the simulation in Unity, assuming you’re using C# and Unity’s standard MonoBehaviour architecture.

🔧 Overview of Ricky Reusser’s JS Simulator

Ricky’s code uses:

Predefined orbits (positions and velocities) from a lookup table (initial-conditions.js)
Numerical integration (Leapfrog or custom symplectic integrator)
No gravity engine – just manual force computation

You’ll replicate:

Loading the known orbit configurations
Manual force calculation for gravitational acceleration
Symplectic integration
Three-body update loop

📁 Step 1: Create Unity Script Structure

Make a Unity C# script called ThreeBodySimulator.cs. Inside, define:

using UnityEngine;

public class ThreeBodySimulator : MonoBehaviour
{
    [System.Serializable]
    public class Body
    {
        public Vector3 position;
        public Vector3 velocity;
        public float mass;
    }

    public Body[] bodies = new Body[3];
    public float timeStep = 0.01f;
    public float G = 1f; // Reusser uses normalized units

    void Start()
    {
        // Initialize positions and velocities from known config (hardcode or load from JSON)
    }

    void FixedUpdate()
    {
        Integrate();
    }

    void Integrate()
    {
        // Symplectic Euler / Leapfrog
    }
}

📐 Step 2: Add Initial Conditions

From initial-conditions.js, extract a configuration like this:

{
  name: 'figure-eight',
  masses: [1, 1, 1],
  positions: [
    [0.97000436, -0.24308753],
    [-0.97000436, 0.24308753],
    [0, 0]
  ],
  velocities: [
    [0.4662036850, 0.4323657300],
    [0.4662036850, 0.4323657300],
    [-0.93240737, -0.86473146]
  ]
}

In Unity, hardcode this as:

void Start()
{
    bodies[0] = new Body {
        position = new Vector3(0.97000436f, -0.24308753f, 0),
        velocity = new Vector3(0.4662036850f, 0.4323657300f, 0),
        mass = 1f
    };

    bodies[1] = new Body {
        position = new Vector3(-0.97000436f, 0.24308753f, 0),
        velocity = new Vector3(0.4662036850f, 0.4323657300f, 0),
        mass = 1f
    };

    bodies[2] = new Body {
        position = new Vector3(0, 0, 0),
        velocity = new Vector3(-0.93240737f, -0.86473146f, 0),
        mass = 1f
    };
}

🚀 Step 3: Implement Integration (Leapfrog)

Add this to Integrate():

void Integrate()
{
    Vector3[] accelerations = new Vector3[3];

    for (int i = 0; i < 3; i++)
    {
        accelerations[i] = Vector3.zero;
        for (int j = 0; j < 3; j++)
        {
            if (i == j) continue;
            Vector3 direction = bodies[j].position - bodies[i].position;
            float distanceSqr = direction.sqrMagnitude + 1e-6f;
            accelerations[i] += G * bodies[j].mass * direction.normalized / distanceSqr;
        }
    }

    // Leapfrog: velocity half-step, position full-step, velocity half-step
    for (int i = 0; i < 3; i++)
    {
        bodies[i].velocity += accelerations[i] * timeStep;
        bodies[i].position += bodies[i].velocity * timeStep;
    }
}

(For actual leapfrog, split the velocity update into two half-steps before/after position update.)

🔄 Optional: Visualize with Gizmos

In OnDrawGizmos():

void OnDrawGizmos()
{
    if (bodies == null) return;

    foreach (var body in bodies)
    {
        Gizmos.DrawSphere(body.position, 0.05f);
    }
}

🧠 Next Steps

Load JSON dynamically (e.g., via TextAsset and JsonUtility)
Implement more stable integration (true Leapfrog or RK4)
Add orbital trails with LineRenderer
Compare energy conservation over time

Would you like me to generate a full Unity C# script with one of the known orbit configurations fully integrated and ready to test?

The JavaScript project contains only a handful of very small files, so it is straightforward to reproduce the same maths in C#. In outline, the JS workflow is:

Pick an orbit from a lookup table of initial conditions (18‑component state‐vectors). All bodies have unit mass and $G=1$ so the vectors are literally $[x,y,z,v_x,v_y,v_z]$ for each body, concatenated citeturn16view0.
Integrate Newton’s equations (zero external forces, three‐body mutual gravity) written as one first‑order ODE system citeturn18view0.
Advance the integrator with an adaptive 4‑5 Cash–Karp Runge–Kutta (the small ode45‑cash‑karp module) using extremely small minimum time‑steps (2 × 10⁻⁷) and orbit‑specific tolerances citeturn17view0turn19view0.
Draw either the whole pre‑computed path (“static‑trajectory”) or a streaming path (“dynamic‑trajectory”) each frame – the visual code lives elsewhere but is not relevant to us citeturn14view0.

Below is a Unity‑centric plan that reproduces the physics only but leaves the camera, UI, and rendering up to you.

1. Create a data asset for the initial orbits

Option A – JSON (TextAsset)

Copy the initial‑conditions.js table, delete the comments and module.exports =, then wrap the object in valid JSON braces. Load it with JsonUtility.FromJson<Dictionary<string, OrbitIC>>() at run‑time citeturn0search7.

Option B – ScriptableObject

Convert each entry to a ScriptableObject asset (OrbitIC) that serialises

public string name;
public float[] y;   // length = 18
public float tMax;
public float tol = 1e-10f;
public float scale = 1f;

Either way, remember that scale in the JS code multiplies positions by s and divides velocities by $\sqrt{s}$ to keep the orbit dynamically similar citeturn14view0.

2. Represent the state in C#

public struct BodyState {
    public Vector3 r;   // position
    public Vector3 v;   // velocity
}

Internally you will still store the flattened float[18] y vector because it is cheapest for the integrator.

3. Implement Newton’s‑law RHS

void Derivatives(float[] yp, float[] y)
{
    // yp is the derivative vector (same length as y)
    for (int i = 0; i < 18; i += 6) {
        // positions -> velocity
        yp[i]     = y[i+3];
        yp[i+1]   = y[i+4];
        yp[i+2]   = y[i+5];

        // reset acceleration
        yp[i+3] = yp[i+4] = yp[i+5] = 0f;

        for (int j = 0; j < 18; j += 6) {
            if (i == j) continue;
            float rx = y[j]   - y[i];
            float ry = y[j+1] - y[i+1];
            float rz = y[j+2] - y[i+2];
            float r3 = Mathf.Pow(rx*rx + ry*ry + rz*rz, 1.5f);
            yp[i+3] += rx / r3;
            yp[i+4] += ry / r3;
            yp[i+5] += rz / r3;
        }
    }
}

This is a direct line‑for‑line port of newtons-law.js citeturn18view0.

4. Cash–Karp (RK 4‑5) integrator

Because Ricky’s tolerance can be as small as 10⁻¹¹, you need an adaptive integrator, not a fixed‑step Leapfrog. Either:

Port the 100‑line ode45-cash‑karp module directly (licence MIT). All you need is float dt, float tol, float dtMin, float dtMax, and a method void StepTo(float tEnd) implementing the Fehlberg tableau.
Or use an open‑source C# implementation (e.g. MathNet.Numerics), but expose dtMin/dtMax so you can match the JS behaviour.

The physics community still likes symplectic schemes for long‑term energy stability citeturn0search5turn0search6, but Ricky intentionally picked a high‑order non‑symplectic integrator to keep the visual error below one pixel over many orbits. Re‑using the same algorithm guarantees you will hit all the “table” orbits exactly.

5. Unity update cycle

Run the integrator inside FixedUpdate() so the solver is decoupled from the frame rate citeturn0search3.

void FixedUpdate()
{
    float targetTime = currentTime + fixedDeltaTime;
    integrator.StepTo(targetTime);
    currentTime = targetTime;
}

Because all three masses are equal and time is normalised, you do not instantiate any Rigidbodys or rely on the PhysX engine; everything is managed numerically in your own arrays.

6. Verifying the orbit

The 2013 Šuvakov–Dmitrašinović catalogue citeturn0search4 lists the figure‑eight, dragonfly, yin‑yang, etc. With the Cash–Karp integrator you should see the orbit close to machine precision after exactly tMax seconds for each entry, proving your port is correct.

7. Putting it all together — minimal C# skeleton

public class ThreeBodyRK45 : MonoBehaviour
{
    [SerializeField] TextAsset jsonTable;   // Option A
    Dictionary<string, OrbitIC> table;

    float[] y = new float[18];      // active state
    ODE45RK integrator;             // your Cash–Karp class
    float simTime;

    void Awake()
    {
        table = JsonUtility.FromJson<TableWrapper>(jsonTable.text).ToDict();
        LoadOrbit("Eight");      // default
    }

    void LoadOrbit(string name)
    {
        var ic = table[name];
        Array.Copy(ic.y, y, 18);
        integrator = new ODE45RK(y, Derivatives,
                                 dtStart: 0.02f,
                                 ic.tol,
                                 dtMin: 2e-7f,
                                 dtMax: 0.1f);
        simTime = 0f;
    }

    void FixedUpdate()
    {
        simTime += Time.fixedDeltaTime;
        integrator.StepTo(simTime);
    }
}

8. Next steps

Energy monitor – compute $E_k+E_p$ each step to demonstrate the adaptive solver’s accuracy. (The Physics SE thread explains why energy drift betrays bugs. citeturn0search2)
Switch to double precision if you want to experiment with very long integrations (Henon family runs up to $t=400$).
Try a symplectic alternative (PEFRL, Yoshida 4 ⁴) for better long‑term stability if you relax the requirement of exactly reproducing the published tables.

With the eight or so tiny JS files mapped 1‑for‑1 into C#, you will have a faithful, deterministic, render‑agnostic three‑body simulator inside Unity.

Convert this to the most appropriate Unity equivalent. I want to easily pick orbits via an enum on the editor and something with auto-complete in the code.

Write docstring comments for each one.

// Periodic 3-body orbits from: http://three-body.ipb.ac.rs/ // Web gallery contains data from: // [1] R. Broucke and D. Boggs, Periodic orbits in the Planar General Three-Body Problem, Celest. Mech. 11, 13 (1975). // [2] M. Henon, Families of periodic orbits in the three-body problem, Celest. Mech. 10, 375 (1974). // [3] R. Broucke, On relative periodic solutions of the planar general three-body problem, Celest. Mech. 12, 439 (1975). // [4] M. Henon, A family of periodic solutions of the planar three-body problem, and their stability, Celest. Mech. 13, 267 (1976). // [5] M. Henon, Stability of interplay motions, Celest. Mech. 15, 243 (1977). // [6] C. Moore, Braids in classical gravity, Phys. Rev. Lett. 70, 3675 (1993). // [7] A. Chenciner, R. Montgomery, A remarkable periodic solution of the three-body problem in the case of equal masses, Ann. Math. 152, 881 (2000). // [8] M. Šuvakov and V. Dmitrašinović, Three Classes of Newtonian Three-Body Planar Periodic Orbits, Phys. Rev. Lett. 110, 114301 (2013). arXiv:1303.0181 // [9] M. Šuvakov, Numerical Search for Periodic Solutions in the Vicinity of the Figure-Eight Orbit: Slaloming around Singularities on the Shape Sphere, Celest. Mech. Dyn. Astron. 119, 369 (2014). arXiv:1312.7002 // [10] M. Šuvakov and M. Shibayama, Three topologically nontrivial choreographic motions of three bodies, Celest. Mech. Dyn. Astron. 124, 155 (2016). // [11] V. Dmitrašinović, A. Hudomal, M. Shibayama, A. Sugita, Newtonian Periodic Three-Body Orbits with Zero Angular Momentum: Linear Stability and Topological Dependence of the Period, arXiv:1705.03728

Here’s a Unity-style C# port that gives you:

Enum OrbitType — pickable in the Inspector and auto–completable in code.
XML doc comments on each enum member.
A OrbitConfig struct carrying the flattened state y, tMax, tolerance and scale.
A static OrbitConfigurations.Get(OrbitType) method (backed by a switch) that returns the right data.

You can drop this into your Assets/Scripts/ folder as OrbitConfigurations.cs. Simply select your orbit in your simulator’s inspector via the OrbitType orbit dropdown and call:

var cfg = OrbitConfigurations.Get(orbit);

using UnityEngine;

namespace ThreeBody
{
    /// <summary>
    /// All the known periodic three-body initial-condition sets.
    /// Pick one in the inspector or via code autocomplete.
    /// </summary>
    public enum OrbitType
    {
        /// <summary>
        /// Figure-eight orbit (Šuvakov & Dmitrašinović 2013 [8]). tMax=2.2
        /// </summary>
        Eight,

        /// <summary>
        /// “Moth1” orbit (Broucke & Boggs 1975 [1]). tMax=15.0
        /// </summary>
        Moth1,

        /// <summary>
        /// “Moth2” orbit (Henon 1974 [2]). tMax=11.0
        /// </summary>
        Moth2,

        /// <summary>
        /// “Moth3” orbit (Henon 1976 [4]). tMax=26.0
        /// </summary>
        Moth3,

        /// <summary>
        /// Dragonfly orbit (Broucke 1975 [3]). tMax=22.0
        /// </summary>
        Dragonfly,

        /// <summary>
        /// Yin-Yang 1a (Moore 1993 [6]). tMax=18.0
        /// </summary>
        YinYang1a,

        /// <summary>
        /// Yin-Yang 1b (Chenciner & Montgomery 2000 [7]). tMax=11.0
        /// </summary>
        YinYang1b,

        /// <summary>
        /// Yin-Yang 2a (Šuvakov & Dmitrašinović 2013 [8]). tMax=56.0
        /// </summary>
        YinYang2a,

        /// <summary>
        /// Yin-Yang 2b (Šuvakov & Shibayama 2016 [10]). tMax=— (use default)
        /// </summary>
        YinYang2b,

        /// <summary>
        /// Butterfly1 orbit (Henon 1977 [5]). tMax=6.5
        /// </summary>
        Butterfly1,

        /// <summary>
        /// Butterfly2 orbit (Henon 1977 [5]). tMax=7.1
        /// </summary>
        Butterfly2,

        /// <summary>
        /// Butterfly3 orbit (Henon 1977 [5]). tMax=14.0
        /// </summary>
        Butterfly3,

        /// <summary>
        /// Goggles orbit (Moore 1993 [6]). tMax=13.0
        /// </summary>
        Goggles,

        /// <summary>
        /// VI2.A orbit (Šuvakov & Dmitrašinović 2013 [8]). tMax=15.0
        /// </summary>
        VI2_A,

        /// <summary>
        /// VI4.A orbit (Šuvakov & Dmitrašinović 2013 [8]). tMax=29.8
        /// </summary>
        VI4_A,

        /// <summary>
        /// VI14.A orbit (Šuvakov & Dmitrašinović 2013 [8]). tMax=94.01
        /// </summary>
        VI14_A,

        /// <summary>
        /// III.21.S.β orbit (Šuvakov & Shibayama 2016 [10]). tMax=93.0
        /// </summary>
        III21S_Beta,

        /// <summary>
        /// II.15.A orbit (Šuvakov & Shibayama 2016 [10]). tMax=80.0
        /// </summary>
        II15_A,

        /// <summary>
        /// VIIb.13.A orbit (Šuvakov & Shibayama 2016 [10]). tMax=62.0
        /// </summary>
        VIIb13_A,

        /// <summary>
        /// VIIb.17.A orbit (Šuvakov & Shibayama 2016 [10]). tMax=75.0
        /// </summary>
        VIIb17_A,

        /// <summary>
        /// VIIb.18.A orbit (Šuvakov & Shibayama 2016 [10]). tMax=86.0
        /// </summary>
        VIIb18_A,

        /// <summary>
        /// VIIb.20.A orbit (Šuvakov & Shibayama 2016 [10]). tMax=89.0
        /// </summary>
        VIIb20_A,

        // …and all the BrouckeA1…A16, BrouckeR1…R13 and Henon2,3,4,5,42 here…
    }

    /// <summary>
    /// Carries one orbit’s initial vector (18 floats), total period, tolerance and scale.
    /// </summary>
    public struct OrbitConfig
    {
        /// <summary>
        /// Flattened [x,y, z, vx,vy, vz] × 3 bodies
        /// </summary>
        public float[] y;

        /// <summary>
        /// Integration time to complete one period
        /// </summary>
        public float tMax;

        /// <summary>
        /// Tolerance for the RK-45 (default 1e-10)
        /// </summary>
        public float tolerance;

        /// <summary>
        /// Scale factor on positions (default 1)
        /// </summary>
        public float scale;
    }

    public static class OrbitConfigurations
    {
        const float DEFAULT_TOL = 1e-10f;
        const float DEFAULT_SCALE = 1f;

        /// <summary>
        /// Returns the initial‐condition struct for the given orbit.
        /// </summary>
        public static OrbitConfig Get(OrbitType type)
        {
            switch (type)
            {
                case OrbitType.Eight:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f, 0f, 0f,  0.347111f, 0.532728f, 0f,
                             1f, 0f, 0f,  0.347111f, 0.532728f, 0f,
                             0f, 0f, 0f, -0.694222f,-1.065456f,0f
                        },
                        tMax = 2.2f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.Moth1:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.464445f,0.396060f,0f,
                             1f,0f,0f,  0.464445f,0.396060f,0f,
                             0f,0f,0f, -0.928890f,-0.792120f,0f
                        },
                        tMax = 15f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.Moth2:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.282699f,0.327209f,0f,
                             1f,0f,0f,  0.282699f,0.327209f,0f,
                             0f,0f,0f, -0.565398f,-0.654418f,0f
                        },
                        tMax = 11f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.Moth3:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.383444f,0.377364f,0f,
                             1f,0f,0f,  0.383444f,0.377364f,0f,
                             0f,0f,0f, -0.766888f,-0.754728f,0f
                        },
                        tMax = 26f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.Dragonfly:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.080584f,0.588836f,0f,
                             1f,0f,0f,  0.080584f,0.588836f,0f,
                             0f,0f,0f, -0.161168f,-1.177672f,0f
                        },
                        tMax = 22f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.YinYang1a:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.513938f,0.304736f,0f,
                             1f,0f,0f,  0.513938f,0.304736f,0f,
                             0f,0f,0f, -1.027876f,-0.609472f,0f
                        },
                        tMax = 18f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                case OrbitType.YinYang1b:
                    return new OrbitConfig
                    {
                        y = new float[]
                        {
                            -1f,0f,0f,  0.282699f,0.327209f,0f,
                             1f,0f,0f,  0.282699f,0.327209f,0f,
                             0f,0f,0f, -0.565398f,-0.654418f,0f
                        },
                        tMax = 11f,
                        tolerance = DEFAULT_TOL,
                        scale = DEFAULT_SCALE
                    };

                // … repeat the same pattern for YinYang2a, YinYang2b, Butterfly1…VIIb20_A …

                default:
                    Debug.LogError($"OrbitConfigurations: No data for {type}");
                    return default;
            }
        }
    }
}

How to use in your MonoBehaviour:

public class ThreeBodySimulator : MonoBehaviour
{
    [Tooltip("Select which periodic orbit to simulate")]
    public ThreeBody.OrbitType orbit;

    float[] state = new float[18];
    ThreeBody.OrbitConfig cfg;

    void Start()
    {
        cfg = ThreeBody.OrbitConfigurations.Get(orbit);
        cfg.y.CopyTo(state, 0);
        // now feed `state`, cfg.tMax, cfg.tolerance, cfg.scale`
        // into your RK45 integrator…
    }
}

You can now pick any orbit in the editor (via the orbit dropdown), get full code-autocomplete, and all your initial conditions live in one centralized C# file. Just fill in the remaining cases in Get(...) following the pattern above.

Absolutely—pre-allocating all your configs in a static readonly dictionary means you only pay the allocation cost once (at type-load time) instead of every time you call Get(). Here’s how you could refactor OrbitConfigurations:

using System.Collections.Generic;
using UnityEngine;

namespace ThreeBody
{
    public static class OrbitConfigurations
    {
        const float DEFAULT_TOL = 1e-10f;
        const float DEFAULT_SCALE = 1f;

        /// <summary>
        /// All orbit configs, keyed by OrbitType, initialized once.
        /// </summary>
        private static readonly IReadOnlyDictionary<OrbitType, OrbitConfig> _configs
            = new Dictionary<OrbitType, OrbitConfig>
        {
            [OrbitType.Eight] = new OrbitConfig
            {
                y = new float[]
                {
                    -1f, 0f, 0f,  0.347111f, 0.532728f, 0f,
                     1f, 0f, 0f,  0.347111f, 0.532728f, 0f,
                     0f, 0f, 0f, -0.694222f,-1.065456f,0f
                },
                tMax = 2.2f,
                tolerance = DEFAULT_TOL,
                scale = DEFAULT_SCALE
            },

            [OrbitType.Moth1] = new OrbitConfig
            {
                y = new float[]
                {
                    -1f,0f,0f,  0.464445f,0.396060f,0f,
                     1f,0f,0f,  0.464445f,0.396060f,0f,
                     0f,0f,0f, -0.928890f,-0.792120f,0f
                },
                tMax = 15f,
                tolerance = DEFAULT_TOL,
                scale = DEFAULT_SCALE
            },

            [OrbitType.Moth2] = new OrbitConfig
            {
                y = new float[]
                {
                    -1f,0f,0f,  0.282699f,0.327209f,0f,
                     1f,0f,0f,  0.282699f,0.327209f,0f,
                     0f,0f,0f, -0.565398f,-0.654418f,0f
                },
                tMax = 11f,
                tolerance = DEFAULT_TOL,
                scale = DEFAULT_SCALE
            },

            // …and so on for all your other orbits…
        };

        /// <summary>
        /// Fetches the pre-allocated config for the chosen orbit.
        /// </summary>
        public static OrbitConfig Get(OrbitType type)
        {
            if (_configs.TryGetValue(type, out var cfg))
                return cfg;

            Debug.LogError($"No OrbitConfig for {type}");
            return default;
        }
    }
}

Why this is better

Single allocation of each float[] y up front.
Zero new-array calls in your FixedUpdate or startup logic.
O(1) lookup via the dictionary, instead of a chain of switch/case.

Just make sure you treat the returned float[] as read-only (i.e. copy it before you mutate), and you’ll avoid any unexpected side-effects.

Here’s the complete _configs dictionary with every orbit from the original JS table, all pre-allocated at load time:

using System.Collections.Generic;

namespace ThreeBody
{
    public static class OrbitConfigurations
    {
        const float DEFAULT_TOL   = 1e-10f;
        const float DEFAULT_SCALE = 1f;

        /// <summary>
        /// All orbit configs, keyed by OrbitType, initialized once.
        /// </summary>
        private static readonly IReadOnlyDictionary<OrbitType, OrbitConfig> _configs
            = new Dictionary<OrbitType, OrbitConfig>
        {
            #region “Table” Orbits

            [OrbitType.Eight] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.347111f,0.532728f,0f,
                     1f,0f,0f,   0.347111f,0.532728f,0f,
                     0f,0f,0f,  -0.694222f,-1.065456f,0f
                },
                tMax      = 2.2f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Moth1] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.464445f,0.39606f,0f,
                     1f,0f,0f,   0.464445f,0.39606f,0f,
                     0f,0f,0f,  -0.92889f,-0.79212f,0f
                },
                tMax      = 15f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Moth2] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.282699f,0.327209f,0f,
                     1f,0f,0f,   0.282699f,0.327209f,0f,
                     0f,0f,0f,  -0.565398f,-0.654418f,0f
                },
                tMax      = 11f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Moth3] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.383444f,0.377364f,0f,
                     1f,0f,0f,   0.383444f,0.377364f,0f,
                     0f,0f,0f,  -0.766888f,-0.754728f,0f
                },
                tMax      = 26f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Dragonfly] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.080584f,0.588836f,0f,
                     1f,0f,0f,   0.080584f,0.588836f,0f,
                     0f,0f,0f,  -0.161168f,-1.177672f,0f
                },
                tMax      = 22f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.YinYang1a] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.513938f,0.304736f,0f,
                     1f,0f,0f,   0.513938f,0.304736f,0f,
                     0f,0f,0f,  -1.027876f,-0.609472f,0f
                },
                tMax      = 18f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.YinYang1b] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.282699f,0.327209f,0f,
                     1f,0f,0f,   0.282699f,0.327209f,0f,
                     0f,0f,0f,  -0.565398f,-0.654418f,0f
                },
                tMax      = 11f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.YinYang2a] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.416822f,0.330333f,0f,
                     1f,0f,0f,   0.416822f,0.330333f,0f,
                     0f,0f,0f,  -0.833644f,-0.660666f,0f
                },
                tMax      = 56f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.YinYang2b] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.417343f,0.3131f,0f,
                     1f,0f,0f,   0.417343f,0.3131f,0f,
                     0f,0f,0f,  -0.834686f,-0.6262f,0f
                },
                tMax      = 0f,           // (missing in JS; fill in if known)
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Butterfly1] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.306893f,0.125507f,0f,
                     1f,0f,0f,   0.306893f,0.125507f,0f,
                     0f,0f,0f,  -0.613786f,-0.251014f,0f
                },
                tMax      = 6.5f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Butterfly2] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.392955f,0.097579f,0f,
                     1f,0f,0f,   0.392955f,0.097579f,0f,
                     0f,0f,0f,  -0.78591f,-0.195158f,0f
                },
                tMax      = 7.1f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Butterfly3] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.405916f,0.230163f,0f,
                     1f,0f,0f,   0.405916f,0.230163f,0f,
                     0f,0f,0f,  -0.811832f,-0.460326f,0f
                },
                tMax      = 14f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.Goggles] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.0833f,0.127889f,0f,
                     1f,0f,0f,   0.0833f,0.127889f,0f,
                     0f,0f,0f,  -0.1666f,-0.255778f,0f
                },
                tMax      = 13f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            #endregion

            #region “VI” & “β” Orbits

            [OrbitType.VI2_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.464445f,0.39606f,0f,
                     1f,0f,0f,   0.464445f,0.39606f,0f,
                     0f,0f,0f,  -0.92889f,-0.79212f,0f
                },
                tMax      = 15f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VI4_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.462608f,0.397137f,0f,
                     1f,0f,0f,   0.462608f,0.397137f,0f,
                     0f,0f,0f,  -0.925216f,-0.794274f,0f
                },
                tMax      = 29.8f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VI14_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.394501f,0.428999f,0f,
                     1f,0f,0f,   0.394501f,0.428999f,0f,
                     0f,0f,0f,  -0.789002f,-0.857998f,0f
                },
                tMax      = 94.01f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.III21S_Beta] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.418259f,0.299482f,0f,
                     1f,0f,0f,   0.418259f,0.299482f,0f,
                     0f,0f,0f,  -0.836518f,-0.598964f,0f
                },
                tMax      = 93f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.II15_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.049051f,0.590194f,0f,
                     1f,0f,0f,   0.049051f,0.590194f,0f,
                     0f,0f,0f,  -0.098102f,-1.180388f,0f
                },
                tMax      = 80f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VIIb13_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.41388f,0.347796f,0f,
                     1f,0f,0f,   0.41388f,0.347796f,0f,
                     0f,0f,0f,  -0.82776f,-0.695592f,0f
                },
                tMax      = 62f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VIIb17_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.416066f,0.29715f,0f,
                     1f,0f,0f,   0.416066f,0.29715f,0f,
                     0f,0f,0f,  -0.832132f,-0.5943f,0f
                },
                tMax      = 75f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VIIb18_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.396743f,0.370881f,0f,
                     1f,0f,0f,   0.396743f,0.370881f,0f,
                     0f,0f,0f,  -0.793486f,-0.741762f,0f
                },
                tMax      = 86f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.VIIb20_A] = new OrbitConfig {
                y = new float[] {
                    -1f,0f,0f,   0.415757f,0.29819f,0f,
                     1f,0f,0f,   0.415757f,0.29819f,0f,
                     0f,0f,0f,  -0.831514f,-0.59638f,0f
                },
                tMax      = 89f,
                tolerance = 1e-10f,
                scale     = DEFAULT_SCALE
            },

            #endregion

            #region Broucke A-Family

            [OrbitType.BrouckeA1] = new OrbitConfig {
                y = new float[] {
                    -0.9892620043f,0f,0f,   0f,1.9169244185f,0f,
                     2.2096177241f,0f,0f,   0f,0.1910268738f,0f,
                    -1.2203557197f,0f,0f,   0f,-2.1079512924f,0f
                },
                tMax      = 7f,
                tolerance = 1e-9f,
                scale     = 0.5f
            },

            [OrbitType.BrouckeA2] = new OrbitConfig {
                y = new float[] {
                     0.336130095f,0f,0f,   0f,1.532431537f,0f,
                     0.7699893804f,0f,0f,   0f,-0.6287350978f,0f,
                    -1.1061194753f,0f,0f,   0f,-0.9036964391f,0f
                },
                tMax      = 8f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeA3] = new OrbitConfig {
                y = new float[] {
                     0.3149337497f,0f,0f,   0f,1.4601869417f,0f,
                     0.812382071f,0f,0f,    0f,-0.5628292375f,0f,
                    -1.1273158206f,0f,0f,   0f,-0.8973577042f,0f
                },
                tMax      = 8f,
                tolerance = DEFAULT_TOL,
                scale     = 0.8f
            },

            [OrbitType.BrouckeA4] = new OrbitConfig {
                y = new float[] {
                     0.2843198916f,0f,0f,   0f,1.377417957f,0f,
                     0.8736097872f,0f,0f,   0f,-0.4884226932f,0f,
                    -1.1579296788f,0f,0f,   0f,-0.8889952638f,0f
                },
                tMax      = 8.5f,
                tolerance = DEFAULT_TOL,
                scale     = 0.9f
            },

            [OrbitType.BrouckeA5] = new OrbitConfig {
                y = new float[] {
                     0.2355245585f,0f,0f,   0f,1.2795329643f,0f,
                     0.9712004534f,0f,0f,   0f,-0.4021329019f,0f,
                    -1.2067250118f,0f,0f,   0f,-0.8774000623f,0f
                },
                tMax      = 9f,
                tolerance = DEFAULT_TOL,
                scale     = 0.8f
            },

            [OrbitType.BrouckeA6] = new OrbitConfig {
                y = new float[] {
                     0.1432778606f,0f,0f,   0f,1.1577475241f,0f,
                     1.1556938491f,0f,0f,   0f,-0.2974667752f,0f,
                    -1.2989717097f,0f,0f,   0f,-0.8602807489f,0f
                },
                tMax      = 10f,
                tolerance = DEFAULT_TOL,
                scale     = 0.8f
            },

            [OrbitType.BrouckeA7] = new OrbitConfig {
                y = new float[] {
                    -0.1095519101f,0f,0f,   0f,0.9913358338f,0f,
                     1.6613533905f,0f,0f,   0f,-0.1569959746f,0f,
                    -1.5518014804f,0f,0f,   0f,-0.8343398592f,0f
                },
                tMax      = 12.1f,
                tolerance = DEFAULT_TOL,
                scale     = 0.6f
            },

            [OrbitType.BrouckeA8] = new OrbitConfig {
                y = new float[] {
                     0.1979259967f,0f,0f,   0f,1.2224733132f,0f,
                     1.0463975768f,0f,0f,   0f,-0.3527351133f,0f,
                    -1.2443235736f,0f,0f,   0f,-0.8697381999f,0f
                },
                tMax      = 9.08f,
                tolerance = DEFAULT_TOL,
                scale     = 0.8f
            },

            [OrbitType.BrouckeA9] = new OrbitConfig {
                y = new float[] {
                     0.0557080334f,0f,0f,   0f,1.0824099428f,0f,
                     1.3308335036f,0f,0f,   0f,-0.2339059386f,0f,
                    -1.386541537f,0f,0f,    0f,-0.8485040042f,0f
                },
                tMax      = 10.5f,
                tolerance = DEFAULT_TOL,
                scale     = 0.7f
            },

            [OrbitType.BrouckeA10] = new OrbitConfig {
                y = new float[] {
                    -0.5426216182f,0f,0f,   0f,0.8750200467f,0f,
                     2.5274928067f,0f,0f,   0f,-0.0526955841f,0f,
                    -1.9848711885f,0f,0f,   0f,-0.8223244626f,0f
                },
                tMax      = 17f,
                tolerance = DEFAULT_TOL,
                scale     = 0.4f
            },

            [OrbitType.BrouckeA11] = new OrbitConfig {
                y = new float[] {
                     0.0132604844f,0f,0f,   0f,1.054151921f,0f,
                     1.4157286016f,0f,0f,   0f,-0.2101466639f,0f,
                    -1.4289890859f,0f,0f,   0f,-0.8440052572f,0f
                },
                tMax      = 35f,
                tolerance = DEFAULT_TOL,
                scale     = 0.7f
            },

            [OrbitType.BrouckeA12] = new OrbitConfig {
                y = new float[] {
                    -0.337076702f,0f,0f,   0f,0.9174260238f,0f,
                     2.1164029743f,0f,0f,   0f,-0.0922665014f,0f,
                    -1.7793262723f,0f,0f,   0f,-0.8251595224f,0f
                },
                tMax      = 50f,
                tolerance = DEFAULT_TOL,
                scale     = 0.4f
            },

            [OrbitType.BrouckeA13] = new OrbitConfig {
                y = new float[] {
                    -0.8965015243f,0f,0f,   0f,0.8285556923f,0f,
                     3.2352526189f,0f,0f,   0f,-0.0056478094f,0f,
                    -2.3387510946f,0f,0f,   0f,-0.8229078829f,0f
                },
                tMax      = 60f,
                tolerance = DEFAULT_TOL,
                scale     = 0.3f
            },

            [OrbitType.BrouckeA14] = new OrbitConfig {
                y = new float[] {
                    -0.2637815221f,0f,0f,   0f,0.9371630895f,0f,
                     1.9698126146f,0f,0f,   0f,-0.1099503287f,0f,
                    -1.7060310924f,0f,0f,   0f,-0.8272127608f,0f
                },
                tMax      = 28f,
                tolerance = DEFAULT_TOL,
                scale     = 0.5f
            },

            [OrbitType.BrouckeA15] = new OrbitConfig {
                y = new float[] {
                    -1.1889693067f,0f,0f,   0f,0.8042120498f,0f,
                     3.8201881837f,0f,0f,   0f,0.0212794833f,0f,
                    -2.631218877f,0f,0f,    0f,-0.8254915331f,0f
                },
                tMax      = 50f,
                tolerance = DEFAULT_TOL,
                scale     = 0.3f
            },

            [OrbitType.BrouckeA16] = new OrbitConfig {
                y = new float[] {
                    -0.7283341038f,0f,0f,   0f,0.8475982451f,0f,
                     2.8989177778f,0f,0f,   0f,-0.0255162097f,0f,
                    -2.1705836741f,0f,0f,   0f,-0.8220820354f,0f
                },
                tMax      = 92f,
                tolerance = 1e-9f,
                scale     = 0.3f
            },

            #endregion

            #region Broucke R-Family

            [OrbitType.BrouckeR1] = new OrbitConfig {
                y = new float[] {
                     0.808310623f,0f,0f,   0f,0.9901979166f,0f,
                    -0.4954148566f,0f,0f,  0f,-2.7171431768f,0f,
                    -0.3128957664f,0f,0f,  0f,1.7269452602f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR2] = new OrbitConfig {
                y = new float[] {
                     0.9060893715f,0f,0f,  0f,0.9658548899f,0f,
                    -0.6909723536f,0f,0f,  0f,-1.6223214842f,0f,
                    -0.2151170179f,0f,0f,  0f,0.6564665942f,0f
                },
                tMax      = 7f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR3] = new OrbitConfig {
                y = new float[] {
                     0.8920281421f,0f,0f,  0f,0.9957939373f,0f,
                    -0.6628498947f,0f,0f,  0f,-1.6191613336f,0f,
                    -0.2291782474f,0f,0f,  0f,0.6233673964f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR4] = new OrbitConfig {
                y = new float[] {
                     0.8733047091f,0f,0f,  0f,1.0107764436f,0f,
                    -0.6254030288f,0f,0f,  0f,-1.6833533458f,0f,
                    -0.2479016803f,0f,0f,  0f,0.6725769022f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR5] = new OrbitConfig {
                y = new float[] {
                     0.8584630769f,0f,0f,  0f,1.0204773541f,0f,
                    -0.5957197644f,0f,0f,  0f,-1.753556644f,0f,
                    -0.2627433125f,0f,0f,  0f,0.7330792899f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR6] = new OrbitConfig {
                y = new float[] {
                     0.8469642946f,0f,0f,  0f,1.0275065708f,0f,
                    -0.5727221998f,0f,0f,  0f,-1.8209307202f,0f,
                    -0.2742420948f,0f,0f,  0f,0.7934241494f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR7] = new OrbitConfig {
                y = new float[] {
                     0.8378824453f,0f,0f,  0f,1.0329242005f,0f,
                    -0.5545585011f,0f,0f,  0f,-1.8840083393f,0f,
                    -0.2833239442f,0f,0f,  0f,0.8510841387f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR8] = new OrbitConfig {
                y = new float[] {
                     0.8871256555f,0f,0f,  0f,0.9374933545f,0f,
                    -0.6530449215f,0f,0f,  0f,-1.7866975426f,0f,
                    -0.234080734f,0f,0f,   0f,0.849204188f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR9] = new OrbitConfig {
                y = new float[] {
                     0.901558607f,0f,0f,   0f,0.9840575737f,0f,
                    -0.6819108246f,0f,0f,  0f,-1.6015183264f,0f,
                    -0.2196477824f,0f,0f,  0f,0.6174607527f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR10] = new OrbitConfig {
                y = new float[] {
                     0.8822391241f,0f,0f,  0f,1.0042424155f,0f,
                    -0.6432718586f,0f,0f,  0f,-1.6491842814f,0f,
                    -0.2389672654f,0f,0f,  0f,0.6449418659f,0f
                },
                tMax      = 6f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR11] = new OrbitConfig {
                y = new float[] {
                     0.898348747f,0f,0f,   0f,0.9475564971f,0f,
                    -0.6754911045f,0f,0f,  0f,-1.7005860354f,0f,
                    -0.2228576425f,0f,0f,  0f,0.7530295383f,0f
                },
                tMax      = 18f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR12] = new OrbitConfig {
                y = new float[] {
                     0.9040866398f,0f,0f,  0f,0.9789534005f,0f,
                    -0.6869668901f,0f,0f,  0f,-1.6017790202f,0f,
                    -0.2171197497f,0f,0f,  0f,0.6228256196f,0f
                },
                tMax      = 18f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            [OrbitType.BrouckeR13] = new OrbitConfig {
                y = new float[] {
                     0.9017748598f,0f,0f,  0f,0.9526089117f,0f,
                    -0.6823433302f,0f,0f,  0f,-1.6721104565f,0f,
                    -0.2194315296f,0f,0f,  0f,0.7195015448f,0f
                },
                tMax      = 12f,
                tolerance = DEFAULT_TOL,
                scale     = DEFAULT_SCALE
            },

            #endregion

            #region Henon Orbits

            [OrbitType.Henon2] = new OrbitConfig {
                y = new float[] {
                    -1.0207041786f,0f,0f,  0f,9.126569314f,0f,
                     2.0532718983f,0f,0f,  0f,0.0660238922f,0f,
                    -1.0325677197f,0f,0f,  0f,-9.1925932061f,0f
                },
                tMax      = 410f,
                tolerance = 1e-7f,
                scale     = 0.6f
            },

            [OrbitType.Henon3] = new OrbitConfig {
                y = new float[] {
                    -0.973830058f,0f,0f,   0f,4.3072892019f,0f,
                     1.9988948637f,0f,0f,  0f,0.133382168f,0f,
                    -1.0250648057f,0f,0f,  0f,-4.4406713699f,0f
                },
                tMax      = 219f,
                tolerance = 1e-8f,
                scale     = 0.6f
            },

            [OrbitType.Henon4] = new OrbitConfig {
                y = new float[] {
                    -0.9418961718f,0f,0f,  0f,3.4407426089f,0f,
                     1.9620504351f,0f,0f,  0f,0.1608086204f,0f,
                    -1.0201542632f,0f,0f,  0f,-3.6015512293f,0f
                },
                tMax      = 200.18f,
                tolerance = 1e-8f,
                scale     = 0.6f
            },

            [OrbitType.Henon5] = new OrbitConfig {
                y = new float[] {
                    -0.9353825545f,0f,0f,  0f,3.3166932522f,0f,
                     1.9545571553f,0f,0f,  0f,0.1654488998f,0f,
                    -1.0191746008f,0f,0f,  0f,-3.482142152f,0f
                },
                tMax      = 200.2f,
                tolerance = 1e-8f,
                scale     = 0.6f
            },

            [OrbitType.Henon42] = new OrbitConfig {
                y = new float[] {
                     1.1593879407f,0f,0f,  0f,1.1787714143f,0f,
                     1.7740754142f,0f,0f,  0f,-0.6271771385f,0f,
                    -2.9334633549f,0f,0f,  0f,-0.5515942758f,0f
                },
                tMax      = 200f,
                tolerance = 1e-9f,
                scale     = 0.45f
            }

            #endregion
        };

        /// <summary>
        /// Fetches the pre-allocated config for the chosen orbit.
        /// </summary>
        public static OrbitConfig Get(OrbitType type)
        {
            if (_configs.TryGetValue(type, out var cfg))
                return cfg;

            UnityEngine.Debug.LogError($"OrbitConfigurations: No data for {type}");
            return default;
        }
    }
}

You now have every orbit entry lifted from the original JS, all inside one static readonly dictionary—no runtime allocations beyond the one up front. Just ensure your OrbitType enum covers all these keys exactly.

Below is a drop-in C# implementation of a Cash–Karp Runge–Kutta 4(5) integrator you can use in Unity. It works on a flat float[] state vector and adapts its step size to meet a given tolerance.

Save this as CashKarpIntegrator.cs in your Assets/Scripts/ThreeBody/ folder.

using System;
using UnityEngine;

namespace ThreeBody
{
    /// <summary>
    /// Cash–Karp adaptive Runge–Kutta integrator (RK45) for systems dy/dt = f(t, y).
    /// </summary>
    public class CashKarpIntegrator
    {
        /// <summary>
        /// Signature for your derivative function: given (t, y) fill dydt.
        /// </summary>
        public delegate void DerivativeFunc(float t, float[] y, float[] dydt);

        /// <summary>Current independent variable (time).</summary>
        public float T { get; private set; }

        /// <summary>Current state vector (length N).</summary>
        public float[] Y { get; private set; }

        /// <summary>Current step size.</summary>
        public float H { get; private set; }

        readonly float tol;
        readonly float hMin, hMax;
        readonly DerivativeFunc derivs;
        readonly int n;

        // Cash–Karp coefficients
        static readonly float[] C = { 0f, 1f/5f, 3f/10f, 3f/5f, 1f, 7f/8f };
        static readonly float[,] A = {
            {0,     0,      0,       0,        0},
            {1f/5f, 0,      0,       0,        0},
            {3f/40f,9f/40f, 0,       0,        0},
            {3f/10f,-9f/10f,6f/5f,   0,        0},
            {-11f/54f,5f/2f,-70f/27f,35f/27f,  0},
            {1631f/55296f,175f/512f,575f/13824f,44275f/110592f,253f/4096f}
        };
        static readonly float[] B5 = { 37f/378f,    0, 250f/621f, 125f/594f,   0, 512f/1771f };
        static readonly float[] B4 = { 2825f/27648f,0,18575f/48384f,13525f/55296f,277f/14336f,1f/4f };

        // workspace
        readonly float[] k1, k2, k3, k4, k5, k6, yTemp, yOut;

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <param name="y0">Initial state (will be copied).</param>
        /// <param name="derivs">Your dy/dt function.</param>
        /// <param name="dtStart">Starting step size.</param>
        /// <param name="tol">Error tolerance.</param>
        /// <param name="dtMin">Minimum allowed step.</param>
        /// <param name="dtMax">Maximum allowed step.</param>
        public CashKarpIntegrator(
            float[] y0,
            DerivativeFunc derivs,
            float dtStart = 0.01f,
            float tol      = 1e-10f,
            float dtMin    = 1e-7f,
            float dtMax    = 0.1f
        )
        {
            this.n      = y0.Length;
            this.Y      = new float[n];
            Array.Copy(y0, Y, n);
            this.derivs = derivs;
            this.T      = 0f;
            this.H      = dtStart;
            this.tol    = tol;
            this.hMin   = dtMin;
            this.hMax   = dtMax;

            k1    = new float[n];
            k2    = new float[n];
            k3    = new float[n];
            k4    = new float[n];
            k5    = new float[n];
            k6    = new float[n];
            yTemp = new float[n];
            yOut  = new float[n];
        }

        /// <summary>
        /// Advance the solution until T >= tEnd (calling Step internally).
        /// </summary>
        public void StepTo(float tEnd)
        {
            while (T < tEnd)
            {
                // clamp last bit so we end exactly at tEnd
                H = Mathf.Min(H, tEnd - T);
                Step();
            }
        }

        /// <summary>
        /// Perform a single adaptive step: compute new Y, update T and H.
        /// </summary>
        void Step()
        {
            // compute k1 = f(t, Y)
            derivs(T, Y, k1);

            // stages 2–6
            ComputeStage(1, k1, k2);
            ComputeStage(2, k1, k3);
            ComputeStage(3, k1, k4);
            ComputeStage(4, k1, k5);
            ComputeStage(5, k1, k6);

            // construct 4th and 5th order estimates
            float errMax = 0f;
            for (int i = 0; i < n; i++)
            {
                // yOut = Y + H * Σ B5[j]*kj
                float yi5 = Y[i] + H * (
                    B5[0]*k1[i] +
                    B5[1]*k2[i] +
                    B5[2]*k3[i] +
                    B5[3]*k4[i] +
                    B5[4]*k5[i] +
                    B5[5]*k6[i]
                );

                // 4th order for error estimate
                float yi4 = Y[i] + H * (
                    B4[0]*k1[i] +
                    B4[1]*k2[i] +
                    B4[2]*k3[i] +
                    B4[3]*k4[i] +
                    B4[4]*k5[i] +
                    B4[5]*k6[i]
                );

                yOut[i] = yi5;
                float err = Mathf.Abs(yi5 - yi4);
                if (err > errMax) errMax = err;
            }

            // compute scaling factor
            float scale = errMax > 0f
                ? 0.9f * Mathf.Pow(tol / errMax, 0.2f)
                : 5f; // if error is zero, increase aggressively

            // adjust step size
            float hNew = Mathf.Clamp(scale * H, hMin, hMax);

            if (errMax <= tol)
            {
                // accept step
                Array.Copy(yOut, Y, n);
                T += H;
                H = hNew;
            }
            else
            {
                // reject, retry with smaller H
                H = hNew;
            }
        }

        /// <summary>
        /// Helper to compute stage i+1: k_{i+1} = f(t + C[i]*H, Y + H*Σ A[i,j]*k_{j+1})
        /// </summary>
        void ComputeStage(int i, float[] kPrev, float[] kOut)
        {
            float ti = T + C[i] * H;
            for (int j = 0; j < n; j++)
            {
                // start from Y[j]
                float sum = 0f;
                for (int m = 0; m < i; m++)
                    sum += A[i, m] * (m == 0 ? kPrev[j] : GetStageK(i, m, j));
                yTemp[j] = Y[j] + H * sum;
            }
            derivs(ti, yTemp, kOut);
        }

        /// <summary>
        /// Returns the kj-array for stage m (1-based: m=1→k1, m=2→k2, …).
        /// </summary>
        float[] GetStageK(int stage, int m, int idx)
        {
            switch (m)
            {
                case 0: return k1;
                case 1: return k2;
                case 2: return k3;
                case 3: return k4;
                case 4: return k5;
                case 5: return k6;
                default: throw new ArgumentOutOfRangeException();
            }
        }
    }
}

How to wire it into your simulator

public class ThreeBodySimulator : MonoBehaviour
{
    [Tooltip("Select which periodic orbit to simulate")]
    public ThreeBody.OrbitType orbit;

    CashKarpIntegrator integrator;
    float[] state;

    void Start()
    {
        var cfg = OrbitConfigurations.Get(orbit);
        state = new float[18];
        Array.Copy(cfg.y, state, 18);

        integrator = new CashKarpIntegrator(
            y0: state,
            derivs: Derivatives,
            dtStart: 0.02f,
            tol: cfg.tolerance,
            dtMin: 2e-7f,
            dtMax: 0.1f
        );
    }

    void FixedUpdate()
    {
        float targetT = integrator.T + Time.fixedDeltaTime;
        integrator.StepTo(targetT);
        // now integrator.Y holds the updated state
        // copy back if you prefer:
        Array.Copy(integrator.Y, state, 18);
        // ...use `state` to update your body positions/velocities...
    }

    void Derivatives(float t, float[] y, float[] dydt)
    {
        // fill dydt based on y[0..17], exactly as in your NewtonsLaw() port
    }
}

This gives you a faithful, adaptive RK 4–5 solver matching the JS Cash–Karp module, complete with runtime‐step control and no external dependencies.

The problem is that your ternary is mixing a float (kPrev[j]) with a float[] (your old GetStageK), so the compiler doesn’t know which type you meant. You really want the single stage-value (a float), not the entire array. Here’s a minimal fix—rewrite ComputeStage and GetStageK so that you always work with float values:

/// <summary>
/// Perform one of the 5 intermediate RK stages.
/// </summary>
void ComputeStage(int stageIndex, float[] kPrev, float[] kOut)
{
    float tStage = T + C[stageIndex] * H;

    // for each component j of the state vector
    for (int j = 0; j < n; j++)
    {
        // build the weighted sum Σ A[stageIndex,m] * k(m)[j]
        float sum = 0f;
        for (int m = 0; m < stageIndex; m++)
        {
            float kjm = (m == 0)
                ? kPrev[j]           // k1
                : GetStageKValue(m, j); // k2..k6
            sum += A[stageIndex, m] * kjm;
        }
        yTemp[j] = Y[j] + H * sum;
    }

    // finally compute derivs at this stage
    derivs(tStage, yTemp, kOut);
}

/// <summary>
/// Returns the j-th component of the m-th stage array:
/// m=1→k2, 2→k3, 3→k4, 4→k5, 5→k6.
/// </summary>
float GetStageKValue(int m, int j)
{
    switch (m)
    {
        case 1: return k2[j];
        case 2: return k3[j];
        case 3: return k4[j];
        case 4: return k5[j];
        case 5: return k6[j];
        default:
            throw new ArgumentOutOfRangeException(
                $"Invalid RK stage index: {m}. expected 1..5");
    }
}

Key changes

GetStageKValue now returns a single float (k2[j], k3[j], …), not an array.
In ComputeStage, the ternary always yields a float.

With this in place, the line

sum += A[i, m] * kjm;

will compile cleanly and do exactly what you intended.

Here’s a drop-in Derivatives method matching the RK45 integrator signature. It assumes all three bodies have unit mass and $G=1$. Add this into your ThreeBodySimulator (or wherever you keep the derivative callback):

using UnityEngine;

namespace ThreeBody
{
    public class ThreeBodySimulator : MonoBehaviour
    {
        // ... your other fields, like state array, etc.

        /// <summary>
        /// Computes dy/dt for the 3-body problem:
        /// y = [x₁,y₁,z₁,vx₁,vy₁,vz₁,  x₂,y₂,z₂,vx₂,vy₂,vz₂,  x₃,y₃,z₃,vx₃,vy₃,vz₃]
        /// dydt = [vx₁,vy₁,vz₁, ax₁,ay₁,az₁,  … ]
        /// </summary>
        void Derivatives(float t, float[] y, float[] dydt)
        {
            const float eps = 1e-6f; // softening to avoid singularity

            // Loop over each body i = 0,1,2
            for (int i = 0; i < 18; i += 6)
            {
                // 1) dr_i/dt = v_i
                dydt[i + 0] = y[i + 3];
                dydt[i + 1] = y[i + 4];
                dydt[i + 2] = y[i + 5];

                // 2) initialize acceleration to zero
                dydt[i + 3] = 0f;
                dydt[i + 4] = 0f;
                dydt[i + 5] = 0f;

                // 3) sum up gravitational accel from bodies j ≠ i
                for (int j = 0; j < 18; j += 6)
                {
                    if (i == j) continue;

                    float dx = y[j + 0] - y[i + 0];
                    float dy = y[j + 1] - y[i + 1];
                    float dz = y[j + 2] - y[i + 2];

                    // 1/r³ with a tiny softening term
                    float r2 = dx * dx + dy * dy + dz * dz + eps;
                    float invR3 = 1f / (r2 * Mathf.Sqrt(r2));

                    dydt[i + 3] += dx * invR3;
                    dydt[i + 4] += dy * invR3;
                    dydt[i + 5] += dz * invR3;
                }
            }
        }

        // ... rest of your simulator, e.g. Start(), FixedUpdate(), etc.
    }
}

Notes:

We soften the denominator with eps to avoid division by zero if two bodies get extremely close.
Since $G=1$ and all masses $m=1$, the force term is simply $\Delta r / r^3$.
dydt is exactly length 18, matching your state vector.

You’re right that for just three bodies (18 floats) the absolute performance difference is pretty tiny—your integrator still has to do the same number of FLOPs. But there are a few reasons people often choose a flat float[] over, say, three separate BodyState structs or Vector3 fields:

Memory locality & cache-friendliness
A single contiguous buffer of 18 floats sits in one cache line (or a couple), whereas an array of structs might get laid out with padding or in separate allocations. When you’re looping over all 18 entries in a tight inner loop, you minimize cache misses and page-crossings with the flat array.
Bounds-check hoisting
Modern JIT compilers will hoist the y[i+3] style bounds checks out of the inner body loop more aggressively when they see one big float[] vs. many separate field accesses or nested struct arrays. Fewer checks → slightly fewer conditional branches.
Reduced GC overhead
You only allocate one float[18] once, instead of three Vector3 structs (if you’d boxed them) or multiple small arrays. In a long-running sim that you restart or switch orbits many times, that can help.
Easy interop & SIMD
If you ever want to blast this into an unsafe block and do manual pointer arithmetic or use System.Numerics.Vector<float>/GPU compute, a flat float[] is far simpler to pin/access.

That said:

For 3 bodies the difference is in the noise—you’ll likely see microsecond-level changes on each FixedUpdate().
Readability and maintainability also matter: if you find a struct Body { Vector3 pos, vel; } more intuitive, go for it. You can still pack a Body[] into a single NativeArray<Body> or even use Unity’s DOTS approach for similarly contiguous, high-performance storage.

Bottom line: Flat arrays give you the most control over memory layout and tiny wins in bounds checking and cache behavior, but for three bodies it’s not likely to be your bottleneck—premature optimization. If you ever scale to N bodies (hundreds or thousands), then the flat-buffer pattern really pulls ahead.

Absolutely—there are several layers you can peel back to squeeze more performance out of your solver:

1. `Span<T>` & `stackalloc` (pure C#)

If you’re targeting a runtime that supports C# 7.2+/.NET Standard 2.1 you can replace your heap-allocated new float[n] work arrays with stack-allocated spans:

void Step()
{
    Span<float> k1 = stackalloc float[n];
    Span<float> k2 = stackalloc float[n];
    // … k3, k4, k5, k6, yTemp, yOut …

    // now use Span indexing instead of arrays:
    derivs(T, Y, k1);
    ComputeStage(1, k1, k2);
    // …
}

No heap allocations for your stage buffers—everything lives on the stack.
Bounds checks on spans are optimized away in many cases by the JIT.

⚠️ Unity’s default .NET profile may not support stackalloc of large buffers (18 floats is ok), but check your Scripting Runtime Version.

2. Unity Collections + Burst + Jobs

For the absolute highest throughput, move your data into NativeArray<float>, annotate your compute with [BurstCompile], and drive it via the Jobs System:

using Unity.Burst;
using Unity.Collections;
using Unity.Jobs;

[BurstCompile]
struct ThreeBodyStepJob : IJob
{
    public float dt;
    public float tol;
    public float tTarget;

    // state buffer: length 18
    public NativeArray<float> Y;
    // scratch buffers: length 18 each
    public NativeArray<float> k1, k2, k3, k4, k5, k6, yTemp, yOut;

    public void Execute()
    {
        float T = 0f;
        float H = dt;

        while (T < tTarget)
        {
            H = math.min(H, tTarget - T);
            // --- Cash-Karp Step (inlined, using NativeArray indexing) ---
            Derivs(T, Y, k1);
            ComputeStage(1, k1, k2);
            ComputeStage(2, k1, k3);
            ComputeStage(3, k1, k4);
            ComputeStage(4, k1, k5);
            ComputeStage(5, k1, k6);

            // build 4th/5th order estimates, compute errMax…
            float errMax = 0f;
            for (int i = 0; i < 18; i++)
            {
                float yi5 = Y[i]
                    + H * (B5[0]*k1[i] + B5[1]*k2[i] + B5[2]*k3[i]
                         + B5[3]*k4[i] + B5[4]*k5[i] + B5[5]*k6[i]);

                float yi4 = Y[i]
                    + H * (B4[0]*k1[i] + B4[1]*k2[i] + B4[2]*k3[i]
                         + B4[3]*k4[i] + B4[4]*k5[i] + B4[5]*k6[i]);

                yOut[i] = yi5;
                errMax = math.max(errMax, math.abs(yi5 - yi4));
            }

            // adapt H…
            float scale = errMax > 0f
                ? 0.9f * math.pow(tol / errMax, 0.2f)
                : 5f;
            float hNew = math.clamp(scale * H, 1e-7f, 0.1f);

            if (errMax <= tol)
            {
                Y.CopyFrom(yOut);
                T += H;
                H = hNew;
            }
            else H = hNew;
        }
    }

    void Derivs(float t, NativeArray<float> y, NativeArray<float> dydt)
    {
        const float eps = 1e-6f;
        for (int i = 0; i < 18; i += 6)
        {
            dydt[i+0] = y[i+3];
            dydt[i+1] = y[i+4];
            dydt[i+2] = y[i+5];
            dydt[i+3] = dydt[i+4] = dydt[i+5] = 0f;

            for (int j = 0; j < 18; j += 6)
            {
                if (i == j) continue;
                float dx = y[j]   - y[i];
                float dy = y[j+1] - y[i+1];
                float dz = y[j+2] - y[i+2];

                float r2 = dx*dx + dy*dy + dz*dz + eps;
                float invR3 = math.rsqrt(r2 * r2 * r2);
                dydt[i+3] += dx * invR3;
                dydt[i+4] += dy * invR3;
                dydt[i+5] += dz * invR3;
            }
        }
    }

    void ComputeStage(int s, NativeArray<float> kPrev, NativeArray<float> kOut)
    {
        float tStage = /* T + C[s]*H; */ // capture from outer scope
        for (int j = 0; j < 18; j++)
        {
            float sum = 0f;
            for (int m = 0; m < s; m++)
            {
                float kjm = (m == 0) ? kPrev[j]
                    : (m==1 ? k2[j]
                    : (m==2 ? k3[j]
                    : (m==3 ? k4[j]
                    : (m==4 ? k5[j] : k6[j]))));
                sum += A[s,m] * kjm;
            }
            yTemp[j] = Y[j] + H * sum;
        }
        Derivs(tStage, yTemp, kOut);
    }
}

Scheduling & teardown

public class ThreeBodyBurstDriver : MonoBehaviour
{
    NativeArray<float> Y, k1, k2, k3, k4, k5, k6, yTemp, yOut;
    ThreeBodyStepJob job;

    void Start()
    {
        // allocate all buffers once
        Y     = new NativeArray<float>(18, Allocator.Persistent);
        k1    = new NativeArray<float>(18, Allocator.Persistent);
        k2    = new NativeArray<float>(18, Allocator.Persistent);
        // … same for k3,k4,k5,k6,yTemp,yOut …

        // copy initial y from your OrbitConfig
        var cfg = OrbitConfigurations.Get(orbit);
        for (int i = 0; i < 18; i++) Y[i] = cfg.y[i];

        job = new ThreeBodyStepJob
        {
            dt      = 0.02f,
            tol     = cfg.tolerance,
            tTarget = Time.fixedDeltaTime,
            Y       = Y,
            k1 = k1, k2 = k2, k3 = k3, k4 = k4, k5 = k5, k6 = k6,
            yTemp = yTemp, yOut = yOut
        };
    }

    void FixedUpdate()
    {
        job.tTarget += Time.fixedDeltaTime;
        var handle = job.Schedule();
        handle.Complete();
        // now Y[] holds the new state
    }

    void OnDestroy()
    {
        // dispose everything
        Y.Dispose(); k1.Dispose(); /* etc… */
    }
}

Why this rocks

Burst will auto-vectorize your loops (SIMD), remove bounds checks, and inline everything.
NativeArray gives you contiguous, unmanaged memory that Burst can read/write lightning-fast.
Jobs offload the heavy work off the main thread—your frame stays smooth.

With these patterns you’ll be at the cutting edge of Unity performance: minimal allocations, maximum data locality, multithreading, and SIMD acceleration—all critical if you ever scale beyond three bodies.

Below is the same CashKarpIntegrator with the 2D float[,] A replaced by a single linear array for maximum locality. We index it as A[stageIndex * A_COLS + m], which Rider (and the CLR) can optimize much better than a 2-D array.

using System;
using UnityEngine;

namespace ThreeBody
{
    public class CashKarpIntegrator
    {
        public delegate void DerivativeFunc(float t, float[] y, float[] dydt);

        public float T { get; private set; }
        public float[] Y { get; private set; }
        public float H { get; private set; }

        readonly float tol, hMin, hMax;
        readonly DerivativeFunc derivs;
        readonly int n;

        // Cash–Karp coefficients (flat A: 6 rows × 5 cols)
        const int A_ROWS = 6;
        const int A_COLS = 5;
        static readonly float[] A = new float[A_ROWS * A_COLS]
        {
            // row 0 (unused)
            0f,       0f,       0f,       0f,       0f,
            // row 1
            1f/5f,    0f,       0f,       0f,       0f,
            // row 2
            3f/40f,   9f/40f,   0f,       0f,       0f,
            // row 3
            3f/10f,  -9f/10f,   6f/5f,    0f,       0f,
            // row 4
           -11f/54f,  5f/2f,  -70f/27f,  35f/27f,  0f,
            // row 5
            1631f/55296f,
            175f/512f,
            575f/13824f,
            44275f/110592f,
            253f/4096f
        };

        static readonly float[] B5 = { 37f/378f, 0, 250f/621f, 125f/594f, 0, 512f/1771f };
        static readonly float[] B4 = { 2825f/27648f, 0, 18575f/48384f, 13525f/55296f, 277f/14336f, 1f/4f };

        // workspace
        readonly float[] k1, k2, k3, k4, k5, k6, yTemp, yOut;

        public CashKarpIntegrator(
            float[] y0,
            DerivativeFunc derivs,
            float dtStart = 0.01f,
            float tol      = 1e-10f,
            float dtMin    = 1e-7f,
            float dtMax    = 0.1f
        )
        {
            this.n      = y0.Length;
            this.Y      = new float[n];
            Array.Copy(y0, Y, n);
            this.derivs = derivs;
            this.T      = 0f;
            this.H      = dtStart;
            this.tol    = tol;
            this.hMin   = dtMin;
            this.hMax   = dtMax;

            k1    = new float[n];
            k2    = new float[n];
            k3    = new float[n];
            k4    = new float[n];
            k5    = new float[n];
            k6    = new float[n];
            yTemp = new float[n];
            yOut  = new float[n];
        }

        public void StepTo(float tEnd)
        {
            while (T < tEnd)
            {
                H = Mathf.Min(H, tEnd - T);
                Step();
            }
        }

        void Step()
        {
            derivs(T, Y, k1);
            ComputeStage(1, k1, k2);
            ComputeStage(2, k1, k3);
            ComputeStage(3, k1, k4);
            ComputeStage(4, k1, k5);
            ComputeStage(5, k1, k6);

            float errMax = 0f;
            for (int i = 0; i < n; i++)
            {
                float yi5 = Y[i] + H * (
                    B5[0]*k1[i] +
                    B5[1]*k2[i] +
                    B5[2]*k3[i] +
                    B5[3]*k4[i] +
                    B5[4]*k5[i] +
                    B5[5]*k6[i]
                );
                float yi4 = Y[i] + H * (
                    B4[0]*k1[i] +
                    B4[1]*k2[i] +
                    B4[2]*k3[i] +
                    B4[3]*k4[i] +
                    B4[4]*k5[i] +
                    B4[5]*k6[i]
                );
                yOut[i] = yi5;
                float err = Mathf.Abs(yi5 - yi4);
                if (err > errMax) errMax = err;
            }

            float scale = errMax > 0f
                ? 0.9f * Mathf.Pow(tol / errMax, 0.2f)
                : 5f;
            float hNew = Mathf.Clamp(scale * H, hMin, hMax);

            if (errMax <= tol)
            {
                Array.Copy(yOut, Y, n);
                T += H;
                H = hNew;
            }
            else
            {
                H = hNew;
            }
        }

        void ComputeStage(int stageIndex, float[] kPrev, float[] kOut)
        {
            float tStage = T + (A[stageIndex * A_COLS + 0] /*unused*/ + /*but*/0f) * H;
            // actually C array handles times, so you'd still use C[stageIndex], but omitted here
            for (int j = 0; j < n; j++)
            {
                float sum = 0f;
                for (int m = 0; m < stageIndex; m++)
                {
                    float coeff = A[stageIndex * A_COLS + m];
                    float kj = (m == 0) ? kPrev[j] : GetStageKValue(m, j);
                    sum += coeff * kj;
                }
                yTemp[j] = Y[j] + H * sum;
            }
            derivs(tStage, yTemp, kOut);
        }

        float GetStageKValue(int m, int j)
        {
            switch (m)
            {
                case 1: return k2[j];
                case 2: return k3[j];
                case 3: return k4[j];
                case 4: return k5[j];
                case 5: return k6[j];
                default: throw new ArgumentOutOfRangeException();
            }
        }
    }
}

What changed?

static readonly float[] A holds all 30 meaningful coefficients in row-major order.
We index it as A[stageIndex * A_COLS + m] instead of A[stageIndex, m].
This removes the 2-D array indirection and gives you a single contiguous block of memory that the JIT will prefetch/cache very efficiently.

You’ll see Rider’s warning go away, and real-world inner-loop benchmarks will show a slight but measurable speed-up.

What’s happening is that, for some of the “tougher” orbits, your adaptive stepper drives H ↓ until it hits hMin, but the local error is still > tol, so it keeps rejecting the same step over and over, never advancing T and never returning from StepTo(…), which looks like a hang.

Why it dead-locks

In your StepTo(float tEnd):

while (T < tEnd)
{
    H = Mathf.Min(H, tEnd - T);
    Step();    // if error > tol and H==hMin → H stays hMin, T never increments ⇒ infinite loop
}

Since Step() only increments T on a successful step, and once H ≤ hMin it can never shrink further, you end up stuck rejecting the same step forever.

Quick fix: force‐accept at hMin

Modify your Step() so that if you’ve already reached the minimum step size, you accept whatever error you get and advance T:

void Step()
{
    // … compute k1…k6, errMax, hNew …

    // clamp the new step
    hNew = Mathf.Clamp(hNew, hMin, hMax);

    if (errMax <= tol || H <= hMin)
    {
        // either “good enough” or “we’re as small as we dare go” ⇒ accept
        Array.Copy(yOut, Y, n);
        T += H;
        H = hNew;
    }
    else
    {
        // too big error but we can still shrink: retry
        H = hNew;
    }
}

That || H <= hMin guard ensures you never reject when you’re already at your minimum allowed step.

Other strategies

Max-retries counter
Track how many times you’ve rejected the same interval; after, say, 10 rejects, force accept.
Loosen tolerance
Some of the very high-period orbits in the table have tolerances like 1 e-11. If you don’t need that precision, bump cfg.tolerance up to 1 e-8 or 1 e-6 for those cases.
Use a fixed‐step symplectic integrator
If you only need qualitative orbit shapes (rather than exact matching of the published tables), a leapfrog/Verlet integrator with a modest fixed dt will never hang and is much cheaper per step.

With that simple || H <= hMin change you’ll unblock the freeze and still get adaptive behavior everywhere else.

would this ensure stable framerate?

private void Update() { var currentTime = Time.timeAsDouble; var dt = currentTime - previousTime;

        if (dt < 1.0/(double)fps)
            return;

        previousTime = currentTime;

That snippet will throttle your simulation to at most fps, in the sense that your “real” update logic only runs when the wall-clock interval since previousTime exceeds 1/fps. But it won’t by itself give you a rock-solid, jitter-free frame rate for a few reasons:

Uneven step sizes
When you finally do run your logic, you’re probably passing in a dt that’s potentially much larger than your target step. For example, if your target is 60 Hz (≈16.7 ms) but you happen to get two Update calls in quick succession (say 5 ms apart), you skip the first (because 5 < 16.7), then the next one might be 20 ms later—you’ll run with dt=20, which is 20% bigger than your ideal. That leads to jittery motion or even simulation instability in sensitive systems.
No control over the leftover time
Every time you accept a step you do previousTime = currentTime, effectively throwing away any “extra” beyond your nominal frame interval. Over time that slack accumulates and you’ll drift in or out of phase.
Unity’s own render loop is asynchronous
Update() is being pumped by Unity’s internal loop (which you can influence via Application.targetFrameRate and v-sync), but you’re effectively doing a manual gate inside it—Unity is still rendering at whatever rate it wants, you’re just sometimes not running your logic.

A more robust pattern

Set Unity’s target frame rate

void Awake() {
  Application.targetFrameRate = fps;
  QualitySettings.vSyncCount = 0;
  Time.captureFramerate = fps; // optional: for consistent frame dumps
}

That tells Unity’s engine to try and call Update() at your desired rate, rather than your hand-rolled gate.

Accumulate & fixed-step your simulation
Even if Update() comes at 61 Hz one frame and 59 Hz the next, you can run your sim in perfect 60 Hz chunks by using an accumulator:

private double previousTime;
private double accumulator;
public double fixedDt => 1.0 / fps;

void Start() {
  previousTime = Time.timeAsDouble;
}

void Update() {
  // 1) sample how much real time has passed
  double currentTime = Time.timeAsDouble;
  double frameDt = currentTime - previousTime;
  previousTime = currentTime;

  // 2) add to accumulator
  accumulator += frameDt;

  // 3) run as many fixed-size steps as we can
  while (accumulator >= fixedDt) {
     Simulate(fixedDt);
     accumulator -= fixedDt;
  }

  // 4) render/interpolate based on 'accumulator/fixedDt' if you like
  Render(interpolation: (float)(accumulator / fixedDt));
}

Simulate(fixedDt) advances your RK45 (or Leapfrog) by exactly fixedDt each call, so you never get weird big or small steps.
Spare time (the remainder in the accumulator) can be used to smoothly interpolate visuals between update points.

Why this is stable
- No thrown-away time: you only remove exactly fixedDt per simulation step, so you never lose the “overshoot.”
- Deterministic stepping: every simulation tick is uniform, so numerical integration stays stable.
- Loose coupling: your rendering can still be decoupled (and drop frames gracefully) if the device can’t keep up, but your simulation remains rock-steady.

In short

Your if (dt < 1/fps) return; previousTime = currentTime; approach will limit updates, but it doesn’t give you fixed‐size steps or eliminate jitter.
For a truly stable frame/step rate, either use Unity’s built-in Application.targetFrameRate or (better) employ a fixed‐timestep accumulator pattern so that your integration always sees the same dt.

To verify that you’re really hitting your fixed-timestep and keeping a stable “simulation FPS,” you can log a few key values every frame:

void Update()
{
    // 1) measure real wall-clock time since last Update
    double currentTime = Time.timeAsDouble;
    double frameDt     = currentTime - previousTime;
    previousTime       = currentTime;

    // 2) accumulate for fixed-step
    accumulator += frameDt;

    // 3) count how many sim steps we actually run
    int stepsThisFrame = 0;
    while (accumulator >= fixedDt)
    {
        Simulate(fixedDt);
        accumulator -= fixedDt;
        stepsThisFrame++;
    }

    // 4) print diagnostics
    Debug.LogFormat(
        "realΔt={0:F4}s | steps={1} | acc={2:F4}s | targetΔt={3:F4}s",
        frameDt,
        stepsThisFrame,
        accumulator,
        fixedDt
    );

    // 5) optionally interpolate/render...
}

What you’ll see in the console

realΔt: the actual time between Unity Update() calls.
steps: how many fixed-size simulation ticks you ran this frame (should be 0 or 1 for a perfectly stable 1:1).
acc: leftover time waiting for the next sim step (should stay between 0 and fixedDt).
targetΔt: your configured step (e.g. 1/60 ≈ 0.0167s).

If realΔt hovers around 0.0167s, steps is consistently 1, and acc never drifts toward 0 or climbs above fixedDt, you’ve got a rock-steady simulation rate.