Previous ToC Up Next

## 2.1. Phase Space

Bob: Now that we have a shared time step code nbody_sh1.rb, I'd like to make sure that it really gives the same results as the previous code nbody_cst1.rb that used constant time steps. Having good energy accuracy is one thing, but it would be much better to check the actual particle positions.

Alice: How about computing the distance between the results of two runs, in 6N-dimensional phase space? A whole N-body system can be viewed as a single point in a space that spans all the degrees of freedom of the positions and velocities of all particles.

More generally, if we model a system in D dimensions, a single particle will have D components for its position as well as its velocity, or 2D in total. N particles will have 2DN such degrees of freedom. For planar N-body systems, we need a 4N-dimensional phase space, and for N-body systems in the real three-dimensional world we are dealing with a 6N-dimensional phase space.

Bob: In practice this means, I presume, that you use Pythagoras to compute the distance in such a large space. In other words you take the square root of the sum of the squares of all the differences, right?

Alice: Indeed.

Bob: So another way of looking at your suggestion is that you compute the average distance between corresponding particles in the usual 6-dimensional phase space. Apart from a normalization factor N, the result would be the same.

Alice: Yes, that's right. I like the more abstract picture, but your description boils down to the same thing.

Bob: I find it a lot easier to think about 2D dimensions than about 2ND dimensions. In either case, the practical procedure is that we have to subtract two N-body systems, point by point, and then we have to compute the norm or `size' of the resulting system, a type of absolute value.

Alice: Ah, so you do like abstract thinking, after all! How about defining a minus operator, that allows us to subtract two N-body systems, nb1 and nb2, just by writing nb1-nb2? That would give us a nicely compact notation. And the 2DN-dimensional distance between the two systems could then be given by (nb1-nb2).abs if we define the correct absolute value operator abs for the NBody class.

Bob: Good idea! However, there is no guarantee that two N-body systems will show their particles in the output in the exact same order. In other words, what we need is a way to identify which particle in the one system corresponds to which particle in the other system.

Alice: I guess we have to number them, by giving each body a unique identifier.

## 2.2. Numbering Bodies

Bob: Given that every object in Ruby already has an identifier called object_id, a natural name for such an ID would be body_id. Now that is easy to implement. How about writing a nifty tool like this, in file nbody_set_id.rb:

``` #!/usr/local/bin/ruby -w

require "nbody.rb"

options_text= <<-END

Description: Takes an N-body system, and gives each body a unique ID
Long description:
This program accepts an N-body system, and assigns a number to each
body consecutively, as an instance variable @body_id which takes on
integer values.

