Previous | ToC | Up | Next |

** Alice**: A job well done! Now every ruby program that starts with
`require "acs"` will automatically receive options for specifying
precision, indentation, and verbosity levels for reporting diagnostics
on the standard error channel as well as in the stories that appear on
the standard output channel.

** Bob**: I'm glad we have all that in place. Now, * finally* let's get back
to our task: letting some real physics come out! We've been talking about
determining the time of first binary formation.

In order to do so, we have to recognize binaries when they are formed.
Well, I added a quick-and-dirty binary finding part to the diagnostics
of our individual time step integrator. I started with `world3.rb`,
but I'm calling it `world4.rb` now, with that addition.

At the top of `world4.rb`, I'm including a file,
`binary.rb`, that contains a new class ` Binary`. Here it is:

require "nbody.rb" class Binary def initialize(body1, body2) @m1 = body1.mass @m2 = body2.mass @mass = @m1 + @m2 @reduced_mass = ( @m1 * @m2 ) / ( @mass ) @pos = (@m1*body1.pos + @m2*body2.pos)/@mass @vel = (@m1*body1.vel + @m2*body2.vel)/@mass @rel_pos = body2.pos - body1.pos @rel_vel = body2.vel - body1.vel end def rel_kinetic_energy 0.5 * @reduced_mass * @rel_vel * @rel_vel end def rel_potential_energy -( @m1 * @m2 / sqrt( @rel_pos * @rel_pos ) ) end def rel_energy rel_kinetic_energy + rel_potential_energy end def angular_momentum_squared r_cross_v = @rel_pos.cross(@rel_vel) @reduced_mass**2 * r_cross_v * r_cross_v end def semi_major_axis -( @m1 * @m2 ) / ( 2 * rel_energy ) end def eccentricity e_sq = 1 - angular_momentum_squared / ( @reduced_mass * @m1 * @m2 * semi_major_axis ) e_sq = 0.0 if e_sq < 0.0 # to avoid round-off to slightly negative numbers sqrt(e_sq) end def period 2*PI/sqrt( @mass / semi_major_axis**3 ) end end

** Alice**: I see. You create a new instance of the ` Binary` class by giving
it two bodies. The initializer then immediately determines the relative
position and velocity as well as the total mass and the reduced mass.
From that point on, you can ask the binary what its eccentricity or
semi-major axis is, or its period, or some other piece of information.
In each case, the answer is provided through a call to a method within the
` Binary` class, where the method computes the answer on the fly.

** Bob**: Indeed. I could have done a more complete job, by adding information
about all six orbital elements, preferably in various coordinate systems.
For now, though, the main information we are interested in is the hardness
of the binary, its internal energy.

** Alice**: Which you call `rel_energy` to indicate that it is the energy
in the relative motion of the two stars, independent of the energy in the
center-of-mass motion of the binary. And you provide the semi-major axis
because that is directly related to the internal energy.

** Bob**: Yes, the semi-major axis * a* is inversely proportional to the binary
binding energy `rel_energy`, for given values of the two masses.
And once I had computed * a* I thought I might as well also compute * e*,
the eccentricity of the orbit, which has a simple relation with the
relative angular momentum, as you can see in the `binary.rb` listing.

** Alice**: Now how did you use this class to report binary diagnostics within
the integrator?

** Bob**: In `world4.rb`o, in the class `WorldEra`, I have added
the following method:

def binary_diagnostics(t) v = 1 acs_log(v, take_snapshot(t).binary_diagnostics) end

** Alice**: So this method takes a snapshot at the requested time, asks the
snapshot to return the binary diagnostics, and then reports that information
with our new `acs_log` function.

** Bob**: I could of course have given the value 1 directly to the first argument
of `acs_log`, but I prefered to keep it a free variable, ` v`, which
is easier to notice and to chance later, if so desired. For the time being,
it will have the value 1, which means a verbosity level of 1.

** Alice**: Let me see whether I got all this correctly now. If the user invokes
the integrator with

kali world4.rb --verbosity 0 --acs_verbosity 2no binary diagnostics will appear on the standard error stream, but it will be written in the story of

kali world4.rb --verbosity 1 --acs_verbosity 0will show binary diagnostics will on the standard error stream, but will not write anything in the stories that appear on the standard output stream. In general,

kali world4.rb --verbosity v1 --acs_verbosity v2will only show diagnostics information on the screen when and will only add the information to the

** Bob**: That is indeed correct. And I must say, I'm really glad to have such
a versatile and general tool for direction information where it is
needed!

** Alice**: Let me look how you have implemented the real work, in the
`binary_diagnostics` method within the ` WorldSnapshot` class:

def binary_diagnostics v = 1 c = Clop.option prec = c.binary_diag_precision s = "" @body.each do |bi| @body.each do |bj| if bj.body_id > bi.body_id b = Binary.new(bi, bj) if b.rel_energy < 0 and b.semi_major_axis <= c.max_semi_major_axis s += " [#{bi.body_id}, #{bj.body_id}] : " s += sprintf("a = %.#{prec}e ; ", b.semi_major_axis) s += sprintf("e = %.#{prec}e ; ", b.eccentricity) s += sprintf("T = %.#{prec}e\n", b.period) end end end end s end

