Previous | ToC | Up | Next |

** Alice**: Hi Bob! Did you recover from our mathematical adventures?

** Bob**: Well, I must say, I'm glad to be back in the normal world.
Let me show you how to deal with a * real* example of Plummer's model,
rather than with a mathematical abstraction.

** Alice**: You mean the code you showed me, which triggered our lengthy
discussion?

** Bob**: Yes, but I made one modification, following your suggestion
that I could factor out the spherical angle treatment, which was so
similar for creating the position and velocity vectors. I've named
it `mkplummer2.rb`. There is now one new method:

def spherical(r) vector = Vector.new theta = acos(frand(-1, 1)) phi = frand(0, 2*PI) vector[0] = r * sin( theta ) * cos( phi ) vector[1] = r * sin( theta ) * sin( phi ) vector[2] = r * cos( theta ) vector end

It takes only one argument, the absolute value of the vector that will
be returned by this method with random values for the two spherical
angles and . You can see how I
invoke this method from within the inner loop of the ` mkplummer` method:

nb.body.each do |b| b.mass = 1.0/n radius = 1.0 / sqrt( rand ** (-2.0/3.0) - 1.0) b.pos = spherical(radius) 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) end

** Alice**: I like this better, yes, it is more modular.

** Bob**: And a lot shorter, so you can see more clearly what is happening now.

** Alice**: just to see the raw output, can you show me what a 3-body realization
looks like, of Plummer's model?

** Bob**: Sure. Here is one:

|gravity> kali mkplummer2.rb -n 3 ==> 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 = 16 Incremental indentation : add_indent = 2 actual seed used : 84602776384174400121148759980534896287 ACS NBody Array body Body body[0] Float mass 3.3333333333333331e-01 Vector pos 1.1107754173294841e+00 8.7598177115218653e-02 1.1024225203258864e+00 Vector vel -2.8803050799931779e-01 4.2827424016947752e-01 -2.4596499224386348e-01 Body body[1] Float mass 3.3333333333333331e-01 Vector pos -2.5020638780233173e-01 4.5646218917292558e-01 -4.7024123561676934e-01 Vector vel -5.4336097846255682e-01 -4.8179216868648034e-01 6.1323589122991451e-01 Body body[2] Float mass 3.3333333333333331e-01 Vector pos 7.0982750725799604e-01 4.5418343016141471e-01 -9.8938292838395003e-01 Vector vel 3.5888119559927067e-01 1.3607894924065120e-01 -3.6610065188231811e-01 SCA

** Bob**: 42.

** Alice**: Of course, that is the answer. But what a mundane question it was!

** Bob**: Different questions can have the same answer.
Here is the first attempt:

|gravity> kali mkplummer2.rb -n 3 -s 42 | tail -2 ==> Plummer's Model Builder <== Number of particles : N = 3 pseudorandom number seed given : 42 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 : 42 1.0707835396652617e-01 3.2726929775243468e-01 2.3909866955874151e-01 SCAAnd now let me repeat the command with the same seed:

|gravity> kali mkplummer2.rb -n 3 -s 42 | tail -2 ==> Plummer's Model Builder <== Number of particles : N = 3 pseudorandom number seed given : 42 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 : 42 1.0707835396652617e-01 3.2726929775243468e-01 2.3909866955874151e-01 SCAAnd now with a different seed:

|gravity> kali mkplummer2.rb -n 3 -s 43 | tail -2 ==> Plummer's Model Builder <== Number of particles : N = 3 pseudorandom number seed given : 43 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 : 43 1.0397916156880645e-01 -1.4513970203331489e-01 1.6854199026403588e-01 SCA

** Bob**: I thought about this checking question, so I cobbled together a tool
to do a quick energy check. Here it is:

#!/usr/local/bin/ruby -w 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 def ekin # kinetic energy 0.5*@mass*(@vel*@vel) end def epot(body_array) # potential energy p = 0 body_array.each do |b| unless b == self r = b.pos - @pos p += -@mass*b.mass/sqrt(r*r) end end p end end class NBody attr_accessor :time, :body def initialize @body = [] end def ekin # kinetic energy e = 0 @body.each{|b| e += b.ekin} e end def epot # potential energy e = 0 @body.each{|b| e += b.epot(@body)} e/2 # pairwise potentials were counted twice end def write_diagnostics etot = ekin + epot print <<END E_kin = #{sprintf("%.3g", ekin)} ,\ E_pot = #{sprintf("%.3g", epot)} ,\ E_tot = #{sprintf("%.3g", etot)} END end end include Math nb = ACS_IO.acs_read(NBody) nb.write_diagnostics

This is all very similar to what we did within our N-body integrators.

** Alice**: Let's start with a few small realizations, to see how large the
fluctuations are from case to case, in the kinetic and potential energies,
as well as the total energy.

|gravity> kali mkplummer2.rb -n 10 | kali energy.rb ==> Plummer's Model Builder <== Number of particles : N = 10 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 : 314084175335806434408394307870262799293 E_kin = 0.108 , E_pot = -0.283 , E_tot = -0.175 |gravity> kali mkplummer2.rb -n 10 | kali energy.rb ==> Plummer's Model Builder <== Number of particles : N = 10 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 : 269773554575491416225742644056394692430 E_kin = 0.133 , E_pot = -0.227 , E_tot = -0.0938 |gravity> kali mkplummer2.rb -n 10 | kali energy.rb ==> Plummer's Model Builder <== Number of particles : N = 10 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 : 237787308285128236747451513485022117506 E_kin = 0.199 , E_pot = -0.371 , E_tot = -0.172

