16
\$\begingroup\$

A nearly massless cat is dropped in space (don't worry, with a space suit and everything) at the point (x, y, z) with velocity (vx, vy, vz). There is an fixed, infinitely dense planet (with volume of 0) at the point (0, 0, 0) and it attracts objects at distance r with acceleration 1/r^2. According to Newtonian gravity, where does the object go after time t?

Nearly massless in this case means you are outputting the value of lim (mass --> 0) <position of cat>. The mass is affected by the planet's gravity, but the planet is not affected by the cat's gravity. In other words, the central body is fixed.

This is somewhat similar to Code Golf: What is the spaceship's fate? [floating point version], but this is different because it is measuring accuracy.

You may implement a solution based on a simulation, which must run in less than 3 seconds, OR you may implement a program that gives an exact value (must also run in less than 3 seconds). See scoring details below. If you implement a simulation, it does not have to be exact, but your score will be lower because of the inaccuracy.

Input: x y z vx vy vz t, not necessarily integers representing x, y, z coordinates, velocity in the x, y, and z directions and time respectively. It is guaranteed that the speed of the cat is strictly less than the escape velocity at that altitude. Input may be taken from anywhere, including parameters to a function. The program must run in less than three seconds on my laptop for t < 2^30, which means, if you are running a simulation, you must adjust your timestep accordingly. If you are planning on hitting the 3 second limit for every test case, make sure there is a tuneable parameter that can make it more accurate/less accurate for speed gains, so that I can make it run in three seconds on my computer.

Output: x y z, the position after time t.

Since the two-body problem can be solved perfectly, it is in theory possible to get a perfect, correct answer.

Scoring: For any test case, the error is defined to be the distance between your output and the "true" output. The true output is defined to be the one that the test case snippet generates. If the error is less than 10^(-8), the error is rounded down to zero. Your score is the average error on 100 (or more) random test cases. If you write a perfectly accurate answer, you should get a score of 0; lowest score wins, and ties will be broken by code length.

Test cases:

1 0 0 0 -1 0 1000000000 --> 0.83789 -0.54584 0

In this case , the orbit is perfectly circular with period 2*pi, so after circling 159154943 times, the cat ends up at approximately (0.83789, -0.54584). This isn't a test case your code will be tested on; if you submit a perfectly accurate answer, however, you may want to test it on this.

The snippet below generates random additional test cases and will be used to judge submissions; let me know if there is a bug with this:

var generateTestCase = function() {
  var periapsis = Math.random() * 10 + 0.5;
  var circularVelocity = Math.sqrt(1 / periapsis);
  var escapeVelocity = Math.sqrt(2) * circularVelocity;
  var v = Math.random() * (escapeVelocity - circularVelocity) + circularVelocity;
  var a = 1 / (2 / periapsis - v * v);
  var ecc = (a - periapsis) * (2 / periapsis - v * v);
  var semiLatus = a * (1 - ecc * ecc);
  var angM = Math.sqrt(semiLatus);
  var period = 2 * Math.PI * a * Math.sqrt(a);
  var getCoordAfterTime = function(t) {
    var meanAnomaly = t / (a * Math.sqrt(a));

    function E(e, M, n) {
      var result = M;
      for (var k = 1; k <= n; ++k) result = M + e * Math.sin(result);
      return result;
    }
    var eAnomaly = E(ecc, meanAnomaly, 10000);
    var radius = a * (1 - ecc * Math.cos(eAnomaly));
    var theta = Math.acos((semiLatus / radius - 1) / ecc);
    if (t > period / 2) theta = -theta;
    return {
      r: radius,
      t: theta,
      x: radius * Math.cos(theta),
      y: radius * Math.sin(theta)
    };
  }
  var t1 = Math.random() * period;
  var p1 = getCoordAfterTime(t1);
  var v1 = Math.sqrt(2 / p1.r - 1 / a);
  var velocityAngle1 = Math.acos(angM / p1.r / v1);
  var horizontal = Math.atan(-p1.x / p1.y);
  if (t1 < period / 2) {
    velocityAngle1 *= -1;
    horizontal += Math.PI;
  }
  var vx = v1 * Math.cos(horizontal + velocityAngle1);
  var vy = v1 * Math.sin(horizontal + velocityAngle1);
  var t2 = Math.random() * period;
  var p2 = getCoordAfterTime(t2);
  var v2 = Math.sqrt(2 / p2.r - 1 / a);
  var dummyPeriods = Math.floor(Math.random() * 100) * period + period;
  var rotate = function(x, y, z, th1, th2, th3) {
    var x1 = x * Math.cos(th1) - y * Math.sin(th1);
    var y1 = x * Math.sin(th1) + y * Math.cos(th1);
    var z1 = z;
    var x2 = x1 * Math.cos(th2) + z1 * Math.sin(th2);
    var y2 = y1;
    var z2 = -x1 * Math.sin(th2) + z1 * Math.cos(th2);
    var x3 = x2;
    var y3 = y2 * Math.cos(th3) - z2 * Math.sin(th3);
    var z3 = y2 * Math.sin(th3) + z2 * Math.cos(th3);
    return [x3, y3, z3];
  }
  var th1 = Math.random() * Math.PI * 2,
    th2 = Math.random() * 2 * Math.PI,
    th3 = Math.random() * 2 * Math.PI;
  var input = rotate(p1.x, p1.y, 0, th1, th2, th3)
    .concat(rotate(vx, vy, 0, th1, th2, th3))
    .concat([t2 - t1 + dummyPeriods]);
  var output = rotate(p2.x, p2.y, 0, th1, th2, th3);
  return input.join(" ") + " --> " + output.join(" ");
}
var $ = function(x) {
  return document.querySelector(x);
}
$("#onetestcase").onclick = function() {
  $("#result1").innerHTML = generateTestCase();
}
$("#huntestcases").onclick = function() {
  var result = [];
  for (var i = 0; i < 100; ++i) result.push(generateTestCase());
  $("#result100").innerHTML = result.join("\n");
}
<div>
  <button id="onetestcase">Generate</button>
  <pre id="result1"></pre>
</div>
<div>
  <button id="huntestcases">Generate 100 test cases</button>
  <pre id="result100"></pre>
</div>

\$\endgroup\$
  • \$\begingroup\$ Is time t given in seconds? If so, would the velocity be given in units per second, or something smaller? \$\endgroup\$ – R. Kap Apr 4 '16 at 18:59
  • \$\begingroup\$ @R. Kap It doesn't matter. t is given in unit time, whatever that is, and velocity will use the same unit. Whether it is in seconds or hours, the answer will be the same. \$\endgroup\$ – soktinpk Apr 4 '16 at 19:04
  • \$\begingroup\$ nearly massless cat Well, what would the exact mass of the cat be? Should we just use 0 as a value for the mass of this cat? \$\endgroup\$ – R. Kap Apr 4 '16 at 19:16
  • \$\begingroup\$ @R. Kap Yes. But it is still affected by gravity (usually, Newton didn't consider objects with no mass to be affected by gravity). So we should consider it to have arbitrarily small mass, and your answer is actually the position as the mass of the cat goes to zero. The main point is that the planet itself is not affected at all by the cat. \$\endgroup\$ – soktinpk Apr 4 '16 at 19:21
  • 2
    \$\begingroup\$ @soktinpk it might be easier to just explicitly say that the central body is fixed. \$\endgroup\$ – Maltysen Apr 4 '16 at 20:30
6
+100
\$\begingroup\$

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])
| improve this answer | |
\$\endgroup\$
  • \$\begingroup\$ What does the @ do? \$\endgroup\$ – R. Kap Apr 15 '16 at 3:42
  • 1
    \$\begingroup\$ It’s a new operator in Python 3.5 that NumPy overloads for numpy.dot (dot product/matrix multiplication). See PEP 465. \$\endgroup\$ – Anders Kaseorg Apr 15 '16 at 3:44
  • \$\begingroup\$ It's great that it's golfed, but this is code-challenge, could you make it a bit clearer, I had some scratchwork in Python done, and can calculate the anomaly, theta, eccentricity, period, etc, but I'm stuck on determining the sign of theta and determining the rotation from x-y reference plane to 3d space. Still, this is really great stuff \$\endgroup\$ – miles Apr 15 '16 at 4:08
  • \$\begingroup\$ @miles Since ties are broken by code length, it makes sense for this to be golfed. \$\endgroup\$ – Mego Apr 15 '16 at 5:01
  • \$\begingroup\$ That's true, as I was working on an exact solution too, since the test case generator only creates elliptical orbits \$\endgroup\$ – miles Apr 15 '16 at 5:22
2
\$\begingroup\$

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.

| improve this answer | |
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.