PERORB(1NEMO) manual page

HTML automatically generated with rman
Table of Contents

Name

perorb - find periodic orbits in symmetric potentials

Description

perorb searches for periodic orbits in two ways:

In the iterative mode it will integrate an orbit for a half or whole period, and iteratively converge towards a periodic orbit assuming regular phase space where the surface of section is a single point. It uses a convergence criterion in the surface of section coordinates. It then can proceed, and attempt to find family members with slightly different energies or spatial extent. The orbits can be stored in a binary file along the way, and also a table with major an minor axis rotation curves, period, energy and average angular momentum can be stored in a table.

In the interactive mode (ncross>1) the selected orbits are integrated for a specified number of surface of section (S.O.S.) crossings, and their S.O.S. coordinates are plotted; from a visual inspection the central periodic orbit(s) can then be estimated at this energy, and used to get good estimates in the iterative mode to find the periodic orbits.

This technique has been used in various forms by various researchers, and is described in its most original form in Henon (1965, Ann.Astr. 28, 499). See also El-Sabaa & Sherief (1990) Ap&SS 167, 305 for some more recent work.

Alternatives approaches are the Particle Swarm Optimization (Skokos et al. 2005). There is also code by T.S.v Albada, but not yet released in NEMO.

Parameters

The following parameters are recognized in any order if the keyword is also given:

