Previous ToC Up Next

4. Measuring Binary Formation Times

4.1. Getting Blunt

Bob: Enough with all these preparations! My fingers are itching to do some lab experiments.

Alice: Fair enough. We have already decided to find the typical time of first binary formation, starting from a Plummer model. We have learned how to let Ruby run a Ruby program, and we have taught our world4.rb work horse how to identify binaries. Time to put it all together.

Bob: You bet. And this time I will make sure that we don't get drawn into niceties and generalities. Let's be blunt and get results first.

Alice: No more mister Nice Guy, you mean?

Bob: I'm happy to be nice again in a little while, but for now, let's get moving. First question: how to do a single run, and determine the time of first binary formation.

Alice: I suppose we can let world4.rb run for a while, then analyse the output, restart world4.rb again if no binary of the right hardness has been formed, and so on.

Bob: Too nice. Too complicated. Later on, fine, but let's make a shortcut. How about letting world4.rb run and run until the first binary appears, and then just kill the program?

Alice: You are getting blunt, aren't you!

Bob: I sure am. Let me try something. How about this, as a lovely short test program? I'll call it test1.rb: 

 #!/usr/local/bin/ruby -w
 
 require "acs.rb"
 
 options_text= <<-END
 
   Description: Find the first reported binary, print its story and exit
   Long description:
     This program accepts a stream of Nbody snapshots, and checks the story of
     each one until it finds the first reported binary.  It then prints that
     story on the standard output, and exits.
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
 
         kali mkplummer.rb -n 4 -s 1 | kali world4.rb -t 1000 -x 0.25 | kali #{$0}
 
 
   END
 
 clop = parse_command_line(options_text)
 
 while nb = ACS_IO.acs_read
   raise "class #{nb.class} is not NBody" unless nb.class == NBody
   nb.story.each_line do |line|
     if line =~ / :  a = /
       print nb.story
       exit
     end
   end
 end

Alice: A lovely program to kill our integrator? Funny choice of phrase.

Bob: But that's what it does. It reads in whatever our integrator spits out, and as soon as it finds a diagnostics line reporting a binary being formed, it exits.

Alice: And by doing so, it breaks the pipe between the integrator and the test program, rudely stopping the integrator in its tracks. This is about the most unelegant way of programming I've every seen!

Bob: So be it. Let me show you in more detail, perhaps you will then appreciate the simplicity, if not the elegance, behind the madness. See, the particle output from world4.rb gets piped into test1.rb which then reachs each story line by line. If the story contains a line that contains a semicolon followed by two spaces and then a =, it assumes that it has found a report of a binary that has been formed.

The first time it finds such a binary, it halts execution and forces the whole pipeline to end its operation. So it is up to the user to give the correct value to the option in world4.rb that governs the maximum semi-major axis that triggers binary reporting.

In other words, if you want to detect when the first binary is formed that has a semi-major axis that is smaller than 0.25, you should invoke the integrator as world4.rb -x 0.25. Wider binaries will then go unreported, and only binaries of interest will trigger the collapse of the pipeline.

4.2. Lots of Output

Alice: Soon I'll insist to clean up this criminal way of killing integrators. But for now, let's see whether it all works.

Bob: It does, and I'll show you. See, I even gave it a suggestion how to get started: 

 |gravity> kali world4.rb ---help
 
     This program evolves an N-body code with a fourth-order Hermite Scheme,
     using individual time steps.  Note that the program can be forced to let
     all particles share the same (variable) time step with the option -a.
 
     This is a test version, for the ACS data format
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
     kali mkplummer.rb -n 4 -s 1 | kali world4.rb -t 1
