diff options
Diffstat (limited to 'Resources/assets')
| -rw-r--r-- | Resources/assets/n-body/script.js | 598 | 
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(); + +  } +}  + | 
