Published
Edited
Jun 14, 2022
Importers
Insert cell
Insert cell
Insert cell
arrowValues = [
0, // top left x, relative to refX
0, // top left y, relative to refY
11, // tip of the arrow x
5, // tip of the arrow Y (relative to refY)
0, // bottom of the arrow x
10 // bottom of the arrow y (full arrow width)
]
Insert cell
{
const svg = d3.create('svg').attr('viewBox', [0, 0, width, 25]);
// Inspired by: https://observablehq.com/@harrystevens/linear-algebra?ui=classic

// do a .attr("marker-end", "url(#arrow)")
let av = arrowValues;
svg
.append('defs')
.append('marker')
.attr('id', 'arrow')
.attr('viewBox', [0, 0, 11, 11])
.attr('refX', "10") // how far "back" from end of the path to start
.attr('refY', "5") // how far "up" from the end of the path to start
.attr('markerWidth', "6")
.attr('markerHeight', "6")
.attr('orient', 'auto-start-reverse')
.append('path')
.attr('d', `M ${av[0]} ${av[1]} L ${av[2]} ${av[3]} L ${av[4]} ${av[5]} z`)
.attr('fill', 'black')
;

svg
.append('path')
.attr('d', d3.line()([[10, 10], [100, 10]]))
.attr('stroke', 'black')
.attr('marker-end', 'url(#arrow)')
.attr('fill', 'none');

return svg.node();
}
Insert cell
Insert cell
Insert cell
Insert cell
{

const width = 600, height = 300,
x = d3.scaleLinear([-10, 30], [0, width]),
y = d3.scaleLinear([-5, 15], [height, 0]);

const drawVec = (v) => {
return d3.line()(v.map((p) => [x(p[0]), y(p[1])]));
};

const svg = d3.create("svg")
.style("overflow", "visible")
.attr("width", width)
.attr("height", height)
.call(d => drawAxes(d, x, y));

const v1 = [[0, 0], [12, 6]];

addArrowHead(svg);

// Draw the line
svg
.append('path')
.attr('d', d3.line()(v1.map((p) => [x(p[0]), y(p[1])])))
.attr('stroke', 'black')
.attr('marker-end', 'url(#arrow)')
.attr('fill', 'none');

return svg.node();
}
Insert cell
drawAxes = (sel, x, y, options) => {
const width = x.range()[1], height = y.range()[0];
const xgrid = sel.append("g")
.attr("class", "grid x")
.call(d3
.axisBottom(x)
.ticks(20) // Roughly how many ticks you want TODO: might want to calculuate this based on x/y range size to keep it square?
.tickSize(height) // How long to make the ticks, e.g. height if you want them to be a full grid
);
xgrid.selectAll("line").style("stroke", "#eee");
xgrid.selectAll("text").remove();
xgrid.select(".domain").remove();
const ygrid = sel.append("g")
.attr("class", "grid y")
.call(d3.axisLeft(y).tickSize(width))
.attr("transform", `translate(${width})`);
ygrid.selectAll("line").style("stroke", "#eee");
ygrid.selectAll("text").remove();
ygrid.select(".domain").remove();

const xaxis = d3.axisBottom(x);
const yaxis = d3.axisLeft(y);
if (options && options.tickFormat) {
xaxis.tickFormat(options.tickFormat);
yaxis.tickFormat(options.tickFormat);
}

// Axes bold line
sel.append("g")
.attr("class", "axis x")
.attr("transform", `translate(0, ${y(0)})`) // TODO, maybe make this center around 0,0 rather than middle of the chart.
.call(xaxis);
sel.append("g")
.attr("class", "axis y")
.attr("transform", `translate(${x(0)})`)
.call(yaxis);

return sel;
}
Insert cell
Insert cell
createStandardGrid = (xrange, yrange, config) => {
xrange = xrange || [-20, 20];
yrange = yrange || [-10, 10];
config = config || {};
let width = config.width || 600;
let height = config.height || 300;

const x = d3.scaleLinear(xrange, [0, width]),
y = d3.scaleLinear(yrange, [height, 0]);

let svg = config.svg || d3.create("svg");
svg
//.style("overflow", "visible")
.attr("width", width)
.attr("height", height)
.call(d => drawAxes(d, x, y, config));

// const g = svg.append("g")
// .attr("cursor", "grab");

// svg.call(d3.zoom()
// .extent([[0, 0], [600, 300]])
// .scaleExtent([1, 8])
// .on("zoom", zoomed));

// function zoomed({transform}) {
// svg.attr("transform", transform);
// }

addArrowHead(svg, 'black');
addArrowHead(svg, 'blue');
addArrowHead(svg, 'green');
addArrowHead(svg, 'red');
addArrowHead(svg, 'gray');


return {svg, x, y}
}
Insert cell
Insert cell
{
let graph = createStandardGrid();

drawVector(graph, [[0, 6], [12, 0]], 'blue');

drawVector(graph, [[-5, -8], [-7, 8]], 'green');

drawVector(graph, [[18, 8], [-15, 2]], 'red');
return graph.svg.node();
}
Insert cell
Type JavaScript, then Shift-Enter. Ctrl-space for more options. Arrow ↑/↓ to switch modes.

