Orbits

Here are some examples of how to compute orbits within NEMO. First we need to select a potential. So, let's take the flat rotation curve potential with a softened core, which fits most disk galaxies, although this represents most of the spherical halo potential. Potentials are implemented via the keywords potname=, potpars= and potfile=, and all carry the same meaning accross programs. Examples of potentials in source code form are in $NEMO/src/orbit/potential/data. For example, let us look at halo.c, which has just two parameters, the peak rotation curve and the core radius. Let us first just plot the rotation curve:

    rotcurves name1=halo pars1=0,1,1 radii=0:8:0.1

(the first parameter in pars1= is the pattern speed, which we will not use here)

Now we need to compute an orbit in this potential. You can give initial conditions (the advantage of this logarithmic potential is that there is no escape energy, so you don't need to worry about making sure the energy is negative to keep it within bounds), integrate them and vizualize the orbit as follows:

 mkorbit - 1 0 0 0 1 0 potname=halo | orbint - - 10000 0.01 | orbplot -

As you can see, this orbit is not closed, actually, it looks closed after three revolutions. For fun, you can perturb it a little bit, for example in the initial X position, and integrate it 10 times longer:

 mkorbit - 1.01 0 0 0 1 0 potname=halo | orbint - - 100000 0.01 | orbplot -

Let us plot how the radius varies with time. Since orbplot cannot handle arbitrary variables yet, we need to convert to a snapshot (where we loose time) and use the index 'i' as time. You can also choose the angular momentum (jtot), but it will be constant and 1.0 in this example.

 mkorbit - 1 0 0 0 1 0 potname=halo | orbint - - 10000 0.01 | otos - - | snapprint - i,r | tabplot -

To find a periodic orbit, you can either solve for circular orbit analytically if the potential allows this (for this potential the relationship is

 v/v0 = r/sqrt(r**2+rc**2)

or you compute them using a surface of section (orbsos) or periodic orbit finder using variational principles (perorb). For example:

 
 mkorbit - 1.01 0 0 0 1 0 potname=halo | orbint - - 100000 0.01 | orbsos - | tabplot - 3 4

or try to check if the orbit integrator can produce a circular orbit at radius 1:

 mkorbit - 1 0 0 0 '1/sqrt(2)' 0 potname=halo | orbint - - 10000 0.01 | orbplot -

Now we need to rewrite halo.c and add a small spiral perturbation.