ooo.org

Object Oriented Orbits: a primer on Newtonian physics

1 What does it take to simulate orbits?

2 What we need

2.1 Before we can simulate orbits, we need to a few things

• a model of space to organize the simulated bodies
• a dynamic rule to update the locations of bodies in space

3 Euclid’s axioms

3.1 The first complete model of space ever recorded was compiled by Euclid in ancient Greece.

4 Axiom 1

4.1 Between any two points ${\color{red} A}$ and ${\color{blue} B}$, a line segment $L$ can be drawn

5 Axiom 2

5.1 A line segment $L$ can be extended indefinitely to a larger line segment $L’$, that contains $L$

6 Axiom 3

6.1 A circle can be drawn at any point with any radius

7 Axiom 4

7.1 All right angles are congruent

8 Axiom 5 (The Parallel Postulate)

8.1 Given a line $L$ and a point $p$ not on the line, there is exactly one line $L’$ through $p$ that doesn’t intersect $L$

9 Theorems

9.1 From these five axioms, we can deduce many useful things, the most useful for our purposes will be the Pythagorean theorem.

9.2 We can use this to compute distance

10 Axioms 1 and 2 and vectors

10.1 Vectors are directed line segments, which can be scaled by real numbers, so axioms 1 and 2 are relevant for vectors

1. Given any two points $A$ and $B$, a vector $\vec{v}$ exists whose tail is $A$ and head is $B$
2. Given any vector $\vec{v}$ and any real number $c$, $c\vec{v}$ extends $\vec{v}$ by a factor of $c$

11 Some terminology

11.1 We call the initial point of a vector its tail

11.2 The final point of the vector is called its head

12 Vectors can be added

12.1 Given any two vectors $\vec{v}$ and $\vec{w}$ with the same tail, their sum $\vec{v} + \vec{w}$ can be visualized using a parallelogram:

12.2 This uses axiom 5, and this operation is commutative

13 Vectors and coordinate systems

13.1 Given a coordinate system, we can represent vectors using pairs (2D) or triples (3D) of real numbers:

13.2 There is a special point, $\vec{0}$ which is just the origin.

14 Vectors have a ‘dot product’

14.1 Given any two vectors $\vec{v} = (v_1,v_2,v_3)$ and $\vec{w} = (w_1,w_2,w_3)$

14.2 their dot product $\vec{v} ⋅ \vec{w} = v_1w_1 + v_2w_2 + v_3w_3$

14.3 Useful fact: $\vec{v} ⋅ \vec{w} = |v||w|cos(θ)$

14.4 That also implies that $\sqrt{\vec{v} ⋅ \vec{v}}$ is the length of the vector

15 Distance between vectors

15.1 We are using vectors to represent points in space, so we will compute the distance between the points $V$ and $W$ by computing $\sqrt{(\vec{v}-\vec{w})⋅ (\vec{v}-\vec{w})}$. This dot product magic just follows from the Pythagorean theorem.

16 Vectors in Ruby (components)

16.1 Now that we have a model of space, we can start writing some ruby code

• a Vector has components (the coordinates)

class Vector
  attr_reader :components

  def initialize(components)
    @components = components
  end
end

17 Vectors in Ruby (algebra)

• a Vector can be added to another vector
• a Vector can be multiplied by a scalar

class Vector
  def +(vector)
    sums = components.zip(vector.components).
                      map {|(vi,wi)| vi+wi }
    Vector.new(sum)
  end

  def *(scalar)
    Vector.new(components.map{|c| scalar*c })
  end
end

18 Vectors in Ruby (equality and dot product)

• we can compare two vectors for equality
• we can take the dot product of two vectors and get the scalar

class Vector
  def ==(vector)
    components == vector.components
  end

  def dot(vector)
    pairs = components.zip(vector.components)
    pairs.map {|(vi,wi)| vi*wi }.
          inject(&:+)
  end
end

19 Time

19.1 Now we have a decent model of space, we can move on to the dynamic rule, it will be a way to update the state of the bodies in space over time.

20 Relation between position, time and velocity

20.1 We can represent the path a body takes using a function $\vec{x}(t)$.

20.2 The velocity is then just the rate of change of position with respect to time

21 Relation between velocity and acceleration

21.1 Similarly, the acceleration is the rate of change of velocity with respect to time

22 Newton’s 1st Law states that

22.1 Bodies travel in straight lines with constant velocity unless a force is acting on it

