partially elastic collisions

Readme.md

@@ -59,3 +59,5 @@ TODO
 - [ ] Enhancement: Track farthest reaches, min/max in each dimension (x, y)
 - [x] Fix: Unpause panning when initiated while sim is paused
 - [ ] Enhancement: Refactor to use viewOrigin as center of display canvas
+- [ ] Feature: Dark Matter
+  - [ ] Enhancement: Clumps of dark matter to help gather stray objects in the universe

system.js

@@ -155,8 +155,7 @@ export class System {
       this.forEachObject((A, i) => {
         this.forEachObject((B, j) => {
           if (this.objectsOverlap(A, B)) {
-            // this.mergeObjects(i, j);
+            this.bounceObjects(i, j, elapsedTime);
-            this.bounceObjects(i, j);
           }
         }, {alive: true, startWith: i + 1});
       });
@@ -397,7 +396,7 @@ export class System {
     return d < A.radius + B.radius;
   }
 
-  bounceObjects(i, j) {
+  bounceObjects(i, j, elapsedTime) {
     const A = this.object(i);
     const B = this.object(j);
     const totalMass = A.mass + B.mass;
@@ -421,6 +420,9 @@ export class System {
       return;
     }
 
+    // TODO: Relative tangential velocity, and initial rotation of each object, should impart angular momentum
+    // Also TODO: In partially elastic case, compute partial angular momentum transfer
+
     // Are these objects sticking together?
     // One way to determine this is,
     // does either object have enough KE to escape their gravitational potential?
@@ -443,11 +445,67 @@ export class System {
     const vA = add(mult(tangent, vAt), mult(normal, vAnNew));
     const vB = div(add(mult(A.velocity, A.mass), mult(B.velocity, B.mass), mult(vA, -A.mass)), B.mass);
 
-    // const Ki = A.kineticEnergy + B.kineticEnergy;
+    // What about position?? The objects are still overlapping.
-    A.velocity = vA;
+    // If we fire and forget, we're hoping the object escapes within the next frame.
-    B.velocity = vB;
+    // But it often doesn't. We catch some of these by verifing a positive relative radial velocity.
-    // const Kf = A.kineticEnergy + B.kineticEnergy;
+    // However, if the objects have sufficient tangential velocity, 
-    // console.log('Collision: Zscaled', Zscaled, 'Energy before', Ki, 'Energy after', Kf, 'Fraction change', (Kf - Ki) / Ki);
+    // or if we're imparting tangential velocity due to initial rotation in partially elastic cases,
+    // then they can remain overlapping in the next frame despite having substantial velocity.
+    // Then, in the next frame, since the masses are so close together, they accelerate toward each other
+    // and may remain overlapping in the next frame; and so they can orbit in this overlapping manner.
+    // However, this is an unstable orbit due to its sensitivity to parameters, including framerate!
+    // Therefore it is desireable not only to mitigate, but to prevent this scenario by enforcing that
+    // the objects never remain overlapping at the end of a frame computation!
+    // The question then becomes, what is the best accuracy we can reasonably hope to achieve?
+    // We are altering the objects gravitational potentials if we move them from their simulated positions.
+    // This might be fine given the relative infrequence of collisions.
+
+    // 1. So a first order solution would be to back the object to its initial position, and leave it there.
+    // 2. An improvement could be to move it from its initial position, to the position it would have after the same duration
+    // as the current frame had elapsed, but with the velocity we've computed for it after collision.
+    // 3. A further improvement would be to linearly interpolate along its computed path, and then update its position
+    // using the remaining time in the current interval.
+    // 4. A yet further improvement would be to interpolate quadratically when estimating the point of collision!
+   
+    // The linear interpolation solution sounds likely to be sufficient, given that the risk of error is not even a change in
+    // the resulting deflection angle, but only of the exact position of the object.
+    // What are likely issues with the simplest solution? An object could perhaps get stuck in place if we fail to update its position.
+    // That brings us to the second solution. It sounds likely fine; the risk of error is just a slight change in the object's position,
+    // which is guaranteed to be less than a frame's worth of distance travelled.
+    // TODO: Specify desired level of accuracy ^_^
+
+    // Well, by experimenting we found the problem with the second solution; because the objects are close together,
+    // giving them an altitude boost greatly diminishes their gravitational potential, increasing their separation.
+    // It turns out that's the problem with the first solution too!
+    // But it depends on how far into the frame the objects collide. Which brings us to interpolation.
+   
+    // Determine time of collision, relative to start of frame.
+    // We know the objects are accelerating toward each other, possibly rapidly.
+    // Therefore we really need to interpolate quadratically.
+
+    // But even that will not be accurate, as the force is increasing substantially during the interval as well.
+    //
+    // What we may really want to do is to run the whole frame computation multiple times, searching for a moment before impact, to whatever specified precision.
+    // Then we re-run the whole frame from that point...
+    // Another pathway here is to approximate the resulting position but then to adjust
+    // the resulting velocity in order to preserve the object's final energy,
+    // by subtracting for the decrease in gravitational potential energy.
+
+    // The objects never should have gotten so close... we should probably do interpolation and re-run whole frame.
+    
+    A.position = add(A.currentPosition, mult(elapsedTime, vA));
+    B.position = add(B.currentPosition, mult(elapsedTime, vB));
+    
+    const G = this.sim.getOption('param.gravity');
+    const dV = 2 * G * (
+      1 / magnitude(sub(A.position, B.position)) - 
+      1 / magnitude(sub(A.currentPosition, B.currentPosition))
+    );
+    const vAmag = Math.sqrt(square(A.velocity) - B.mass * dV);
+    const vBmag = Math.sqrt(square(B.velocity) - A.mass * dV);
+
+    A.velocity = mult(vA, div(vAmag, magnitude(vA)));
+    B.velocity = mult(vB, div(vBmag, magnitude(vB)));
   }
 
   mergeObjects(i, j) {