Moving some moons

This week, let’s continue exploring some concepts that the Advent of Code (AoC) challenges bring back to my memory. During Day 12 this year, Eric Wastl offered to examine how we could simulate a very basic newtonian system of 4 moons that attract each other and move in a periodic pattern over time. Just like for Day 15, I also made a little visualization to make the explanation a bit more catchy to the eye 🙂

Day 12 was a physics simulation of sorts that made us study how a system of 4 objects interacts. Eric decided upon some basic gravity- and velocity-simulating rules that can be used to model the evolution of 4 small moons in space. The puzzle was about using those rules to move moons. The gist is as follows.

Understanding the puzzle

Let’s say those 4 moons are initially positioned somewhere in space. Then, you know that they attract each other “in pairs”:

  • two moons attract each other in the x, y and z axis by “pulling” the other towards them; for example, just considering the x axis, if moon A has a smaller x coordinate than moon B, then the x velocity of moon A will increase by 1 and the x velocity of moon B will decrease by one
  • this happens for every possible pair of moons: moon A – moon B, moon A – moon C, moon A – moon D, moon B – moon C, moon B – moon D, moon C – D

Here is a rough idea of how two moons would attract on the x and y axis – we assume that they have the same z coordinate and therefore neither has an impact on the other’s z velocity:

 

Once it has been updated, the velocity of each moon acts on its position and moves it a bit around. All this process (gravity and velocity computation) counts as one iteration.

Note: this physical model is very simplified compared to the reality – in particular, the fact that the objects’ masses don’t impact the way their velocity changes is quite unrealistic. Still, this gives a fun and easy-to-get model with already some cool results.

The point of the puzzle is to simulate the moons’ movement for a given number of iterations in order to get the final total energy of the system (meaning the sum of the potential and kinetic energy of all the moons).

Besides calculating the movements per se, Day 12 challenge was also about finding out the periodicity of this movement; we are told that there is a finite number of iterations after which the moons will all be in a state that was encountered previously (i.e. each moon has the exact same position and velocity for the 3 axis as it did during a previous iteration). This means that from that point on, history will repeat itself and we will get the same evolution.

One additional goal I had was to visualize the moons movement with a little 3D animation.

The problem can thus be decomposed in 4 big parts:

  1. how do we find all the possible moon pairs?
  2. how do we apply gravity and velocity at each iteration?
  3. how do find the period of the movement?
  4. how do visualize our result in a neat way to actually see the moons move around?

Finding the moon pairs using combinatorics

To find all the possible pairings, we can use combinatorics to extract combinations from our list. Combinations are selections of items from a list where the order does not matter. In our case, our list contains 4 items: A, B, C and D; this gives us 6 combinations: A – B, A – C, A – D, B – C, B – D, C – D.

The important thing is that, since order does not matter, two pairs that link the same objects only in a different order are considered to be the same combination. For example, A – D and D – A are the same combination. This is why combinations are interesting for us here: by keeping only the unique pairs, we avoid applying the gravity pull twice on the same pair of objects.

Note: on the other hand, we would have 12 permutations because for those, order matters and thus A – D and D – A are different permutations.

Doing some physics

Once we have our list of moon pairs, we can easily apply our “gravity rules” iteration per iteration. The process is quite straight-forward now that we have detailed the process:

  • first, go through the list of pairs; for each pair, check if the two moons have a different coordinate and thus attract each other on this axis for the x, y and z axis – this may update the vx, vy and vz values of those two moons
  • then, go through the list of moons; for each moon, update its position by its velocity on the x, y and z axis

The corresponding Python code looks like that (it is a bit more commented with docstrings in the Github repo):

from itertools import combinations

def move_moons(timesteps, moons, objs):
   # prepare moons with their velocities
   moons = [ [ x, y, z, 0., 0., 0. ] for x, y, z in moons ]
   # prepare all the unique moon pairs
   moon_pairs = list(combinations(range(len(moons)), 2))

   for time in range(timesteps):
      # apply gravity
      for i1, i2 in moon_pairs:
         x1, y1, z1, vx1, vy1, vz1 = moons[i1]
         x2, y2, z2, vx2, vy2, vz2 = moons[i2]
         if x1 > x2:
            moons[i1][3] -= 1; moons[i2][3] += 1
         elif x1 < x2:
            moons[i1][3] += 1; moons[i2][3] -= 1
         if y1 > y2:
            moons[i1][4] -= 1; moons[i2][4] += 1
         elif y1 < y2:
            moons[i1][4] += 1; moons[i2][4] -= 1
         if z1 > z2:
            moons[i1][5] -= 1; moons[i2][5] += 1
         elif z1 < z2:
            moons[i1][5] += 1; moons[i2][5] -= 1

       # apply velocity
       for i, moon in enumerate(moons):
          moon[0] += moon[3]
          moon[1] += moon[4]
          moon[2] += moon[5]

Python has a built-in module called itertools that gives us access to lots of efficient container generators and combinatoric tools: permutations, combinations, products… it’s really useful when you need to compute such things quickly!

Searching for the period of the movement

We want to investigate what this period is, but we know it can be quite long (more than a million iterations). Just simulating the process over and over will take ages! However there is a trick to drastically reduce the computation time: taking advantage of the fact that in this problem, each axis is independent. In other words, we can find out the periodicity of the x positions and velocities (we’ll call it px), the one of the y positions and velocities (py) and the one of the z positions and velocities (pz).

