MEGNO

docs/chaos.md

@@ -0,0 +1,68 @@
+# Chaos indicators
+REBOUND supports different chaos indicators.
+All of these make use of variational equations, but most of the complexity is hidden.
+
+## Initialization
+If you want to use a chaos indicator in REBOUND, first add all the particles to your simulations.
+Then, initialize the variational particles and MEGNO variables with
+=== "C"
+    ```c
+    struct reb_simulation* r = reb_create_simulation();
+    // ... add particles ...
+    reb_tools_megno_init(r);
+    ```
+=== "Python"
+    ```python
+    sim = rebound.Simulation()
+    # ... add particles ...
+    sim.init_megno()
+    ```
+
+REBOUND uses random numbers to initialize the variational particles.
+The initial seed is chosen based on the current time and the process id. 
+This ensures the seed if different every time you run the simulation.
+See the discussion on [random sampling](randomsampling.md) for more details.
+
+If you want to have reproducible result, you can specify the seed manually:
+=== "C"
+    ```c
+    struct reb_simulation* r = reb_create_simulation();
+    // ... add particles ...
+    reb_tools_megno_init_seed(r, 0); // 0 is the initial seed
+    ```
+=== "Python"
+    ```python
+    sim = rebound.Simulation()
+    # ... add particles ...
+    sim.init_megno(seed=0) # 0 is the initial seed
+    ```
+## Accessing chaos indicators
+Once you've initialized the chaos indicators, you can integrate the simulation normally.
+To print out the MEGNO value or the largest Lyapunov characteristic number (LCN), use the following syntax:
+
+=== "C"
+    ```c
+    struct reb_simulation* r = reb_create_simulation();
+    // ... add particles ...
+    reb_tools_megno_init_seed(r);
+    // ... integrate ...
+    printf("MEGNO = %f\n", reb_tools_calculate_megno(r));
+    printf("LCN = %f\n", reb_tools_calculate_lyapunov(r));
+    ```
+=== "Python"
+    ```python
+    sim = rebound.Simulation()
+    # ... add particles ...
+    sim.init_megno()
+    # ... integrate ...
+    print("MEGNO", sim.calculate_megno())
+    print("LCN", sim.calculate_lyapunov())
+    ```
+!!! Important
+    Using chaos indicators is not always straightforward. 
+    It can be particularly tricky to figure out how long to integrate for.
+    If the integration time is too short, you might not capture the Lyapunov timescale accurately. 
+    On the other hand, if the integration time is too long, you can run into problems as well because quantities tend to grow exponentially with time in chaotic systems.
+    In the worst case, the magnitude of the variational particles' coordinates will reach the maximum range of double precision floating point numbers. 
+    Ideally, you would want to rescale the variational particles whenever that happens. 
+    However, this is currently not implemented in these convenience methods and you will have to use variational particles directly.

mkdocs.yml

@@ -30,6 +30,7 @@ nav:
     - units.md
     - orbitalelements.md
     - simulationarchive.md
+    - chaos.md
     - c_randomsamplingfunctions.md
     - c_outputfunctions.md
     - miscellaneous.md

src/rebound.h

@@ -52,6 +52,7 @@ extern volatile sig_atomic_t reb_sigint;  ///< Graceful global interrupt handler
 struct reb_simulation;
 struct reb_display_data;
 struct reb_treecell;
