Previous | ToC | Up | Next |

** Alice**: Did you start building a star cluster, using Plummer's model as a
blueprint?

** Bob**: Yes, here is a quick and dirty version that I cobbled together.
It was a lot easier than I expected, since I stumbled upon a nice
recipe, which must be a classic in stellar dynamics. I followed it
line by line, in my implementation. It was published in the paper
*A comparison of Numerical Methods for the Study of Star Cluster
Dynamics*, by Sverre Aarseth, Michel Henon, and Roland Wielen,
which appeared ages ago, in 1974, in Astron. Astroph. **37**, 183.

** Alice**: I remember reading that paper as a student.

** Bob**: You are * that* ancient?

** Alice**: Not quite. I must have just started high school at that time.
I mean that I read it as an undergraduate when I wanted to do a small
N-body project for my final thesis. My adviser recommended it. And
indeed, it was an influential paper by some of the masters in the
field, around that time.

I remember it well: the paper made the first detailed quantitative comparison between the Monte-Carlo Fokker-Planck method and direct N-body integration. They even showed the results for runs with a few hundred particles. Hard to believe that they got so much computing done, given the slow speed of computers at that time.

** Bob**: That is impressive, given that computers have increased their speed
by a factor of ten every five years or so. That means that in 2004 computers
are a million times faster than what they had available. I bet that my
camera is a lot more powerful than the computers in their university centers.

** Alice**: And I'm sure your camera has way more storage as well. But
what was actually most impressive was how much information they were
able to squeeze out of their calculations, given their limited
resources. Anyway, why did you bring up this paper?

** Bob**: Because it contains a nifty recipe for constructing a Plummer's model.

** Alice**: Hey, I didn't remember that. Really?

** Bob**: They wrote it as an appendix. Everything was done analytically, except
for the fact that they used a rejection technique at the last step, to
get the velocities.

** Alice**: I must not have looked at the appendix. Well, that is convenient.
Can you show me your code?

** Bob**: Here it is, and I warned you, I haven't commented it or cleaned it up:
I just translated their recipe straight from the appendix into Ruby.

There are three parts to the file `mkplummer1.rb`. The first part
contains the class definitions. It is rather short, since I realized I
could do everything with stripped-down versions of the ` Body` class and the
` Nbody` class:

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(n = 0) @body = [] for i in 0...n @body[i] = Body.new end end end

** Alice**: That all looks familiar. So you only need to create an N-body
system, and then print it out.

** Bob**: Indeed. The actual work is done here:

include Math def frand(low, high) low + rand * (high - low) end def mkplummer(n, seed) if seed == 0 srand else srand seed end nb = NBody.new(n) nb.body.each do |b| b.mass = 1.0/n radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0) theta = acos(frand(-1, 1)) phi = frand(0, 2*PI) b.pos[0] = radius * sin( theta ) * cos( phi ) b.pos[1] = radius * sin( theta ) * sin( phi ) b.pos[2] = radius * cos( theta ) x = 0.0 y = 0.1 while y > x*x*(1.0-x*x)**3.5 x = frand(0,1) y = frand(0,0.1) end velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25) theta = acos(frand(-1, 1)) phi = frand(0, 2*PI) b.vel[0] = velocity * sin( theta ) * cos( phi ) b.vel[1] = velocity * sin( theta ) * sin( phi ) b.vel[2] = velocity * cos( theta ) end STDERR.print " actual seed used\t: ", srand, "\n" nb.acs_write end

The ` mkplummer` method takes two arguments: ` n`, which is the number
of stars in our star cluster, and ` seed`, which is the seed for the
random number generator.

If you want to create a different model, each time you invoke
` mkplummer`, you can choose `seed = 0`. That will invoke
` srand` without any argument, as you can see from the first two lines.
Ruby guarantees that such a call to ` srand` will return a different
seed each time.

However, if you do want to reproduce the same result, you can give a
specific seed by choosing `seed = k`, where ` k` is a positive
integer. Calling ` srand` with that number will actually use that
seed. This will be useful when you are debugging the ` mkplummer` code
itself, or when you are debugging another code, that takes the ` mkplummer`
results as input.

I'm a bit puzzled about one thing, though: the manual tells me that
the call `srand 0` would have the same effect as calling ` srand`
with nu arguments. This would suggest that I could replace the first
four lines of ` mkplummer` by

srand seedwithout having to include the

After seeding the random number generator, I create a new N-body system
with ` n` particles, and then I enter a loop, in which I initialize each
of those particles in turn.

