diff --git a/index.html b/index.html
index 0e3fca4..3295298 100644
--- a/index.html
+++ b/index.html
@@ -86,16 +86,16 @@ aside > container > div {
 <button id="button-surface-earth" onclick="set_surface(SURFACE_EARTH)">Earth</button>
 </div>
 <div>
-<button onclick="make_charges(0);faces=[]">0</button>
-<button onclick="make_charges(1);faces=[]">1</button>
-<button onclick="make_charges(2);faces=[]">2</button>
-<button onclick="make_charges(3);faces=[]">3</button>
-<button onclick="make_charges(4);faces=[]">4</button>
-<button onclick="make_charges(5);faces=[]">5</button>
-<button onclick="make_charges(6);faces=[]">6</button>
-<button onclick="make_charges(7);faces=[]">7</button>
-<button onclick="make_charges(8);faces=[]">8</button>
-<button onclick="make_charges(9);faces=[]">9</button>
+<button onclick="make_particles(0);faces=[]">0</button>
+<button onclick="make_particles(1);faces=[]">1</button>
+<button onclick="make_particles(2);faces=[]">2</button>
+<button onclick="make_particles(3);faces=[]">3</button>
+<button onclick="make_particles(4);faces=[]">4</button>
+<button onclick="make_particles(5);faces=[]">5</button>
+<button onclick="make_particles(6);faces=[]">6</button>
+<button onclick="make_particles(7);faces=[]">7</button>
+<button onclick="make_particles(8);faces=[]">8</button>
+<button onclick="make_particles(9);faces=[]">9</button>
 </div>
 </container>
 </aside>
diff --git a/sketch.js b/sketch.js
index d0a1b95..b8f9c3d 100644
--- a/sketch.js
+++ b/sketch.js
@@ -1,7 +1,7 @@
 let camera;
 let red;
 
-let charges = [];
+let particles = [];
 let faces = [];
 
 let sphere_radius;
@@ -14,6 +14,7 @@ let surface = SURFACE_CIRCLES;
 let physics = false;
 let skeleton = false;
 let polytope = false;
+let charge = -1;
 
 let buttons_surface;
 let checkbox_physics;
@@ -70,9 +71,9 @@ function draw() {
   camera.centerZ = 0;
 
   make_lights();
-  if (physics) move_charges(charges, 8e-4);
+  if (physics) move_particles(particles, 8e-4);
 
-  draw_charges(sphere_radius);
+  draw_particles(sphere_radius);
   if (skeleton) draw_skeleton(sphere_radius);
   if (polytope) {
     if (physics || faces.length === 0) find_faces();
@@ -88,27 +89,27 @@ function face_dist_sq([v1, v2, v3]) {
 
 function find_faces() {
   faces = [];
-  for (let i = 2; i < charges.length; i += 1) {
+  for (let i = 2; i < particles.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 polytope
         // enclosing all vertices ...
-        const p1 = charges[i].position;
-        const p2 = charges[j].position;
-        const p3 = charges[k].position;
+        const p1 = particles[i].position;
+        const p2 = particles[j].position;
+        const p3 = particles[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 r = 1; r < particles.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;
+            const q1 = particles[r].position;
+            const q2 = particles[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
@@ -133,7 +134,7 @@ function find_faces() {
             );
             plane_separates_vertices ||= t >= 0 && t <= 1;
             if (plane_separates_vertices) break;
-            euler_formula ||= charges.length * 2 - faces.length == 4;
+            euler_formula ||= particles.length * 2 - faces.length == 4;
             if (euler_formula) break;
           }
           if (plane_separates_vertices || euler_formula) break;
@@ -170,21 +171,21 @@ function draw_skeleton(radius) {
   fill(0xff);
   sphere(4);
   stroke(0xbf);
-  for (let charge of charges) {
+  for (let particle of particles) {
     line(
       0,
       0,
       0,
-      charge.position.x * radius,
-      charge.position.y * radius,
-      charge.position.z * radius,
+      particle.position.x * radius,
+      particle.position.y * radius,
+      particle.position.z * radius,
     );
   }
   pop();
 }
 
-function make_charges(n) {
-  charges = [];
+function make_particles(n) {
+  particles = [];
   for (let i = 0; i < n; i += 1) {
     let position;
     if (i === 0) {
@@ -192,21 +193,22 @@ function make_charges(n) {
     } else {
       position = p5.Vector.random3D();
     }
-    charges.push({
+    particles.push({
       position: position,
       velocity: createVector(),
       acceleration: createVector(),
+      charge: charge,
       color: red,
     });
   }
 }
 
-function draw_charges(radius) {
+function draw_particles(radius) {
   push();
   noStroke();
-  for (let charge of charges.values()) {
-    ambientMaterial(charge.color);
-    let position = charge.position.copy();
+  for (let particle of particles) {
+    ambientMaterial(particle.color);
+    let position = particle.position.copy();
     position.mult(radius);
     push();
     translate(position.x, position.y, position.z);
@@ -292,7 +294,7 @@ function keyPressed() {
   } else if (key == 'g') {
     toggle_polytope();
   } else if (key >= '0' && key <= '9') {
-    make_charges(int(key));
+    make_particles(int(key));
     faces = [];
   }
 }
diff --git a/thomson-problem.js b/thomson-problem.js
index d13cc31..3843784 100644
--- a/thomson-problem.js
+++ b/thomson-problem.js
@@ -1,14 +1,18 @@
-function move_charges(charges, force_constant) {
-  for (let charge of charges) {
-    charge.acceleration.setMag(0);
+function move_particles(particles, force_constant) {
+  for (let particle of particles) {
+    particle.acceleration.setMag(0);
   }
-  for (let i = 0; i < charges.length; i += 1) {
+  for (let i = 0; i < particles.length; i += 1) {
     for (let j = 0; j < i; j += 1) {
       const displacement = p5.Vector.sub(
-        charges[i].position,
-        charges[j].position,
+        particles[i].position,
+        particles[j].position,
+      );
+      const force_mag = (
+        particles[i].charge * particles[j].charge
+        / displacement.magSq()
+        * force_constant
       );
-      const force_mag = 1 / displacement.magSq() * force_constant;
       // XXX possible extension: divide by charge's mass
       const ai_mag = force_mag;
       const aj_mag = force_mag;
@@ -23,16 +27,16 @@ function move_charges(charges, force_constant) {
       }
       ai.mult(ai_mag);
       aj.mult(aj_mag);
-      project_onto_plane(ai, charges[i].position);
-      project_onto_plane(aj, charges[j].position);
-      charges[i].acceleration.add(ai);
-      charges[j].acceleration.add(aj);
+      project_onto_plane(ai, particles[i].position);
+      project_onto_plane(aj, particles[j].position);
+      particles[i].acceleration.add(ai);
+      particles[j].acceleration.add(aj);
     }
   }
-  for (let charge of charges) {
-    charge.velocity = charge.velocity.add(charge.acceleration);
-    charge.position = charge.position.add(charge.velocity);
-    charge.position.normalize();
+  for (let particle of particles) {
+    particle.velocity = particle.velocity.add(particle.acceleration);
+    particle.position = particle.position.add(particle.velocity);
+    particle.position.normalize();
   }
 }