+struct reb_variational_configuration;
 
 struct reb_particle {
     double x;
@@ -72,7 +73,7 @@ struct reb_particle {
     struct reb_simulation* sim; // Pointer to the parent simulation.
 };
 
-// Generic 3d vector, for internal use only.
+// Generic 3d vector
 struct reb_vec3d {
     double x;
     double y;
@@ -278,26 +279,6 @@ enum REB_STATUS {
     REB_EXIT_COLLISION = 7,     // The integration ends early because two particles collided. 
 };
 
-
-/**
- * @brief Struct describing the properties of a set of variational equations.
- * @details One struct describes one or more sets of variational equations.
- * If testparticle is set to -1, then it is assumed that all particles are massive
- * and all particles influence all other particles. If testparticle is >=0 then 
- * the particle with that index is assumed to be a testparticle, i.e. it does not 
- * influence other particles. For second order variational equation, index_1st_order_a/b 
- * is the index in the particle array that corresponds to the 1st order variational 
- * equations.
- */
-struct reb_variational_configuration{
-    struct reb_simulation* sim; // Reference to the simulation.
-    int order;                  // Order of the variational equation. 1 or 2. 
-    int index;                  // Index of the first variational particle in the particles array.
-    int testparticle;           // Is this variational configuration describe a test particle? -1 if not.
-    int index_1st_order_a;      // Used for 2nd order variational particles only: Index of the first order variational particle in the particles array.
-    int index_1st_order_b;      // Used for 2nd order variational particles only: Index of the first order variational particle in the particles array.
-};
-
 // IDs for content of a binary field. Used to read and write binary files.
 enum REB_BINARY_FIELD_TYPE {
     REB_BINARY_FIELD_TYPE_T = 0,
@@ -747,6 +728,7 @@ void reb_tools_init_plummer(struct reb_simulation* r, int _N, double M, double R
 void reb_run_heartbeat(struct reb_simulation* const r);  // used internally
 
 // Functions to add and initialize particles
+struct reb_particle reb_particle_nan(void); // Returns a reb_particle structure with fields/hash/ptrs initialized to nan/0/NULL. 
 void reb_add(struct reb_simulation* const r, struct reb_particle pt);
 void reb_add_fmt(struct reb_simulation* r, const char* fmt, ...);
 struct reb_particle reb_particle_new(struct reb_simulation* r, const char* fmt, ...);    // Same as reb_add_fmt() but returns the particle instead of adding it to the simualtion.
@@ -779,6 +761,45 @@ enum reb_input_binary_messages {
     REB_INPUT_BINARY_WARNING_CORRUPTFILE = 512,
 };
 
+// Chaos indicators
+void reb_tools_megno_init(struct reb_simulation* const r);
+void reb_tools_megno_init_seed(struct reb_simulation* const r, unsigned int seed);
+double reb_tools_calculate_megno(struct reb_simulation* r);
+double reb_tools_calculate_lyapunov(struct reb_simulation* r);
+
+// Variational equations
+
+// Struct describing the properties of a set of variational equations.
+// If testparticle is set to -1, then it is assumed that all particles are massive
+// and all particles influence all other particles. If testparticle is >=0 then 
+// the particle with that index is assumed to be a testparticle, i.e. it does not 
+// influence other particles. For second order variational equation, index_1st_order_a/b 
+// is the index in the particle array that corresponds to the 1st order variational 
+// equations.
+struct reb_variational_configuration{
+    struct reb_simulation* sim; // Reference to the simulation.
+    int order;                  // Order of the variational equation. 1 or 2. 
+    int index;                  // Index of the first variational particle in the particles array.
+    int testparticle;           // Is this variational configuration describe a test particle? -1 if not.
+    int index_1st_order_a;      // Used for 2nd order variational particles only: Index of the first order variational particle in the particles array.
+    int index_1st_order_b;      // Used for 2nd order variational particles only: Index of the first order variational particle in the particles array.
+};
+
+// Add and initialize a set of first order variational particles
+// If testparticle is >= 0, then only one variational particle (the test particle) will be added.
+// If testparticle is -1, one variational particle for each real particle will be added.
+// Returns the index of the first variational particle added
+int reb_add_var_1st_order(struct reb_simulation* const r, int testparticle);
+
+// Add and initialize a set of second order variational particles
+// Note: a set of second order variational particles requires two sets of first order variational equations.
+// If testparticle is >= 0, then only one variational particle (the test particle) will be added.
+// If testparticle is -1, one variational particle for each real particle will be added.
+// index_1st_order_a is the index of the corresponding first variational particles.
+// index_1st_order_b is the index of the corresponding first variational particles.
+// Returns the index of the first variational particle added
+int reb_add_var_2nd_order(struct reb_simulation* const r, int testparticle, int index_1st_order_a, int index_1st_order_b);
+
 
 /**
  * @brief This function calculates the first/second derivative of a Keplerian orbit. 
@@ -939,63 +960,6 @@ void reb_transformations_whds_to_inertial_posvel(struct reb_particle* const part
 
 
 
-/**
- * @brief Add and initialize a set of first order variational particles
- * @param r The rebound simulation to be considered
- * @param testparticle This flag determines if the set of variational particles is for a testparticle or not.
- * If testparticle is >= 0, then only one variational particle (the test particle) will be added.
- * If testparticle is -1, one variational particle for each real particle will be added.
- * @return Returns the index of the first variational particle added
- **/
-int reb_add_var_1st_order(struct reb_simulation* const r, int testparticle);
-
-/**
- * @brief Add and initialize a set of second order variational particles
- * @details Note that a set of second order variational particles requires two sets of first order variational equations.
- * @param r The rebound simulation to be considered
- * @param testparticle This flag determines if the set of variational particles is for a testparticle or not.
- * If testparticle is >= 0, then only one variational particle (the test particle) will be added.
- * If testparticle is -1, one variational particle for each real particle will be added.
- * @param index_1st_order_a The index of the corresponding first variational particles.
- * @param index_1st_order_b The index of the corresponding first variational particles.
- * @return Returns the index of the first variational particle added
- **/
-int reb_add_var_2nd_order(struct reb_simulation* const r, int testparticle, int index_1st_order_a, int index_1st_order_b);
-
-/** 
- * @brief Init the MEGNO particles, enable MEGNO calculation
- * @param r The rebound simulation to be considered
- */
-void reb_tools_megno_init(struct reb_simulation* const r);
-
-/** 
- * @brief Init the MEGNO particles, enable MEGNO calculation, and specify a seed for the random number generation.
- * @param r The rebound simulation to be considered
- * @param seed The seed to use for the random number generator
- */
-void reb_tools_megno_init_seed(struct reb_simulation* const r, unsigned int seed);
-
-/**
- * @brief Get the current MEGNO value
- * @param r The rebound simulation to be considered
- * @return Returns the current value of the MEGNO
- */
-double reb_tools_calculate_megno(struct reb_simulation* r);
-
-/**
- * @brief Returns the largest Lyapunov characteristic number (LCN), or maximal Lyapunov exponent
- * @details MEGNO needs to be enabled to calculate this value.
- * @param r The rebound simulation to be considered
- * @return Returns the current CN
- */
-double reb_tools_calculate_lyapunov(struct reb_simulation* r);
-
-/**
- * @brief Returns a reb_particle structure with fields/hash/ptrs initialized to nan/0/NULL. 
- * @return reb_particle with fields initialized to nan.
- */
-struct reb_particle reb_particle_nan(void);
-
 /**
  * @brief Print out an error message, then exit in a semi-nice way.
  */