From 48d3baf9803f673a101f6360a3d40b8211ec3e10 Mon Sep 17 00:00:00 2001 From: Navan Chauhan Date: Mon, 16 Oct 2023 16:04:29 -0600 Subject: update dependencies --- docs/posts/2023-04-30-n-body-simulation.html | 1065 +++++++++++++------------- 1 file changed, 536 insertions(+), 529 deletions(-) (limited to 'docs/posts/2023-04-30-n-body-simulation.html') diff --git a/docs/posts/2023-04-30-n-body-simulation.html b/docs/posts/2023-04-30-n-body-simulation.html index cb0bf6c..e12380a 100644 --- a/docs/posts/2023-04-30-n-body-simulation.html +++ b/docs/posts/2023-04-30-n-body-simulation.html @@ -89,46 +89,46 @@ 
/*
  * Earth - Sun Orbit Plot
  * Taken from Numerics tutorial
- */
-
-const G = 6.67e-11;
-const Msun = 2e30;
-const AU = 1.5e11;
-const v0 = Math.sqrt(G * Msun / AU); // SI
-
-function dR(r, v) {
-  const dv = [-G * Msun / Math.pow(r[0] ** 2 + r[1] ** 2, 3 / 2) * r[0], -G * Msun / Math.pow(r[0] ** 2 + r[1] ** 2, 3 / 2) * r[1]];
-  const dr = [...v];
-  return [dr, dv];
-}
-
-// initialize system
-let r = [-AU, 0];
-const theta = Math.atan2(r[1], r[0]);
-let v = [-v0 * Math.sin(theta), v0 * Math.cos(theta)];
-
-const t = Array.from({ length: 1001 }, (_, i) => i / 100 + 0.0); // years
-const yearSec = 365 * 24 * 3600;
-const dt = (t[1] - t[0]) * yearSec; // s
-const x4Plot = Array.from({ length: t.length }, () => 0);
-const y4Plot = Array.from({ length: t.length }, () => 0);
-
-// integrate using RK4!
-for (let i = 0; i < t.length; i++) {
-  const k1 = dR(r, v).map(x => x.map(y => y * dt));
-  const k2 = dR(r.map((ri, j) => ri + k1[0][j] / 2), v.map((vi, j) => vi + k1[1][j] / 2)).map(x => x.map(y => y * dt));
-  const k3 = dR(r.map((ri, j) => ri + k2[0][j] / 2), v.map((vi, j) => vi + k2[1][j] / 2)).map(x => x.map(y => y * dt));
-  const k4 = dR(r.map((ri, j) => ri + k3[0][j]), v.map((vi, j) => vi + k3[1][j])).map(x => x.map(y => y * dt));
-  r = r.map((ri, j) => ri + (k1[0][j] + 2 * k2[0][j] + 2 * k3[0][j] + k4[0][j]) / 6);
-  v = v.map((vi, j) => vi + (k1[1][j] + 2 * k2[1][j] + 2 * k3[1][j] + k4[1][j]) / 6);
-  x4Plot[i] = r[0];
-  y4Plot[i] = r[1];
-}
-
-// make data for plot
-var sun = { x: 0, y: 0 };
-const earth = { x: x4Plot.map(x => x / AU), y: y4Plot.map(y => y / AU) };
-const circle = { x: Array.from({ length: 1001 }, (_, i) => Math.cos(i / 100 * 2 * Math.PI)), y: Array.from({ length: 1001 }, (_, i) => Math.sin(i / 100 * 2 * Math.PI)) };
+ */
+
+const G = 6.67e-11;
+const Msun = 2e30;
+const AU = 1.5e11;
+const v0 = Math.sqrt(G * Msun / AU); // SI
+
+function dR(r, v) {
+  const dv = [-G * Msun / Math.pow(r[0] ** 2 + r[1] ** 2, 3 / 2) * r[0], -G * Msun / Math.pow(r[0] ** 2 + r[1] ** 2, 3 / 2) * r[1]];
+  const dr = [...v];
+  return [dr, dv];
+}
+
+// initialize system
+let r = [-AU, 0];
+const theta = Math.atan2(r[1], r[0]);
+let v = [-v0 * Math.sin(theta), v0 * Math.cos(theta)];
+
+const t = Array.from({ length: 1001 }, (_, i) => i / 100 + 0.0); // years
+const yearSec = 365 * 24 * 3600;
+const dt = (t[1] - t[0]) * yearSec; // s
+const x4Plot = Array.from({ length: t.length }, () => 0);
+const y4Plot = Array.from({ length: t.length }, () => 0);
+
+// integrate using RK4!
+for (let i = 0; i < t.length; i++) {
+  const k1 = dR(r, v).map(x => x.map(y => y * dt));
+  const k2 = dR(r.map((ri, j) => ri + k1[0][j] / 2), v.map((vi, j) => vi + k1[1][j] / 2)).map(x => x.map(y => y * dt));
+  const k3 = dR(r.map((ri, j) => ri + k2[0][j] / 2), v.map((vi, j) => vi + k2[1][j] / 2)).map(x => x.map(y => y * dt));
+  const k4 = dR(r.map((ri, j) => ri + k3[0][j]), v.map((vi, j) => vi + k3[1][j])).map(x => x.map(y => y * dt));
+  r = r.map((ri, j) => ri + (k1[0][j] + 2 * k2[0][j] + 2 * k3[0][j] + k4[0][j]) / 6);
+  v = v.map((vi, j) => vi + (k1[1][j] + 2 * k2[1][j] + 2 * k3[1][j] + k4[1][j]) / 6);
+  x4Plot[i] = r[0];
+  y4Plot[i] = r[1];
+}
+
+// make data for plot
+var sun = { x: 0, y: 0 };
+const earth = { x: x4Plot.map(x => x / AU), y: y4Plot.map(y => y / AU) };
+const circle = { x: Array.from({ length: 1001 }, (_, i) => Math.cos(i / 100 * 2 * Math.PI)), y: Array.from({ length: 1001 }, (_, i) => Math.sin(i / 100 * 2 * Math.PI)) };
@@ -150,76 +150,77 @@ Next, the function dR takes the position r and velocity v of Earth as input and

