Previous ToC Up Next

10. Two More Code Versions

10.1. Standard Units

Alice: In order to adapt your code mkplummer2.rb to standard units, we have to change scales, both in position space and in velocity space. Let us start with positions. In our old system, we had a structural length . In that system, the virial radius was:

Let us call this quantity our scale factor. In standard units, the virial radius . This means that we have to divide the virial radius by the scale factor, when we switch from the old units to standard units.

Bob: And this means that all distances have to be scaled this way, since that is the definition of a scale factor after all. I'll copy the previous code in a file called mkplummer3.rb, and add this scaling. That will take only two lines. Here, right at the top of the mkplummer method, I'm adding this line: 

   scalefactor = 16.0 / (3.0 * PI)

Then at the place where we initialize the positions for each body, we have to divide by that scale factor, just as we had to divide the virial radius by that factor: 

     b.pos = spherical(radius) / scalefactor

Alice: That must be right. Now the velocity scaling is a bit more tricky. How shall we approach that?

Bob: Velocities come into the kinetic energy. Since the masses of the particles are not affected, the velocities must scale like the square root of the energy.

Alice: Good point! The total mass of a cluster with stars is unity, and each star therefore has a mass of in both systems, or old system of units, and the standard system of units. It is only a change in the magnitude of their velocities that can make up for a change in kinetic energy -- or in total or potential energies, for that matter, since all these quantities are related. And we know that the potential energy can only scale with the distances, since in both the old system and the new one.

Bob: In other words, we have for the potential energy of the cluster as a whole:

and

Since they must remain proportional to each other, with the virial factor of two between them, we know that

Alice: That is an extremely sloppy notation, using and as arbitrary lengths and velocities, but I see what you mean, and yes, obviously the square of the velocities transform inversely proportional to the distances.

Bob: This means that I can add the following conversion factor to the code, at the point where the velocities are calculated: 

     b.vel = spherical(velocity) * sqrt(scalefactor)

10.2. Checking Quartiles

Alice: That must be correct, but let's make sure that we get the right energy and the right quartiles again -- and this time in standard units.

Bob: Let's start with the quartiles: 

 |gravity> kali mkplummer3.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  : 259442867243987632366378221437785255631
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.478
   r(1/2) = 0.762
   r(3/4) = 1.28
 |gravity> kali mkplummer3.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  : 106794077952737035657087400769366713948
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4805
   r(1/2) = 0.7715
   r(3/4) = 1.266
 |gravity> kali mkplummer3.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  : 96818168461728251751749907043384566601
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4772
   r(1/2) = 0.7701
   r(3/4) = 1.29
Alice: These are indeed what we had derived a bit earlier, in eqs. (126) and (127):

Well done!

10.3. Checking Energy

Bob: Now let's see whether the energies come out standardized as well. That should be easy to verify: in the standard units, the total energy should be -1/4. Let's check: 

 |gravity> kali mkplummer3.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 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  : 79118596476922924266900277323025514556
   E_kin = 0.248 , E_pot =  -0.488 , E_tot = -0.24
 |gravity> kali mkplummer3.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 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  : 18939911909289684081959561570789443691
   E_kin = 0.246 , E_pot =  -0.492 , E_tot = -0.246
 |gravity> kali mkplummer3.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 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  : 257729252563068529020018340246685814098
   E_kin = 0.253 , E_pot =  -0.482 , E_tot = -0.228
Alice: Just as it should be, within the expected statistical errors of a few percent. Great! Now we can generate a truly standard realization of Plummer's model.

10.4. Quiet Start

Bob: There are a couple other features I'd like to add to the code: a quiet start, and center-of-mass adjustment. Let's start with the first one.

In my original code, each star is given a position in the star cluster independently of all other stars. If you make a small star cluster, it is possible that you wind up with a significant excess of stars, in the core, say, or equally likely a significant lack of stars, there or somewhere else.

I have found it helpful in some of my simulations to start in a more quiet way, in which the stars are initially layered, each occupying a position somewhere in their proper mass shell. In other words, you divide the N-body system as an onion into N different concentric shells, centered on the center of the cluster, and each having the same amount of mass. You then sprinkle one star into each different shell.