Then, the period of the full system is the least common multiple of px, py and pz which is the smallest integer that is divisible by these three numbers; indeed, the global period is the “first time” that the 3 periods happen at the same time.

Visualizing the result with Blender

The final part is about visualization: all those calculations are really nice but you don’t actually picture the moons movement just by looking at these numbers. Wouldn’t it be cool if we could truly make some movie of the simulation?

To do that, I’ve used the famous 3D software called Blender. This open-source and free tool is really powerful and can be used for lots of things: small scenes modeling, complex character creation and animation, or even full 3D movie production. And something that is very interesting with Blender is that, among other features, it allows us to do actions through Python scripts. In particular, we can quite easily instantiate, rotate, scale and translate objects in our 3D scene. And since I already had a Python solution for the AoC challenge, I basically just needed to transfer it to Blender and adapt some of it to touch upon my scene.

Note: I also made a small change in the moon’s initial coordinates to “recenter” the system – I simply computed the barycenter of the set of moons and did the difference to update their initial positions and “move” the barycenter to the origin of the scene.

The script consists in several blocks:

  • a few imports to find combinations easily, to do some random computation and to access the Blender Python API:
import bpy
from random import random
from itertools import combinations
  • some global variables for the scene (re)creation:
materials = bpy.data.materials
N_STEPS = 152
RATE = 12
scene = bpy.data.scenes["Scene"]
scene.frame_start = 1
scene.frame_end = (N_STEPS + 1) * RATE
  • a little clean-up of the scene in case there already are moons in it:
for o in bpy.data.objects:
   if 'Moon' in o.name:
     o.select = True
     bpy.ops.object.delete()
  • a function to instantiate the moons at their initial positions as spheres of different colors:
def instantiate_moons(moons):
   # find barycenter of the moons to recenter the system
   xs, ys, zs = zip(*moons)
   bary_x = sum(xs) / len(moons)
   bary_y = sum(ys) / len(moons)
   bary_z = sum(zs) / len(moons)
   # set initial animation frame
   bpy.context.scene.frame_current = 1
   # prepare objects in the scene
   objs = []
   for i, (x,y,z) in enumerate(moons):
      # (update the moon's position to realign with the barycenter)
      moons[i] = (x-bary_x,y-bary_y,z-bary_z)
      # instantiate a sphere
      bpy.ops.mesh.primitive_uv_sphere_add(
         location=(x-bary_x,y-bary_y,z-bary_z)
      )
      bpy.ops.object.shade_smooth()
      # select object and set its properties
      obj = bpy.context.active_object
      obj.name = 'Moon%d' % i
      s = 0.5 + random()
      bpy.ops.transform.resize(value=(s,s,s))
      obj.data.materials.append(materials[i])
      bpy.ops.anim.keyframe_insert_menu(type='Location')
      obj.keyframe_insert(data_path='location')
      objs.append(obj)
   return objs
  • finally, the function to simulate the moons’ movement for a given number of steps (in blue, the differences with the plain Python function from above that did not have any Blender in it):
def move_moons(timesteps, moons, objs):
   # prepare moons with their velocities
   moons = [ [ x, y, z, 0., 0., 0. ] for x, y, z in moons ]
   # prepare all the unique moon pairs
   moon_pairs = list(combinations(range(len(moons)), 2))

   for time in range(timesteps):
      bpy.context.scene.frame_current = (time + 1) * RATE
      # apply gravity
      for i1, i2 in moon_pairs:
         x1, y1, z1, vx1, vy1, vz1 = moons[i1]
         x2, y2, z2, vx2, vy2, vz2 = moons[i2]
         if x1 > x2:
            moons[i1][3] -= 1; moons[i2][3] += 1
         elif x1 < x2:
            moons[i1][3] += 1; moons[i2][3] -= 1
         if y1 > y2:
            moons[i1][4] -= 1; moons[i2][4] += 1
         elif y1 < y2:
            moons[i1][4] += 1; moons[i2][4] -= 1
         if z1 > z2:
            moons[i1][5] -= 1; moons[i2][5] += 1
         elif z1 < z2:
            moons[i1][5] += 1; moons[i2][5] -= 1
      # apply velocity
      for i, moon in enumerate(moons):
         moon[0] += moon[3]
         moon[1] += moon[4]
         moon[2] += moon[5]
         obj = objs[i]
         obj.location = (moon[0], moon[1], moon[2])
         obj.keyframe_insert(data_path='location')

We end up with a series of PNG images that we can know transform into a movie, for example with the common audio/video conversion tool ffmpeg (note: if you want more details on how to do this conversion, check out the README in the Github repository). Here is an example of final video:

 

It was really cool to mix Python and Blender on this. It’s really interesting to see how those powerful tools can communicate and how we can use them to build a neat movie from a simple puzzle. Once again, big thanks to Eric for all his Advent of Code puzzles (if he ever reads this…) and stay tuned for more articles about science and tech!

References
  1. Advent of Code’s about: https://adventofcode.com/2019/about
  2. Blender’s website: https://www.blender.org/
  3. ffmpeg’s website: https://ffmpeg.org/
  4. My Github’s repository for AoC: https://github.com/MinaPecheux/Advent-Of-Code
  5. I. Wikimedia Foundation, “Least common multiple” (https://en.wikipedia.org/wiki/Least_common_multiple), September 2019. [Online; last access 23-December-2019].
  6. I. Wikimedia Foundation, “Barycenter” (https://en.wikipedia.org/wiki/Barycenter), November 2019. [Online; last access 23-December-2019].

Leave a Reply

Your email address will not be published. Required fields are marked *