Plotting the orbit

-

Now that we have the data for the Sun's position, Earth's orbit, and the reference circle, we can plot them using Plotly.js. -

- 
let traceSun = {
-      x: [sun.x],
-      y: [sun.y],
-      mode: "markers",
-      marker: {
-        symbol: "star",
-        size: 10,
-        color: "gold",
-      },
-      name: "Sun",
-    };
- 

-    const traceEarth = {
-      x: earth.x,
-      y: earth.y,
-      mode: "lines",
-      line: {
-        color: "white"
-      },
-      name: "Earth",
-    };
- 

-    const traceOrbit = {
-      x: circle.x,
-      y:circle.y,
-      mode: "lines",
-      line: {
-        color: "crimson",
-        dash: "dash"
-      },
-      name: "1 AU Circle",
-    };
- 

-    const earthSunLayout = {
-      title: "Earth-Sun Orbit",
-      xaxis: {
-        title: "x [AU]",
-        range: [-1.1,1.1],
-        showgrid: true,
-        gridcolor: "rgba(255,255,255,0.5)",
-        gridwidth: 1,
-        zeroline: true,
-        tickmode: "auto",
-        nticks: 5,
-      },
-      yaxis: {
-        title: "y [AU]",
-        range: [-1.1,1.1],
-        showgrid: true,
-        gridcolor: "rgba(255,255,255,0.5)",
-        gridwidth: 1,
-        zeroline: false,
-        tickmode: "auto",
-        nticks: 5,
-      },
-      margin: {
-        l: 50,
-        r: 50,
-        b: 50,
-        t: 50,
-        pad: 4,
-      },
-      paperbgcolor: "black",
-      plotbgcolor: "black",
-    };
-    Plotly.newPlot("plot",[traceSun,traceEarth,traceOrbit], earthSunLayout);
-
-

+

Now that we have the data for the Sun's position, Earth's orbit, and the reference circle, we can plot them using Plotly.js.

+ +
+
    let traceSun = {
+      x: [sun.x],
+      y: [sun.y],
+      mode: "markers",
+      marker: {
+        symbol: "star",
+        size: 10,
+        color: "gold",
+      },
+      name: "Sun",
+    };
+
+    const traceEarth = {
+      x: earth.x,
+      y: earth.y,
+      mode: "lines",
+      line: {
+        color: "white"
+      },
+      name: "Earth",
+    };
+
+    const traceOrbit = {
+      x: circle.x,
+      y:circle.y,
+      mode: "lines",
+      line: {
+        color: "crimson",
+        dash: "dash"
+      },
+      name: "1 AU Circle",
+    };
+
+    const earthSunLayout = {
+      title: "Earth-Sun Orbit",
+      xaxis: {
+        title: "x [AU]",
+        range: [-1.1,1.1],
+        showgrid: true,
+        gridcolor: "rgba(255,255,255,0.5)",
+        gridwidth: 1,
+        zeroline: true,
+        tickmode: "auto",
+        nticks: 5,
+      },
+      yaxis: {
+        title: "y [AU]",
+        range: [-1.1,1.1],
+        showgrid: true,
+        gridcolor: "rgba(255,255,255,0.5)",
+        gridwidth: 1,
+        zeroline: false,
+        tickmode: "auto",
+        nticks: 5,
+      },
+      margin: {
+        l: 50,
+        r: 50,
+        b: 50,
+        t: 50,
+        pad: 4,
+      },
+      paper_bgcolor: "black",
+      plot_bgcolor: "black",
+    };
+    Plotly.newPlot("plot",[traceSun,traceEarth,traceOrbit], earthSunLayout);
+
+

Figure of 8 orbit

@@ -240,78 +241,78 @@ Next, the function dR takes the position r and velocity v of Earth as input and

The code for this simulation is very similar to the Earth-Sun orbit simulation, except that we now have three bodies instead of two. We also use a different set of initial conditions to generate the figure of 8 orbit.

