Celestial dynamics

dynamics.py

class dynamics.Planet(position, velocity, mass=1.0)[source]

A planet (or some other point-like object).

A planet has a position (a 3-vector), a velocity (3-vector) and a mass (a scalar).

Parameters:

  • position (array) – coordinates \([x, y, z]\)

  • velocity (array) – velocity components \([v_x, v_y, v_z]\)

  • mass (float) – mass \(m\)

accelerate(force, dt)[source]

Set a new velocity for the particle as

\[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{2m}\vec{F} \Delta t\]
Parameters:

  • force (array) – force acting on the planet \([F_x, F_y, F_z]\)

  • dt (float) – time step \(\Delta t\)

move(force, dt)[source]

Set a new position for the particle as

\[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v} \Delta t + \frac{1}{2m}\vec{F} (\Delta t)^2\]
Parameters:

  • force (array) – force acting on the planet \([F_x, F_y, F_z]\)

  • dt (float) – time step \(\Delta t\)

save_position()[source]

Save the current position of the particle in the list ‘trajectory’.

Note: in a real large-scale simulation one would never save trajectories in memory. Instead, these would be written to a file for later analysis.

dynamics.animate(particles, bounds=[[- 1, 3], [- 2, 2]])[source]

Animates the simulation.

Parameters:

  • particles (list) – list of Planet objects

  • bounds (array) – list of lower and upper bounds for the plot as [[xmin, xmax], [ymin, ymax]]

dynamics.calculate_forces(particles, g=0.1)[source]

Calculates the total force applied on each body.

All bodies interact via an inverse square force (Newton’s gravitation). Specifically, the force on body \(j\) due to body \(i\) is

\[\vec{F}_{i \to j} = - \frac{g m_i m_j}{r^2_{ij}} \hat{r}_{i \to j}.\]

The forces are returned as a numpy array where each row contains the force on a body and the columns contain the x, y and z components.

[ [ fx0, fy0, fz0 ],
  [ fx1, fy1, fz1 ],
  [ fx2, fy2, fz2 ],
  ...               ]
Parameters:

  • particles (list) – a list of Planet objects

  • g (float) – gravitational constant

Returns:

forces

Return type:

array

dynamics.calculate_kinetic_energy(particles)[source]

Calculates the total kinetic energy of the system.

\[K_\text{total} = \sum_i \frac{1}{2} m_i v_i^2.\]
Parameters:

particles (list) – a list of Planet objects

Returns:

kinetic energy \(K\)

Return type:

float

dynamics.calculate_momentum(particles)[source]

Calculates the total momentum of the system.

\[\vec{p}_\text{total} = \sum_i \vec{p}_i = \sum_i m_i \vec{v}_i\]
Parameters:

particles (list) – a list of Planet objects

Returns:

momentum components \([p_x, p_y, p_z]\)

Return type:

array

dynamics.calculate_potential_energy(particles, g=0.1)[source]

Calculates the total potential energy of the system.

All bodies interact via an inverse square force (Newton’s gravitation), which has an inverse distance potential energy.

\[U = \sum_{i < j} - \frac{g m_i m_j}{r_{ij}}.\]
Parameters:

  • particles (list) – a list of Planet objects

  • g (float) – gravitational constant

Returns:

potential energy \(U\)

Return type:

float

dynamics.distance_squared(particle1, particle2)[source]

Calculates the squared distance between two atoms,

\[r^2_{ij} = (x_i-x_j)^2 + (y_i-y_j)^2 + (z_i-z_j)^2.\]
Parameters:

  • particle1 (Planet) – the first body

  • particle2 (Planet) – the second body

Returns:

the squared distance \(r^2_{ij}\)

Return type:

float

dynamics.draw(frame, xtraj, ytraj, ztraj, bounds)[source]

Draws a representation of the particle system as a scatter plot.

Used for animation.

Parameters:

  • frame (int) – index of the frame to be drawn

  • xtraj (array) – x-coordinates of all particles at different animation frames

  • ytraj (array) – y-coordinates at all particles at different animation frames

  • ztraj (array) – z-coordinates at all particles at different animation frames

  • bounds (array) – list of lower and upper bounds for the plot as [[xmin, xmax], [ymin, ymax]]

dynamics.euler(particles, dt, time, trajectory_dt=1.0)[source]

Euler algorithm for integrating the equations of motion, i.e., advancing time.

The algorithm works as follows:

  • Forces are calculated at the new positions, \(\vec{F}(t)\).

  • Positions are updated using velocities and forces,
    \[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v}(t) \Delta t + \frac{1}{2m}\vec{F}(t) (\Delta t)^2.\]
  • Velocities are updated using the forces
    \[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{m}\vec{F}(t) \Delta t\]

These operations are done using the methods calculate_forces(), update_velocities() and update_positions().

Unfortunately, the Euler algorithm is not stable and should not be used for proper simulations!

Parameters:

  • particles (list) – a list of Planet objects

  • dt (float) – time step \(\Delta t\)

  • time (float) – the total system time to be simulated

  • trajectory_dt (float) – the positions of particles are saved at these these time intervals - does not affect the dynamics in any way

Returns:

kinetic energy time series, potential energy time series

