Previous | ToC | Up | Next |

** 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)

** 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

** 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

** 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

** 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

** 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.273Indeed, 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

Previous | ToC | Up | Next |