Insert cell
## Plot a function
Insert cell
{
let graph = createStandardGrid();

let fn = (x) => 0.1*x**2;

let pts = [];

for (let i = -20; i <= 20; i += 0.1) {
pts.push([i, fn(i)]);
}

graph.svg.append("path")
.attr("d", d3.line()(pts.map(p => [graph.x(p[0]), graph.y(p[1])])))
.attr('stroke', 'black')
.attr('fill', 'none');

return graph.svg.node();
}
Insert cell
vectorAdd = (v1, v2) => {
return v1.map((v, i) => v+v2[i])
}
Insert cell
vectorScale = (v, k) => v.map(v0 => v0*k)
Insert cell
makeStateVec = (nodes) => {
return nodes.map(n => n.x)
.concat(nodes.map(n => n.y))
.concat(nodes.map(n => n.vx))
.concat(nodes.map(n => n.vy));
}
Insert cell
Insert cell
// Returns a function that will calculate new accelerations given an input vector of positions and velocities.
// It calculates them based on the positions and the gravitational forces
// It uses masses from the initial 'nodes' list passed in to calculate acceleration due to gravity.
gravityField = (nodes, G) => {
// nodes are the nodes exhibiting gravitational forces, generally celestial bodies

// the s input below is supposed to be the values affecting each other, but since we're just doing fields,
// s will always be the initial state.
G = G || 1;
const n = nodes.length;
// Borrowed from https://observablehq.com/@mcmcclur/gravitation-and-the-n-body-problem
let x_accel = (s, i, j) => {
if (i == j) { return 0 };
let m2 = nodes[j].m;

let x1 = s[i];
let y1 = s[i + n];
let x2 = s[j];
let y2 = s[j + n];
return (m2 * (x2 - x1)) / ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** (3 / 2);
};

let y_accel = (s, i, j) => {
if (i == j) { return 0 };

let m2 = nodes[j].m;
let x1 = s[i];
let y1 = s[i + n];
let x2 = s[j];
let y2 = s[j + n];
return (m2 * (y2 - y1)) / ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** (3 / 2);
};

let gr = (s, t) => {
return s
// These two copy the stored velocities in the state to return them. By lining it up this way,
// we can rely on generic vector addition in the Runge-Kutta function to compute the next steps
// approximations. So we line velocity up with postion, and acceleration up with velocity, then
// scale and add in RK.
.slice(2 * n, 3 * n) // copy vx's [0n: px, 1n: py, 2n: vx, 3n: vy]
.concat(s.slice(3 * n, 4 * n)) // copy vy's ^
.concat(d3.range(n).map(i => {
return G*d3.sum(d3.range(nodes.length).map((j) => x_accel(s, i, j)))
}))
.concat(d3.range(n).map(i => {
return G*d3.sum(d3.range(nodes.length).map((j) => y_accel(s, i, j)))
}));
};

return gr
}
Insert cell
circularizeOrbitsSimple = (nodes, G) => {
// Assumes the first node is the central body with mass.
// Also assumes they're on the x axis, but that probably shouldn't be a limit.

G = G || 1;

const n0 = nodes[0];
for (let i = 1; i < nodes.length; i += 1) {
let n = nodes[i];

// This assumes the mass of each lesser node is
// 10x smaller than the central node. If that's not
// the case, use the circularizeOrbits2B function
let v = Math.sqrt((G*n0.m)/(n.x - n0.x)) ;

n.vy = v;
}

return nodes;
}
Insert cell
circularizeOrbits2B = (nodes, G) => {
G = G || 1;

let n1 = nodes[0];
let n2 = nodes[1];

const calcCom = (p1, p2) => {
let Rx = p2.x - p1.x;
let m1 = p1.m;
let m2 = p2.m;

let r1x = -1*(m2*Rx) /(m1 + m2);
let r2x = (m1*Rx) /(m1 + m2);

let Ry = p2.y - p1.y;

let r1y = -1*(m2*Ry) /(m1 + m2);
let r2y = (m1*Ry) /(m1 + m2);
return {r1: {x: r1x, y: r1y}, r2: {x: r2x, y: r2y}};
}

const com = calcCom(nodes[0], nodes[1]);

const R = {
x: com.r2.x - com.r1.x,
y: com.r2.y - com.r1.y,
}

//return R;

// Setting the origin at the center of mass makes these calculations much easier.
let calcV = (m, r, R) => {
return {
// || 1 is a hack to avoid divide by 0.
// if the y vector for position is 0, then there's no y velocity.
// This is broken, I'm calculating the total velocity, I need to put it perpendicular
// I'm not sure where I got this formula...
vy: Math.sqrt((G*m*Math.abs(r.x))/(R.x*R.x || 1)),
vx: Math.sqrt((G*m*Math.abs(r.y))/(R.y*R.y || 1))

}
}
let v1 = calcV(n2.m, com.r1, R);
let v2 = calcV(n1.m, com.r2, R);

return [
{...n1, x: com.r1.x, y: com.r1.y, vx: v1.vx, vy: -v1.vy},
{...n2, x: com.r2.x, y: com.r2.y, vx: v2.vx, vy: v2.vy},
]
}
Insert cell
Insert cell
rk_step = (f, y, t, dt) => {
let mdt = dt/2;

let k1 = f(y, t);
let k2 = f(vectorAdd(y, vectorScale(k1, mdt)), t + mdt);
let k3 = f(vectorAdd(y, vectorScale(k2, mdt)), t + mdt);
let k4 = f(vectorAdd(y, vectorScale(k3, dt)), t + mdt);

return vectorAdd(
y, vectorScale(
vectorAdd(vectorAdd(vectorAdd(k1, vectorScale(k2, 2)), vectorScale(k3, 2)), k4)
, dt/6));
}
Insert cell
{
const fns = [(x) => x + 0.1, (x) => x * 2, (x) => x - 1];

// const fns = [];

const cur = 5;

const expected = (cur + .1)*2 -1;


const got = fns.reduce((x, f) => (f(x)), cur)
return {expected, got}
}
Insert cell
nBodyOrbits = (nodes, steps, dt, G, mutations) => {
const t0 = 0;

mutations = mutations || [];

dt = dt || 1;
G = G || 1;

let n = nodes.length;
let s0 = makeStateVec(nodes);
let F = gravityField(nodes, G);

let states = [s0];
let cur = s0;
for (let i = 0; i < steps; i++) {
cur = mutations.reduce((c, m) => m(c, i), cur);
cur = rk_step(F, cur, t0 + i * dt, dt);
states.push(cur);
}

return {states, nodes, F, dt}
}
Insert cell
Insert cell
visVivaEq = tex`v^2 = GM(\frac{2}{r} - \frac{1}{a})`
Insert cell
visViva = (mu, r, a) => {
return Math.sqrt(mu*((2/r) - 1/a))
}
Insert cell
calcOrbPeriod = (a, mu) => 2*Math.PI*Math.sqrt((a**3)/mu)
Insert cell
orbitTools = {
return {visViva, calcOrbPeriod, nBodyOrbits, circularizeOrbitsSimple, circularizeOrbits2B}
}
Insert cell
Insert cell
d3 = require('d3@7.3.0')
Insert cell
import {Scrubber} from "@mbostock/scrubber"
Insert cell
# Some links

- Principia, Kerbal Space Program mod to give N-body functionality: https://github.com/mockingbirdnest/Principia/wiki/A-guide-to-going-to-the-Mun-with-Principia

- Orbital Maneuver: https://en.wikipedia.org/wiki/Orbital_maneuver
- Interplanetary Transport Network : https://en.wikipedia.org/wiki/Interplanetary_Transport_Network
- https://observablehq.com/@mcmcclur/gravitation-and-the-n-body-problem
- https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods
-
Insert cell

Purpose-built for displays of data

Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
Learn more