./ch12.html

ToC

12.1. Units Adjustment

Bob: We now have a code that produces particle positions and velocities, drawn from a Plummer distribution function represented in proper standard units of length, time, and mass, and properly centered at the center of mass, in position as well as velocity. What more could we possibly want? I think we can call it a day.

Alice: Before doing so, there is just one thing that is still bothering me. Even though our Plummer realizations are now perfectly centered, their units are not quite right.

Bob: But I thought we had checked that? Starting with mkplummer3.rb we have made sure to use standard units.

Alice: Well, let's see what happens for really low N values:

 |gravity> kali mkplummer5.rb -n 3 -s 1 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 1
 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  : 1
   E_kin = 0.0661 , E_pot =  -0.319 , E_tot = -0.252
 |gravity> kali mkplummer5.rb -n 3 -s 2 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 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  : 2
   E_kin = 0.121 , E_pot =  -0.273 , E_tot = -0.152
 |gravity> kali mkplummer5.rb -n 3 -s 3 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 3
 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  : 3
   E_kin = 0.204 , E_pot =  -0.386 , E_tot = -0.182
 |gravity> kali mkplummer5.rb -n 3 -s 4 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 4
 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  : 4
   E_kin = 0.124 , E_pot =  -0.184 , E_tot = -0.0595

Bob: Hmm, that doesn't look like the desired total energy value of minus one quarter. And I certainly don't like the total energy to become positive! That means that some of our realizations are actually unbound!

Alice: The problem is that, even though our underlying distribution function has been scaled correctly, any small-number realization will introduce fluctuations in the actual numbers that are picked out.

Bob: Ah, of course, this is just what happened with the center of mass. Even though an ensemble of many realizations will shown the average center of mass position and velocity values to be almost zero, individual realizations need to be adjusted, as we just did. Similarly, we will have to rescale the positions and velocities, to make any single realization come out with the right energy.

Alice: Effectively, we have to rescale length and time units; length for the positions, and in addition separately time for the velocities.

12.2. Implementation

Bob: Let's create a new file. How about calling it mkplummer.rb, without a number attached to it now, to show that this will be our final version, that we can use as a work horse?

Alice: Hope springs eternal. But go ahead, we can always rename it to mkplummer.rb6, like in a few minutes.

Bob: No, I think this will be really it. Well, we will see.

Alice: We have to do two separate things. Just getting the total energy to be -1/4 is not good enough: we want to make sure that we have an accurate virial equilibrium, with a kinetic energy of +1/4 and a potential energy of -1/2. This is the reason that we have to scale length and time units separately.

Bob: In analogy with our way to adjust the center of mass, we can now adjust the units. If I understand correctly what you just said, it should be this, right?

   def adjust_units
     alpha = -epot / 0.5
     beta = ekin / 0.25
     @body.each do |b|
       b.pos *= alpha
       b.vel /= sqrt(beta)
     end
   end

Alice: That indeed looks correct. And you invoke that function after you invoke the center-of-mass adjustment.

Bob: Ah, yes, the order is important. Not for the potential energy, since that only depends on relative distances. But the kinetic energy depends on the square of the velocity differences between each particle and the center of mass. Good point.

You know, let me modularize the mkplummer method a bit further. I'm sure you would love that!

Alice: Any move to more modularity is welcome, within reason, but something tells me you won't unreasonably overshoot toward modularity.

Bob: You bet I won't. But with all the adjustments we are now making, it would seem more clear to isolate the actual sampling procedure in a separate methods plummer_sample, and to leave the rest of the administrative details to the higher-level function mkplummer, which can deal with picking the right seed, invoking the sampling function, pushing the sampling results on the stack, and doing all the final adjustments.

Alice: You're becoming a true modularity spokesman!

Bob: I'll ignore that. This is just common sense.

Alice: Let's hope it becomes more common.

12.3. Treating Even the Vacuum

Bob: Here is the top level method:

 def mkplummer(c)
   if c.seed == 0
     srand
   else
     srand c.seed
   end
   nb = NBody.new
   c.n.times do |i|
     b = plummer_sample
     b.mass = 1.0/c.n
     b.body_id = i
     nb.body.push(b)
   end
   nb.adjust_center_of_mass if c.n > 0
   nb.adjust_units if c.n > 1
   nb.acs_log(1, "             actual seed used\t: #{srand}\n")
   nb.acs_write($stdout, false, c.precision, c.add_indent)
 end

