5. Examples (*)¶

Now that we have a reasonable idea how NEMO is structured and used, we should be ready to go through some real examples Some of the examples below are short versions of scripts available online in one of the directories (check $NEMO/scripts and $NEMOBIN). The manual pages programs(8NEMO) and intro(1NEMO) are useful to find (and cross-reference) programs if you’re a bit lost. Each program manual should also have some references to closely related programs, in the “SEE ALSO” section.

5.1. N-body experiments¶

In this section we will describe how to set up an N-body experiment, run, display and analyze it. In the first example, we shall set up a head-on collision between two spherical “galaxies” and do some simple analysis.

5.1.1. Setting it up¶

In Filestructure (*) we already used mkplummer to create a Plummer model; so here we shall use the program mkommod (“MaKe an Osipkov-Merritt MODel”) to make two random N-body realizations of a King model with dimensionless central potential $W_c = 7$ and 100 particles each. The small number of particles is solely for the purpose of getting results within a reasonable time. Adjust it to whatever you can afford on your CPU and test your patience and integrator.

1% mkommod in=$NEMODAT/k7isot.dat out=tmp1 nbody=100 seed=123

These models are produced in so-called RMS-units in which the gravitational constant $G=1$, the total mass $M=1$, and binding energy $E=-1/2$. In case you would like virial units (see also: Heggie & Mathieu, $E=-1/4$, in: The use of supercomputers in stellar dynamics ed. Hut & McMillan Springer 1987, pp.233) the models have to be rescaled using snapscale:

2% snapscale in=tmp1 out=tmp1s rscale=2 "vscale=1/sqrt(2.0)"

Note the use of the quotes in the expression, to prevent the shell to give special meaning to the parenthesis, which are shell meta characters.

The second galaxy is made in a similar way, with a different seed (the default seed=0 takes the time of the day):

3% mkommod in=$NEMODAT/k7isot.dat out=tmp2 nbody=100 seed=987

This second galaxy needs to be rescaled too, if you want virial units:

4% snapscale in=tmp2 out=tmp2s rscale=2 "vscale=1/sqrt(2.0)"

We then set up the collision by stacking the two snapshots, albeit with a relative displacement in phase space. The program snapstack was exactly written for this purpose:

5% snapstack in1=tmp1s in2=tmp2s out=i001.dat deltar=4,0,0 deltav=-1,0,0

The galaxies are initially separated by 4 unit length and approaching each other with a velocity consistent with infall from infinity (parabolic encounter). The particles assembled in the data file i001.dat are now ready to be integrated.

To look at the initials conditions we could use:

6% snapplot i001.dat xrange=-5:5 yrange=-5:5

which is displayed in Figure X

5.1.2. Integration using hackcode1¶

We then run the collision for 20 time units, with the standard N-body integrator based on the Barnes “hierarchical tree” algorithm:

7% hackcode1 in=i001.dat out=r001.dat tstop=20 freqout=2 \
   freq=40 eps=0.05 tol=0.7 options=mass,phase,phi > r001.log

The integration frequency relates to the integration timestep as freq=40, the softening length eps=0.05, and opening angle or tolerance tol=$\theta$. A major output of masses, positions and potentials of all particles is done every 1/freqout = 0.5 time units, which corresponds to about 1/5 of a crossing time. The standard output of the calculation is diverted to a file r001.log for convenience. This is an (ASCII) listing, containing useful statistics of the run, such as the efficiency of the force calculation, conserved quantities etc. Some of this information is also stored in diagnostic sets in the structured binary output file r001.dat.

As an exercise, compare the output of the following two commands:

8% more r001.log
9% tsf r001.dat | more

5.1.3. Display and Initial Analysis¶

As described in the previous subsection, hackcode1 writes various diagnostics in the output file. A summary of conservation of energy and center-of-mass motion can be graphically displayed using snapdiagplot:

10% snapdiagplot in=r001.dat

The program snapplot displays the evolution of the particle distribution, in projection (in the yt package this is called a phase plot):

11% snapplot in=r001.dat

