About Ncorpi\(\mathcal{O}\,\)N


The development of Ncorpi\(\mathcal{O}\)N was initiated in 2023 by Jérémy Couturier, as a tool to simulate collisional and fragmenting disks made up of many particles. Jérémy is a postdoctoral associate at University of Rochester , under the supervision of professors Miki Nakajima and Alice Quillen. His decision to write a new \(N\)-body software was motivated by their work on the formation of the Moon from a protolunar disk, subsequent to a giant collision between the young Earth and a Mars-sized object.


Previous \(N\)-body integrators did not implement the possibily that, upon a violent collision, particles may fragment instead of just bouncing or merging. Ncorpi\(\mathcal{O}\)N solves this issue with a buil-in fragmentation model that relies on crater scaling and ejecta models to implement realistic fragmentations.

The most difficult challenge in a \(N\)-body code with a large number \(N\) of particles is to time-efficiently treat mutual interactions between them (collisions and gravity). Ncorpi\(\mathcal{O}\)N comes with four different built-in modules for mutual interactions computation, one of them being the fast multipole method Falcon of Walter Dehnen, that Ncorpi\(\mathcal{O}\)N adapts so it can detect collisions as well as compute mutual gravity.


Ncorpi\(\mathcal{O}\)N is entirely written in the C language. Since it was primarily developed to study the formation of the Moon, the particles are referred to as moonlets, even though other types of systems can be efficiently simulated with Ncorpi\(\mathcal{O}\)N, such as protoplanetary disks or planetary rings.

Moonlet structure

The particles in Ncorpi\(\mathcal{O}\)N are stored in an array of struct moonlet. This array of size N_max is allocated once and for all at the beginning of the simulation. The integer N_max is chosen by the user and denotes the maximum number of moonlets that the simulation can handle. An error is returned if the number of particles ever exceeds N_max. The moonlet structure contains the cartesian coordinates of a particle and is defined as

struct moonlet {
      typ x;
      typ y;
      typ z;
      typ vx;
      typ vy;
      typ vz;
      typ mass;
      typ radius;
where the reserved keyword typ can denote double, long double or float, depending on the required numerical precision. This choice is such that sizeof(struct moonlet) is equal to a line of processor cache (generally 64 bytes), when typ = double.

Mutual interactions

Ncorpi\(\mathcal{O}\)N can handle mutual interactions between the particles in four different ways :