in=in-file: An optional input file in standard orbit(5NEMO) format. If given, the input conditions are taken from the first phase point in this orbit (using the keyword last= you can select the last orbit). Default: none. out=out-file output file, into which the orbits are accumulated. The file is in standard orbit(5NEMO) format. If no filename is given, the iterative mode is followed. [Default: none].
freqout=rate: Frequency for output of the timeslots for the periodic orbits. Only used in iterative mode. [default: 1].
maxsteps=number: maximum number of reserved integration steps. [default: 5000].
dt=time-step: time step (constant) to be taken [default: 0.05].
phase=trial-phase: Initial trial phase space coodinates for orbit. The full 6 phase-space coordinates can specified, although you can also specify the (x0,v0) for dirint=0 or (y0,u0) for dirint=1. If a 7th (or 3rd) number is given (optional), this is taken as the trial energy for the first orbit. The trial launch velocity is then calculated from the available potential at the specified launch position.
step=step,perturb: Step in phase space to take to get next orbit (the step is taken in the direction given by the dir keyword. The second number is the perturbation in direction dir used to start the first iteration that should result in finding the periodic orbit. Default: 0.1,0.01
dir=x|y: Direction in which to step for next orbit. This also determines in which surface of section the searching for the periodic orbit is done. Only x and y are allowed for now, since only planar XY-orbits are handled properly.
norbit=n: Number of orbits to investigate. [Default: 1]
accuracy=eps: A relative accuracy parameter to which the final surface of section coordinate is to be accurate. Only used in iterative mode. [Default: 0.001].
maxiter=: Maximum iterations allowed to find a periodic orbit. Default: 50.
ncross=count: Number of crossings through surface of section used to plot/iterate. This count does include the half-period crossings if the keyword period=2 is choosen. If set larger than 1, interactive mode is used. [default: 1]
period=2|1: The number of surface of section coordinates to generate per period. Must be either 1 or 2. A value of 2 make search faster, but a bit more difficult for irregular orbits or near resonances. Also non-symmetric orbits cannot be found this way. [default: 2].
potname=name: name of file of potential(5NEMO) descriptor [default: plummer].
potpars=par-list: List of parameters to the potential descriptor. The first parameter must be the pattern speed in the x-y plane, although rotating frames of reference are not yet supported. The remaining parameters are used by the _inipotential() routine in the potential descriptor. [default: none - let them be defined by routine itself].
potfile=name: name of an optional datafile to the potential descriptor. This might be an N-body snapshot or list of spline fit coefficients etc. [default: none].
tab=table_file: If set to a filename, a table, which normally appears on stdout is now written to a file. It contains the two axes rotation curves, number of steps in the orbit, numer of iterations in surface of section, the full period and the energy: pos,vel,pos,vel,nsteps,niter,period,energy. [Default: none].
mode=int_mode: Specify the integration mode. Any one of euler, leapfrog, rk2 or rk4 can be given. [Default: rk4].
last=t|f: Controls which orbit from the input file is used for the initial conditions. By default, the first one, otherwise the last. Default: f.
headline=: Random user verbiage added to the output file. Default:none

Caveats

The potential (potname= etc) needs to have a proper potential defined for this routine to work.

Examples

By default, perorb will find a circular orbit in a plummer potential near R=1:

perorb
INIPOTENTIAL Plummer: [3d version]
### Warning [perorb]: DRY RUN: No output orbit file created
#  x0       v0       y1        x1    NPT NITER PERIOD   ETOT    LZ_MEAN
  ETOT_ERR
1.007320 0.595679 1.007323 -0.595674 533 3 10.625226 -0.527107 -0.600039 4.08912e-06

The printed columns are pos1(t=0), vel2(t=0), pos2(t=P/2), vel1(t=P/2), l, iter, Period, Energy, meanAngMom. See also orbstat(1NEMO) for another view of the statistics of these orbits, showing the curious property of the 1:1 resonance orbit in a logarithmic potential

% perorb potname=log potpars=0,1,1,0.9 norbit=10 out=log.orb
#  x0      v0        y1        x1    NPT NITER PERIOD   ETOT   LZ_MEAN
   ETOT_ERR
0.956311 0.675829 0.629557 -1.140210 280 5 5.560441 0.984048 -0.679534 -1.4864e-05
1.020431 0.731892 0.706384 -1.174306 290 5 5.761082 1.094431 -0.785260 -1.55914e-06
1.086113 0.781263 0.781720 -1.205717 300 5 5.972464 1.203861 -0.892278 -1.50248e-05
1.153376 0.825038 0.856210 -1.234575 311 5 6.194336 1.312000 -1.000653 -7.58749e-06
1.222193 0.864064 0.930308 -1.260978 323 4 6.426479 1.418565 -1.110458 -1.37877e-05
1.292542 0.898976 1.004360 -1.285151 335 4 6.668597 1.523326 -1.221855 -7.93163e-06
1.364403 0.930272 1.078516 -1.307186 348 4 6.920329 1.626100 -1.334623 -1.17395e-05
1.437707 0.958411 1.153092 -1.327410 361 4 7.181385 1.726745 -1.448891 -8.47802e-06
1.512414 0.983745 1.228047 -1.345782 374 4 7.451368 1.825161 -1.564586 -9.25693e-06
1.588487 1.006573 1.303499 -1.362555 388 4 7.729841 1.921275 -1.681411 -8.9171e-06
% orbstat log.orb
# T    E    x_max    y_max    u_max    v_max    j_mean    j_sigma
0 0.984048 0.629601 0.956311 0.675829 1.14031 -0.679593 0.0252821
0 1.09443 0.706429 1.02043 0.731892 1.1744 -0.785327 0.0292126
0 1.20386 0.781738 1.08611 0.781263 1.20575 -0.892351 0.0331657
0 1.312 0.85621 1.15338 0.825038 1.23457 -1.00073 0.037188
0 1.41857 0.930308 1.22219 0.864064 1.26098 -1.11054 0.0412889
0 1.52333 1.00436 1.29254 0.898976 1.28515 -1.22194 0.0454228
0 1.6261 1.07857 1.3644 0.930272 1.30726 -1.33472 0.0496375
0 1.72675 1.15309 1.43771 0.958411 1.32741 -1.44899 0.0538812
0 1.82516 1.22807 1.51241 0.983745 1.34581 -1.56469 0.0581483
0 1.92128 1.30352 1.58849 1.00657 1.36259 -1.68152 0.0624934
% orbstat log.orb | tabmath - - ’(%4*%5)/(%3*%6)’ all
0.900218 
0.900214 
0.900228 
0.900224 
0.900222 
0.90022 
0.900205 
0.900233 
0.900213 
0.900214

The last command shows the ratio of the angular momentum along the major and minor axis, and that they equal the axis ratio (q) of the potential.

In interactive mode the Time,Ycross,VYcross,Energy for each surface of section are listed, and can be piped into tabplot to get an idea of phase space and help you getter better initial conditions.

1% perorb potname=log potpars=0.1,1,1,0.9 phase=0.1,-0.2 ncross=50
#  Time    Ycross    VYcross    Energy
2.36021 0.0947884 -0.0301556 0.0322215
4.70673 0.0817796 -0.0484523 0.0322859
2% perorb potname=log potpars=0.1,1,1,0.9 phase=0.1,-0.2 ncross=50 | tabplot
- 2 3
3% perorb potname=log potpars=0.1,1,1,0.9 phase=0.1,-0.2 ncross=50 | tablsqfit
- 2 3 fit=ellipse

 x-center:        0.0732815
 y-center:        1.0291e-06

4% perorb potname=log potpars=0.1,1,1,0.9 phase=0.1,-0.2,0.0322215
# y0      u0        x1        v1     NPT   NITER PERIOD   ETOT  LZ_MEAN
  ETOT_ERR
0.074508 -0.225473 -0.166434 -0.100333 234 6 4.650688 0.032222 0.029217 -3.71096e-06

In the 1st example we establish that energy is roughly conserved (and the value is 0.0322215). The 2nd example plots the SOS. The 3rd example takes column 2 and 3 and fits an ellipse and reports the center. The 4th example tries to reproduce the periodic orbit with that specified energy.

Author

Peter Teuben

Files

pjt/orbit    original sources

Update History

1980s           V0.x Teuben’s thesis work (written in SHELTRAN)    PJT
22-may-90    V1.0 created from old Cyber CDC-7600 program ’PERORB’  PJT
24-may-91    V1.1 rotating XY-frames, fixed energy option  PJT
19-apr-95    V1.5 various, rk4 is now default integrator    PJT
1-mar-03    minor code cleanup, less lies in the man page, added etot_err    PJT
19-aug-04    V1.6 added last= and fixed an allocation problem    PJT
29-oct-2019    V1.7 allow phase= with 2 or 3 values    PJT

Table of Contents

Name
Synopsis
Description
Parameters
Caveats
Examples
See Also
Author
Files
Update History