Depending on the actual graphics (yapp) interface snapplot was compiled with, you may have to hit the RETURN key, push a MOUSE BUTTON or just WAIT to advance from one to the next frame.

The snapplot program has a very powerful tool built into it which makes it possible to display any projection the user wants.

As an example consider:

12% snapplot in=r001.dat xvar=r yvar="x*vy-y*vx" xrange=0:10 yrange=-2:2 \
             "visib=-0.2<z&&z<0.2&&i%2==0"

plots the angular momentum of the particles along the z axis, $J_z = x.v_y - y.v_x$ , against their radius, $r$, but only for the even numbered particles, (i%2==0) within a distance of 0.2 of the X-Y plane ($-0.2<z && z<0.2$). Again note that some of the expressions are within quotes, to prevent the shell of giving them a special meaning.

The xvar, yvar and visib expressions are fed to the C compiler (during runtime!) and the resulting object file is then dynamically loaded into the program for execution. The expressions must contain legal C expressions and depending on their nature must return a value in the context of the program. {it E.g.} {tt xvar} and {tt yvar} must return a real value, whereas {tt visib} must return a boolean (false/true or 0/non-0) value. This should be explained in the manual page of the corresponding programs.

In the context of snapshots, the expression can contain basic body variables which are understood to the bodytrans(3NEMO) routine. The real variables {tt x, y, z, vx, vy, vz} are the cartesian phase-space coordinates, {tt t} the time, {tt m} the mass, {tt phi} the potential, {tt ax,ay,az} the cartesian acceleration and {tt aux} some auxiliary information. The integer variables are {tt i}, the index of the particle in the snapshot (0 being the first one in the usual C tradition) and {tt key}, another spare slot.

For convenience a number of expressions have already been pre-compiled (see also Table ref{t:bodytrans}), e.g. the radius r= $\sqrt{x^2+y^2+z^2}$ = sqrt(x*x+y*y+z*z), and velocity v = $\sqrt{v_x^2+v_y^2+v_z^2}$ = sqrt(vx*vx+vy*vy+vz*vz). Note that {tt r} and {tt v} themselves cannot be used in expressions, only the basic body variables listed above can be used in an expression.

When you need a complex expression that has be used over and over again, it is handy to be able to store these expression under an alias for later retrieval. With the program {tt bodytrans} it is possible to save such compiled expressions object files under a new name.

Some precompiled *bodytrans* expressions¶
name	type	expression
0	int	`0`
1	int	`1`
i	int	`i`
key	int	`key` (see also real version below)
0	real	`0.0`
1	real	`1.0`
r	real	`sqrt(xx+yy+z*z)`
		or: $\|\rvec\|$
ar	real	`sqrt(xax+yay+zaz)/sqrt(xx+yy+zz)`
ar	real	`(xax+yay+zaz)/sqrt(xx+yy+zz)` or: (rvec$cdot$avec)/$\|$rvec$\|$ \
aux	real	`aux`
ax	real	`ax`
ay	real	`ay`
az	real	`az`
etot	real	`phi+0.5(vxvx+vyvy+vzvz)` or: $phi$ + vvec$^2$/2 \
i	real	`i`
jtot	real	`sqrt(sqr(xvy-yvx)+sqr(yvz-zvy)+sqr(zvx-xvz))` or: $\|$rvec$times$vvec$\|$ \
key	real	`key` (see also int} version above)
m	real	`m`
phi	real	`phi`
r	real	`sqrt(xx+yy+z*z)`
		or: $\|$rvec$\|$ \
t	real	`t`
v	real	`sqrt(vxvx+vyvy+vz*vz)`
		or: $\|$vvec$\|$
vr	real	`(xvx+yvy+zvz)/sqrt(xx+yy+zz)`
		or: or: (rvec$cdot$vvec)/$\|$rvec$\|$
vt	real	`sqrt((vxvx+vyvy+vzvz)-sqr(xvx+yvy+zvz)/(xx+yy+z*z))`
		or: $sqrt{}$(vvec$^2$-(rvec$cdot$vvec)$^2$/$\|$rvec$\|^2$)\