** Alice**: Let's look at that loop later. I would like to see the overall
structure of the code first.

** Bob**: Okay: after the loop finishes for the last particle, I print out
the random number seed that was used, as a check to make sure that we can
later recreate the same model realization.

** Alice**: Why would you do that? You already gave the seed value as an
argument.

** Bob**: If I give a positive value, then indeed this is superfluous. But
if I give a seed value of zero, the system chooses a pseudo-random seed
for me. And if I would not echo the seed, I would not be able to reproduce
the run later on.

** Alice**: I see. But why do you call ` srand` again in the print statement?

** Bob**: Because Ruby tells me that the value that ` srand` returns is the
value of the * previous* seed. Don't ask why. It is defined that way.
So by calling ` srand` again, I get the value of the seed that I have
to give the next time, in order to reproduce the previous run, even if
at that time I gave a value zero for the second argument of ` mkplummer`.
It works: I tested it out.

** Alice**: I take your word for it. Can you show me how you invoke ` mkplummer`?

** Bob**: Here it is:

options_text= <<-END Description: Plummer's Model Builder Long description: This program creates an N-body realization of Plummer's Model. (c) 2004, Piet Hut and Jun Makino; see ACS at www.artcompsi.org The algorithm used is described in Aarseth, S., Henon, M., & Wielen, R., Astron. Astroph. 37, 183 (1974). Short name: -n Long name: --n_particles Value type: int Default value: 1 Variable name: n Print name: N Description: Number of particles Long description: Number of particles in a realization of Plummer's Model. Each particles is drawn at random from the Plummer distribution, and therefore there are no correlations between the particles. Physical Units are used in which G = M = a = 1, where G is the gravitational constant M is the total mass of the N-body system a is the structural length, with potential U(r) = GM/(r^2 + a^2)^{1/2} Short name: -s Long name: --seed Value type: int Default value: 0 Description: pseudorandom number seed given Print name: Variable name: seed Long description: Seed for the pseudorandom number generator. If a seed is given with value zero, a preudorandom number is chosen as the value of the seed. The seed value used is echoed separately from the seed value given, to allow the possibility to repeat the creation of an N-body realization. Example: |gravity> kali mkplummer1.rb -n 42 -s 0 . . . pseudorandom number seed given : 0 actual seed used : 1087616341 . . . |gravity> kali mkplummer1.rb -n 42 -s 1087616341 . . . pseudorandom number seed given : 1087616341 actual seed used : 1087616341 . . . END c = parse_command_line(options_text) mkplummer(c.n, c.seed)

It is mostly command line argument parser, using our nifty new ` Clop`
class based parser. It contains the `here document' that defines the
options, followed first by the command `parse_command_line`
that does what it says it does, and then by the command ` mkplummer`.

** Alice**: Let us now look at the inner loop of ` mkplummer`, where each
particle is given its initial values for its mass, position, and velocity.
The first line is simple:

b.mass = 1.0/n

Each particle acquires a mass that is : you are creating an equal-mass system, where the total mass is normalized to be unity.

** Bob**: Indeed. It would be possible to give a mass spectrum, of course,
but that can be done later, when we are ready for that. For now, I
just wanted to construct a minimal model.

** Alice**: Fair enough. Then you pick values for position and velocity
components. What is the basic idea here?

** Bob**: The idea is to proceed in two steps. You start by sprinkling
particles in space, as a random realization of the mass density
distribution of Plummer's model. This means that you have to be
careful to determine the radial position of the particles with the
right statistical weight. Because of spherical symmetry, the angles
can be randomly chosen from a two-dimensional spherical surface.

The second step is to give each particle a velocity with a random direction and a magnitude that is also random, but drawn from the appropriate velocity distribution at that point in space.

** Alice**: In other words, you really are sprinkling particles into
phase space, the six-dimensional space that is the direct product
of configuration space and velocity space.

** Bob**: Yes, but even though you can look at phase space as a single
six-dimensional space, that still leaves the fact that there is the
three-by-three structure left for the two quite different subspaces.

What I mean is: you have to start by picking an appropriately random point in what you call configuration space, and what is normally just called `space', containing all possible positions in three dimensions. Only after you know the position of a particle in that subspace, can you determine the potential energy of that particle. And only when you know the potential energy, can you know what type of velocities are admissible, in order to keep the particle bound and to give each velocity the correct statistical weight.