(c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org

example:
kali #{\$0} -n 0

Short name:          -n
Long name:            --starting_number
Value type:           int
Default value:        1
Description:          value of @body_id for 1st body
Variable name:        n
Long description:
This option allows the user to start numbering the particles at a
value different from the default value (0, say, or 10, or whatever).

END

c = parse_command_line(options_text)

class Body
attr_accessor :body_id
end

i = c.n
nb.body.each do |b|
b.body_id = i
i += 1
end
nb.acs_write(\$stdout, false, c.precision)
```

## 2.3. A Test Count

You see, I even allowed for taste: the default is to start numbering at one, but you can start at zero as well, or at any other place you like. Let's test it with our Plummer Model:

``` |gravity> kali nbody_set_id.rb -h
Takes an N-body system, and gives each body a unique ID
-n  --starting_number  : value of @body_id for 1st body [default: 1]
--verbosity            : Screen Output Verbosity Level  [default: 1]
--acs_verbosity                : ACS Output Verbosity Level     [default: 1]
--precision            : Floating point precision       [default: 16]
--indentation          : Incremental indentation        [default: 2]
-h  --help             : Help facility
---help                        : Program description (the header part of --help)
|gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5
==> Plummer's Model Builder <==
Number of particles            : N = 2
pseudorandom number seed given : 0
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
==> Takes an N-body system, and gives each body a unique ID <==
value of @body_id for 1st body : n = 1
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 5
Incremental indentation                : add_indent = 2
actual seed used  : 252462866881473732163958989257426001914
ACS
NBody
Array body
Body body[0]
Fixnum body_id
1
Float mass
5.00000e-01
Vector pos
1.94755e-01   9.73956e-02  -1.22820e-01
Vector vel
1.78002e-01   5.49335e-01  -4.08101e-01
Body body[1]
Fixnum body_id
2
Float mass
5.00000e-01
Vector pos
-1.94755e-01  -9.73956e-02   1.22820e-01
Vector vel
-1.78002e-01  -5.49335e-01   4.08101e-01
Float time
0.00000e+00
SCA
|gravity> kali mkplummer.rb -n 2 | kali nbody_set_id.rb --precision 5 -n 42
==> Plummer's Model Builder <==
Number of particles            : N = 2
pseudorandom number seed given : 0
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
actual seed used  : 210733041374937108452252448370138453664
==> Takes an N-body system, and gives each body a unique ID <==
value of @body_id for 1st body : n = 42
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 5
Incremental indentation                : add_indent = 2
ACS
NBody
Array body
Body body[0]
Fixnum body_id
42
Float mass
5.00000e-01
Vector pos
1.16550e-01  -4.85085e-02   2.15785e-01
Vector vel
-5.72883e-02   3.45794e-01   6.14121e-01
Body body[1]
Fixnum body_id
43
Float mass
5.00000e-01
Vector pos
-1.16550e-01   4.85085e-02  -2.15785e-01
Vector vel
5.72883e-02  -3.45794e-01  -6.14121e-01
Float time
0.00000e+00
SCA
```

## 2.4. Defensive Programming

Bob: Now, with identifiers in place, we can extend our nbody.rb file by adding a subtraction operator:

```   def -(other)
if other.class != self.class
raise "other.class.name = #{other.class.name} != #{self.class.name}"
end
if (n = other.body.size) != body.size
raise "other.body.size = #{other.body.size} != #{body.size}"
end
nb = NBody.new
body.each_index do |i|
if (id = body[i].body_id) == nil
raise "body[#{i}].body_id == nil"
end
ob = other.body.find{|oi| oi.body_id == id}
if ob == nil
raise "body_id = #{body[i].body_id} not found in other N-body system"
end
nb.body[i] = Body.new
nb.body[i].pos = body[i].pos - ob.pos
nb.body[i].vel = body[i].vel - ob.vel
end
nb
end
```

This includes proper testing of particle identities.

Alice: And I see that you are returning the newly constructed nbody system nb at the end of the method. Also, you're getting quite defensive in your programming!

Bob: This is a situation where defensive programming is called for. In general, I don't like to litter my code with many error checking lines all over the place, but in this particular case, it is quite likely that the difference operator will be invoked for incompatible N-body systems. So I'm checking whether the other object is really of NBody class, whether it has the same number of particles as the calling NBody object, and most importantly, I am then checking whether each particle in the calling system has a counter part in the other system, with the same body_id.

Alice: Applause, applause! As you know, I like to check things, as a matter of habit, in order to prevent future surprises that may be hard to debug. But I agree, you can overdo with anything, and there is a lot to say for keeping code small. In this case, though, the whole code is so small that it doesn't distract to have these extra safeguards built in.

## 2.5. Putting it together

Bob: With our new tools in hand, it now becomes really trivial to check for the actual size of the difference between two N-body systems, by defining the absolute value method, also within the file nbody.rb:

```   def abs
a = 0
body.each{ |b| a += b.pos*b.pos + b.vel*b.vel }
sqrt a
end
```

Alice: Sometimes we may want to look only at the positions, restricting ourselves to configuration space, rather than to phase space.

Bob: That's easy. Here is a position-only version of abs:

```   def abs_pos
a = 0
body.each{ |b| a += b.pos*b.pos }
sqrt a
end
```

Alice: Let me try to implement the module that gives the difference between two N-body systems, as the 6N-dimensional distance in phase space, or the 2DN-dimensional distance for the general case of D spatial dimensions. It's time to practice my Ruby skills again. Let's call the file nbody_diff.rb:

``` #!/usr/local/bin/ruby -w

require "nbody.rb"

options_text= <<-END

Description: 6N-dimensional phase space distance between two N-body systems
Long description:
This program accepts two N-body systems, and computes the distance between
them in the 2DN-dimensional phase space of all positions and velocities,
if they are given in D spatial dimensions.  The values of the masses of
the particles are not used, and each particle is given equal weight.  In
order to compare corresponding particles in the two systems, the body_id
instance variable of each body is used as an identifier.

(c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org

example:
cat run1.out run2.out | kali #{\$0}

Short name:          -r
Long name:            --positions_only
Value type:           bool
Description:          Using only positions, not velocities
Variable name:        position_only_flag
Long description:
This option allows the user to compute the distance between two N-body
systems using only particle information; velocity information is not used.

END

c = parse_command_line(options_text)

class Body
attr_accessor :body_id
end

nb = nb1 - nb2
d = nb.body[0].pos.size
n = nb.body.size
print "#{2*d}N-dim. phase space dist. for two #{n}-body systems:  "
if c.position_only_flag == false
printf(" %.#{c.precision}e\n", nb.abs)
else
printf(" %.#{c.precision}e\n", nb.abs_pos)
end
```

## 2.6. Checking it out

Bob: Looks fine, as far as I can see. You put in options to allow for more precision, without burdening the reader with all those extra digits in the general case. And you also gave an option for positions-only comparisons.

Time to test it! We still have the initial conditions lying around for a circular binary, I believe. Ah yes, here it is, in circular_binary.in:

``` ACS
NBody
Array body
Body body[0]
Float mass
4
Vector pos
1  0
Vector vel
0  1
Body body[1]
Float mass
4
Vector pos
-1  0
Vector vel
0  -1
SCA
```
To start with, let's see whether two identical systems are assigned a distance of zero correctly:

``` |gravity> kali nbody_set_id.rb < circular_binary.in > tmp0.in
==> Takes an N-body system, and gives each body a unique ID <==
value of @body_id for 1st body : n = 1
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
|gravity> cat tmp0.in tmp0.in | kali nbody_diff.rb
==> 6N-dimensional phase space distance between two N-body systems <==
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
4N-dim. phase space dist. for two 2-body systems:   0.0000000000000000e+00
```
Alice: Good! Let's perturb our binary slightly. If we use the most popular Pythagorean triangle, we can immediately check whether our Cartesian distances in phase space come out correctly. Here is a file pert_circ_binary.in:

``` ACS
NBody
Array body
Body body[0]
Float mass
4
Vector pos
1  0.03
Vector vel
0  1
Body body[1]
Float mass
4
Vector pos
-1  0
Vector vel
0  -1.04
SCA
```
And here is our first real result:

``` |gravity> kali nbody_set_id.rb < pert_circ_binary.in > tmp1.in
==> Takes an N-body system, and gives each body a unique ID <==
value of @body_id for 1st body : n = 1
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
|gravity> cat tmp0.in tmp1.in | kali nbody_diff.rb
==> 6N-dimensional phase space distance between two N-body systems <==
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 16
Incremental indentation                : add_indent = 2
4N-dim. phase space dist. for two 2-body systems:   5.0000000000000031e-02
```
Just what it should be! Now for positions only, and a higher precision:

``` |gravity> cat tmp0.in tmp1.in | kali nbody_diff.rb -r --precision 5
==> 6N-dimensional phase space distance between two N-body systems <==
Using only positions, not velocities
Screen Output Verbosity Level  : verbosity = 1
ACS Output Verbosity Level     : acs_verbosity = 1
Floating point precision       : precision = 5
Incremental indentation                : add_indent = 2
4N-dim. phase space dist. for two 2-body systems:   3.00000e-02
```
Perfect. I think we're ready to move on.
 Previous ToC Up Next