Previous | ToC | Up | Next |

** Bob**: Hi Alice, what's next? I can think of a few further improvements
to my Plummer code, but do you have something specific, you'd like to add?

** Alice**: I'm happy that we have done this energy checking, which seems to
indicate that your code probably does produce Plummer's model correctly.
However, I'd like to be even more sure, since it would really be terrible
if we start our future simulations with the wrong initial
conditions. That would make debugging our evolution codes quite confusing,
to say the least.

** Bob**: There are many ways to check things. The question is, what would
be the simplest way. Ah, how about just looking at mass shells,
shells with a certain mass fraction, and check that they have the
correct radius? After all, we computed the cumulative mass
as a function of the radius already.

** Alice**: That is a good idea. For one thing, the half-mass radius should
come out with the correct value. And for good measure, we may as well
compute the quarter-mass radius and the three-quarter-mass radius. In
other words: the radii of the three quartiles, within which 1/4, 1/2,
and 3/4 of the total mass are enclosed.

** Bob**: That we should be able to do from scratch. We can start again with
minimal versions of the ` Body` and ` Nbody` classes, just enough to
read a snapshot in, and initialize everything properly. The only thing to
do is add a method ` quartiles` that prints out the three quartiles.

Fortunately, Ruby has a handy ` sort` method, that comes with the ` Array`
class. As usual, it has an associated `sort!` method that affects
the array itself, unlike ` sort` that returns a new array. Let's see.
We have to sort the particles in radial order, using ,
the distance from the center. It will be easier to use the square instead,
. This can be done in three lines:

def order_squared_radii a = [] @body.each{|b| a.push b.pos*b.pos} a.sort! end

The first line defines ` a` to be an empty array. The second line fills
the array with the squares of the radial positions of all particles.
Then the third line sorts that array. Couldn't be simpler!

** Alice**: Can I write the ` quartiles` method? I'd like to get some more
experience with Ruby.

** Bob**: Here is the keyboard!

** Alice**: Since you have defined `order_squared_radii` as a
member function, we may as well define ` quartiles` as a member function
too. We can then immediately invoke `order_squared_radii` to
create our order list of particle positions. All we have to do then
is to cut the list in four equal parts, and document the places where
we cut the list. Hey, that is too simple to learn much about Ruby!
This should do the job:

def quartiles a = order_squared_radii n = a.size r_1 = a[(n/4.0).round - 1] r_2 = a[(n/2.0).round - 1] r_3 = a[(n*3/4.0).round - 1] print "The values of the three quartiles for r(M) are:\n" print " r(1/4) = " printf("%.4g\n", r_1) print " r(1/2) = " printf("%.4g\n", r_2) print " r(3/4) = " printf("%.4g\n", r_3) end

** Bob**: I bet it does. I saw you looking up ` round`: indeed, that is
a predefined method that rounds a floating point number to an integer.
In case of a number exactly in between, such as 2.5, it rounds up, to 3.0.
You could have used integer division directly, instead of floating
point division, but I agree that what you wrote might be easier to understand.

** Alice**: We could quibble about exactly where to choose the boundaries
for the quartiles. For example, for a four-particle system, things come
out just right, but for a five particle system, the half mass radius
encloses 3/5 of the mass. We could have take the average of the
two-particle and three-particle distance, but that would have been
overkill, in my opinion. Anyway, the fluctuations in the positions of
the particles, in an N-body system will be of order
, so we don't have to worry about enclosing
one particle more or less, at least for .

** Bob**: Let me print out the whole program, and write it out to a file called
`quartiles1.rb`:

require "acs" class Body attr_accessor :mass, :pos, :vel def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0]) @mass, @pos, @vel = mass, pos, vel end end class NBody attr_accessor :time, :body def initialize @body = [] end def order_squared_radii a = [] @body.each{|b| a.push b.pos*b.pos} a.sort! end def quartiles a = order_squared_radii n = a.size r_1 = a[(n/4.0).round - 1] r_2 = a[(n/2.0).round - 1] r_3 = a[(n*3/4.0).round - 1] print "The values of the three quartiles for r(M) are:\n" print " r(1/4) = " printf("%.4g\n", r_1) print " r(1/2) = " printf("%.4g\n", r_2) print " r(3/4) = " printf("%.4g\n", r_3) end end include Math nb = ACS_IO.acs_read(NBody) nb.quartiles

** Alice**: Let's try it out:

|gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 100 pseudorandom number seed given : 0 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 : 157807557255272290078302474298427494808 The values of the three quartiles for r(M) are: r(1/4) = 0.642 r(1/2) = 1.737 r(3/4) = 4.359 |gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 100 pseudorandom number seed given : 0 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 : 251979634788701045369825442993623037445 The values of the three quartiles for r(M) are: r(1/4) = 0.5245 r(1/2) = 1.511 r(3/4) = 3.524 |gravity> kali mkplummer2.rb -n 100 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 100 pseudorandom number seed given : 0 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 : 116741143007870386439078079644181238717 The values of the three quartiles for r(M) are: r(1/4) = 0.5138 r(1/2) = 1.192 r(3/4) = 3.349Quite noisy. Not too surprising: for the inner quartile we are dealing with only 25 particles, so we should expect noise. Indeed, there are fluctuations of about twenty percent in the first quartile.