  • Brute-forcely : Running time proportional to \(N^2\)
  • With the standard Barnes-Hut tree code : Running time proportional to \(N \ln N\)
  • With a mesh-grid based algorithm : Running time proportional to \(N\)
  • With the fast multipole method Falcon : Running time proportional to \(N\)

The brute-force method refers to the straightforward algorithm that goes over all \(N(N-1)/2\) pairs of particles. It does not require any particular structure and is fastest for \(N < 100\). It however becomes prohibitively slow for large values of \(N\).

Mesh-based algorithm

This algorithm uses a mesh, or 3D-grid, to treat mutual interactions between the moonlets. Collisions are searched, and mutual gravity is computed, only between particles that are neighbours (a neighbours of a particle is a particle in the same cell or in an adjacent cell). The user must choose the expected number of neighbours per particle, and the mesh-size (sidelength of a cell) is updated at every timestep to meet the user's requirement.

The cells of the mesh are stored in an array. With \(1000\) cells per dimension, this method requires an array of size \(1000^3\). In order not to allocate a too large chunk of memory, the array contains pointers towards the cells, instead of the cells themselves, as pointers are much lighter than cells. A cell is just a chain (also called unrolled linked list) of the ids (in the moonlet array defined above) of the moonlets it contains. It is defined as the structure

struct chain {
      int ids[max_ids_per_node]; //The array of ids in that node of the chain
      struct chain * queue;      //A pointer towards the next node of the chain
      int how_many;              //How many ids are in this node

The mesh is a hash table that remembers what ids are in a given cell. It is an array of pointers towards chains of ids, and has the type struct chain **. Most cells of the mesh are empty, and in that case, the hash table contains a NULL pointer.

FalcON method and Barnes-Hut tree code

Both the Barnes-Hut code and falcON fast multipole method uses an octree (renamed boxdot tree). It is defined in C as the structure

struct boxdot {
      struct boxdot * oct[8];    //The eight octants, or children of that cell.
      struct chain * dots;       //A chain containing the ids of the moonlets in that cell.
      typ corner[3];             //Top-left-front corner of the node.
      typ sidelength;            //Sidelength of the node (half that of its parent).
      int id;                    //The unique id of that cell, determined by its Hilbert-Peano order : index in the flattree
      int how_many;              //The number of moonlets in that cell.
      int rotation;              //A rotation used for determination of the Hilbert-Peano digit
      int level;                 //The level of that node in the boxdot tree

This data structure is convenient for the construction of the octree, but its poor locality (children are not next to each other in memory) leads to frequent cache misses, and travelling in the tree requires multiple time-consuming pointer dereferences. Therefore, it is used only for the construction, and falcON and the Barnes-Hut code are performed on a regular array that I name flattree. The flattree is merely an array of nodes (or cells) of the tree. A node is defined by the data structure

struct node {
      typ com[3];            //When treating self-gravity : Center of mass and expansion center. When treating collisions : Average position of the moonlets in the node
      typ C1[3];             //Acceleration, or first order tensor field of interactions
      typ center[3];         //Center of the node
      typ sidelength;        //Sidelength of the node
      typ M0;                //When treating self-gravity : Mass of the cell or zeroth multipole moment. When treating collisions : max(R_i + v_i*timestep)
      typ r_max;             //Convergence radius of the Taylor expansion
      typ r_crit;            //When treating self-gravity : r_max/theta. When treating collisions : r_max + M0
      int * dots;            //Array of ids of moonlets contained in that node
      int idParent;          //Index of the parent node
      int idFirstChild;      //Index of the first child (Child with smallest unique id)
      int how_many_children; //The number of children that node has
      int how_many_dots;     //The number of moonlets that node contains

and contains relevant informations to compute mutual gravity and to search for collisions. The multipole moments \(\mathbf{M}^{(2)}\) to \(\mathbf{M}^{(5)}\) and the tensor fields of interactions \(\mathbf{C}^{(2)}\) to \(\mathbf{C}^{(6)}\) are stored in separate arrays, indexed identically to the flattree. The index of a node in the flattree is given by its Hilbert-Peano order, derived from a 3D version of Hilbert's 1891 space filling curve. This ensures maximum locality. See Ncorpi\(\mathcal{O}\)N's paper for details on the Hilbert-Peano order and explanations on how these two methods work.

When building the octree, a two to three speed-up factor is achieved by adding the moonlets in the Hilbert order of the previous timestep. Initially, only the root cell exists and all of its children are NULL. It is centered on the central mass and has a sidelength root_sidelength. The user defines the quantity subdivision_threshold and I proceed in the following manner to add a particle to a cell struct boxdot * Boxdot:

if (Boxdot -> how_many < subdivision_threshold){
      if (Boxdot -> how_many == 0){
            Allocate memory to the cell;
      Boxdot -> how_many++;
      Add the id of the moonlet recursively to the cell's children;
      if (Boxdot -> how_many == subdivision_threshold){
            Add all the moonlet's ids of the cell to its children;
      Boxdot -> how_many++;
Add the id of the moonlet to the chain Boxdot -> dots;

In order to build the octree, this sequence of instructions is repeated as long as there are moonlets left to add. It has a time-complexity \(\mathcal{O}(\ln N)\), hence an overall complexity \(\mathcal{O}(N \ln N)\) for the whole building procedure.

Set Up


Ncorpi\(\mathcal{O}\)N is an open-source software distributed with the GNU general public license. It can be obtained from the following github repository. Navigate to the directory where you would like Ncorpi\(\mathcal{O}\)N to be installed and run the command

git clone git@github.com:Jeremycouturier/NcorpiON.git

Later on, you can download updates of Ncorpi\(\mathcal{O}\)N by running git pull in the same directory. Alternatively, if you don't have git installed and don't want to install it, you can simply download the repository by hand from github.

Parameter file

The parameter file located at src/parameters.h contains all the data needed to run a particular simulation. The user uninterested in modifying Ncorpi\(\mathcal{O}\)N won't need to touch any other file than this one. Decimal points should not be removed or added to this file. That is, integers must stay that way, as well as floating point numbers.

Physical constants

The physical constants relative to the simulation are defined by these lines

#define Rearth 1.0                   //Radius of the Earth (or central mass) is the unit of length
#define Mearth 1.0                   //Mass of the Earth (or central mass) is the unit of mass
#define G 39.47841760435743          //Gravitational constant is set to 4*pi^2, so that the unit of time is the surface orbital period
#define density 0.1448               //The density of the moonlets (3344 kg/m^3) in Mearth/Rearth^3.
#define Tearth 3.4076                //Earth's sideral period in units of the surface orbital period. Must be larger than 1. Today value is 17.038

By default in Ncorpi\(\mathcal{O}\)N, the units of length, mass and time are the radius of the central body, its mass, and the orbital period at its surface, respectively.

Initial conditions

The parameter file is also used to determine the initial conditions of the simulation.

#define N_max 10000                  //The maximum number of moonlets the simulation can handle
#define N_0 1000                     //The initial number of moonlets. Must be less than or equal to N_max.
#define M_0 0.02214                  //Total (expected) moonlet mass at t=0
#define t_end 512.0                  //Total simulation time (in surface orbital period)
#define time_step 0.0078125          //Timestep of the simulation.
#define output_step 128              //Output occurs every output_step timestep. Matters only if write_to_files_bool is 1
#define radius_stddev 0.57           //Standard deviation of moonlet's radii at t=0 (drawn uniformly), in units of the mean radius. Must be less than 1/sqrt(3) to prevent negative radius
#define eccentricity_min 0.0         //Minimal eccentricity             for a moonlet at t=0
#define eccentricity_max 0.2         //Maximal eccentricity             for a moonlet at t=0
#define sma_min 2.9                  //Minimal semi-major axis          for a moonlet at t=0
#define sma_max 14.0                 //Maximal semi-major axis          for a moonlet at t=0
#define inclination_min 0.0          //Minimal inclination (in radians) for a moonlet at t=0
#define inclination_max 0.174533     //Maximal inclination (in radians) for a moonlet at t=0
#define low_dumping_threshold 2.0    //Moonlets falling below this threshold (in Earth radii) are dumped from the simulation (collision with the Earth or disruption by tidal forces)
#define high_dumping_threshold 200.0 //Moonlets going  beyond this threshold (in Earth radii) are dumped from the simulation (assumed unbounded)
#define max_ids_per_node 173         //The maximum number of ids in each node of the unrolled linked lists (chains). Choose such that sizeof(struct chain) be a multiple of the cache line
#define softening_parameter 0.0      //The softening parameter for mutual gravitational interations, in units of the sum of the radii.
#define seed 778345128               //The seed used for random number generation. Unimportant if seed_bool is 0.
#define switch_to_brute_force 512    //Threshold for N below which the program switches to the brute-force method for mutual interactions. Unimportant if brute_force_bool is 1

The orbital elements are drawn randomly from a uniform distribution whose bounds are defined by these lines. The moonlet's radii are drawn uniformy from the range \(\left[\bar{R}\left(1-\sqrt{3}\sigma\right),\bar{R}\left(1+\sqrt{3}\sigma\right)\right]\), where

  • \(\bar{R}^3=\frac{3M_0}{4\pi\rho\left(1+3\sigma^2\right)}\)
  • \(\sigma\) is radius_stddev
  • \(M_0\) is M_0/N_0
  • \(\rho\) is density

Note that if you require non-random initial conditions, you can achieve that by modifying the function populate in src/structure.c


Boolean values are set to determine how Ncorpi\(\mathcal{O}\)N should behave

#define write_to_files_bool      1   //Determines if the simulation writes to output files. Set to 0 to run speed tests, or if you are satisfied with what is displayed in the terminal   
#define make_animation_bool      1   //Determines if animations of the simulation are produced. write_to_files_bool must be 1
#define seed_bool                0   //Determines if the seed for random number generation is chosen by the user. If seed_bool is 0, then the seed is the number of seconds since 01/01/1970
#define J2_bool                  1   //Determines if the contribution from the symmetrical equatorial bulge is taken into account in the simulation
#define Sun_bool                 0   //Determines if the perturbations from the Sun are taken into account in the simulation
#define collision_bool           1   //Determines if the moonlets are able to collide
#define mutual_bool              1   //Determines if there are mutual gravitational interactions between the moonlets.
/******** Boolean relative to conservation of the total momentum or total angular momentum ********/
#define tam_bool                 0   //Determines if the total angular momentum should be conserved upon merging or fragmenting impact. If tam_bool is 0 then the total momentum is conserved

/******** Booleans relative to collision resolution. Exactly one of them must be 1 ********/
#define elastic_collision_bool   0   //Determines if the collisions are all elastic.
#define inelastic_collision_bool 0   //Determines if the collisions are all inelastic with inelasticity parameter collision_parameter.
#define instant_merger_bool      0   //Determines if the collisions all result in a merger.
#define fragmentation_bool       1   //Determines if the outcome of the collision should be either a merge, a partial fragmentation, a full fragmentation, or a catastrophic disruption, depending on the relative velocity and according to the fragmentation model of Ncorpi\(\mathcal{O}\)N. 

/******** Booleans relative to the treatment of mutual interactions (collisions and self-gravity). Exactly one of them must be 1. ********/
#define brute_force_bool         0   //Determines if a brute-force method should be used for mutual interactions.
#define falcON_bool              1   //Determines if falcON algorithm     should be used for mutual interactions. (Often the best speed/accuracy compromise for large N)
#define standard_tree_bool       0   //Determines if a standard tree code should be used for mutual interactions. (Significantly slower than falcON for the same accuracy).
#define mesh_bool                0   //Determines if the  mesh algorithm  should be used for mutual interactions. (Fastest but gravity with neighbours and three largest moonlets only).

You can choose, among other things,how collisions should be resolved or how mutual interactions should be handled. Here we choose to handle mutual interactions with FalcON (by far best) and to resolve collisions with the built-in fragmentation model of Ncorpi\(\mathcal{O}\)N.


If the mesh \(\mathcal{O}(N)\) algorithm is used to handle mutual interactions, then theses lines of the parameter file must be updated

#define collision_cube_min 80.0      //Minimal sidelength of the mesh-grid centered on the Earth. The mesh-size will never be less than collision_cube_min/collision_cube_cells
#define collision_cube_cells 1024    //The number of mesh cells per dimension of the collision cube. Must be even. For 16+ GiB of RAM (resp. 8 or 4 GiB), choose ~1000 (resp. ~800 or ~500)
#define how_many_neighbours 16.0     //The desired expected number of neighbours for a moonlet

If a tree algorithm is used to handle mutual interactions (that is, if either falcON_bool or standard_tree_bool is 1), then these lines of the parameter file must be updated

#define expansion_order 3            //The order p of the Taylor expansions. This is p = 3 in Dehnen (2002). NcorpiON allows up to p = 6. Minimum is 1 since order 0 yields no acceleration
#define theta_min 0.5                //Minimal value of the tolerance parameter theta. Must be strictly less than 1. Sensible values are 0.2 < theta_min < 0.7
                                     //The precision (and computational time) of the mutual gravity computed by Ncorpion increases with increasing p and decreasing theta_min.
#define subdivision_threshold 17     //A cubic cell is not divided as long as it contains at most that many moonlets. Called s in Dehnen (2002). Must be > 0
                                     //The computational time depends a lot on this threshold. Suggested values are s ~ 26 if p = 3 and s ~ 80 if p = 6, but should be tweaked by the user.
#define root_sidelength 80.0         //Sidelength of the root cell. Must be machine representable. Particles outside of the root cell only feel Earth's gravity and do not affect others
#define level_max 25                 //The maximum allowed number of levels in the tree. Root is at level 0 and a cell at level level_max - 1 is never divided.
#define child_multipole_threshold 1  //If number of particles/number of children is at most this threshold then the multipole moments of a cell are computed directly from the particles. Otherwise, they are computed from that of the children

The following lines are specifically needed when falcON is used

#define N_cc_pre 8                   //For two cells with N1 and N2 moonlets, if N1N2 < N_cc_pre, then the interaction is computed brute-forcely, no matter if they are well-separated or not
#define N_cc_post 64                 //For two cells with N1 and N2 moonlets, if N1N2 < N_cc_post and they are not well-separated, then the interaction is computed brute-forcely
#define N_cs 64                      //If N1 < N_cs, the self-interaction of a cell containing N1 moonlets is computed brute-forcely
#define N_cc_collision 16            //This is N_cc_post. For two non well-separated nodes with Na and Nb moonlets, if NaNb < N_cc_collision, then collisions are searched brute-forcely, otherwise, the largest node is subdivised. N_cc_pre is 0 when falcON is used for collision detection.
#define N_cs_collision 12            //For a node with Na moonlets, if Na < N_cs_collision, then collisions are searched brute-forcely, otherwise, the node is subdivised

Whereas these lines are related to the standard Barnes-Hut tree code

#define N_cb_pre 6                   //If N1 < N_cb_pre, the interaction between a moonlet and a cell with N1 moonlets is performed brute-forcely, no matter if they are well-separated or not
#define N_cb_post 16                 //If N1 < N_cb_post and they are not well-separated, the interaction between a moonlet and a cell with N1 moonlets is performed brute-forcely
#define N_cb_collision 16            //This is N_cb_post. For a moonlet not well-separated from a node with Na moonlets, if Na < N_cb_collision, then collisions are searched brute-forcely, otherwise, the node is subdivised. N_cb_pre is 0 when the standard tree code is used for collision detection.

Some of theses values are not straightforward to understand, but details are given in functions TreeWalk (for falcON) and StandardTree of Sects. 3.3.4 and 3.3.5 of Ncorpi\(\mathcal{O}\)N's paper.


The behaviour of Ncorpi\(\mathcal{O}\)N's fragmentation model is determined by these lines of the parameter file

#define beta_slope 2.6470588235294118//Slope of the power law for fragment size distribution in Leinhardt and Stewart (2012). Must be such that N_tilde is integer.
#define N_tilde 15                   //2*beta_slope/(3 - beta_slope). This is the ratio between the ejected mass and the mass of the second largest fragment -> N° of fragments in the tail
#define mu_parameter 0.55            //The exponent of the impact velocity in the coupling parameter. See Table 3 of Housen & Holsapple (2011)
#define C1_parameter 1.5             //A dimensionless parameter of impact theories. See Table 3 of Housen & Holsapple (2011)
#define k_parameter 0.2              //A dimensionless parameter of impact theories. See Table 3 of Housen & Holsapple (2011)
#define frag_threshold 0.000000005   //Ejected mass threshold below which the collision results in a merger. If the ejected mass is above that threshold but the mass of the second largest fragment is below, then the tail is reunited into a single moonlet. If the mass of the second largest fragment is above that threshold then the collision results in a full fragmentation.                             
#define pq_min_max {-1,3,-1,1}       //Extremal integer values for p_k and q_k to determine the position of the tail fragments with respect to the largest fragment. Must define a rectangle containing exactly N_tilde points with integer coordinates. More precisely, if pq_min_max = {a,b,c,d}, then we must have N_tilde = (b-a+1)(d-c+1). See NcorpiON's paper (Sect. 5.3.5) for details

Running a simulation

If write_to_files_bool is 1, then Ncorpi\(\mathcal{O}\)N writes to output files, at a frequency determined by output_step. The path towards the simulation files is determined by the string path of the parameter file, for example

#define path "./ncorpiON/"

The file ./ncorpiON/e.txt would then contain the eccentricities. Each line corresponds to an output step and has as many columns as the number \(N\) of particles the simulation contained when the output occured.

Compilation + Run + Animation

Once the file src/parameters.h is ready and you are in the src folder, you can compile Ncorpi\(\mathcal{O}\)N, run the simulation and produce an animation all with the unique command


The C code will compile and run. At the end of the simulation, if make_animation_bool is 1, then the same C code will call python3 to make images of each output of the simulation. Finally, ffmpeg will be called to assemble the images into an animated gif and a mp4 video of the simulation.


Alternatively, if you only want to try to compile the code, you can achieve that by running

make clean && make && make clean

The command make needs to be installed for a successfull compilation.


If you decided to compile independently, you can now run Ncorpi\(\mathcal{O}\)N with the command


Making animations

Ncorpi\(\mathcal{O}\)N comes with an add-on implemented in python3 that can makes images out of the simulations files. Then, ffmpeg is used to animate these images. This add-on is only called by the C code if make_animation_bool has been set to 1 in the parameters file.

Image creation

For the add-on to work, you will need to have python3 installed with the libraries matplotlib and numpy. By default, image creation is parallelized on 8 processor cores, but this can be changed by updating the file src/image_creation.sh.

Image assembly

Once the images are produced, the C code will automatically call ffmpeg to assemble the animation. If you need to do that independently, simply run

./ncorpion_animation.sh ./ncorpion/

Make sure that the path given in argument is the same as indicated in the parameters file. The files ./ncorpion/gif/ncorpion.gif and ./ncorpion/gif/ncorpion.mp4 will be created.

Here is the kind of satellite-forming simulation that Ncorpi\(\mathcal{O}\)N is able to carry out in a few minutes


Coming Soon