And here we go: 

 |gravity> kali mkplummer.rb -n 4 -s 1 | kali world4.rb -t 1000 -x 0.25 | kali test1.rb
 ==> Find the first reported binary, print its story and exit <==
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 ==> Plummer's Model Builder <==
 Number of particles            : N = 4
 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
 ==> Individual Time Step, Individual Integration Scheme Code <==
 Choice of integration method   : integration_method = hermite
 Determines the time step size  : dt_param = 0.01
 Initial timescale factor       : init_timescale_factor = 0.01
 Duration of an era             : dt_era = 0.0078125
 Maximum time step (units dt_era): dt_max_param = 1.0
 Diagnostics output interval    : dt_dia = 1.0
 Snapshot output interval       : dt_out = 1.0
 Prune Factor                   : prune_factor = 0
 Duration of the integration    : t = 1000.0
 Max cputime diagnost. interval : cpu_time_max = 60
 Maximum semi-major axis                : max_semi_major_axis = 0.25
 Binary Diagnostics Precision   : binary_diag_precision = 4
 Screen Output Verbosity Level  : verbosity = 1
 ACS Output Verbosity Level     : acs_verbosity = 1
 Floating point precision       : precision = 16
 Incremental indentation                : add_indent = 2
 t = 0 (after 0, 0, 0 steps <,=,> t)
           E_kin = 0.25 ,     E_pot =  -0.5 ,      E_tot = -0.25
           E_tot - E_init = 0
           (E_tot - E_init) / E_init = -0
   [1, 2] :  a = 1.8007e-01 ; e = 8.8222e-01 ; T = 6.7900e-01
 t = 1 (after 4329, 0, 6 steps <,=,> t)
           E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
           E_tot - E_init = 1.64e-10
           (E_tot - E_init) / E_init = -6.58e-10
   [1, 2] :  a = 1.8175e-01 ; e = 9.4611e-01 ; T = 6.8852e-01
 t = 2 (after 7067, 0, 6 steps <,=,> t)
           E_kin = 0.171 ,     E_pot =  -0.421 ,      E_tot = -0.25
           E_tot - E_init = -9.53e-10
           (E_tot - E_init) / E_init = 3.81e-09
   [1, 2] :  a = 1.8109e-01 ; e = 9.2257e-01 ; T = 6.8479e-01
 t = 0 (after 0, 0, 0 steps <,=,> t)
           E_kin = 0.25 ,     E_pot =  -0.5 ,      E_tot = -0.25
           E_tot - E_init = 0
           (E_tot - E_init) / E_init = -0
   [1, 2] :  a = 1.8007e-01 ; e = 8.8222e-01 ; T = 6.7900e-01
 t = 1 (after 4329, 0, 6 steps <,=,> t)
           E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
           E_tot - E_init = 1.64e-10
           (E_tot - E_init) / E_init = -6.58e-10
   [1, 2] :  a = 1.8175e-01 ; e = 9.4611e-01 ; T = 6.8852e-01
Alice: It works all right! That much I have to credit you for. But I can't say I like it.

Bob: We can do a mopping up operation later. Let's move on.

4.3. Time of Binary Formation

Alice: Okay. So you have managed to wrestle the story with the binary formation information out of the integrator. What we really need is not the story, but the time of first binary formation, given the specified limit on the size of the binary. How do you propose to do that?

Bob: Simple, just ask for it, that's the beauty of our ACS data format. Here it is, as test2.rb 

 #!/usr/local/bin/ruby -w
 
 require "acs.rb"
 
 options_text= <<-END
 
   Description: Find when the first binary froms, print that time and exit
   Long description:
     This program accepts a stream of Nbody snapshots, and checks the story of
     each one until it finds the first reported binary.  It then prints the
     time at which that binary was first found, on the standard output, and
     exits.
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
 
         kali mkplummer.rb -n 4 -s 1 | kali world4.rb -t 1000 -x 0.25 | kali #{$0}
 
 
   END
 
 class NBody
   attr_reader :time
 end
 
 clop = parse_command_line(options_text)
 
 while nb = ACS_IO.acs_read
   raise "class #{nb.class} is not NBody" unless nb.class == NBody
   nb.story.each_line do |line|
     if line =~ / :  a = /
       print nb.time, "\n"
       exit
     end
   end
 end

Alice: You see the payoff that we got for not being blunt, and taking the time to set up a well crafted ACS IO mechanism.

Bob: Noted. Let's see whether this works, but let me spare you the lots and lots of output that we got in the previous run. I'll just report the end of the story. Taking the same initial conditions, here we go once again: 

 |gravity> kali mkplummer.rb -n 4 -s 1 > tmp.in
 ==> Plummer's Model Builder <==
 Number of particles            : N = 4
 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
 |gravity> (kali world4.rb -t 1000 -x 0.25 < tmp.in | kali test2.rb > tmp.out) >& tmp.err