Alice: I'm not sure whether that really helps. After you let the system evolve for a while, you will soon develop fluctuations that obey Poissonian statistics. You won't keep your neat layering for long. As a matter of fact, within a small fraction of the crossing time, your particles will visit shells other than the ones they were born in.

Bob: That is true, but even so, I like to start with a more orderly bunch of stars. For one thing, the quartiles will come out better initially.

Alice: I'm still not convinced that it will help, but it won't hurt either. Go right ahead!

Bob: I'll copy the previous version in a new file mkplummer4.rb Now I have to be careful. Given me a minute. . . . . Ah, this should do it. Here is the new version of the mkplummer method: 

 def mkplummer(n, seed)
   if seed == 0
     srand
   else
     srand seed
   end
   scalefactor = 16.0 / (3.0 * PI)
   nb = NBody.new(n)
   cumulative_mass_min = 0                                                    
   cumulative_mass_max = 1.0/n                                                
   nb.body.each do |b|
     b.mass = 1.0/n
     cumulative_mass = frand(cumulative_mass_min, cumulative_mass_max)        
     cumulative_mass_min = cumulative_mass_max                                
     cumulative_mass_max += 1.0/n                                             
     radius = 1.0 / sqrt( cumulative_mass ** (-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)
   end
   STDERR.print "             actual seed used\t: ", srand, "\n"
   nb.acs_write
 end

10.5. Quiet Indeed

Alice: Can you point out what you changed?

Bob: There are only two places where I added something. Before entering the nb.body.each loop, I define the inner mass shell -- in fact, an inner mass sphere; the central layer is a sphere, and all the subsequent ones are shells -- as follows: 

   cumulative_mass_min = 0                                                    
   cumulative_mass_max = 1.0/n

Then, at the beginning of the loop, I pick a random value for the cumulative mass within the constraint that the value for the cumulative mass has to lie within the range present within the current mass shell: 

     cumulative_mass = frand(cumulative_mass_min, cumulative_mass_max)

As soon as I have done that, I shift the boundaries of the mass shell up by one shell, to make them ready for the next traversal of the loop: 

     cumulative_mass_min = cumulative_mass_max                                
     cumulative_mass_max += 1.0/n

The last difference is that the radius is now determined from the value for the cumulative mass that I had just found, rather than from using the random number generator here directly: 

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

Alice: How about testing the new version?

Bob: The total energy should remain unchanged, still -1/4. Here it is: 

 |gravity> kali mkplummer4.rb -n 1000 | kali energy.rb
 ==> Plummer's Model Builder <==
 Number of particles            : N = 1000
 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  : 65204294661663778940336406294820512989
   E_kin = 0.252 , E_pot =  -0.499 , E_tot = -0.246
Alice: Fair enough. And your quartiles should come out wonderfully accurate, even for, say, 100 particles, by construction.

Bob: I hope so! Let's try: 

 |gravity> kali mkplummer4.rb -n 100 | kali quartiles.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  : 290538234552319837019298862196807084315
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4762
   r(1/2) = 0.7566
   r(3/4) = 1.268
 |gravity> kali mkplummer4.rb -n 100 | kali quartiles.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  : 144809416188173600253986381061280801268
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4712
   r(1/2) = 0.7675
   r(3/4) = 1.258
 |gravity> kali mkplummer4.rb -n 100 | kali quartiles.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  : 46461356975499902229965866264015114536
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4722
   r(1/2) = 0.7668
   r(3/4) = 1.273
Indeed, as expected. And of course it will be even more accurate for 10,000 particles: 

 |gravity> kali mkplummer4.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  : 206965737790539113565393631360305885539
 The values of the three quartiles for r(M) are:
   r(1/4) = 0.4777
   r(1/2) = 0.7685
   r(3/4) = 1.281
Alice: Quiet indeed. You have managed to shove the statistical noise under the rug.
Previous ToC Up Next