and here is where the actual sampling is done, the place where a single particle receives its initial position and velocity values, before they later will be adjusted:

 def plummer_sample
   b = Body.new
   scalefactor = 16.0 / (3.0 * PI)
   radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0)
   b.pos = spherical(radius) / scalefactor
   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)
   b.vel = spherical(velocity) * sqrt(scalefactor)
   b
 end

Alice: Ah, I see that you only invoke the adjust_units method if you have two or more particles. Of course, if you have only a single particle, then its velocity, after adjusting it to the center of mass frame, will be zero.

Bob: Yes, I didn't want to risk dividing by zero, and in any case, there is no need to scale anything for a one-particle system already at rest.

Alice: But what is the meaning of the if statement for the center-of-mass adjustment?

Bob: Oh, I thought I might as well make it work for any reasonable value, including the vacuum, with zero particles. Without that clause, for c.n = 0 the method adjust_center_of_mass would try to set the length of the position and velocity vectors of the center of mass to the same length as that of the position vector of the first particle, which would be non-existent, and so you would get an error message -- even though the rest of the body of adjust_center_of_mass would work fine for zero particles.

Alice: I like building in such careful statements. Who knows, at some point we may decide to build slews of N-body models with different particle numbers, and it is nice to have them all behave well, for any reasonable and even not so reasonable numbers, including 0. However, what will happen for negative numbers?

Bob: I haven't tested that. But when I look at mkplummer above, it would seem that the do loop never gets executed, so effectively if will still give a zero-particle system, without running into any more serious error.

Alice: That may or may not be what we want. In any case, we can discuss more careful exception handling some other time.

12.4. Bells and Whistles

Bob: I agree. What next. Ah, let me take out this quiet start business. It was an interesting idea, but as we noticed above, with center of mass adjustment it will get partly screwed up anyway. Actually, I just realized a much more important reason not to do a radial layering of particles. Pretty soon we may want to give some particles extra properties. We could introduce a mass spectrum, or promordial binaries, or what not. If we make sure that the original distribution of particles is truly random, without any layering bias, we are less likely to wind up with unrealistic distributions . . .

Alice: . . . such as having the lighter particles all in the center and the heavier ones outside. Yes, I see what you mean. Let's keep it simple and forget about being too quiet.

Bob: Finally, we can use some of the bells and whistles that we have added before, when we wrote the ACS I/O routines. There we allowed the user to specify the number of digits accuracy and indentation prefered. If you're dealing with only a thousand particles, you may not want double precision, if that will cut down your inital file length by a factor half. Let us test the options we have provided:

 |gravity> kali mkplummer.rb -h
 Plummer's Model Builder
 -n  --n_particles      : Number of particles            [default: 1]
 -s  --seed             : pseudorandom number seed given [default: 0]
 --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 3 --precision 7 --indentation 4
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 7
 Incremental indentation                : add_indent = 4
              actual seed used  : 57666404336617107274195290498234884809
 ACS
     NBody 
         Array body
             Body body[0]
                 Fixnum body_id
                     0
                 Float mass
                        3.3333333e-01
                 Vector pos
                       -1.8861446e-02   1.3663158e-01  -1.8421583e-01
                 Vector vel
                        6.9844162e-01  -3.2014228e-01   4.5152305e-01
             Body body[1]
                 Fixnum body_id
                     1
                 Float mass
                        3.3333333e-01
                 Vector pos
                        3.9659039e-01  -3.8162750e-01   7.1921257e-02
                 Vector vel
                       -7.9737169e-02   2.3731016e-01   4.9509386e-02
             Body body[2]
                 Fixnum body_id
                     2
                 Float mass
                        3.3333333e-01
                 Vector pos
                       -3.7772894e-01   2.4499591e-01   1.1229457e-01
                 Vector vel
                       -6.1870445e-01   8.2832115e-02  -5.0103244e-01
         String story 4
                          actual seed used      : 57666404336617107274195290498234884809
             
 SCA