Alice: Quite a relief, to see such a silent run! What did you catch in the standard output?

Bob: Hopefully the time of first binary formation: 

 |gravity> cat tmp.out
 1.0
Alice: Congratulations! Of course, it is not yet the actual time of binary formation. It is the first time the binary has been spotted in a story, as part of the particle output. So if we do this type of output once every time unit, a binary forming at, say, time 3.14 will be noticed only at time 4.0.

Bob: Correct. And there is no reason not to do an output more often, as soon as we find a good way to suppress all the unnecessary visible output.

Alice: did the standard error channel give the write diagnostics?

Bob: Presumably, but let's check: 

 |gravity> tail tmp.err
 t = 1 (after 4329, 0, 6 steps <,=,> t)
           E_kin = 0.313 ,     E_pot =  -0.563 ,      E_tot = -0.25
           E_tot - E_init = 1.64e-10
           (E_tot - E_init) / E_init = -6.58e-10
   [1, 2] :  a = 1.8175e-01 ; e = 9.4611e-01 ; T = 6.8852e-01
 t = 2 (after 7067, 0, 6 steps <,=,> t)
           E_kin = 0.171 ,     E_pot =  -0.421 ,      E_tot = -0.25
           E_tot - E_init = -9.53e-10
           (E_tot - E_init) / E_init = 3.81e-09
   [1, 2] :  a = 1.8109e-01 ; e = 9.2257e-01 ; T = 6.8479e-01
Looks like it!

4.4. A Test with One Run

Alice: What is next? I presume you're going to generate these times of first binary formation user a higher-level Ruby program?

Bob: Exactly. And we know how to do that, as we have demonstrated with our test program test0.rb. We can apply the same procedure. Let me call the new program test3.rb.

Just to stay on the safe side, let me ask test3.rb to orchestrate only a single run, for now, using test2.rb to produce the time in which the first binary forms in that run. Once we know how to do that, we can ask it to perform a whole series of runs.

How about this: 

 #!/usr/local/bin/ruby -w
 
 require "acs.rb"
 
 options_text= <<-END
 
   Description: Gathers statistics for binary formation times
   Long description:
     This program runs a number of simulations, starting with a Plummer model,
     and reports some statistics concerning the times of first binary formation.
     NOTE: TEST VERSION: prints only the result from one run
 
     (c) 2005, Piet Hut, Jun Makino; see ACS at www.artcompsi.org
 
     example:
 
         kali #{$0} -n 4 -s 1 -x 0.25
 
 
   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.
 
     Standard Units are used in which G = M = 1 and E = -1/4, where
       G is the gravitational constant
       M is the total mass of the N-body system
       E is the total energy of the N-body system
 
 
   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
         . . .
 
 
   Short name:           -x
   Long name:            --max_semi_major_axis
   Value type:           float
   Default value:        #{VERY_LARGE_NUMBER}
   Description:          Maximum value of semi major axis
   Variable name:        max_semi_major_axis
   Long description:
     This option allows the user to limit the number of binaries detected
     by discarding binaries with a semi-major axis larger than the specified
     number.  This is useful in situation such as the initial state for a cold
     collapse situation, where every star is formally bound to every other star.
 
 
   END
 
 c = parse_command_line(options_text)
 
 print `kali mkplummer.rb -n #{c.n} -s #{c.seed} --verbosity 0 | kali world4.rb \
 -t 1000 -x #{c.max_semi_major_axis} --verbosity 0 | kali test2.rb --verbosity 0`

Alice: Yes, that should do just the type of thing we've seen before. And by setting the normal verbosity level to 0 you surpress the screen output. Good! That will clean up things. And since the acs_verbosity will retain its default value of 1, the binary diagnostics will still be passed on to the story in the particle output, so that test2.rb can use it to produce the final number.

Bob: Let's see whether it indeed does what it should do, once more with the same numbers as before: 

 |gravity> test3.rb -n 4 -s 1 -x 0.25
 test3.rb: 許可がありません.
Alice: Wonderful. Our first automatically generated lab result number.

Bob: No stopping us now!
Previous ToC Up Next