API Reference

Main module

These are the 3 main classes of the module. PoincareMapper is used for the computations, Tomography for visualization and PoincareCollection for exporting and importing.

class poincarepy.PoincareMapper(pot, crossing_function=<function event_yplanecross>, max_integ_time=200, dx=1e-08, dvx=1e-08)

Bases: object

Class that handles computations (orbit integration, map generation,..) based on a physical potential

This is the base class that, provided with a Potential object, performs dynamics computation in this potential. These include: orbit integration, mapping, generation of Poincaré maps as well as collections of those at multiple energies, and periodic orbit searching.

In the following, y is understood as a 4-vector containing the 4 dof of the system. By default, the map is drawn in the y[1] = 0 section of the phase space, and the y[3] variable is eliminated by the Hamiltonian integral. Thus, the map space is M = {y[0], y[2]}. In a cartesian frame, we would have y = [x,y,vx,vy] and M = {x,vx} with y=0 and vy calculated via H. In a meridonal plane frame, we would have y = [r,z,vr,vz] and M = {r,vr}, with z=0 and vz calculated via H.

Other coordinate frames are not yet implemented and may lead to inconsistencies (if for example the kinetic energy is not in a cartesian form K = 0.5(y[2]^2 + y[3]^2)).

Parameters:

  • pot (Potential) – The physical potential that describes the dynamical system

  • crossing_function (callable) – Condition function used to compute Poincaré maps, has the signature f(t,y) = 0 when a point is to be drawn. One sets the crossing direction by defining an attribute on the function object as f.direction = +1, -1, where +/- indicate the direction of crossing.

  • max_integ_time (float) – The maximum allowed integration time (in physical/internal units). If a computation exceeds this limit, an error is raised. In general, the effective integration time in most computations will be determined by a condition (e.g. number of crossing events) rather than a total time, so this value only serves as an upper limit.

  • dx (float) – Finite difference step used for the first dynamical variable (x,r,.) in the jacobian computation

  • dvx (float) – Finite difference step used for the second dynamical variable (vx,vr,.) in the jacobian computation

map(q, E, N=1)

Map a point q to its Poincare map after N crossings (in 2D)

Parameters:

  • q (array-like size (2,)) – Starting point in phase space

  • E (float) – Energy of the orbit

  • N (int) – Number of crossings (order of the map)

Returns:

q’ – Result of the mapping. Returns None if the starting point was outside of the zero-velocity curve of the potential at energy E

Return type:

array (2,)

jac(q, E, N=1)

2D-Jacobian matrix of the map() function

Parameters:

  • q (array-like size (2,)) – Starting point in phase space

  • E (float) – Energy of the orbit

  • N (int) – Number of crossings (order of the map)

Returns:

J – Jacobian matrix of the mapping. Returns None if the starting point was outside of the zero-velocity curve of the potential at energy E

Return type:

array (2,2)

find_periodic_orbit(q0, E, N=1, print_result=False, print_progress=False, maxiter=100, eps=1e-05)

Starting from q0 at energy E, find an N-periodic orbit

Parameters:

  • q0 (array-like size (2,)) – Initial guess of the periodic orbit search

  • E (float) – Energy of the orbit

  • N (int) – Number of crossings (order of the map)

  • print_result (bool) – Print the found periodic orbit

  • print_progress (bool) – Print the steps performed during the search

  • maxiter (int) – Maximum number of iterations for the orbit search

  • eps (float) – Precision |q_n-q_n-1| < eps required to find an orbit

Returns:

q* – Phase-space location of the found periodic orbit. Returns None if either the starting point was outside of the zero-velocity curve of the potential at energy E or if the search reached its maximum number of allowed iterations.

Return type:

array (2,) or None

integrate_orbit(q, E, N=1, N_points_orbit=1000)

Integrate an orbit with given starting point

Same principle as map() but the whole history of crossings, as well as the corresponding configuration-space orbit is returned, with a resolution provided as argument.

Parameters:

  • q (array-like size (2,)) – Starting point in phase space

  • E (float) – Energy of the orbit

  • N (int) – Number of crossings of the section plane

  • N_points_orbit (int or None) – Number of points in the returned configuration space orbit. Increase for smoother orbits, decrease to reduce memory usage. Does not affect the computed sections in any way. Set to None to use integrator timesteps, which is faster, but cannot be used to generate sections since the number of points is then variable.

Returns:

  • sec (array (2,N)) – Phase space section of the integrated orbit, i.e. the N points at which the orbit has crossed the section plane in a given direction

  • orb (array(2,)*) – Configuration space orbit, length depending on integration time required to reach N crossings

integrate_orbit_full(y0, tf, N_points_orbit=1000)

Integrate an orbit with given starting point in the full phase space

Given a point y0 in the 4D phase space, integrate until tf. This method does not calculate crossing points and outputs only a phase space orbit.