-
function deltaR(coords, masses, nBodies, G) {
-    let x = coords[0];
-    let y = coords[1];
-    let vx = coords[2];
-    let vy = coords[3];
-
-    let delta = math.clone(coords);
-
-    for (let n = 0; n < nBodies; n++) {
-        let xn = x[n];
-        let yn = y[n];
-        let deltaVx = 0.0;
-        let deltaVy = 0.0;
-
-        for (let i = 0; i < nBodies; i++) {
-            if (i !== n) {
-                let sep = Math.sqrt(Math.pow(xn - x[i], 2) + Math.pow(yn - y[i], 2)); // Euclidean distance
-                deltaVx -= G * masses[i] * (xn - x[i]) / Math.pow(sep, 3);
-                deltaVy -= G * masses[i] * (yn - y[i]) / Math.pow(sep, 3);
-            }
-        }
-
-        delta[2][n] = deltaVx;
-        delta[3][n] = deltaVy;
-    }
-
-    delta[0] = vx;
-    delta[1] = vy;
-
-    return delta;
-}
-
-function step(coords, masses, deltaT, nBodies = 3, G = 6.67408313131313e-11) {
-    let k1 = math.multiply(deltaT, deltaR(coords, masses, nBodies, G));
-    let k2 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k1, 0.5)), masses, nBodies, G));
-    let k3 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k2, 0.5)), masses, nBodies, G));
-    let k4 = math.multiply(deltaT, deltaR(math.add(coords, k3), masses, nBodies, G));
-
-    coords = math.add(coords, math.multiply(math.add(k1, math.multiply(2.0, k2), math.multiply(2.0, k3), k4), 1/6));
-
-    return coords;
-}
-
-    // Initial conditions setup
-    let M = [1, 1, 1];
-    let x = [-0.97000436, 0.0, 0.97000436];
-    let y = [0.24208753, 0.0, -0.24208753];
-    let vx = [0.4662036850, -0.933240737, 0.4662036850];
-    let vy = [0.4323657300, -0.86473146, 0.4323657300];
-    let Ei = -1 / Math.sqrt(Math.pow(2 * 0.97000436, 2) + Math.pow(2 * 0.24208753, 2)) - 2 / Math.sqrt(Math.pow(0.97000436, 2) + Math.pow(0.24208753, 2)) + 0.5 * (math.sum(math.add(math.dotPow(vx, 2), math.dotPow(vy, 2))));
-
-    function linspace(start, stop, num) {
-        const step = (stop - start) / (num - 1);
-        return Array.from({length: num}, (_, i) => start + (step * i));
-    }
-
-    let coords = [x, y, vx, vy];
-    const time = linspace(0, 6.3259, 1001);
-    let deltaT = time[1] - time[0];
-
-    let X = math.zeros(3, time.length).toArray();
-    let Y = math.zeros(3, time.length).toArray();
-    let VX = math.zeros(3, time.length).toArray();
-    let VY = math.zeros(3, time.length).toArray();
-
-    for (let i = 0; i < time.length; i++) {
-        coords = step(coords, M, deltaT, 3, 1);
-        X.forEach((_, idx) => X[idx][i] = coords[0][idx]);
-        Y.forEach((_, idx) => Y[idx][i] = coords[1][idx]);
-        VX.forEach((_, idx) => VX[idx][i] = coords[2][idx]);
-        VY.forEach((_, idx) => VY[idx][i] = coords[3][idx]);
-    }
+
function deltaR(coords, masses, nBodies, G) {
+    let x = coords[0];
+    let y = coords[1];
+    let vx = coords[2];
+    let vy = coords[3];
+
+    let delta = math.clone(coords);
+
+    for (let n = 0; n < nBodies; n++) {
+        let xn = x[n];
+        let yn = y[n];
+        let deltaVx = 0.0;
+        let deltaVy = 0.0;
+
+        for (let i = 0; i < nBodies; i++) {
+            if (i !== n) {
+                let sep = Math.sqrt(Math.pow(xn - x[i], 2) + Math.pow(yn - y[i], 2)); // Euclidean distance
+                deltaVx -= G * masses[i] * (xn - x[i]) / Math.pow(sep, 3);
+                deltaVy -= G * masses[i] * (yn - y[i]) / Math.pow(sep, 3);
+            }
+        }
+
+        delta[2][n] = deltaVx;
+        delta[3][n] = deltaVy;
+    }
+
+    delta[0] = vx;
+    delta[1] = vy;
+
+    return delta;
+}
+
+function step(coords, masses, deltaT, nBodies = 3, G = 6.67408313131313e-11) {
+    let k1 = math.multiply(deltaT, deltaR(coords, masses, nBodies, G));
+    let k2 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k1, 0.5)), masses, nBodies, G));
+    let k3 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k2, 0.5)), masses, nBodies, G));
+    let k4 = math.multiply(deltaT, deltaR(math.add(coords, k3), masses, nBodies, G));
+
+    coords = math.add(coords, math.multiply(math.add(k1, math.multiply(2.0, k2), math.multiply(2.0, k3), k4), 1/6));
+
+    return coords;
+}
+
+    // Initial conditions setup
+    let M = [1, 1, 1];
+    let x = [-0.97000436, 0.0, 0.97000436];
+    let y = [0.24208753, 0.0, -0.24208753];
+    let vx = [0.4662036850, -0.933240737, 0.4662036850];
+    let vy = [0.4323657300, -0.86473146, 0.4323657300];
+    let Ei = -1 / Math.sqrt(Math.pow(2 * 0.97000436, 2) + Math.pow(2 * 0.24208753, 2)) - 2 / Math.sqrt(Math.pow(0.97000436, 2) + Math.pow(0.24208753, 2)) + 0.5 * (math.sum(math.add(math.dotPow(vx, 2), math.dotPow(vy, 2))));
+
+    function linspace(start, stop, num) {
+        const step = (stop - start) / (num - 1);
+        return Array.from({length: num}, (_, i) => start + (step * i));
+    }
+
+    let coords = [x, y, vx, vy];
+    const time = linspace(0, 6.3259, 1001);
+    let deltaT = time[1] - time[0];
+
+    let X = math.zeros(3, time.length).toArray();
+    let Y = math.zeros(3, time.length).toArray();
+    let VX = math.zeros(3, time.length).toArray();
+    let VY = math.zeros(3, time.length).toArray();
+
+    for (let i = 0; i < time.length; i++) {
+        coords = step(coords, M, deltaT, 3, 1);
+        X.forEach((_, idx) => X[idx][i] = coords[0][idx]);
+        Y.forEach((_, idx) => Y[idx][i] = coords[1][idx]);
+        VX.forEach((_, idx) => VX[idx][i] = coords[2][idx]);
+        VY.forEach((_, idx) => VY[idx][i] = coords[3][idx]);
+    }
@@ -329,365 +330,369 @@ Next, the function dR takes the position r and velocity v of Earth as input and

