add js file for sim

1 files changed, 598 insertions, 0 deletions

diff --git a/Resources/assets/n-body/script.js b/Resources/assets/n-body/script.js
new file mode 100644
index 0000000..e87d92c
--- /dev/null
+++ b/Resources/assets/n-body/script.js
@@ -0,0 +1,598 @@
+// Global variables
+var plotOrbit = true 
+var plotClassic = false
+var plotRandom = false
+var plotIndex = 0
+var delay = 10
+var tailLength = 1;
+var nBodies = 256;
+
+/*
+ * 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)) };
+
+/*
+ * Generic Functions 
+ * 
+ *
+ */
+
+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 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 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 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)];
+}
+
+
+function initCondGen(nBodies, vRange = [-7e3, 7e3], posRange = [-35, 35]) {
+    const m = Array.from({length: nBodies}, () => Math.random() * 1500 / 10);
+    
+    const rad = m.map(x => Math.pow(x, 0.8));
+    const minV = vRange[0], maxV = vRange[1];
+    const minPos = posRange[0], maxPos = posRange[1];
+    const posList = [];
+
+    function checkPos(randPos, n, posList, rad) {
+        for (let i = 0; i < n - 1; i++) {
+            const dist = Math.sqrt(Math.pow(posList[i][0] - randPos[0], 2) + Math.pow(posList[i][1] - randPos[1], 2));
+            if (dist * 1.5e11 < (rad[n] + rad[i])) {
+                return false;
+            }
+        }
+        return true;
+    }
+
+    function genPos(nBodies, posList, rad, minPos, maxPos) {
+        posList.push([Math.random() * (maxPos - minPos) + minPos, Math.random() * (maxPos - minPos) + minPos]);
+        for (let n = 1; n < nBodies; n++) {
+            let acceptPos = false;
+            while (acceptPos === false) {
+                const randPos = [Math.random() * (maxPos - minPos) + minPos, Math.random() * (maxPos - minPos) + minPos];
+                acceptPos = checkPos(randPos, n, posList, rad);
+                if (acceptPos === true) {
+                    posList.push(randPos);
+                }
+            }
+        }
+        return posList;
+    }
+
+    const pos = genPos(nBodies, posList, rad, minPos, maxPos).map(x => x.map(y => y * 1.5e11));
+    const coords = [new Array(nBodies).fill(0), new Array(nBodies).fill(0), new Array(nBodies).fill(0), new Array(nBodies).fill(0)];
+    const v = [];
+
+    for (let i = 0; i < nBodies; i++) {
+        coords[0][i] = pos[i][0];
+        coords[1][i] = pos[i][1];
+        const V = [Math.random() * (maxV - minV) + minV, Math.random() * (maxV - minV) + minV];
+        v.push(V);
+        coords[2][i] = V[0];
+        coords[3][i] = V[1];
+    }
+
+    return {m: m.map(x => x * 2e30), rad: rad.map(x => x * 7e8), coords: coords};
+}
+
+
+function calculateAndPlot() {
+  try {
+    Plotly.purge("plot");
+  } catch (e) {
+    console.log(e);
+  }
+
+  if (plotOrbit===true) {
+    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);
+  } else if (plotRandom==true) {
+
+    var [coordsRecordR, _, tR] = genNBodyResults(nBodies,1,nBodies*2);
+    //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`,//= ${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();
+
+
+  } else if (plotClassic==true) {
+    // 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 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)
+      }
+    }
+      
+    }
+
+    plotClassicFunc();
+
+  }
+} 
+