Parameters:

  • y0 (1D-array of size 4) – The starting point (x,y,vx,vy) in phase space

  • tf (float) – Integration time

  • N_points_orbit (int) – Number of samples in (0,tf) at which to store the solution. The sampling times are uniformly distributed in (0,tf)

Returns:

yt – The trajectory (x(t),y(t),vx(t),vy(t)) in phase space

Return type:

2D-array of size (4,N_points_orbit)

section(E, xlim, N_orbits, N_points, xdot=0.0, auto_lim=True, N_points_orbit=10000, Nsteps_lim=20, print_progress=False)

Calculate a surface of section at given energy

A surface of section is a collection of N_orbits orbits that fill a given region of the phase space. Each orbit is integrated until N_points crossings of the section plane occur, and at each (oriented) crossing, a point is added to the surface. Phase space is filled by considering ICs distributed linearly in a given or calculated x-range, with xdot_0 a constant value.

Parameters:

  • E (float) – Energy of the section

  • xlim (array-like or tuple) – If auto_lim is True, lower and upper bounds used for the automatic xlim search. If auto_lim is False, these are used as-is to compute the section.

  • N_orbits (int) – Number of orbits in the surface of section

  • N_points (int) – Number of points per orbit in the map (e.g number of crossings of the section plane)

  • xdot (float) – Constant scalar value for the xdot (y) initial condition

  • auto_lim (bool) – Whether or not to automatically determine xlims of the section

  • N_points_orbit (int or None) – Number of points per orbit in the configuration space. Increase for smoother orbits, decrease to reduce memory usage. If set to None, no orbit output will be produced and internally the integrator won’t need to do any interpolation, which makes it run moderately faster.

  • Nsteps_lim (int) – Number of subdivisions of the interval to use when auto-searching. Use higher value for larger intervals.

Returns:

  • secs (array (N_orbits,2,N_points)) – The surface of section. First dim indexes the orbits, second the variable (x,xdot), third are the points

  • orbs (list (N_orbits,) of arrays (2,)*) – Configuration-space orbits corresponding to the sections. List because the number of points in the orbit is not constant (deps on the integration time that was required to reach N_points crossings of the plane). If N_points_orbit = None, orbs = None

  • zvc (array (2,800)) – Zero-velocity curve (limit of the section). Consists of (x,xdot) pairs that contour the surface of section once.

section_collection(E, xlim, N_orbits, N_points, xdot=0.0, N_points_orbit=10000, auto_lim=True, Nsteps_lim=200)

Calculate a set of sections at different energies

Parameters:

  • E (array) – Energies of the sections

  • N_orbits (int) – Number of orbits in each surface of section

  • N_points (int) – Number of points per orbit in the maps (e.g number of crossings of the section plane)

  • xlim (array-like or tuple) – If auto_lim is True, lower and upper bounds used for the automatic xlim search. If auto_lim is False, these are used as-is to compute the section.

  • xdot (array) – Constant scalar values for the xdot (y) initial condition, used if xlim is explicitely provided

  • x0 (array of tuples) – Starting point (initial guess) for the automatic xlim computation

  • N_points_orbit (int) – Number of points for the orbit output. If None, the values at each time step of the integrator will be used (faster)

  • auto_lim (bool) – Whether or not to automatically determine xlims of the section

Returns:

  • sections (array (N_energies,N_orbits,2,N_points)) – The surface of section. First dim indexes the orbits, second the variable (x,xdot), third are the points

  • orbits (array (N_energies,N_orbits,2,N_points_orbit)) – Configuration-space orbits corresponding to the sections. List because the number of points in the orbit is not constant (deps on the integration time that was required to reach N_points crossings of the plane)

  • zvcs (array (N_energies,2,800)) – Zero-velocity curves (limits of the sections). Each entry consists of (x,xdot) pairs that contour that surface of section once.

xlim(E, xdot, a, b, Nsteps=200)

Helper function to calculate physical x-limits at given energy

Use the condition E=phi with y=0, ydot=0 to compute the maximum allowed values for x at given xdot.

Parameters:

  • E (float) – Energy of the map

  • xdot (float) – Conjugate variable to x

  • a (float) – Lower bound for the search

  • b (float) – Upper bound for the search

  • Nsteps (int) – Number of subdivisions of [a,b] to use to find the roots. Default is 200 (Use larger value if a large interval is provided)

zvc(E, x)

Helper function to calculate zero velocity curve

class poincarepy.PoincareCollection(E_list, orbits_list, sections_list: ndarray, zvc_list: ndarray, mapper: PoincareMapper)

Bases: object

Container class for a collection of surfaces of section and related data, used for pickling

This class is essentially a dictionnary of the elements making up a Poincaré Map collection, that are required by the Tomography class. It is used only for convenience when saving (pickling) calculated data for reuse.

