Updating version to 4.2.0

README.md

@@ -1,4 +1,4 @@
-[![Version](https://img.shields.io/badge/rebound-v4.1.1-green.svg?style=flat)](https://rebound.readthedocs.org)
+[![Version](https://img.shields.io/badge/rebound-v4.2.0-green.svg?style=flat)](https://rebound.readthedocs.org)
 [![PyPI](https://badge.fury.io/py/rebound.svg)](https://badge.fury.io/py/rebound)
 [![GPL](https://img.shields.io/badge/license-GPL-green.svg?style=flat)](https://github.com/hannorein/rebound/blob/main/LICENSE)
 [![Paper](https://img.shields.io/badge/arXiv-1110.4876-green.svg?style=flat)](https://arxiv.org/abs/1110.4876)

changelog.md

@@ -4,6 +4,11 @@ This changelog only includes the most important changes in recent updates. For a
 
 ## Version 4.x
 
+### Version 4.2.0
+* It is now possible to programmatically change all aspects of a REBOUND visualization. This can be used to set up default viewing options or to render animations. See the C examples in `animation_solar_system` and `animation_saturn_rings`.
+* Reworked matrix operations in visualization routines to follow the Model-View-Projection paradigm.
+* Fixed an issues where unit tests would fail because a binary file was not deleted.
+
 ### Version 4.1.1
 * Fixed python wheels for windows.
 

examples/animation_saturns_rings/Makefile

@@ -0,0 +1,35 @@
+export OPENGL=0# Set this to 1 to enable OpenGL
+export SERVER=1# Set this to 1 to enable the visualization web server
+include ../../src/Makefile.defs
+
+# CCPROBLEM is defined in Makefile.defs to allow for
+# a compact cross platform Makefile
+.PHONY: all librebound
+all: problem.c librebound
+	@echo "Compiling $< ..."
+	$(CCPROBLEM)
+	@echo ""
+	@echo "Compilation successful. To run REBOUND, execute the file '$(EXEREBOUND)'."
+	@echo ""
+
+librebound:
+	@echo "Compiling shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/
+	@-$(RM) $(LIBREBOUND)
+	@$(LINKORCOPYLIBREBOUND)
+	@echo ""
+
+clean:
+	@echo "Cleaning up shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/ clean
+	@echo "Cleaning up local directory ..."
+	@-$(RM) $(LIBREBOUND)
+	@-$(RM) $(EXEREBOUND)
+
+rebound_webgl.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL enabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -DOPENGL=1 -sSTACK_SIZE=655360 -s USE_GLFW=3 -s FULL_ES3=1 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_webgl.html -o rebound_webgl.html
+
+rebound_console.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL disabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -sSTACK_SIZE=655360 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_console.html -o rebound_console.html

examples/animation_saturns_rings/problem.c

@@ -0,0 +1,97 @@
+/**
+ * Animation of the Saturn's Rings
+ * 
+ * This examples show how to use display_settings to
+ * programmatically change the visualization of a 
+ * REBOUND simulation. Here, we visualize a simulation of
+ * Saturn's rings and rotate the viewwing angle programatically.
+ * To understand what a 4x4 view matrix is, you can read 
+ * up on linear algebra for computer graphics, specifically
+ * the Model-View-Projection (MVP) paradigm. 
+ *
+ */
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "rebound.h"
+
+double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v);
+
+// The heartbeat function is called once a timestep
+void heartbeat(struct reb_simulation* const r){
+    // Construct a rotation. 
+    struct reb_vec3d axis = {.y=1., .z=0.2}; 
+    struct reb_rotation rot = reb_rotation_init_angle_axis(0.003, axis); // small increment every timestep
+
+    // Convert quaternion to rotation matrix
+    struct reb_mat4df rm = reb_rotation_to_mat4df(rot);
+
+    // Apply incremental rotation to view matrix.
+    r->display_settings->view = reb_mat4df_multiply(rm, r->display_settings->view); 
+}
+
+
+int main(int argc, char* argv[]) {
+    struct reb_simulation* r = reb_simulation_create();
+
+    // The heartbeat function handles the visualization in this example.
+    r->heartbeat         = heartbeat;           
+
+    // Setup problem. For more details, see the shearing sheet example.
+    r->opening_angle2    = .5;
+    r->integrator        = REB_INTEGRATOR_SEI;
+    r->boundary          = REB_BOUNDARY_SHEAR;
+    r->gravity           = REB_GRAVITY_TREE;
+    r->collision         = REB_COLLISION_TREE;
+    r->collision_resolve = reb_collision_resolve_hardsphere;
+    r->coefficient_of_restitution = coefficient_of_restitution_bridges;
+    r->ri_sei.OMEGA      = 0.00013143527;       // 1/s
+    r->minimum_collision_velocity = r->ri_sei.OMEGA*0.001;
+    r->G                 = 6.67428e-11;         // N / (1e-5 kg)^2 m^2
+    r->softening         = 0.1;                 // m
+    r->dt                = 1e-3*2.*M_PI/r->ri_sei.OMEGA;  // s
+
+    reb_simulation_configure_box(r, 100, 2, 2, 1); // 100m box
+    r->N_ghost_x = 2; r->N_ghost_y = 2; r->N_ghost_z = 0;
+
+    // Add all ring paricles
+    double mass = 0;
+    while(mass < 400*r->boxsize.x*r->boxsize.y){ // 400kg/m^2 surface density
+        struct reb_particle pt = {0};
+        pt.x         = reb_random_uniform(r, -r->boxsize.x/2.,r->boxsize.x/2.);
+        pt.y         = reb_random_uniform(r, -r->boxsize.y/2.,r->boxsize.y/2.);
+        pt.z         = reb_random_normal(r, 1.);                    // m
+        pt.vy         = -1.5*pt.x*r->ri_sei.OMEGA;
+        double radius     = reb_random_powerlaw(r, 1., 4.,-3);
+        pt.r         = radius;                        // m
+        double        particle_mass = 400.0*4./3.*M_PI*radius*radius*radius; // 400kg/m^3 particle density
+        pt.m         = particle_mass;     // kg
+        reb_simulation_add(r, pt);
+        mass += particle_mass;
+    }
+
+    // Normally the visualization settings are determined by the 
+    // user interface. If we add the display_settings struct to 
+    // the simulation itself, it will overwrite any change the 
+    // user has made and allows us to programatically change any 
+    // settings such as the orientation, zoom, etc. 
+    reb_simulation_add_display_settings(r);
+
+    // This allows you to connect to the simulation using
+    // a web browser. Simply go to http://localhost:1234
+    reb_simulation_start_server(r, 1234);
+
+    // Integrate forever
+    reb_simulation_integrate(r, INFINITY);
+    reb_simulation_free(r);
+}
+
+// This example is using a custom velocity dependend coefficient of restitution
+double coefficient_of_restitution_bridges(const struct reb_simulation* const r, double v){
+    // assumes v in units of [m/s]
+    double eps = 0.32*pow(fabs(v)*100.,-0.234);
+    if (eps>1) eps=1;
+    if (eps<0) eps=0;
+    return eps;
+}
+

examples/animation_solar_system/Makefile

@@ -0,0 +1,35 @@
+export OPENGL=0# Set this to 1 to enable OpenGL
+export SERVER=1# Set this to 1 to enable the visualization web server
+include ../../src/Makefile.defs
+
+# CCPROBLEM is defined in Makefile.defs to allow for
+# a compact cross platform Makefile
+.PHONY: all librebound
+all: problem.c librebound
+	@echo "Compiling $< ..."
+	$(CCPROBLEM)
+	@echo ""
+	@echo "Compilation successful. To run REBOUND, execute the file '$(EXEREBOUND)'."
+	@echo ""
+
+librebound:
+	@echo "Compiling shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/
+	@-$(RM) $(LIBREBOUND)
+	@$(LINKORCOPYLIBREBOUND)
+	@echo ""
+
+clean:
+	@echo "Cleaning up shared library $(LIBREBOUND) ..."
+	$(MAKE) -C ../../src/ clean
+	@echo "Cleaning up local directory ..."
+	@-$(RM) $(LIBREBOUND)
+	@-$(RM) $(EXEREBOUND)
+
+rebound_webgl.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL enabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -DOPENGL=1 -sSTACK_SIZE=655360 -s USE_GLFW=3 -s FULL_ES3=1 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_webgl.html -o rebound_webgl.html
+
+rebound_console.html: problem.c
+	@echo "Compiling problem.c with emscripten (WebGL disabled)..."
+	emcc -O3 -I../../src/ ../../src/*.c problem.c -DSERVERHIDEWARNING -sSTACK_SIZE=655360 -sASYNCIFY -sALLOW_MEMORY_GROWTH -sSINGLE_FILE -sEXPORTED_RUNTIME_METHODS="callMain" --shell-file ../../web_client/shell_rebound_console.html -o rebound_console.html

examples/animation_solar_system/problem.c

@@ -0,0 +1,210 @@
+/**
+ * Animation of the Solar System
+ * 
+ * This examples show how to use display_settings to
+ * programmatically change the visualization of a 
+ * REBOUND simulation. This can be used to render movies.
+ * To understand what a 4x4 view matrix is, you can read 
+ * up on linear algebra for computer graphics, specifically
+ * the Model-View-Projection (MVP) paradigm. 
+ *
+ */
+#include "rebound.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+struct reb_mat4df view0; // Initial view matrix
+double dt0;              // Initial timestep
+
+double inflate_size = 5; // Inflate particle sizes
+
+// Quadratic ease in/out function for smooth animations
+double ease_in_out(double x){
+    if (x<0.0) x=0.0;
+    if (x>1.0) x=1.0;
+    return x < 0.5 ? 2.0*x*x : 1.0-(2.0-2.0* x)*(1.0-1.0*x);
+}
+
+// Heartbeat function changes the display settings every timestep
+// The following animations play in this order:
+// - Zoom in
+// - Rotation around the x axis
+// - Slowing down of simulation
+// - Combined zoom in on Earth and rotation 
+void heartbeat(struct reb_simulation* const r){
+    // Construct a 90 rotation around the x axis
+    struct reb_vec3d a = {.x=1.};                               // x axis
+    struct reb_rotation rot_x = reb_rotation_init_angle_axis(M_PI_2, a); // pi/2 = 90 degrees
+
+    // We start from the original view matrix view0, but we could also 
+    // apply consequitive changes to the current view matrix stored in 
+    // r->display_settings->view
+    struct reb_mat4df view = view0;
+
+    if (r->t < 2.*2.*M_PI){                                         // First 2 years
+
+        // Slerp (interpolate) between no rotation (identity) and a rotation around the x axis
+        float t = ease_in_out(r->t/(2.*2.*M_PI));                   // runs from 0 to 1
+        struct reb_rotation rot_slerp = reb_rotation_slerp(reb_rotation_identity(), rot_x, t);
+
+        // Convert rotation quaternion to 4d rotation matrix and then
+        // operate the rotation matrix on the view matrix.
+        struct reb_mat4df rm = reb_rotation_to_mat4df(rot_slerp);
+        view = reb_mat4df_multiply(rm, view); 
+
+    }else if (r->t > 2.*2.*M_PI && r->t < 4.*2.*M_PI){              // Year 2 to 4
+
+        // Increase length of trail to 64
+        r->display_settings->past_N = 64;
+
+        // Apply same rotation matrix as before
+        struct reb_mat4df rm = reb_rotation_to_mat4df(rot_x);
+        view = reb_mat4df_multiply(rm, view); 
+
+        // Also apply a zoom operation
+        float t = ease_in_out(((r->t/(2.*M_PI) - 2.0)/2.0));        // runs from 0 to 1
+        float s = 1.+30.0*t;                                        // zoom factor 
+        view = reb_mat4df_scale(view, s,s,s);                       // zoom in
+
+    }else if (r->t > 4.*2.*M_PI && r->t < 5.*2.*M_PI){             // Year 4 to 5
+
+
+        // Reduce timestep
+        r->dt = dt0/10.0; 
+        // Hide wires
+        r->display_settings->wire = 0;
+        if (r->N==9){
+            // Add moon (these are not exact parameters, just for illustration)
+            struct reb_particle e = r->particles[3];
+            reb_simulation_add_fmt(r, "a m r primary", 0.0025695553, 3.6943033e-08, 1.1617812e-05*inflate_size, e);
+        }
+
+        // Get Earth-Moon Barycenter
+        struct reb_particle em_com = reb_particle_com_of_pair(r->particles[3], r->particles[9]); 
+
+        // Apply same rotation and zoom as before (we could cache this)
+        struct reb_mat4df rm = reb_rotation_to_mat4df(rot_x);
+        view = reb_mat4df_multiply(rm, view); 
+        view = reb_mat4df_scale(view, 31.,31.,31.);
+
+        // slowly zoom in
+        float ts = ease_in_out((r->t/(2.*M_PI) - 4.0));     // runs from 0 to 1
+        float s = 1.+500.0*ts*ts;                                   // zoom factor 
+        view = reb_mat4df_scale(view, s,s,s);                   // second zoom operation
+
+        // slowly move earth to center
+        float tm = ease_in_out((r->t/(2.*M_PI) - 4.0)*3.);    // runs from 0 to 1
+        view = reb_mat4df_translate(view, -tm*em_com.x, -tm*em_com.y, -tm*em_com.z);    // translate view
+
+    }else{  // Continue until forever
+
+        // Get Earth-Moon Barycenter
+        struct reb_particle em_com = reb_particle_com_of_pair(r->particles[3], r->particles[9]); 
+
+        // Apply same rotation, zoom, and keep centered on earth
+        struct reb_mat4df rm = reb_rotation_to_mat4df(rot_x);
+        view = reb_mat4df_multiply(rm, view); 
+        view = reb_mat4df_scale(view, 31.,31.,31.);             // We could combine the zoom operations
+        view = reb_mat4df_scale(view, 500.,500.,500.);
+        view = reb_mat4df_translate(view, -em_com.x, -em_com.y, -em_com.z);    // translate view
+
+        // Show real size of particles
+        r->display_settings->spheres = 1;
+        // Hide planet trails
+        r->display_settings->past = 0;
+
+    }
+
+    // Store the view matrix
+    r->display_settings->view = view;
+}
+
+int main(int argc, char* argv[]) {
+    struct reb_simulation* r = reb_simulation_create();
+
+    // Solar System initial conditions from NASA Horizons
+    // Units are solar mass, AU, years/2pi.
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 0.9999999999950272, 0.0046524726,
+            -0.007784066163300598, -0.003160235321847074, 0.0002085198739177632, 
+            0.00029808520873668176, -0.0003933583120820104, -3.0307129383172144e-06);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 1.6601208254808336e-07, 1.6313735e-05,
+            -0.0197573020117554, -0.4659490539727674, -0.036512748618975056, 
+            1.3072205503890428, 0.04106840511679157, -0.11649038132810724);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 2.447838287784771e-06, 4.0453784e-05,
+            -0.31310669579085715, -0.6609240473612075, 0.008792509815427053, 
+            1.058788984690766, -0.5006529730168531, -0.06794916369487035);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 3.0404326489511185e-06, 4.2587571e-05,
+            -0.7304336184787301, 0.6677833148611707, 0.0001754893039781267, 
+            -0.6964579069368259, -0.7370893850360778, 4.193835863036079e-05);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 3.2271560828978514e-07, 2.2702195e-05,
+            0.22987379117805698, -1.4195661773490975, -0.035304801938727, 
+            0.833220942267478, 0.2040504545094634, -0.01614911011819359);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 0.0009547919099366768, 0.0004778945,
+            3.278533805383686, 3.7537759447476775, -0.08892301493418556, 
+            -0.3352619068991809, 0.3093717438204439, 0.006217554149365762);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 0.0002858856700231729, 0.0004028667,
+            9.053130670591704, -3.529112154855013, -0.2990863565225466, 
+            0.09969170490350149, 0.30152611855823136, -0.009212249609912534);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 4.366249613200406e-05, 0.00017085136,
+            12.148791396671397, 15.380702941916812, -0.10026566241603418, 
+            -0.18110086740935336, 0.1310657378884879, 0.00283290820462697);
+    reb_simulation_add_fmt(r, "m r x y z vx vy vz", 5.151383772628957e-05, 0.00016553712,
+            29.840954544603573, -1.6778892885804566, -0.6531621601938832, 
+            0.00903845442721642, 0.18328001746914682, -0.003982606081774178);
+
+    reb_simulation_move_to_com(r);
+
+    // Inflate sized for illustration
+    for(int i=0;i<r->N;i++){
+        r->particles[i].r *= inflate_size;
+    }
+
+    // We use the WHFast integrator and a fixed timestep of 2 days
+    r->integrator = REB_INTEGRATOR_WHFAST;
+    r->dt = 2./365.25 *2.0*M_PI;
+
+
+    // Starting the REBOUND visualization server. This
+    // allows you to visualize the simulation by pointing 
+    // your web browser to http://localhost:1234
+    reb_simulation_start_server(r, 1234);
+
+    // Normally the visualization settings are determined by the 
+    // user interface. If we add the display_settings struct to 
+    // the simulation itself, it will overwrite any change the 
+    // user has made and allows us to programatically change any 
+    // settings such as the orientation, zoom, etc. 
+    reb_simulation_add_display_settings(r);
+
+    // Initially, rotate to an edge on view
+    struct reb_vec3d a = {.x=1.};                               // x axis
+    struct reb_rotation rot_x = reb_rotation_init_angle_axis(M_PI_2, a); // pi/2 = 90 degrees
+    struct reb_mat4df rm = reb_rotation_to_mat4df(rot_x);
+    r->display_settings->view = reb_mat4df_multiply(rm, r->display_settings->view);
+
+    // Store initial view matrix and timestep so we can use it as a reference
+    view0 = r->display_settings->view;
+    dt0 = r->dt;
+
+    // The heartbeat function is called once per timestep and handles the 
+    // scripted view changes
+    r->heartbeat = heartbeat;
+
+    // Show particles as points (not as spheres with their real size)
+    r->display_settings->spheres = 0;
+    // Show orbits
+    r->display_settings->wire = 1;
+    // Show 4 past particle positions
+    r->display_settings->past = 1;
+    r->display_settings->past_N = 4;
+
+    // Slow the simulation down to at most 120 timesteps per second.
+    r->usleep = 8333;
+
+    // Then keep running forever.
+    reb_simulation_integrate(r, INFINITY);
+
+    // Cleanup 
+    reb_simulation_free(r);
+}
+