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.
name |
type |
expression |
---|---|---|
0 |
int |
|
1 |
int |
|
i |
int |
|
key |
int |
|
0 |
real |
|
1 |
real |
|
r |
real |
|
or: \(|\rvec|\) |
||
ar |
real |
|
ar |
real |
|
aux |
real |
|
ax |
real |
|
ay |
real |
|
az |
real |
|
etot |
real |
|
i |
real |
|
jtot |
real |
|
key |
real |
|
m |
real |
|
phi |
real |
|
r |
real |
|
or: $|$rvec$|$ \ |
||
t |
real |
|
v |
real |
|
or: $|$vvec$|$ |
||
vr |
real |
|
or: or: (rvec$cdot$vvec)/$|$rvec$|$ |
||
vt |
real |
|
or: $sqrt{}$(vvec$^2$-(rvec$cdot$vvec)$^2$/$|$rvec$|^2$)\ |
||
vx |
real |
|
vy |
real |
|
vz |
real |
|
x |
real |
|
y |
real |
|
z |
real |
|
glon |
real |
\(l\) , |
glat |
real |
\(b\) , |
mul |
real |
|
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:
1 2 3 4 5 6 7 | % 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)
\(\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