As before, you introduce the variable * v* right at the start and set it
equal to 1, which is the normal verbosity level. And then in the third
line you introduce a measure of precision, ` prec`.

** Bob**: Yes. I could have used the normal precision specified in our
standard ACS precision, which by default is 16 digits long for floating
point variables. However, it would be destracting to list the semimajor
axis of a binary in full double precision. Therefore I introduced a
special command line option to describe the binary diagnostics precision.
Here, you can find it among the ever growing list of options for
`world4.rb` :

|gravity> kali world4.rb -h Individual Time Step, Individual Integration Scheme Code -g --integration_method: Choice of integration method [default: hermite] -c --step_size_control : Determines the time step size [default: 0.01] -f --init_timescale_factor: Initial timescale factor [default: 0.01] -e --era_length : Duration of an era [default: 0.0078125] -m --max_timestep_param: Maximum time step (units dt_era) [default: 1] -d --diagnostics_interval: Diagnostics output interval [default: 1] -o --output_interval : Snapshot output interval [default: 1] -y --pruned_dump : Prune Factor [default: 0] -t --time_period : Duration of the integration [default: 10] -u --cpu_time_max : Max cputime diagnost. interval [default: 60] -i --init_out : Output the initial snapshot -r --world_output : World output format, instead of snapshot -a --shared_timesteps : All particles share the same time step -x --max_semi_major_axis: Maximum semi-major axis [default: 1.0e+30] --binary_diag_precision : Binary Diagnostics Precision [default: 4] --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)

** Bob**: This is of course rather wastful of computer time, and not very
efficient. If we would be dealing with tens of thousands of
particles, we may want to think of a more clever scheme. But for now,
this seems like the simplest way of doing things.

** Alice**: I agree. Premature optimization is the root of all evil, as we
already said! So for each pair you first check whether it is bound. Only
if the relative energy is negative do you continue and check the semi-major
axis. And only if that value is smaller than the maximum value allowed,
do you print diagnostics information.

** Bob**: Well, ` sprint` diagnostics information, you mean; the information is
printed on a string, and the string is passed back to the calling function.

** Alice**: And that calling function, also named `binary_diagnostics`
within the ` WorldEra` class, can use that string both for printing on the
standard error stream and for adding it to the total ACS wrapped output on
the standard output stream. Got ya! And in the list of options above
I already saw how the user can set the maximum semi-major axis value.

Great! It all seems to be quite transparent. Last question: how and where
does the `WorldEra#binary_diagnostics` method get invoked? Most
likely on the ` World` level, and close to the other diagnostics. But I
don't see the word binary in the listing of clas ` World`.

** Bob**: Almost correct. You are right in that the calls indeed follow the
calls to the `WorldEra#write_diagnostics` methods, literally
so, but they occur in the module ` Output`. But then the module
` Output` gets mixed into the class ` World`, so a language lawyer might
argue that the calls occur at the ` World` level. Hera are the two
places where it happens, first for normal diagnostics output:

def diagnostics(t_target, dt_dia) dia_output = false while @t_dia <= t_target @era.write_diagnostics(@t_dia, @initial_energy) @era.binary_diagnostics(@t_dia) @t_dia += dt_dia dia_output = true end dia_output end

and then for unscheduled diagnostics output, when a CPU limit is reached before a whole era is finished:

def unscheduled_diagnostics(dt_dia) t = @era.wordline_with_minimum_interpolation.worldpoint.last.time unless diagnostics(t, dt_dia) @era.write_diagnostics(t, @initial_energy, true) @era.binary_diagnostics(t) end end

** Alice**: No need for lawyers here; I think I was close enough. Now I'd
like to see your binary diagnostics in actions. How shall I begin?

** Bob**: Easy enough, why not just do what the programs suggest you do,
when you use the `---help` option.

** Alice**: I keep forgetting that we have such a nice way to do the
hand-holding, for users and developers alike! Here we are:

|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 1So you suggest I literally do that? Let me try:

|gravity> kali mkplummer.rb -n 4 -s 1 | kali world4.rb -t 1 ==> 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 = 1.0 Max cputime diagnost. interval : cpu_time_max = 60 Maximum semi-major axis : max_semi_major_axis = 1.0e+30 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 [0, 1] : a = 2.1938e+00 ; e = 6.9252e-01 ; T = 2.8872e+01 [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 [0, 1] : a = 4.9021e-01 ; e = 8.4660e-01 ; T = 3.0498e+00 [1, 2] : a = 1.8175e-01 ; e = 9.4611e-01 ; T = 6.8852e-01 ACS NBody Array body Body body[0] Vector acc 4.3680007628230150e-01 -7.7920066607005689e-01 -3.4447512499978329e-02 Array acsmixins Module acsmixins[0] Integrator_hermite Fixnum body_id 0 Float dt_param 1.0000000000000000e-02 Vector jerk 3.2148293860481963e-01 1.5455084837998969e+00 5.0309934483466201e-01 Float mass 2.5000000000000000e-01 Float maxstep 5.3298829935937153e-03 Float minstep 2.8993571645200446e-05 Float next_time 1.0026061296123183e+00 Fixnum nsteps 332 Vector pos 7.9772604070178846e-03 8.2919239954665480e-01 8.6949862920921966e-02 Float time 1.0000000000000000e+00 Vector vel -6.0338461103000307e-01 -2.5365567416813917e-01 -4.2283983072465980e-01 Body body[1] Vector acc 2.9407213696345624e+00 4.3737103517008107e+00 2.8696484413452512e+00 Array acsmixins Module acsmixins[0] Integrator_hermite Fixnum body_id 1 Float dt_param 1.0000000000000000e-02 Vector jerk -3.0679123445456916e+01 -6.0941785672989397e+01 -4.1544144235617914e+01 Float mass 2.5000000000000000e-01 Float maxstep 2.8967962731903940e-03 Float minstep 1.0126966959234096e-05 Float next_time 1.0007045202601113e+00 Fixnum nsteps 1933 Vector pos 3.4274367211420625e-01 1.1210653405826151e-01 6.9410494803967479e-03 Float time 1.0000000000000000e+00 Vector vel -1.7471717282768337e-01 -4.6103585065601543e-01 -3.2046114392320907e-01 Body body[2] Vector acc -3.4927708485553790e+00 -3.7533254304484691e+00 -2.8624660563443318e+00 Array acsmixins Module acsmixins[0] Integrator_hermite Fixnum body_id 2 Float dt_param 1.0000000000000000e-02 Vector jerk 3.0274773693388109e+01 5.9364365363757656e+01 4.1104361307013662e+01 Float mass 2.5000000000000000e-01 Float maxstep 2.8967962777788347e-03 Float minstep 1.0126966976886642e-05 Float next_time 1.0007045842393165e+00 Fixnum nsteps 1933 Vector pos 4.5401351630391895e-01 2.5523386444086688e-01 1.0661490776678197e-01 Float time 1.0000000000000000e+00 Vector vel 9.5626266521217429e-01 2.6335235690235648e-01 2.0752875999010723e-01 Body body[3] Vector acc 1.1435586922134454e-01 1.5721237097257557e-01 2.7242540252499236e-02 Array acsmixins Module acsmixins[0] Integrator_hermite Fixnum body_id 3 Float dt_param 1.0000000000000000e-02 Vector jerk 8.8699613819200801e-02 3.4954588596844444e-02 -6.0072687819262158e-02 Float mass 2.5000000000000000e-01 Float maxstep 7.8125000000000278e-03 Float minstep 1.6386823913682258e-04 Float next_time 1.0002615373733836e+00 Fixnum nsteps 135 Vector pos -8.0473444786322912e-01 -1.1965327941527619e+00 -2.0050581825457850e-01 Float time 1.0000000000000000e+00 Vector vel -1.7816088524364576e-01 4.5133917639638632e-01 5.3577221678447662e-01 String story 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 [0, 1] : a = 2.1938e+00 ; e = 6.9252e-01 ; T = 2.8872e+01 [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 [0, 1] : a = 4.9021e-01 ; e = 8.4660e-01 ; T = 3.0498e+00 [1, 2] : a = 1.8175e-01 ; e = 9.4611e-01 ; T = 6.8852e-01 Float time 1.0000000000000000e+00 SCA

** Bob**: I bet that's more than you wanted!

** Alice**: Yeah, I have to learn to tame our programs. But it is nice to
see the same diagnostics appearing on the screen and within the story
that now appears within the ACS wrapped particle output.

Let me exercise your new options. I will ask only for binaries that
have a semimajor axis smaller than 0.5, and I want to know the results
only to the first two significant digits. And I certainly don't want
all that output. I could redirect it to `/dev/null`, but let
me instead set the time of next output to, say, 1000, to make sure that
we don't see it within our run. Here goes:

|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 1 -o 1000 -x 0.5 --binary_diag_precision 2 < tmp.in ==> 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 = 1000.0 Prune Factor : prune_factor = 0 Duration of the integration : t = 1.0 Max cputime diagnost. interval : cpu_time_max = 60 Maximum semi-major axis : max_semi_major_axis = 0.5 Binary Diagnostics Precision : binary_diag_precision = 2 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.80e-01 ; e = 8.82e-01 ; T = 6.79e-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 [0, 1] : a = 4.90e-01 ; e = 8.47e-01 ; T = 3.05e+00 [1, 2] : a = 1.82e-01 ; e = 9.46e-01 ; T = 6.89e-01

Previous | ToC | Up | Next |