Return type:

array, array

dynamics.main(filename, algo, dt, simulation_time, recording_dt)[source]

The main program.

Read the initial configuration from a file, run the simulation, and print results.

Parameters:

  • filename (str) – the file with system information

  • algo (str) – either ‘euler’ or ‘verlet’

  • dt (float) – time step \(\Delta t\)

  • time (float) – the total system time to be simulated

  • recording_dt (float) – the positions of particles are saved at these these time intervals - does not affect the dynamics in any way

dynamics.print_progress(step, total)[source]

Prints a progress bar.

Parameters:

  • step (int) – progress counter

  • total (int) – counter at completion

dynamics.read_particles_from_file(filename)[source]

Reads the properties of planets from a file.

Each line should define a single Planet listing its position in cartesian coordinates, velocity components and mass, separated by whitespace:

x0 y0 z0 vx0 vy0 vz0 m0
x1 y1 z1 vx1 vy1 vz1 m1
x2 y2 z2 vx2 vy2 vz2 m2
x3 y3 z3 vx3 vy3 vz3 m3
...
Parameters:

filename (str) – name of the file to read

Returns:

list of Planet objects

Return type:

list

dynamics.show_trajectories(particles, tail=10, skip=10)[source]

Plot a 2D-projection of the trajectory of the system.

The function creates a plot showing the current and past positions of particles.

Parameters:

  • particles (list) – list of Planet objects

  • tail (int) – the number of past positions to include in the plot

  • skip (int) – only every nth past position is plotted - skip is the number n, specifying how far apart the plotted positions are in time

dynamics.update_positions(particles, forces, dt)[source]

Update the positions of all particles using Planet.move() according to

\[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v} \Delta t + \frac{1}{2m}\vec{F} (\Delta t)^2\]
Parameters:

  • particles (list) – a list of Planet objects

  • force (array) – array of forces on all bodies

  • dt (float) – time step \(\Delta t\)

dynamics.update_positions_no_force(particles, dt)[source]

Update the positions of all particles using Planet.move() according to

\[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v} \Delta t\]
Parameters:

  • particles (list) – a list of Planet objects

  • dt (float) – time step \(\Delta t\)

dynamics.update_velocities(particles, forces, dt)[source]

Update the positions of all particles using Planet.accelerate() according to

\[\vec{v}(t+\Delta t) = \vec{v}(t) + \frac{1}{m}\vec{F} \Delta t\]
Parameters:

  • particles (list) – a list of Planet objects

  • force (array) – array of forces on all bodies

  • dt (float) – time step \(\Delta t\)

dynamics.vector(particle1, particle2)[source]

The vector pointing from one planet, \(i\), to another, \(j\),

\[\vec{r}_{i \to j} = \vec{r}_j - \vec{r}_i\]
Parameters:

  • particle1 (Planet) – the first body

  • particle2 (Planet) – the second body

Returns:

components of \(\vec{r}_{i \to j}\), \([x_{i \to j}, y_{i \to j}, z_{i \to j}]\)

Return type:

array

dynamics.velocity_verlet(particles, dt, time, trajectory_dt=1.0)[source]

Verlet algorithm for integrating the equations of motion, i.e., advancing time.

There are a few ways to implement Verlet. The leapfrog version works as follows: First, forces are calculated for all particles and velocities are updated by half a time step, \(\vec{v}(t+\frac{1}{2}\Delta t) = \vec{v}(t) + \frac{1}{2m}\vec{F} \Delta t\). Then, these steps are repeated:

  • Positions are updated by a full time step using velocities but not forces,
    \[\vec{r}(t+\Delta t) = \vec{r}(t) + \vec{v}(t+\frac{1}{2}\Delta t) \Delta t.\]

  • Forces are calculated at the new positions, \(\vec{F}(t + \Delta t)\).

  • Velocities are updated by a full time step using the forces
    \[\vec{v}(t+\frac{3}{2}\Delta t) = \vec{v}(t+\frac{1}{2}\Delta t) + \frac{1}{m}\vec{F}(t+\Delta t) \Delta t\]

These operations are done using the methods calculate_forces(), update_velocities() and update_positions_no_force().

Because velocities were updated by half a time step in the beginning of the simulation, positions and velocities are always offset by half a timestep. You always use the one that has advanced further to update the other and this results in a stable algorithm.

Note

This function is incomplete!

Parameters:

  • particles (list) – a list of Planet objects

  • dt (float) – time step \(\Delta t\)

  • time (float) – the total system time to be simulated

  • trajectory_dt (float) – the positions of particles are saved at these these time intervals - does not affect the dynamics in any way

Returns:

kinetic energy time series, potential energy time series

Return type:

array, array

dynamics.write_particles_to_file(particles, filename)[source]

Write the configuration of the system in a file.

The format is the same as that specified in read_particles_from_file().

Parameters:

  • particles (list) – list of Planet objects

  • filename (str) – name of the file to write

dynamics.write_xyz_file(particles, filename)[source]

Write the configuration of the system in a file.

The information is written in so called xyz format, which many programs can parse.

Parameters:

  • particles (list) – list of Planet objects

  • filename (str) – name of the file to write