vx	real	`vx`
vy	real	`vy`
vz	real	`vz`
x	real	`x`
y	real	`y`
z	real	`z`
glon	real	$l$ , `atan2(y,x)*180/PI` [-180,180]
glat	real	$b$ , `atan2(z,sqrt(xx+yy))*180/PI` [-90,90]
mul	real	`(-vx\sin{l} + vx\cos{l})/r`
mub	real	$(-vx\cos{l}\sin{b} - vy\sin{l}\sin{b}+vz\cos{b})/r$
xait	real	Aitoff projection x [-2,2] T.B.A.
yait	real	Aitoff projection y [-1,1] T.B.A.
ra	real	Right Ascension, converted from radians to degrees (-x)
dec	real	Declination, converted from radians to degrees (y)

As usual an example:

13% bodytrans expr="x*vy-y*vz" type=real file=jz

saves the expression for the angular momentum in a real valued bodytrans expression file, {tt btr_jz.o} which can in future programs be referenced as {tt expr=jz} (whenever a real-valued bodytrans expression is required), {it e.g.}

14% snapplot i001.dat xvar=r yvar=jz xrange=0:5

Alternatively, one can handcode a {it bodytrans} function, compile it, and reference it locally. This is useful when you have truly complicated expressions that do not easily write themselves down on the commandline. The $(x,y)$ AITOFF projection are an example of this. For example, consider the following code in a (local working directory) file {tt btr_r2.c}:

#include <bodytrans.h>

real btr_r2(b,t,i)
Body *b;
real t;
int  i;
{
    return sqrt(x*x + y*y);
}

By compiling this:

15% cc -c btr_r2.c

an object file {tt btr_r2.o} is created in the local directory, which could be used in any real-valued bodytrans expression:

16% snapplot i001.dat xvar=r2 yvar=jz xrange=0:5

For this your environment variable {bf BTRPATH} must have been set to include the local working directory, designated by a dot. Normally your NEMO system manager will have set the search path such that the local working directory is searched before the system one (in {tt $NEMOOBJ/bodytrans}).

5.1.4. Advanced Analysis¶

5.1.5. Generating models¶

5.1.6. Using Unix pipes¶

In most cases a NEMO file can be used in a pipe (usually via in=- and out=-), therefore limiting the need to write files. Here is an example of plotting the measured and expected surface brightness of a homogeneous sphere of 1,000,000 particles with unit mass and unit radius:

 % mkconfig - 1000000 ball seed=0 |\
     snapgrid - - nx=800 ny=800  |\
     ccdellint - 0:1.1:0.01 inc=0 out=- |\
     ccdprint - x=  newline=t label=x |\
     tabmath - - '1.5/pi*sqrt(1-%1**2)' |\
     tabplot -  1 2,3 color=2,3 line=0,0,1,1 point=2,0.1,0,0 \
       xlab="Radius" ylab="Surface Brightness" headline="mkconfig shape=ball" yapp=ball.png/png

A few comments on the highlighted lines: In line 3 the out= keyword is not the second keyword, hence the explicit way it was written with the out=-. In line 5 the expected surface brightness expression is added as the 3rd column to the table in the pipe, then passed on for a quick and dirty plot (shown below).

5.1.7. Handling large datasets¶

One of NEMOs weaknesses is also it’s strong point: programs must generally be able to fit all their data in (virtual) memory. Although programs usually free memory associated with data that is not needed anymore, there is a very clear maximum to the number of particles it can handle in a snapshot. By defaultfootnote{one can recompile NEMO in single precision and define {tt body.h} with less wastefull members} a particle takes up about 100 bytes, which limits the size of a snapshots on workstations.

It may happen that your data was generated on a machine which had a lot more memory then the machine you want to analyze your data on. As long as you have the diskspace, and as long as you don’t need programs that cannot operate on data in serial mode, there is a solution to this problem. Instead of keeping all particles in one snapshot, they are stored in several snapshots of (equal number of) bodies, and as long as all snapshots have the same time and are stored back to back, most programs that can operate serially, will do this properly and know about it. Of course it’s best to split the snapshots on the machine with more memory

% snapsplit in=run1.out out=run1s.out nbody=10000

