diff --git a/index.html b/index.html index cf15339..ff66bc9 100644 --- a/index.html +++ b/index.html @@ -10,12 +10,14 @@ body { padding: 0; margin: 0; background-color: #1b1b1b; + overflow-y: hidden; }
+ diff --git a/sketch.js b/sketch.js index 362e59f..1d3cb0a 100644 --- a/sketch.js +++ b/sketch.js @@ -1,16 +1,33 @@ -let charges = []; +let camera; +let red; + +let charges = []; +let faces = []; -let sphere_mode = 'circles'; let sphere_radius = 200; +const SURFACE_NONE = 0; +const SURFACE_CIRCLES = 1; +const SURFACE_EARTH = 2; + +let surface = SURFACE_CIRCLES; let physics = false; +let skeleton = false; +let polyhedron = false; + function preload() { - earth = loadImage("atlas1.jpg"); + earth_image = loadImage("atlas1.jpg"); } function setup() { - createCanvas(600, 600, WEBGL); + createCanvas(windowWidth, windowHeight, WEBGL); + camera = createCamera(); + red = color(0xbf, 0x00, 0x00); +} + +function windowResized() { + resizeCanvas(windowWidth, windowHeight); } function draw() { @@ -18,51 +35,117 @@ function draw() { background(50); make_lights(); - if (physics) { - move_charges(); - } + if (physics) move_charges(charges); + draw_charges(sphere_radius); + if (skeleton) draw_skeleton(sphere_radius); + if (polyhedron) { + if (physics || faces.length === 0) find_faces(); + draw_faces(sphere_radius); + } draw_sphere(sphere_radius, 25); } -function move_charges() { - for (charge of charges) { - charge.acceleration.setMag(0); - } - for (let i = 0; i < charges.length; i += 1) { - for (let j = 0; j < i; j += 1) { - const displacement = p5.Vector.sub( - charges[i].position, - charges[j].position, - ); - const acceleration_mag = 1 / displacement.mag() * 0.005; - let ai = displacement.copy().normalize().mult(acceleration_mag); - let aj = p5.Vector.mult(ai, -1); - project_onto_plane(ai, charges[i].position); - project_onto_plane(aj, charges[j].position); - charges[i].acceleration.add(ai); - charges[j].acceleration.add(aj); +function face_dist_sq([v1, v2, v3]) { + const center = p5.Vector.add(v1, v2).add(v3).mult(1 / 3); + return createVector(camera.eyeX, camera.eyeY, camera.eyeZ).sub(center).magSq(); +}; + +function find_faces() { + faces = []; + for (let i = 2; i < charges.length; i += 1) { + for (let j = 1; j < i; j += 1) { + for (let k = 0; k < j; k += 1) { + // Check if p1 p2 p3 form a face of the convex polyhedron + // enclosing all vertices ... + const p1 = charges[i].position; + const p2 = charges[j].position; + const p3 = charges[k].position; + const normal = p5.Vector.sub(p2, p1).cross(p5.Vector.sub(p3, p1)); + // ... by checking if the other vertices are on the same + // side of the plane generated by p1 p2 p3 + let plane_separates_vertices = false; + let euler_formula = false; + for (let r = 1; r < charges.length; r += 1) { + for (let s = 0; s < r; s += 1) { + if ( + r === i || r === j || r === k || + s === i || s === j || s === k + ) continue; + const q1 = charges[r].position; + const q2 = charges[s].position; + // Let l(t) := q1 + (q2 - q1) * t. + // L := { l(t) : 0 <= t <= 1 } is the line segment + // between q1 and q2. L intersects the plane + // generated by p1 p2 p3 iff l(t) intersects with + // the plane for some 0 <= t <= 1. + // + // If L intersects the plane, q1 and q2 are on + // opposite sides of the plane generated by p1 p2 p3, + // so p1 p2 p3 can't be a face. If we want the + // polyhedron to be convex. Which we do. + // + // A point k is on the plane generated by p1 p2 p3 iff + // dot(k - p1, normal) = 0. Let n := normal. + // + // dot(l(t) - p1, n) = 0 + // iff dot(q1 + (q2 - q1) * t - p1, n) = 0 + // iff dot(q1 - p1, n) + dot(q2 - q2, n) * t = 0 + // iff t = dot(p1 - q1, n) / dot(q2 - q2, n) + const t = ( + p5.Vector.dot(p5.Vector.sub(p1, q1), normal) / + p5.Vector.dot(p5.Vector.sub(q2, q1), normal) + ); + plane_separates_vertices ||= t >= 0 && t <= 1; + if (plane_separates_vertices) break; + euler_formula ||= charges.length * 2 - faces.length == 4; + if (euler_formula) break; + } + if (plane_separates_vertices || euler_formula) break; + } + if (euler_formula) return; + if (!plane_separates_vertices) { + faces.push([p1, p2, p3]); + } + } } } - for (let i = 0; i < charges.length; i += 1) { - let charge = charges[i]; - charge.velocity = charge.velocity.add(charge.acceleration); - //project_onto_plane(charge.velocity, charge.position); - charge.position = charge.position.add(charge.velocity); - charge.position.normalize(); +} + +function draw_faces(radius) { + // fix OpenGL stacking alpha behaviour + faces.sort((a, b) => face_dist_sq(b) - face_dist_sq(a)); + push(); + strokeWeight(2); + stroke(0x00); + for ([p1, p2, p3] of faces) { + fill(0xbf, 0x7f); + beginShape(TRIANGLES); + vertex(p1.x * radius, p1.y * radius, p1.z * radius); + vertex(p2.x * radius, p2.y * radius, p2.z * radius); + vertex(p3.x * radius, p3.y * radius, p3.z * radius); + endShape(); } + pop(); } -function project_onto_unit_vector(v, unit_vector) { - let size = p5.Vector.dot(v, unit_vector); - return p5.Vector.mult(unit_vector, size); -} - -/// Project `v` onto the plane normal to `unit_vector` -/// Mutates `v` -function project_onto_plane(v, unit_vector) { - let v_proj = project_onto_unit_vector(v, unit_vector); - v.sub(v_proj) +function draw_skeleton(radius) { + push(); + noStroke(); + fill(0xff); + sphere(4); + stroke(0xbf); + for (let charge of charges) { + line( + 0, + 0, + 0, + charge.position.x * radius, + charge.position.y * radius, + charge.position.z * radius, + ); + } + pop(); } function make_charges(n) { @@ -84,6 +167,7 @@ function make_charges(n) { position: position, velocity: createVector(), acceleration: createVector(), + color: red, }); } } @@ -91,8 +175,8 @@ function make_charges(n) { function draw_charges(radius) { push(); noStroke(); - ambientMaterial(0xbf, 0x00, 0x00); for (let charge of charges.values()) { + ambientMaterial(charge.color); let position = charge.position.copy(); position.mult(radius); push(); @@ -107,11 +191,14 @@ function draw_sphere(radius, n_axis_circles) { stroke(0x3f); noFill(); + if (surface === SURFACE_NONE) n_axis_circles = 0; + else if (surface === SURFACE_EARTH) n_axis_circles = 2; + push(); rotateX(TAU / 4); draw_circles( radius, - sphere_mode === 'earth' ? 2 : n_axis_circles, + n_axis_circles, color(0x00, 0x9f, 0xff), color(0xff, 0x9f, 0x00), ); @@ -120,18 +207,18 @@ function draw_sphere(radius, n_axis_circles) { rotateY(TAU / 4); draw_circles( radius, - sphere_mode === 'earth' ? 2 : n_axis_circles, + n_axis_circles, color(0xff, 0x00, 0xff), color(0x00, 0xff, 0x00), ); pop(); - if (sphere_mode === 'earth') { + if (surface === SURFACE_EARTH) { + push(); noStroke(); noFill(); tint(0xff, 0x9f); - texture(earth); - push(); + texture(earth_image); rotateY(TAU / 4); sphere(radius); pop(); @@ -167,11 +254,16 @@ function make_lights() { } function keyPressed() { - if (key == 'd') { - sphere_mode = sphere_mode === 'earth' ? 'circles' : 'earth'; - } else if (key == ' ') { + if (key == ' ') { physics = !physics; + } else if (key == 'd') { + surface = (surface + 1) % 3; + } else if (key == 'f') { + skeleton = !skeleton; + } else if (key == 'g') { + polyhedron = !polyhedron; } else if (key >= '0' && key <= '9') { make_charges(int(key)); + faces = []; } } diff --git a/thomson-problem.js b/thomson-problem.js new file mode 100644 index 0000000..0cbbdaa --- /dev/null +++ b/thomson-problem.js @@ -0,0 +1,44 @@ +function move_charges(charges) { + for (charge of charges) { + charge.acceleration.setMag(0); + } + for (let i = 0; i < charges.length; i += 1) { + for (let j = 0; j < i; j += 1) { + const displacement = p5.Vector.sub( + charges[i].position, + charges[j].position, + ); + let acceleration_mag = 1 / displacement.mag() * 0.001; + let ai; + if (acceleration_mag === Infinity) { + ai = createVector(random(-1, 1), random(-1, 1), random(-1, 1)); + } else { + ai = displacement.copy(); + } + ai = ai.normalize().mult(acceleration_mag); + let aj = p5.Vector.mult(ai, -1); + project_onto_plane(ai, charges[i].position); + project_onto_plane(aj, charges[j].position); + charges[i].acceleration.add(ai); + charges[j].acceleration.add(aj); + } + } + for (let i = 0; i < charges.length; i += 1) { + let charge = charges[i]; + charge.velocity = charge.velocity.add(charge.acceleration); + charge.position = charge.position.add(charge.velocity); + charge.position.normalize(); + } +} + +function project_onto_unit_vector(v, unit_vector) { + let size = p5.Vector.dot(v, unit_vector); + return p5.Vector.mult(unit_vector, size); +} + +/// Project `v` onto the plane normal to `unit_vector` +/// Mutates `v` +function project_onto_plane(v, unit_vector) { + let v_proj = project_onto_unit_vector(v, unit_vector); + v.sub(v_proj) +}