** Bob**: Let's try 10,000 particles instead:

|gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 195924807145294240096344655914317919224 The values of the three quartiles for r(M) are: r(1/4) = 0.6504 r(1/2) = 1.665 r(3/4) = 4.717 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 42611728057165554375136236111523795370 The values of the three quartiles for r(M) are: r(1/4) = 0.6226 r(1/2) = 1.655 r(3/4) = 4.727 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles1.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 323527003508540608555251954323353747057 The values of the three quartiles for r(M) are: r(1/4) = 0.6753 r(1/2) = 1.724 r(3/4) = 4.761Not bad! Fluctuations of a few percent, at most, as it should be. Encouraging.

** Alice**: Now let's check the numbers.
We have used already a few times eq. (26):

** Alice**: Excellent point! And there is an easy check: Let us insert
eq. (96) into the right-hand side of eq. (95),
to check whether we get the original amount of mass back. Let us call
the left-hand side of eq. (95) ,
so that we can distinguish it from the mass we start with,
in eq. (96). Eq. (95) then reads:

** Bob**: Yes, many bugs turn out to be `almost trivial' as well, and yet
they can waste many hours of your time while you're chasing them.
By the way, wouldn't it have been much faster to substitute
eq. (96) into the right-hand side of eq. (94)?
Calling the used in eq. (94)
, we would have gotten:

** Bob**: Let's compute the expected quartile radii, from eq. (96):

** Bob**: What a pity! I was more or less convinced that my Plummer code
was correct. Obviously there is something seriously wrong with the way
particles are sprinkled into space.

** Alice**: What is really strange is the fact that we got the energy correct.
First of all, the virial theorem did hold, something that is typically
violated if you make a random mistake in the energies, since it is unlikely
that you will make the same mistake in the kinetic as in the potential energy.
But secondly, we actually computed the energy, and we found that your code
produced it correctly.

** Bob**: Yes, that is puzzling. Unless . . .

** Alice**: . . . unless my quartile code is wrong.
But that is also very unlikely.
How much can there be wrong in just a few lines of code?

** Bob**: It wouldn't be the first time, to find a nasty bug in a small code.
I guess that's why they are called bugs! Bugs are small, and can hide even
in a code of a single line. Let's have a look.

** Alice**: Yes, the good thing about a small code is that it makes sense to go
through each line, one by one. And there are less than ten lines here
where actual calculations are done.

** Bob**: Let's be systematic. The driver part calls ` quartiles`, and that
method in turn starts by calling `order_squared_radii`. You like
long names, don't you? You could have called it
`place_an_order_to_order_squared_radii`.

** Alice**: Yes, I do like long names, since that way I will remember later
what was doing what. I hate that old style in which it probably would
have been called ` srtrsq` for `sorting r squared'. Forty years ago there
was a reason, when single memory locations were expensive, but that excuse
has long since gone.

** Bob**: Well, I can't for the life of me see anything wrong with
`order_squared_radii`. Back to ` quartiles`. In the second line
you find the number of particles. Nothing wrong with that either.
In the third line you cut the list at the place of the first quarter.
Indeed, in that way you are bound to find the mass shell that contains
one quarter of the mass, and precisely the inner quarter.

** Alice**: Iron logic, and I, too, cannot see anything wrong with that.

** Bob**: So we both agree that `a[(n/4.0).round - 1]` returns
what you have ordered: the squared radius of the inner quarter of . . .

** Alice**: . . . the *squared* radius!!

** Bob**: Ah, yes, that's the problem! You had forgotten that you ordered
the *squares* of the radial distances, for convenience. Instead,
you just assigned the values to the quartile radii themselves.

** Alice**: That must have been the problem. And that immediately explains
why the errors grew bigger for larger values of * r*. Okay, that is easy
to fix. Let us rewrite these three lines, and put the result in a file
called `quartiles.rb`:

r_1 = sqrt(a[(n/4.0).round - 1]) r_2 = sqrt(a[(n/2.0).round - 1]) r_3 = sqrt(a[(n*3/4.0).round - 1])

And let's test it right away:

|gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 201299017348521780587188254399647312211 The values of the three quartiles for r(M) are: r(1/4) = 0.8137 r(1/2) = 1.292 r(3/4) = 2.176 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 42982039013636173889503342183695809112 The values of the three quartiles for r(M) are: r(1/4) = 0.8216 r(1/2) = 1.315 r(3/4) = 2.181 |gravity> kali mkplummer2.rb -n 10000 | kali quartiles.rb ==> Plummer's Model Builder <== Number of particles : N = 10000 pseudorandom number seed given : 0 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 : 198688716986344613646322407980102362651 The values of the three quartiles for r(M) are: r(1/4) = 0.8141 r(1/2) = 1.307 r(3/4) = 2.186

** Alice**: You may laugh, and perhaps I'm getting to much addicted to
testing, but let me do one final little check. Having just made a
really silly mistake determining the quartile radii, I feel
I cannot take anything for granted. Let me see whether the quartile
radii are * really* quartile radii. Starting with:

We have

Similarly with the half-mass radius , and:

** Bob**: Now that is what I would call overkill. But you're right, it doesn't
hurt, and better safe than sorry.

** Alice**: Indeed!

Previous | ToC | Up | Next |