Where will the cat go? (orbital mechanics)
Python 3.5 + NumPy, exact, 186 bytes
from math import*
def o(r,v,t):
d=(r@r)**.5;W=2/d-v@v;U=W**1.5;b=[0,t*U+9]
while 1:
a=sum(b)/2;x=1-cos(a);y=sin(a)/U;k=r@v*x/W+d*y*W
if a in b:return k*v-r*x/W/d+r
b[k+a/U-y>t]=a
This is an exact solution, using a formula I engineered based on Jesper Göranssonhis, “Symmetries of the Kepler problem”, 2015. It uses a binary search to solve the transcendental equation Ax + B cos x + C sin x = D, which has no closed-form solution.
The function expects the position and velocity to be passed in as NumPy arrays:
>>> from numpy import array
>>> o(array([1,0,0]),array([0,-1,0]),1000000000)
array([ 0.83788718, -0.54584345, 0. ])
>>> o(array([-1.1740058273269156,8.413493259550673,0.41996042044140003]),array([0.150014367067652,-0.09438816345868332,0.37294941703455975]),7999.348650387233)
array([-4.45269544, 6.93224929, -9.27292488])
Javascript
This is just to get the ball rolling, since no one seems to be posting answers. Here is a very naive, simple way that can be improved a lot:
function simulate(x, y, z, vx, vy, vz, t) {
var loops = 1884955; // tune this parameter
var timestep = t / loops;
for (var i = 0; i < t; i += timestep) {
var distanceSq = x*x + y*y + z*z; // distance squared from origin
var distance = Math.sqrt(distanceSq);
var forceMag = 1/distanceSq; // get the force of gravity
var forceX = -x / distance * forceMag;
var forceY = -y / distance * forceMag;
var forceZ = -z / distance * forceMag;
vx += forceX * timestep;
vy += forceY * timestep;
vz += forceZ * timestep;
x += vx * timestep;
y += vy * timestep;
z += vz * timestep;
}
return [x, y, z];
}
Testing:
simulate(1, 0, 0, 0, -1, 0, Math.PI*2) --> [0.9999999999889703, -0.0000033332840909716455, 0]
Hey, that's pretty good. It has an error of about 3.333*10^(-6) which isn't enough for it to be rounded down... it's close.
Just for fun:
console.log(simulate(1, 0, 0, 0, -1, 0, 1000000000))
--> [-530516643639.4616, -1000000000.0066016, 0]
Oh well; so it's not the best.
And on a random test case from the generator:
simulate(-1.1740058273269156,8.413493259550673,0.41996042044140003,0.150014367067652,-0.09438816345868332,0.37294941703455975,7999.348650387233)
--> [-4.528366392498373, 6.780385554803544, -9.547824236472668]
Actual:[-4.452695438880813, 6.932249293597744, -9.272924876103785]
With an error of only approximately 0.32305!
This can be improved a lot by using something like Verlet integration or some fancy algorithm. In fact, those algorithms may even get perfect scores, despite being simulations.