23 Newton’s 2nd Law states that

23.1 The vector sum of forces acting on a body is its acceleration times its mass

23.2 Note that forces are represented as vectors

24 Newton’s Law of Universal Gravitation

25 Bodies in Ruby

25.1 the Body class should have a read-only mass

25.2 along with a position and a velocity

class Body
  attr_reader :mass
  attr_accessor :position, :velocity

  def initialize(mass:, position:, velocity:)
    @mass = mass
    @position = Vector.new(position)
    @velocity = Vector.new(velocity)
  end
end

26 Forces on Bodies in Ruby

26.1 Bodies have a method to compute the gravitational force acting on it from another Body.

class Body
  def force_from(body)
    rvec = body.position - position
    r = rvec.norm
    rhat = rvec * (1/r)
    rhat * (Newtonian.G * mass * body.mass / r**2)
  end
end

27 the Universe

28 the Universe in Ruby

28.1 The final class will be Universe, it organizes all the bodies

class Universe
  attr_reader :dimensions, :bodies

  def initialize(dimensions:, bodies:)
    @dimensions = dimensions
    @bodies = bodies
  end
end

28.2 it also has a number of **dimensions**, we can use this to make sure the bodies are all in the same kind of space

29 the Enumerable Universe

29.1 Since force is computed pairwise, we create an iterator for pairs of distinct objects

class Universe
  def each_pair_with_index
    bodies.each_with_index do |body_i, i|
      bodies.each_with_index do |body_j, j|
        next if i == j
        yield [body_i, body_j, i, j]
      end
    end
  end
end

30 The main simulation loop

class Universe
  def evolve(dt)
    forces = bodies.map{ |_| zero_vector }
    each_pair_with_index do |(body_i, body_j, i, j)|
      forces[i] += body_i.force_from(body_j)
    end
    bodies.each_with_index do |_, i|
      a = forces[i] * (1.0 / bodies[i].mass)
      v = bodies[i].velocity 
      bodies[i].velocity += a * dt
      bodies[i].position += v * dt
    end
  end
end

31 The server

31.1 We can serve this up to a browser using

• WEBrick for HTTP
• websocketd for piping STDOUT to a WebSocket server

32 Fork off an HTTP server

rd, wt = IO.pipe
pid = fork do
  rd.close
  server = WEBrick::HTTPServer.new({
    :Port => PORT,
    :BindAddress => "localhost",
    :StartCallback => Proc.new {
      wt.write(1)  # write "1", signal start
      wt.close
    }
  })
  trap('INT') { server.stop }
  server.mount("/", WEBrick::HTTPServlet::FileHandler, './examples')
  server.start
end
# ...

33 Shell out to websocketd

33.1 websocketd converts standard input and output into a fully functioning websocket server, so we can just puts out the universe state

examples = ["binary.rb", "ternary.rb", "random.rb", "figure_eight.rb"]
index = ARGV.last.to_i
# Shell out to websocketd, block until program finishes
system("bin/websocketd \
   -port=8080 \
   ruby #{examples[index]}")

Process.kill('INT', pid)   # kill HTTP server in child process

34 Binary Star system

34.1 Our first application is going to be simulating a binary star system, with two equal-mass stars

35 Find initial conditions

35.1 The two bodies will be traveling in uniform circular motion, so the following relation holds:

36 Given the masses and the distance $r$, we can figure out $a$:

37 Substituting a back in to get $v$

38 Run simulated binary star system

38.1 Pause to run simulations

39 The Three Body Problem

39.1 With only two bodies, it turns out to be possible to solve the equations of motion for all time, exactly.

39.2 With three or more bodies, it is in general impossible

40 However

40.1 The three body problem has been studied since 1747, and there are some well known examples

41 The “Figure Eight” Three Body Orbit

41.1 The paper “A remarkable periodic solution of the three-body problem in the case of equal masses” by Alain Chenciner and Richard Montgomery works out an orbit that looks like this:

42 The initial conditions

42.1 $\vec{x_1} = [0.97000436, -0.24308753]; \vec{x_2} = -\vec{x_1}$

42.2 $\vec{x_3} = \vec{0}$

42.3 $\vec{v_3} = [-0.93240737, -0.86473146]$

42.4 $-2\vec{v_1} = -2\vec{v_2} = \vec{v_3}$

43 Run simulated three-body orbit

43.1 Pause to run simulations