|gravity> kali mkplummer2.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 : 120675979305433593447309484955153184810 E_kin = 0.147 , E_pot = -0.279 , E_tot = -0.131 |gravity> kali mkplummer2.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 : 229279117753660015758434507167504350638 E_kin = 0.139 , E_pot = -0.286 , E_tot = -0.147 |gravity> kali mkplummer2.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 : 219014694910382831760580860722482279459 E_kin = 0.14 , E_pot = -0.284 , E_tot = -0.144

** Alice**: Yes, that is encouraging. The question is: are these magnitudes
right?

** Bob**: Well, we can check by looking up what the kinetic and potential
energy should be for a Plummer model. But something tells me that you'd
rather compute it for yourself.

** Alice**: How did you guess! But we have to compute only one quantity.
Since we start with a model that is in equilibrium, we can assume the
virial theorem to hold, so as soon as we know one of the three energy
diagnostics that you are printing out here, we can determine the other
two immediately.

** Bob**: In that case, the total potential energy of the system would seem
to be the best candidate.

** Alice**: I agree. We already know the density and the potential, so it
is just a matter of integration.

** Bob**: Let's see. The total potential energy must be the sum of the
potential energy that each mass element feels, with respect to the rest
of the star cluster:

** Alice**: Ah, but now you're double counting all pairwise interactions.
You have to divide this by a factor two:

** Alice**: I'm sure it is wrong to leave out the factor two. in the above
equation. But I like the notion of building a cluster up, bit by bit.
Given that we have spherical symmetry, we can make life easiest by letting
each radial shell fall in, as if it was prefabricated. Like building
a pre-fab house: you order the parts and put it together.

** Bob**: So every shell with destination radius and
destination thickness then has a mass
. Even if the shell has a much lower
density at first, and perhaps a much larger thickness, and certainly a
much larger radius, when it settles in it must have that amount of mass,
and since mass is conserved, this must be the mass that the pre-fab shell
had before it got compressed into place. The term pre-fab is not a very
good metaphor, perhaps, when we have to crunch the components, but
let's not worry about words for now.

At the moment that the shell reaches its proper place, the material inside that shell has already been put there, while the material outside that shell is still waiting to fall in. So during the trip from an infinite distance with zero binding energy, down to its final destination, the infalling shell acquires an amount of potential energy of

Now if we look at it this way, and repeat the procedure for all mass shells from the inside outward, we must have for the total potential energy:

And I'm sure that there is no factor two here, no matter what you may argue about pairwise interactions.
** Alice**: Why do you call this energy * top* rather than * pot*?

** Bob**: Because I want to stress the difference with the expression you wrote
down above, and you already claimed the label * pot*. So I just
reversed the letters. Besides, my way of constructing the star
cluster by letting mass rain down from infinity is surely a top-down method.

** Alice**: As you wish. And I agree that there should not be a factor two
in your case, probably because on average, each mass element sees only half
of the rest of the mass. But I must admit, I'm not totally sure that both
expressions are correct. Why don't we work them both out, and see whether
they boil down to the same result.

** Bob**: You go first, since you wrote down the first expression. I'll hand
you the expressions for the density and the potential, from way back when,
eq. (12). Here they are again:

** Bob**: I'll try. I first need to recover the expression for the accumulative
mass, eq. (26):

** Alice**: It is more complicated than my expression, alright, but an
integration by parts may help out. May I try?

** Bob**: Thanks! Now it is clear that and
are identically the same.

** Alice**: All that is left to do is to solve it. Let me try another integration
by parts:

Now how shall we tackle this integral? I wonder how we can simplify this further.

** Bob**: Now this is really simple enough for my taste. Yesterday I let you
get away with going back to square one, but if we keep doing that, we'll
never build an N-body environment.

We should be able to find this integral easily in a book with a table of integrals, or by using a symbolic package. Ah, you see, here it is, in good old Abramowitz and Stegun, right at the beginning of their list of most common integrals:

** Bob**: Substituting eq. (92) into eq. (91), we get:

** Bob**: I'm really curious to see whether my code gives the potential energy
that we have now calculated, at least within the statistical noise.
In my code, following Aarseth *et al.*, I have ,
so I should find a value of

For the kinetic and total energy, we should similarly expect:

** Alice**: And look, those are indeed the type of numbers we got.
Why don't you create a few more 1000-particle realizations.

** Bob**: Here are a few more.

|gravity> kali mkplummer2.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 : 302041592721379705309481006390526637954 E_kin = 0.142 , E_pot = -0.275 , E_tot = -0.132 |gravity> kali mkplummer2.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 : 8717335803285052460172785737126512016 E_kin = 0.145 , E_pot = -0.287 , E_tot = -0.142 |gravity> kali mkplummer2.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 : 283452213331896683662695500815454999402 E_kin = 0.144 , E_pot = -0.298 , E_tot = -0.153You're right: the numbers work out beautifully.

** Alice**: Congratulations!

Previous | ToC | Up | Next |