Alice: Fine to have the extra freedom, but I doubt we'll ever use it; as soon as you start doing a run, you'll probably want to keep the output to full precision.

Bob: In any case, at least for display purposes it all fits within the good old VT100 80-column wide screen.

Alice: Which actually came from the 80-column punched card format. But you're too young to remember that.

Bob: Next time we go to a museum you can point out to me what the world was like when you grew up.

12.5. Checking the Output

Alice: Let's check again, starting with the zero-body problem, while working our way up. And we might as well use your precision cap to keep it punched card printable:

 |gravity> kali mkplummer.rb -n 0 --precision 10
 ==> Plummer's Model Builder <==
 Number of particles            : N = 0
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 10
 Incremental indentation                : add_indent = 2
              actual seed used  : 72266540839129608651576854498790047952
 ACS
   NBody 
     Array body
     String story 2
                    actual seed used    : 72266540839129608651576854498790047952
       
 SCA
 |gravity> kali mkplummer.rb -n 1 --precision 10
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 10
 Incremental indentation                : add_indent = 2
              actual seed used  : 286646769464357533434848278335949129210
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           0
         Float mass
              1.0000000000e+00
         Vector pos
              0.0000000000e+00   0.0000000000e+00   0.0000000000e+00
         Vector vel
              0.0000000000e+00   0.0000000000e+00   0.0000000000e+00
     String story 2
                    actual seed used    : 286646769464357533434848278335949129210
       
 SCA
 |gravity> kali mkplummer.rb -n 2 --precision 10
 ==> Plummer's Model Builder <==
 Number of particles            : N = 2
 pseudorandom number seed given : 0
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 10
 Incremental indentation                : add_indent = 2
              actual seed used  : 241142602265268226534038162525383782985
 ACS
   NBody 
     Array body
       Body body[0]
         Fixnum body_id
           0
         Float mass
              5.0000000000e-01
         Vector pos
              3.1472715644e-02  -6.2992319248e-02   2.3987796040e-01
         Vector vel
             -4.8829226418e-01   5.0651686309e-01  -7.0790763148e-02
       Body body[1]
         Fixnum body_id
           1
         Float mass
              5.0000000000e-01
         Vector pos
             -3.1472715644e-02   6.2992319248e-02  -2.3987796040e-01
         Vector vel
              4.8829226418e-01  -5.0651686309e-01   7.0790763148e-02
     String story 2
                    actual seed used    : 241142602265268226534038162525383782985
       
 SCA

12.6. Checking the Energy

That all works as it should. Now an energy check:

 |gravity> kali mkplummer.rb -n 0 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 0
 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  : 300258510109887149931592557709397573146
   E_kin = 0 , E_pot =  0 , E_tot = 0
 |gravity> kali mkplummer.rb -n 1 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1
 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  : 22162567356998985900389770970377908349
   E_kin = 0 , E_pot =  0 , E_tot = 0
 |gravity> kali mkplummer.rb -n 2 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 2
 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  : 242244929093852391396607793127211875559
   E_kin = 0.25 , E_pot =  -0.5 , E_tot = -0.25

Looks good so far. Even the zero-body behaves! A single particle at rest in the center has neither potential nor kinetic energy. And the two-particle realization is nicely scaled to virial equilibrium.

Bob: Let's try the same four realizations you made earlier, with the very same seeds:

 |gravity> kali mkplummer.rb -n 3 -s 1 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 1
 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  : 1
   E_kin = 0.25 , E_pot =  -0.5 , E_tot = -0.25
 |gravity> kali mkplummer.rb -n 3 -s 2 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 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  : 2
   E_kin = 0.25 , E_pot =  -0.5 , E_tot = -0.25
 |gravity> kali mkplummer.rb -n 3 -s 3 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 3
 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  : 3
   E_kin = 0.25 , E_pot =  -0.5 , E_tot = -0.25
 |gravity> kali mkplummer.rb -n 3 -s 4 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 3
 pseudorandom number seed given : 4
 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  : 4
   E_kin = 0.25 , E_pot =  -0.5 , E_tot = -0.25

Alice: Good! Believe it or not, we may not need to call this version mkplummer6.rb, after all. It may remain mkplummer.rb and become our standard initial conditions generator.