Parameters:

  • E_list (list or array) – Array of energy values corresponding to the energies of each set of orbits

  • orbits_list (list) – List of shape (N_energies,N_orbits) where each element is again a list corresponding to a set of orbits at a given energy. An orbit is an array of shape (2,N) with the first axis giving x/y and N the number of points in the orbit.

  • sections_list (list) – Same idea as for the orbits, except that a surface of section is an array of shape (2,nb_points) with nb_points the number of points in the Poincare map per orbit (number of crossings of the x/y plane)

  • potential (Potential) – Potential that generated the collection (used if the collection is imported and exported)

energylist

E_list

Type:

list

orbitslist

orbits_list

Type:

list

sectionslist

sections_list

Type:

list

nb_energies

Number of energy levels in the collection (= len(energylist))

Type:

int

nb_orbits_per_E

Number of orbits per energy

Type:

int

(The total number of orbits in a collection is nb_energies x nb_orbits_per_E)
save(fname='pcollection.pkl')

Save PoincareCollection to a file

Parameters:

fname (str) – Name of the file to write

classmethod load(fname)

Create PoincareCollection from saved file

Usage: collection = PoincareCollection.load([filename])

Parameters:

fname (str) – Name of the file to read

class poincarepy.Tomography(collection: PoincareCollection, figsize=(12.5, 6), title=None, redraw_orbit: bool = True, axlabels=['$x$', '$\\dot{x}$', '$x$', '$y$'])

Bases: object

Tomographic visualisation of an ensemble of surfaces of section at different energies

Visualization class. Using precalculated data (Poincaré sections at different energies), allows the user to pan through the energy dimension in a tomography-like fashion and see how the phase space of the system (represented with a Poincaré section) changes in this dimension. One can click on a given orbit in the reduced phase space to see how it looks like in configuration space.

In addition, if the user wishes to analyze some part of the phase space in more detail, he or she can redraw orbits in real time, using one of 3 functions: - Single orbit redrawing: This mode is toggled using the “z” key and, by clicking on a point in the reduced phase space, a new orbit is calculated using the initial conditions derived from the clicked point. - Rectangular selection redrawing: This mode is toggled using the “t” key and allows the user to select a rectangular region in which a desired N number of orbits, with ICs lying uniformly on the middle horizontal of the rectangle, are redrawn. - Full view redraw: Is toggled with a button and simply redraws the section at the current energy level with N orbits, where N can be set in a text box.

The class methods are essentially all internal and are thus not described in detail in the API.

Parameters:

  • collection (PoincareCollection) – The collection (i.e. set of surfaces of section at different energies) to visualize.

  • figsize (tuple of floats, optional) – Size of the figure in inches

  • redraw_orbit (bool, optional) – Whether to redraw a selected orbit when switching to a different energy. Default is true

show(idx)

Function that updates the plot when energy is changed :param idx: Index of the energy level to switch to :type idx: int

Included Potentials

These are the potentials that have been implemented. All potentials inherit from the Potential class.

poincarepy.potentials.G_grav = 1

Gravitational constant, set to one by default. Changing this value allows to use different unit systems.

class poincarepy.potentials.Potential

Bases: object

Base class of a physical potential

Each potential must redefine at least 3 methods: phi: The potential at given phase space position accel: The acceleration at given phase space position info: String containing information about the specific potential

accel is used to compute the RHS of the equation dy/dt = RHS(t,y), which is solved by the integrator

phi(y)
accel(y)
RHS(t, y)
info()
plot_x(x0, x1, y=0, Npoints=100, ax=None, label=None)
plotcontour(x0, x1, y0, y1, Npoints=100, levels=20, cmap='viridis', ax=None)
class poincarepy.potentials.LogarithmicPotential(v0=10.0, rc=1.0, q=0.8)

Bases: Potential

phi(y)
accel(y)
info()
class poincarepy.potentials.HomospherePotential(a=1.0, M=1)

Bases: Potential

Potential of a homogeneous sphere

phi(y)
accel(y)
info()
class poincarepy.potentials.PlummerPotential(a=5.0, M=3000.0)

Bases: Potential

phi(y)
accel(y)
info()
class poincarepy.potentials.zRotation(omega)

Bases: Potential

Coriolis: -2*Omega x V Centrifugal: -Omega x (Omega x X)

accel(y)
phi(y)
info()
class poincarepy.potentials.CombinedPotential(*pots: Potential)

Bases: Potential

manage_pots(given_potentials)
phi(y)
accel(y)
info()
class poincarepy.potentials.EffectiveLogarithmic_cylindrical(v0=10.0, rc=1.0, q=0.8, Lz=0.0)

Bases: Potential

phi(y)
accel(y)
info()
class poincarepy.potentials.PointMassPotential(M=1000.0)

Bases: Potential

phi(y)
accel(y)
maxval_x(E)