If it is just one particular program (e.g. snapgrid that needs a lot of extra memory, the following may work:

% snapsplit in=run1.out out=- nbody=1000 times=3.5 |\
    snapgrid in=- out=run1.ccd nx=1000 ny=1000 stack=t

Using tcppipe(1NEMO) one can also pipe data from the machine with large memory (A) to your workstation with smaller memory (B), as can be demonstrate with the following code snippet:

A% mkplummer - 1000 | tcppipe
B% tcppipe A | tsf -

5.2. Images¶

Todo

examples/Images

5.3. Tables¶

Todo

examples/Tables

5.4. Potential¶

Todo

examples/Potential

5.5. Orbits¶

In this section we will describe how to integrate individual stellar orbits, display and analyze them. Be aware that although 3D orbits can be computed the number of utilities to analyze them is rather limited.

Orbits are normally stored in datafile (see also {it orbit(5NEMO)}), and a close conceptual relationship exists between a (single-particle type) {bf snapshot} and an {bf orbit}: an orbit is an ordered series of phase-space coordinates whereas a snapshot is a series of particles with no particular order, but all at the same time.

Since orbits will be computed in an analytical potential, we assume for the remainder of this section that you have familiarized yourself with how to supply potentials to orbit integrator programs. They all share the same triple potname=, potpars=, potfile= keyword interface, as described in Section ref{s:potential}. Many examples of the tricky {tt potpars=} keyword are given in Appendix ref{a:potential}.

5.5.1. Initializing¶

There are a few programs with which orbits can be initialized:

mkorbit is the most straightforward program. You can give simply give it all 6 phase space coordinates, and an orbit file consisting of this one point is generated. It is also possible to give the potential in which the particle is to move, and 5 phase space coordinates plus the energy, or even 4 phase space coordinates and the energy plus the total angular momentum or angular momentum along the Z axis (for axisymmetric systems).

Let’s start with an example of creating a simple orbit by itself with no associated potential.

% mkorbit out=orb1 x=1 y=0 z=0 vx=0 vy=0.2 vz=0
### Warning [mkorbit]: Potential potname= not used; set etot=0.0
pos: 1.000000 0.000000 0.000000
vel: 0.000000 0.200000 0.000000
etot: 0.000000
lz=0.200000

% tsf orb1
char History[59] "mkorbit out=orb1 x=1 y=0 z=0 vx=0 vy=0.2 vz=0 VERSION=3.2b"
set Orbit
  set Parameters
    int Ndim 03
    double Mass 1.00000
    double IOM[3] 0.00000 0.200000 0.00000
    int Nsteps 01
  tes
  set Potential
  tes
  set Path
    double TimePath[1] 0.00000
    double PhasePath[1][2][3] 1.00000 0.00000 0.00000 0.00000 0.200000 0.00000
  tes
tes

perorb is a program that for given initial conditions (similar to the ones described in {tt mkorbit} above) attempts to calculate periodic orbits in that potential. The output file will be a file with one (or more) orbits. This is a bit of an advanced program, and will not be covered here.
stoo is a program that can take a particle position from a snapshot, and turn it into an orbit. For example, sampling some initial conditions from the positions of stars in a Plummer sphere, we could use the following small C-shell code to find some statistical properties of this selected set of orbitsfootnote{For the careful reader: {tt mkplummer} and {tt potname=plummer} actually have different units, and as such this experiment is not properly set up.}

mkplummer out=p100 nbody=p100
foreach i (`nemoinp 0:100:10`)
    stoo in=p100 out=orb$i ibody=$i
    orbint orb$i orb$i.out 10000 0.01 10000 potname=plummer
    orbstat orb$i.out
end

The reverse program, {tt otos} turns an orbit into a snapshot, and may come in handy since the snapshot package has far more advanced analysis programs.

5.5.2. Integration¶

orbint integrates orbits from given initial conditions. If the input orbit has more than 1 step, the last step is taken as the initial conditions. Although the {tt potname=, potpars=, potfile=} keywords can be given, if the input orbit contains…

perorb finds periodic orbits, and stores a full period which should close the orbit. This program finds periodic orbits in the XY plane (i.e. currently it will only find 2D orbits) by searching for the centers of invariant curves in the surface of section.
henyey also finds periodic orbits, but uses Henyey’s methodfootnote{see also van Albada & Sanders, (1982, MNRAS, 201, 303)}. This program has however not been released to the public version of NEMO, and in fact it seem the source code was lost.

5.5.3. Display¶

orbplot is the only orbit plotting program we currently have. For more sophisticated display {tt tabplot} and/or {tt snapplot} would have to be used after transforming the data. Also {tt snapplot} uses the powerful {it bodytrans} expression parser to plot arbitrary body related expressions, although {tt orbplot} can handle both {tt x, y, z} and {tt vx, vy, vz} for the {tt xvar=} and {tt yvar=} keywords. An example of the output of {tt orbplot} is given in Figure ref{f:orbit1}.

5.5.4. Analysis¶

orbstat is an example of a simple program that reads orbits, and displays statistics of it’s 2D (x-y-) coordinates: maximum extent, as well as statistics of the angylar momentum. This program is not suited for 3D orbits yet.

% orbint orb1 orb1.long
% orbstat orb1.out
# T     E       x_max   y_max   u_max   v_max   j_mean  j_sigma
1000 -0.687107 1 0.999958 0.746764 0.746611 0.2 3.83111e-09

orbfour performs a variety of fourier analysis on the coordinates

% orbint orb1 orb1.long 100000 0.01 10000 10 plummer
INIPOTENTIAL Plummer: [3d version]
Pattern speed=0
0.000000 0.020000 -0.707107     -0.6871067811865
100.000000 0.277794 -0.964901     -0.6871067811856
200.010000 0.020912 -0.708019     -0.6871067812165
300.020000 0.271222 -0.958329     -0.6871067812194
400.030000 0.023376 -0.710483     -0.6871067812465
500.040000 0.259253 -0.946360     -0.6871067812551
600.050000 0.027415 -0.714522     -0.6871067812765
700.060000 0.242979 -0.930086     -0.6871067812904
800.070000 0.033056 -0.720163     -0.6871067813065
900.080000 0.223694 -0.910801     -0.6871067813241
Energy conservation: 2.00138e-10
% orbfour orb1.long amode=t
<R> N A0 A1 A2 A3 A4 B1 B2 B3 B4
1 10001 0.000360461 0.334714     0.000150399 -0.000472581 -0.000158864
                   -0.000667155  0.000228086 -0.000725406  0.000103029

% orbfour orb1.long amode=f
<R> N C0 C1 P1 C2 P2 C3 P3 C4 P4
1 10001 0.000360461
        0.334715      -0.114202
        0.000273209   56.5992
        0.000865763 -123.083
        0.000189349  147.035

orbsos computes surface of section coordinates. Since this program does not plot, but produces a simple ascii table, you can pipe the output into tabplot:

% orbsos orb1.long y | tabplot - 3 4  xlab=Y ylab=VY
% orbsos orb1.long x | tabplot - 3 4  xlab=X ylab=VX

will plot either a Y-VY or X-VX surface of section.

orbdim computes the dimensionality of an orbit, i.e. how many integrals of motions it has. Although it requires very long integration times to accurately compute this, it is completely automatic, and does not require an analysis like that for a surface of section (which is also graphic). It is based on an interesting paper by Carnevali & Santangelo (1984, ApJ 281 473-476).
otos transforms an orbit back into a snapshot, thereby giving you the much richer set of analysis tools that are available for {it snapshot}’s.

5.6. Exchanging data¶

Todo

examples/Exchanging Data

5.7. Potential(tmp)¶

Here we list some of the standard potentials available in NEMO, in a variety of units, so not always $G=1$!

Recall that most NEMO program use the keywords potname= for the identifying name, potpars= for an optional list of parameters and potfile= for an optional text string,for example for potentials that need some kind of text file. The parameters listed in potpars= will always have as first parameter the pattern speed in cases where rotating potentials are used. A Plummer potential with mass 10 and core radius 5 would be hence be supplied as: potname=plummer potpars=0,10,5. The plummer potential ignored the potfile= keyword.

plummer Plummer potential (BT, pp.42, eq. 2.47)

\[\Phi = - { M \over { {(r_c^2 + r^2)}^{1/2} } }\]

$\Omega_p$ : Pattern Speed (always the first parameter in potpars=)
$M$ : Total mass (clearly G=1 here)
$r_c$ : Core radius

potname=plummer potpars= $\Omega_p$, $M$, $r_c$

5.8. Images¶

5.8.1. Initializing Images¶

There are a few programs with which images can be initialized:

ccdmath is the most straightforward program. Here is an example of creating an image from scratch:

% ccdmath out=ccd1 fie=%x+%y size=2,4
Generating a map from scratch

% tsf ccd1
set Image
  set Parameters
    int Nx 2
    int Ny 4
    int Nz 1
    double Xmin 0.00000
    double Ymin 0.00000
    double Zmin 0.00000
    double Dx 1.00000
    double Dy 1.00000
    double Dz 1.00000
    double MapMin -4.00000
    double MapMax 0.00000
    int BeamType 0
    double Beamx 0.00000
    double Beamy 0.00000
    double Beamz 0.00000
    double Time 0.00000
    char Storage[5] "CDef"
  tes
  set Map
    double MapValues[2][4] -4.00000 -3.00000 -2.00000 -1.00000
      -3.00000 -2.00000 -1.00000 0.00000
  tes
tes

% ccdprint ccd1 x= y= label=x,y
 Y\X 0 1

3  -1 0
2  -2 -1
1  -3 -2
0  -4 -3

snapgrid converts a snapshot to an image.
fitsccd converts a FITS file to an image. The inverse of this, ccdfits also exists.

nx,ny ->    data[nx][ny]

e.g.  ccdmath out=ccd1   nx=10 ny=5
      gives       double MapValues[10][5]


ccdmath "" - %x 3,2 | tsf - margin=100 | grep MapVal
      MapValues[3][2] -2.00000 -2.00000 -1.00000 -1.00000 0.00000 0.00000
ccdmath "" - %y 3,2 | tsf - margin=100 | grep MapVal
      MapValues[3][2] -2.00000 -1.00000 -2.00000 -1.00000 -2.00000 -1.00000

5.8.2. Galactic Velocity Fields¶

As an example, a special section is devoted here to the analysis of galactic velocity fields.footnote{In this example shell variables such as r=$(nemoinp 0:60) have been replaced with the more portable macro files like @tmp.r. Although the example uses 0:60 and works fine in the shell the example was used under, increasing the number to 256 would fail because of overflowing the maximum characters allowed on the commandline}

The following programs are available:

    ccdvel          create a model velocity field, from scratch
    rotcur          tilted ring model velocity field fitting
    rotcurshape     annulus rotation curve shape fitting to a velocity field
    ccdmath         perform math on images, or use math to create images
    ccdplot         plot (contour/greyscale) an image
    ccdprint        print out pixel values in an imamge


% nemoinp 0:60 > tmp.r
% tabmath tmp.r - "100*%1/(20+%1)" all > tmp.v
% ccdvel out=map1.vel rad=@tmp.r vrot=@tmp.v pa=30 inc=60
% rotcurshape in=map1.vel radii=0,60 pa=30 inc=60 vsys=0 units=arcsec,1 \
              rotcur1=core1,100,20,1,1 tab=-


% ccdmath out=map0.vel fie=0 size=128,128
% rotcurshape map0.vel 0,40 30 45 0 blank=-999 resid=map2.vel \
              rotcur1=plummer,200,10,0,0 fixed=all units=arcsec,1

Since rotcurshape computes a residual velocity field, one can easily create nice model velocity fields from any selected shape by fitting a rotation curve shape to a velocity field of all 0s and keeping all parameters fixed to the requested values:

% ccdmath out=map0.vel fie=0 size=128,128
% rotcurshape map0.vel 0,40 30 45 0 blank=-999 resid=map.vel \
           rotcur1=plummer,200,10,0,0 fixed=all units=arcsec,1
% ccdplot map.vel -100:100:10 blankval=0 cmode=1