2. A Minimal Code

2.1. A Classic Recipe

2.2. Classes

 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.

2.3. Where the Work is Done

 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 seed

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?

2.4. The Driver

 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.

2.5. The Basic Idea

     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.

2.6. Sprinkling Particles in Space

     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!

2.7. Populating Velocity Space

     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