Animation?

-

Now that we have time-series data, we need to animate it. We can use Plotly's animate function, as this does not force a full re-render, saving us some precious GPU and CPU cycles when we are trying to run this in the browser itself -

- 
function plotClassicFunc() {
-      var tailLength = 1;
-      if (plotIndex < tailLength) {
-      tailLength = 0;
-      } else if (plotIndex > time.length) {
-      plotIndex = 0;
-      } else {
-        tailLength = plotIndex - tailLength;
-      }
- 

-      var currentIndex = plotIndex;
- 

-     try {
-         time[currentIndex].toFixed(3);
-      } catch (e) {
-        currentIndex = 0;
-      }
- 

-       const data = [
-        {
-            x: X[0].slice(tailLength, currentIndex),
-            y: Y[0].slice(tailLength, currentIndex),
-            mode: 'lines+markers',
-            marker: {
-                symbol: 'star',
-                size: 8,
-                line: { width: 0 },
-            },
-            line: {
-                width: 2,
-            },
-            name: '',
-        },
-        {
-            x: X[1].slice(tailLength, currentIndex),
-            y: Y[1].slice(tailLength, currentIndex),
-            mode: 'lines+markers',
-            marker: {
-                symbol: 'star',
-                size: 8,
-                line: { width: 0 },
-            },
-            line: {
-                width: 2,
-            },
-            name: '',
-        },
-        {
-            x: X[2].slice(tailLength, currentIndex),
-            y: Y[2].slice(tailLength, currentIndex),
-            mode: 'lines+markers',
-            marker: {
-                symbol: 'star',
-                size: 8,
-                line: { width: 0 },
-            },
-            line: {
-                width: 2,
-            },
-            name: '',
-        },
-    ];
- 

-    // width: 1000, height: 400
-    const layout = {
-        title: '∞ Three-Body Problem: t = ' + time[currentIndex].toFixed(3),
-        xaxis: {
-            title: 'x',
-            range: [-1.1,1.1]
-        },
-        yaxis: {
-            title: 'y',
-            scaleanchor: 'x',
-            scaleratio: 1,
-            range: [-0.5,0.5]
-        },
-        plotbgcolor: 'black',
-        paperbgcolor: 'black',
-        font: {
-            color: 'white',
-        },
-    };
- 

-    try {
-    Plotly.animate("plot", {
-        data: data, layout: layout
-      }, {
-        staticPlot: true,
-        transition: {
-          duration: 0,
-        },
-        frame: {
-          duration: 0,
-          redraw: false,
-        }
-      });
-      } catch (err) {
-        Plotly.newPlot('plot', data, layout);
-      }
- 

- 

-    plotIndex += delay;
-    if (plotClassic===true) {
-      try {
-        requestAnimationFrame(plotClassicFunc);
-        }
-      catch (err) {
-        console.log(err)
-      }
-    }
- 

-    }
-
-

+

Now that we have time-series data, we need to animate it. We can use Plotly's animate function, as this does not force a full re-render, saving us some precious GPU and CPU cycles when we are trying to run this in the browser itself