To sum up, the order of picking a point in configuration space, and then picking a point in velocity space, is important: you couldn't do it the other way around.

** Alice**: Ah, that must correspond to what mathematicians mean when they
call the velocity space the tangent bundle to the configuration space.
Each point in configuration space has its own tangent space, and it is
only because we work in Newtonian flat space that we can pretend that
phase space is a single six-dimensional space, shared by all particles.

** Bob**: Whatever. Enough mathematical terms! Let me show you how the
particles get sprinkled into normal space first, and how we then
populate the velocity space. In practice, rather than doing all the
configuration space sprinkling first, followed by a second loop in
which we populate velocity space, I find it easier to do everything in
one loop. For each particle I determine the three position coordinates,
and with that information I can then immediately give the three velocity
coordinates.

** Alice**: Okay, lets look at configuration space first. I see that you
invoke some magical expression to determine the value of the radius:

radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0)

** Bob**: Yes, this is part of the recipe provided by Aarseth *et al.*,
for choosing a correctly randomized ` radius`, the distance of the star
from the center of the star cluster.

** Alice**: Let's look at that in a moment, after we've gone through the body
of the loop first. Given the value for ` radius`, I see that you assign
the three Cartesian coordinates of the star in the usual way from spherical
coordinates ` radius`, ` theta`, and ` phi`:

b.pos[0] = radius * sin( theta ) * cos( phi ) b.pos[1] = radius * sin( theta ) * sin( phi ) b.pos[2] = radius * cos( theta )

** Bob**: Yes, and randomizing the two angular coordinates was relatively simple.
The value for , for example, is a random number
between and , since every azimuthal angle
is equally likely:

phi = frand(0, 2*PI)

By the way, here I have adapted the Ruby defined random number call
` rand` in the following way, by defining a general floating point
version ` frand`:

def frand(low, high) low + rand * (high - low) end

Since ` rand` returns a value uniformly distributed throughout the
range , `frand(a,b)` returns a value
uniformly distributed throughout the range .

** Alice**: And this is how you initialize the angle ` theta` between the
positive ` z` axis and the position vector of your star:

theta = acos(frand(-1, 1))

** Bob**: To pick a random value for , we have to make
sure that the spherical integration element
gets an equal weight for any value. In other words,
any value for should be equally likely.

Now the highest and lowest values occur for ,
along the positive * z* axis, and for , along
the negative * z* axis. So runs from
to .

All we have to do is to pick a floating point number at random, somewhere in the interval . Let's call the number . This determines a number defined as . By construction, the value of is uniformly random, as was desired for . Hence has the right distribution of values for , and we can simply take .

** Alice**: It sounds complicated when you say it in words, but I see what
you mean. Indeed, that must be correct.

** Bob**: I sometimes think that anything to do with probability gets more
confusing when you try to talk about it. Easy to fool others, and to
fool yourself for that matter. No wonder people talk about lies, damn
lies, and statistics!

** Alice**: Moving right along, we now come across something really mysterious:

x = 0.0 y = 0.1 while y > x*x*(1.0-x*x)**3.5 x = frand(0,1) y = frand(0,0.1) end velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)

** Bob**: Yeah, and fun too. But since you postponed a discussion of the
proper weighting function for choosing the radial position, we should
also postpone a discussion of what goes on here. Briefly, I am using
a rejection technique, following Aarseth *et al.* in order to
determine the magnitude ` velocity` of the velocity vector, in the last
line above.

** Alice**: Okay, I'm happy to wait till later. Now in the next lines you
repeat the same spherical distribution trick you applied in order to
find the position coordinates.

** Bob**: Indeed, and this gives me the velocity coordinates.

** Alice**: You could have put those five lines into a separate method.
Remember the DRY principle: don't repeat yourself.

** Bob**: Good point. Let me do that in the next version of this code.
This is a rather minimal one, and I can think already of several
improvements. For example, we may want to recenter the center of mass
of the star cluster we have created onto the center of the coordinates.
We can also dampen some of the fluctuations we are introducing by
layering the particles more evenly in radial bins, rather than
sprinkling every particle in space, independently of any other particle.

** Alice**: That's fine: I also like to start with a minimal script, so that
we can really test and understand its behavior, before we start adding
bells and whistles. Testing the first version of a new code is half the
work. Once you have one well-tested version, you can use that to make
comparisons, to check whether each new addition still
reproduces the old result. The hardest part is to get an initial result,
and make sure that that one is correct.

Previous | ToC | Up | Next |