Intro

For an upcoming space game I've scripted a basic custom universal gravity simulator. I've set the workspace's gravity value to 0 and decided to use BodyForces to simulate gravitational attraction.

The Problem

I've gotten pretty far on the grav sim, but I've encountered an issue regarding the trajectories of the parts. Namely, when a part starts orbiting another part, its orbit inexplicably gets larger each revolution. An example: https://imgur.com/gallery/G4h1sd4

As you can see, the orbits get progressively wider, and yes, there's only two bodies involved.

Code

Force Updater

01runService.Stepped:Connect(function()
02 
03    local parts = {}
04 
05    for _, part in pairs(workspace:GetDescendants()) do
06        local mass = getMass(part)
07        if part ~= workspace.Terrain and not part:FindFirstChild("DoNotBeAttracted") and mass > 0 then
08            parts[part] = {["pos"] = getPos(part), ["mass"] = mass}
09        end
10    end
11 
12    for part, attr in pairs(parts) do
13 
14        local pos = attr.pos
15        local mass = attr.mass
16 
17        -- Calculate the net gravitational forces
18        local netGForce = Vector3.new(0, 0, 0)          -- Possible bottleneck
19 
20        for _, otherpart in pairs(workspace:GetDescendants()) do
21 
22            if part ~= otherpart and otherpart ~= workspace.Terrain and not part:FindFirstChild("DoNotAttract") then
23 
24                local otherpartMass = getMass(otherpart)
25                if otherpartMass > 0 then
26 
27                    local otherpartPos = getPos(otherpart)
28 
29                    local attraction = gravityModule.calculateGravitationalAttraction(pos, mass, otherpartPos, otherpartMass)
30                    local directionVector = (otherpartPos - pos).Unit
31                    netGForce = netGForce + attraction * directionVector
32                end
33            end
34        end
35 
36        -- Make/update BodyForce
37        updateForces(part, netGForce)
38    end
39 
40    -- Loop through all parts and apply velocity change
41    for _, part in pairs(workspace:GetDescendants()) do
42        if part ~= workspace.Terrain and not part:FindFirstChild("DoNotBeAttracted") then
43            if part:FindFirstChild("GravityForce") and part:FindFirstChild("GravityForceValue") then
44                part.GravityForce.Force = part.GravityForceValue.Value
45            end
46        end
47    end
48end)
 View all 48 lines...

Gravitation Calculation Function

1function module.calculateGravitationalAttraction(pointA, massA, pointB, massB)
2    return module.gravitationalConstant*massA*massB/(pointB-pointA).Magnitude^2
3end

Theories

Could the forces not be updating frequently enough, leading to a delayed gravity, in a sense, which would allow the bodies to drift more than they should? If so, how should I fix this?

Could this be similar to what's happening in real life to the Earth and Moon (very very doubtful)?

Any help is appreciated

Custom planetary gravity = slowly enlarging orbits?