+ +
+
    function plotClassicFunc() {
+      var tailLength = 1;
+      if (plotIndex < tailLength) {
+      tailLength = 0;
+      } else if (plotIndex > time.length) {
+      plotIndex = 0;
+      } else {
+        tailLength = plotIndex - tailLength;
+      }
+
+      var currentIndex = plotIndex;
+
+     try {
+         time[currentIndex].toFixed(3);
+      } catch (e) {
+        currentIndex = 0;
+      }
+
+       const data = [
+        {
+            x: X[0].slice(tailLength, currentIndex),
+            y: Y[0].slice(tailLength, currentIndex),
+            mode: 'lines+markers',
+            marker: {
+                symbol: 'star',
+                size: 8,
+                line: { width: 0 },
+            },
+            line: {
+                width: 2,
+            },
+            name: '',
+        },
+        {
+            x: X[1].slice(tailLength, currentIndex),
+            y: Y[1].slice(tailLength, currentIndex),
+            mode: 'lines+markers',
+            marker: {
+                symbol: 'star',
+                size: 8,
+                line: { width: 0 },
+            },
+            line: {
+                width: 2,
+            },
+            name: '',
+        },
+        {
+            x: X[2].slice(tailLength, currentIndex),
+            y: Y[2].slice(tailLength, currentIndex),
+            mode: 'lines+markers',
+            marker: {
+                symbol: 'star',
+                size: 8,
+                line: { width: 0 },
+            },
+            line: {
+                width: 2,
+            },
+            name: '',
+        },
+    ];
+
+    // width: 1000, height: 400
+    const layout = {
+        title: '∞ Three-Body Problem: t = ' + time[currentIndex].toFixed(3),
+        xaxis: {
+            title: 'x',
+            range: [-1.1,1.1]
+        },
+        yaxis: {
+            title: 'y',
+            scaleanchor: 'x',
+            scaleratio: 1,
+            range: [-0.5,0.5]
+        },
+        plot_bgcolor: 'black',
+        paper_bgcolor: 'black',
+        font: {
+            color: 'white',
+        },
+    };
+
+    try {
+    Plotly.animate("plot", {
+        data: data, layout: layout
+      }, {
+        staticPlot: true,
+        transition: {
+          duration: 0,
+        },
+        frame: {
+          duration: 0,
+          redraw: false,
+        }
+      });
+      } catch (err) {
+        Plotly.newPlot('plot', data, layout);
+      }
+
+
+    plotIndex += delay;
+    if (plotClassic===true) {
+      try {
+        requestAnimationFrame(plotClassicFunc);
+        }
+      catch (err) {
+        console.log(err)
+      }
+    }
+
+    }
+
+

"General" N-Body Solver

Show me the code!

-
function step(coords, masses, deltaT, nBodies = 3, G = 6.67408313131313e-11) {
-    let k1 = math.multiply(deltaT, deltaR(coords, masses, nBodies, G));
-    let k2 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k1, 0.5)), masses, nBodies, G));
-    let k3 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k2, 0.5)), masses, nBodies, G));
-    let k4 = math.multiply(deltaT, deltaR(math.add(coords, k3), masses, nBodies, G));
-
-    coords = math.add(coords, math.multiply(math.add(k1, math.multiply(2.0, k2), math.multiply(2.0, k3), k4), 1/6));
-
-    return coords;
-}
-
-function detectCollisionsEscape(coords, deltaT, maxSep) {
-  const [x, y, vx, vy] = coords;
-  const V = vx.map((v, i) => Math.sqrt(v ** 2 + vy[i] ** 2));
-  const R = V.map(v => v * deltaT);
-  let collision = false, collisionInds = null, escape = false, escapeInd = null;
-
-  for (let n = 0; n < x.length; n++) {
-      const rn = R[n], xn = x[n], yn = y[n];
-      for (let i = 0; i < x.length; i++) {
-          if (i !== n) {
-              const minSep = rn + R[i];
-              const sep = Math.sqrt((xn - x[i]) ** 2 + (yn - y[i]) ** 2);
-              if (sep < minSep) {
-                  collision = true;
-                  collisionInds = [n, i];
-              } else if (sep > maxSep) {
-                  escape = true;
-                  escapeInd = n;
-                  return [collision, collisionInds, escape, escapeInd];
-              }
-          }
-      }
-  }
-  return [collision, collisionInds, escape, escapeInd];
-}
-
-function nBodyStep(coords, masses, deltaT, maxSep, nBodies, G = 6.67408313131313e-11) { // Similar to our step function before, but keeping track of collisions
-  coords = step(coords, masses, deltaT, nBodies, G); // Update the positions as we did before
-  //console.log(detectCollisionsEscape(coords, deltaT, maxSep));
-  let [collision, collisionInds, escape, escapeInd] = detectCollisionsEscape(coords, deltaT, maxSep); // Detect collisions/escapes
-
-
-  if (collision) { // Do inelastic collision and delete extra body (2 -> 1)
-    const [i1, i2] = collisionInds;
-      const [x1, x2] = [coords[0][i1], coords[0][i2]];
-      const [y1, y2] = [coords[1][i1], coords[1][i2]];
-      const [vx1, vx2] = [coords[2][i1], coords[2][i2]];
-      const [vy1, vy2] = [coords[3][i1], coords[3][i2]];
-      const [px1, px2] = [masses[i1] * vx1, masses[i2] * vx2];
-      const [py1, py2] = [masses[i1] * vy1, masses[i2] * vy2];
-      const px = px1 + px2;
-      const py = py1 + py2;
-      const newM = masses[i1] + masses[i2];
-      const vfx = px / newM;
-      const vfy = py / newM;
-      coords[0][i1] = (x1 * masses[i1] + x2 * masses[i2]) / (masses[i1] + masses[i2]); // Center of mass
-      coords[1][i1] = (y1 * masses[i1] + y2 * masses[i2]) / (masses[i1] + masses[i2]);
-      coords[2][i1] = vfx;
-      coords[3][i1] = vfy;
-      coords[0].splice(i2, 1);
-      coords[1].splice(i2, 1);
-      coords[2].splice(i2, 1);
-      coords[3].splice(i2, 1);
-      masses[i1] = newM;
-      masses.splice(i2, 1);
-      nBodies--;
-  }
-  // Could also implement condition for escape where we stop calculating forces but I'm too lazy for now
-  return [coords, masses, nBodies, collision, collisionInds, escape, escapeInd];
-}
-
-function uniform(min, max) {
-  return Math.random() * (max - min) + min;
-}
-
-function deepCopyCoordsArray(arr) {
-  return arr.map(innerArr => innerArr.slice());
-}
-
-function genNBodyResults(nBodies, tStop, nTPts, nBodiesStop = 10, G = 6.67408313131313e-11) {
-
-  var btn = document.getElementById("startSim3");
-  // Set button text to Solving
-  var prevText = btn.innerHTML;
-  btn.innerHTML = "Solving...";
-
-  let coords = [Array(nBodies).fill(0), Array(nBodies).fill(0), Array(nBodies).fill(0), Array(nBodies).fill(0)];
-  const Mstar = 2e30;
-  const Mp = Mstar / 1e4;
-
-  for (let i = 0; i < nBodies; i++) { // Initialize coordinates on ~Keplerian orbits
-      let accept = false;
-      let r = null;
-      while (!accept) { // Prevent a particle from spawning within 0.2 AU too close to "star"
-          r = Math.random() * 2 * 1.5e11; // Say radius of 2 AU
-          if (r / 1.5e11 > 0.2) {
-              accept = true;
-          }
-      }
-      const theta = uniform(0, 2 * Math.PI);
-      const x = r * Math.cos(theta);
-      const y = r * Math.sin(theta);
-      const v = Math.sqrt(G * Mstar / r);
-      const perturbedV = v + v / 1000 * uniform(-1, 1); // Perturb the velocities ever so slightly
-      const vTheta = Math.atan2(y, x);
-      coords[0][i] = x;
-      coords[1][i] = y;
-      coords[2][i] = -perturbedV * Math.sin(vTheta);
-      coords[3][i] = perturbedV * Math.cos(vTheta);
-  }
-
-  //console.log('Initial coords:', coords);
-
-
-  let masses = Array(nBodies).fill(Mp); // Initialize masses
-  masses[0] = Mstar; // Make index one special as the central star
-  coords[0][0] = 0;
-  coords[1][0] = 0;
-  coords[2][0] = 0;
-  coords[3][0] = 0; // Initialize central star at origin with no velocity
-  const yearSec = 365 * 24 * 3600;
-  const time = Array.from({ length: nTPts }, (_, i) => i * tStop / (nTPts - 1) * yearSec); // Years -> s
-  let t = time[0];
-  const deltaT = time[1] - time[0];
-  let tInd = 0;
-  const coordsRecord = [deepCopyCoordsArray(coords)];
-  const massRecord = [masses.slice()]; // Initialize records with initial conditions
-
-
-  while (tInd < nTPts && nBodies > nBodiesStop) {
-    //console.log('Initial coords:', coords);
-    [coords, masses, nBodies] = nBodyStep(coords, masses, deltaT, 10 * 1.5e11, nBodies, G); // Update
-    coordsRecord.push(deepCopyCoordsArray(coords));
-    massRecord.push(masses.slice()); // Add to records
-    tInd++;
-    t = time[tInd];
-    //console.log(`currently at t = ${(t / yearSec).toFixed(2)} years\r`);
-  }
-  console.log(`final time = ${time[tInd] / yearSec} years with ${nBodies} bodies remaining`);
-
-  // Set button text to Start Simulation
-  btn.innerHTML = prevText;
-
-  return [coordsRecord, massRecord, time.slice(0, tInd + 1)];
-}
-
-
- var [coordsRecordR, _, tR] = genNBodyResults(256,1,1001);
-    //console.log(coordsRecordR);
-    const yearSec = 365 * 24 * 3600;
-
-    function createFrame(coordsR) {
-      if (!coordsR || !coordsR[0] || !coordsR[1]) {
-          return [];
-      }
-
-      const traceCentralStar = {
-          x: [coordsR[0][0] / 1.5e11],
-          y: [coordsR[1][0] / 1.5e11],
-          mode: 'markers',
-          type: 'scatter',
-          name: 'Central star',
-          marker: { color: 'gold', symbol: 'star', size: 10 },
-      };
-
-      const xCoords = coordsR[0].slice(1).map(x => x / 1.5e11);
-      const yCoords = coordsR[1].slice(1).map(y => y / 1.5e11);
-
-      const traceOtherBodies = {
-          x: xCoords,
-          y: yCoords,
-          mode: 'markers',
-          type: 'scatter',
-          name: '',
-          marker: { color: 'dodgerblue', symbol: 'circle', size: 2 },
-      };
-
-      return [traceCentralStar, traceOtherBodies];
-  }
-
-
-  function createLayout(i) {
-    return {
-        title: {
-            text: `N-Body Problem: t = ${Number(t[i] / yearSec).toFixed(3)} years`,
-            x: 0.03,
-            y: 0.97,
-            xanchor: 'left',
-            yanchor: 'top',
-            font: { size: 14 },
-        },
-        xaxis: { title: 'x [AU]', range: [-2.1, 2.1] },
-        yaxis: { title: 'y [AU]', range: [-2.1, 2.1], scaleanchor: 'x', scaleratio: 1 },
-        showlegend: false,
-        margin: { l: 60, r: 40, t: 40, b: 40 },
-        width: 800,
-        height: 800,
-        plot_bgcolor: 'black',
-    };
-}
-
-  function animateNBodyProblem() {
-  const nFrames = tR.length;
-
-  for (let i = 0; i < nFrames; i++) {
-      const frameData = createFrame(coordsRecordR[i]);
-      const layout = createLayout(i);
-      //Plotly.newPlot(plotDiv, frameData, layout);
-      try {
-        Plotly.animate("plot", {
-        data: frameData, layout: layout
-      }, {
-        staticPlot: true,
-        transition: {
-          duration: 0,
-        },
-        frame: {
-          duration: 0,
-          redraw: false,
-        }
-      });
-    } catch (err) {
-      Plotly.newPlot('plot', frameData, layout);
-    }
-  }
-}
-
-animateNBodyProblem();
+
function step(coords, masses, deltaT, nBodies = 3, G = 6.67408313131313e-11) {
+    let k1 = math.multiply(deltaT, deltaR(coords, masses, nBodies, G));
+    let k2 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k1, 0.5)), masses, nBodies, G));
+    let k3 = math.multiply(deltaT, deltaR(math.add(coords, math.multiply(k2, 0.5)), masses, nBodies, G));
+    let k4 = math.multiply(deltaT, deltaR(math.add(coords, k3), masses, nBodies, G));
+
+    coords = math.add(coords, math.multiply(math.add(k1, math.multiply(2.0, k2), math.multiply(2.0, k3), k4), 1/6));
+
+    return coords;
+}
+
+function detectCollisionsEscape(coords, deltaT, maxSep) {
+  const [x, y, vx, vy] = coords;
+  const V = vx.map((v, i) => Math.sqrt(v ** 2 + vy[i] ** 2));
+  const R = V.map(v => v * deltaT);
+  let collision = false, collisionInds = null, escape = false, escapeInd = null;
+
+  for (let n = 0; n < x.length; n++) {
+      const rn = R[n], xn = x[n], yn = y[n];
+      for (let i = 0; i < x.length; i++) {
+          if (i !== n) {
+              const minSep = rn + R[i];
+              const sep = Math.sqrt((xn - x[i]) ** 2 + (yn - y[i]) ** 2);
+              if (sep < minSep) {
+                  collision = true;
+                  collisionInds = [n, i];
+              } else if (sep > maxSep) {
+                  escape = true;
+                  escapeInd = n;
+                  return [collision, collisionInds, escape, escapeInd];
+              }
+          }
+      }
+  }
+  return [collision, collisionInds, escape, escapeInd];
+}
+
+function nBodyStep(coords, masses, deltaT, maxSep, nBodies, G = 6.67408313131313e-11) { // Similar to our step function before, but keeping track of collisions
+  coords = step(coords, masses, deltaT, nBodies, G); // Update the positions as we did before
+  //console.log(detectCollisionsEscape(coords, deltaT, maxSep));
+  let [collision, collisionInds, escape, escapeInd] = detectCollisionsEscape(coords, deltaT, maxSep); // Detect collisions/escapes
+
+
+  if (collision) { // Do inelastic collision and delete extra body (2 -> 1)
+    const [i1, i2] = collisionInds;
+      const [x1, x2] = [coords[0][i1], coords[0][i2]];
+      const [y1, y2] = [coords[1][i1], coords[1][i2]];
+      const [vx1, vx2] = [coords[2][i1], coords[2][i2]];
+      const [vy1, vy2] = [coords[3][i1], coords[3][i2]];
+      const [px1, px2] = [masses[i1] * vx1, masses[i2] * vx2];
+      const [py1, py2] = [masses[i1] * vy1, masses[i2] * vy2];
+      const px = px1 + px2;
+      const py = py1 + py2;
+      const newM = masses[i1] + masses[i2];
+      const vfx = px / newM;
+      const vfy = py / newM;
+      coords[0][i1] = (x1 * masses[i1] + x2 * masses[i2]) / (masses[i1] + masses[i2]); // Center of mass
+      coords[1][i1] = (y1 * masses[i1] + y2 * masses[i2]) / (masses[i1] + masses[i2]);
+      coords[2][i1] = vfx;
+      coords[3][i1] = vfy;
+      coords[0].splice(i2, 1);
+      coords[1].splice(i2, 1);
+      coords[2].splice(i2, 1);
+      coords[3].splice(i2, 1);
+      masses[i1] = newM;
+      masses.splice(i2, 1);
+      nBodies--;
+  }
+  // Could also implement condition for escape where we stop calculating forces but I'm too lazy for now
+  return [coords, masses, nBodies, collision, collisionInds, escape, escapeInd];
+}
+
+function uniform(min, max) {
+  return Math.random() * (max - min) + min;
+}
+
+function deepCopyCoordsArray(arr) {
+  return arr.map(innerArr => innerArr.slice());
+}
+
+function genNBodyResults(nBodies, tStop, nTPts, nBodiesStop = 10, G = 6.67408313131313e-11) {
+
+  var btn = document.getElementById("startSim3");
+  // Set button text to Solving
+  var prevText = btn.innerHTML;
+  btn.innerHTML = "Solving...";
+
+  let coords = [Array(nBodies).fill(0), Array(nBodies).fill(0), Array(nBodies).fill(0), Array(nBodies).fill(0)];
+  const Mstar = 2e30;
+  const Mp = Mstar / 1e4;
+
+  for (let i = 0; i < nBodies; i++) { // Initialize coordinates on ~Keplerian orbits
+      let accept = false;
+      let r = null;
+      while (!accept) { // Prevent a particle from spawning within 0.2 AU too close to "star"
+          r = Math.random() * 2 * 1.5e11; // Say radius of 2 AU
+          if (r / 1.5e11 > 0.2) {
+              accept = true;
+          }
+      }
+      const theta = uniform(0, 2 * Math.PI);
+      const x = r * Math.cos(theta);
+      const y = r * Math.sin(theta);
+      const v = Math.sqrt(G * Mstar / r);
+      const perturbedV = v + v / 1000 * uniform(-1, 1); // Perturb the velocities ever so slightly
+      const vTheta = Math.atan2(y, x);
+      coords[0][i] = x;
+      coords[1][i] = y;
+      coords[2][i] = -perturbedV * Math.sin(vTheta);
+      coords[3][i] = perturbedV * Math.cos(vTheta);
+  }
+
+  //console.log('Initial coords:', coords);
+
+
+  let masses = Array(nBodies).fill(Mp); // Initialize masses
+  masses[0] = Mstar; // Make index one special as the central star
+  coords[0][0] = 0;
+  coords[1][0] = 0;
+  coords[2][0] = 0;
+  coords[3][0] = 0; // Initialize central star at origin with no velocity
+  const yearSec = 365 * 24 * 3600;
+  const time = Array.from({ length: nTPts }, (_, i) => i * tStop / (nTPts - 1) * yearSec); // Years -> s
+  let t = time[0];
+  const deltaT = time[1] - time[0];
+  let tInd = 0;
+  const coordsRecord = [deepCopyCoordsArray(coords)];
+  const massRecord = [masses.slice()]; // Initialize records with initial conditions
+
+
+  while (tInd < nTPts && nBodies > nBodiesStop) {
+    //console.log('Initial coords:', coords);
+    [coords, masses, nBodies] = nBodyStep(coords, masses, deltaT, 10 * 1.5e11, nBodies, G); // Update
+    coordsRecord.push(deepCopyCoordsArray(coords));
+    massRecord.push(masses.slice()); // Add to records
+    tInd++;
+    t = time[tInd];
+    //console.log(`currently at t = ${(t / yearSec).toFixed(2)} years\r`);
+  }
+  console.log(`final time = ${time[tInd] / yearSec} years with ${nBodies} bodies remaining`);
+
+  // Set button text to Start Simulation
+  btn.innerHTML = prevText;
+
+  return [coordsRecord, massRecord, time.slice(0, tInd + 1)];
+}
+
+
+ var [coordsRecordR, _, tR] = genNBodyResults(256,1,1001);
+    //console.log(coordsRecordR);
+    const yearSec = 365 * 24 * 3600;
+
+    function createFrame(coordsR) {
+      if (!coordsR || !coordsR[0] || !coordsR[1]) {
+          return [];
+      }
+
+      const traceCentralStar = {
+          x: [coordsR[0][0] / 1.5e11],
+          y: [coordsR[1][0] / 1.5e11],
+          mode: 'markers',
+          type: 'scatter',
+          name: 'Central star',
+          marker: { color: 'gold', symbol: 'star', size: 10 },
+      };
+
+      const xCoords = coordsR[0].slice(1).map(x => x / 1.5e11);
+      const yCoords = coordsR[1].slice(1).map(y => y / 1.5e11);
+
+      const traceOtherBodies = {
+          x: xCoords,
+          y: yCoords,
+          mode: 'markers',
+          type: 'scatter',
+          name: '',
+          marker: { color: 'dodgerblue', symbol: 'circle', size: 2 },
+      };
+
+      return [traceCentralStar, traceOtherBodies];
+  }
+
+
+  function createLayout(i) {
+    return {
+        title: {
+            text: `N-Body Problem: t = ${Number(t[i] / yearSec).toFixed(3)} years`,
+            x: 0.03,
+            y: 0.97,
+            xanchor: 'left',
+            yanchor: 'top',
+            font: { size: 14 },
+        },
+        xaxis: { title: 'x [AU]', range: [-2.1, 2.1] },
+        yaxis: { title: 'y [AU]', range: [-2.1, 2.1], scaleanchor: 'x', scaleratio: 1 },
+        showlegend: false,
+        margin: { l: 60, r: 40, t: 40, b: 40 },
+        width: 800,
+        height: 800,
+        plot_bgcolor: 'black',
+    };
+}
+
+  function animateNBodyProblem() {
+  const nFrames = tR.length;
+
+  for (let i = 0; i < nFrames; i++) {
+      const frameData = createFrame(coordsRecordR[i]);
+      const layout = createLayout(i);
+      //Plotly.newPlot(plotDiv, frameData, layout);
+      try {
+        Plotly.animate("plot", {
+        data: frameData, layout: layout
+      }, {
+        staticPlot: true,
+        transition: {
+          duration: 0,
+        },
+        frame: {
+          duration: 0,
+          redraw: false,
+        }
+      });
+    } catch (err) {
+      Plotly.newPlot('plot', frameData, layout);
+    }
+  }
+}
+
+animateNBodyProblem();

Playground

+ + +
- + +

- +

+