SAIntegrtor: Add example sim for RKF45 called AsteroidFlyBy. #1114

trick_source/trick_utils/SAIntegrator/examples/AsteroidFlyBy/Flyby.cpp

@@ -0,0 +1,59 @@
+#include <math.h>
+#include <stdio.h>
+#include "SAIntegrator.hh"
+
+#define GRAVITATIONAL_CONSTANT 6.674e-11
+#define EARTH_MASS 5.9723e24
+#define EARTH_RADIUS 6367500.0
+
+struct Flyby {
+    double pos[2];
+    double vel[2];
+    double planet_mass;
+    Flyby(double px, double py, double vx, double vy, double m);
+};
+Flyby::Flyby(double px, double py, double vx, double vy, double m) {
+    pos[0] = px;
+    pos[1] = py;
+    vel[0] = vx;
+    vel[1] = vy;
+    planet_mass = m;
+}
+void print_header() {
+    printf ("time, dt, flyby.pos[0], flyby.pos[1], flyby.vel[0], flyby.vel[1]\n");
+}
+void print_state( double t, double dt, Flyby& flyby ) {
+    printf ("%10.10f, %10.10f, %10.10f, %10.10f, %10.10f, %10.10f\n",
+             t, dt, flyby.pos[0], flyby.pos[1], flyby.vel[0], flyby.vel[1]);
+}
+void G( double t, double* state, double derivs[], void* udata) {
+    Flyby* flyby = (Flyby*)udata;
+    double d = sqrt( state[0]*state[0] + state[1]*state[1]);
+    derivs[0] = state[2];
+    derivs[1] = state[3];
+    derivs[2] = -state[0] * GRAVITATIONAL_CONSTANT * flyby->planet_mass / (d*d*d);
+    derivs[3] = -state[1] * GRAVITATIONAL_CONSTANT * flyby->planet_mass / (d*d*d);
+}
+int main ( int argc, char* argv[]) {
+
+    double sim_duration = 25000.0; // s
+    double dt = 60.0; // s
+    double epsilon = 0.000000001;
+    Flyby flyby(-20.0 * EARTH_RADIUS, 2.0 * EARTH_RADIUS, 10000.0 , 0.0, EARTH_MASS);
+    double* state_p[4] = { &flyby.pos[0], &flyby.pos[1], &flyby.vel[0], &flyby.vel[1] };
+
+    double time = 0.0; // s
+    print_header();
+    print_state(time, dt, flyby);
+
+    SA::RKF45Integrator integ(epsilon, dt, 4, state_p, G, &flyby);
+    while (time < sim_duration) {
+        integ.integrate();
+        double last_h = integ.getLastStepSize();
+        time = integ.getIndyVar();
+
+        double r = sqrt( flyby.pos[0]*flyby.pos[0] + flyby.pos[1]*flyby.pos[1]);
+        if (r < 500000.0) { printf("Collision\n"); }
+        print_state(time, last_h, flyby);
+    }
+}

trick_source/trick_utils/SAIntegrator/examples/AsteroidFlyBy/README.md

@@ -0,0 +1,48 @@
+# Flyby
+
+The Flyby program uses the **SA::RKF45Integrator** class to simulate
+an asteroid passing near Earth.
+
+The RKF45Integrator is an adaptive step-size integrator. It adapts the
+integration step-size to maintain a specified accuracy. If a particular step-size 
+doesn't produce the needed accuracy then the step-size is reduced and the integration step is performed again. If the needed accuracy is being produced then the step-size can be increased. There is some over-head in the extra calculations, that estimate the local-error. But, this can be more than made up for by the fact that the step-size is small **only** when necessary. 
+
+For each numerical integration time-step, the simulation program prints:
+
+1. time (s)
+2. the size of the last time step
+2. 2D position vector (m)
+3. 2D velocity vector (m/s)
+
+to ```stdout```, in Comma Separated Values (CSV) format.
+
+### Building & Running the Simulation Program
+
+Generate the results as follows:
+
+```
+$ make
+$ ./Flyby > flyby.csv
+```
+### Plotting the Results
+The Python script, ```plot_position.py``` is provided to plot the results
+in ```flyby.csv ``` using (Python) matplotlib.
+
+Plot the asteroid path as follows:
+
+```
+$ python plot_position.py
+```
+The following shows the path of the asteroid for 25000 seconds (about 7 hours).
+The asteroid starts about 20 Earth-radii from the Earth, traveling at 10000 meters per second ( about 22000 miles per hour). The Earth is at 0,0.
+![Orbit](images/Figure_1.png)
+
+The normal (maximum) step-size (dt) for this simulation is 60 seconds. As the asteroid approaches Earth, and gravitational acceleration increases, the RKF45Integrator decreases its step-size to maintain accurancy. The step-size reaches a minimum of about 3 seconds when closest to Earth. As the asteroid retreats, the step-size returns to normal.
+
+With RKF45, a max step-size of 60 seconds, and epsilon = 0.000000001, this 25000 second simulation requires 1513 steps. With RK4 and a step-size of 3 seconds (to maintain the required accuracy), this simulation would require about 8300 steps. So, it would appear that the overhead of RKF45 can be a worthwhile investment in time.
+
+In a simulation where the asteroid were mostly flying through open space, and rarely encountering another planet, the payoff would be much bigger.
+
+![Orbit](images/Figure_2.png)
+
+

trick_source/trick_utils/SAIntegrator/examples/AsteroidFlyBy/makefile

@@ -0,0 +1,20 @@
+
+RM = rm -rf
+CC = cc
+CPP = c++
+
+CXXFLAGS = -g -Wall
+INCLUDE_DIRS = -I../../include
+LIBDIR = ../../lib
+
+all: Flyby
+
+Flyby: Flyby.cpp
+	$(CPP) $(CXXFLAGS) Flyby.cpp ${INCLUDE_DIRS} -L${LIBDIR} -lSAInteg -o Flyby
+
+clean:
+	${RM} Flyby.dSYM
+
+spotless: clean
+	${RM} Flyby
+	${RM} flyby.csv

trick_source/trick_utils/SAIntegrator/examples/AsteroidFlyBy/plot_position.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy as np
+data = np.genfromtxt('flyby.csv',
+                     delimiter=',',
+                     skip_header=1,
+                     skip_footer=1,
+                     names=['t', 'dt', 'posx','posy','velx','vely'],
+                     dtype=(float, float, float, float, float, float)
+                    )
+
+curve1 = plt.plot(data['posx'], data['posy'], 'C1-')
+plt.title('Flyby')
+plt.xlabel('position')
+plt.ylabel('position')
+plt.grid(True)
+plt.show()

trick_source/trick_utils/SAIntegrator/examples/AsteroidFlyBy/plot_stepsize.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy as np
+data = np.genfromtxt('flyby.csv',
+                     delimiter=',',
+                     skip_header=1,
+                     skip_footer=1,
+                     names=['t', 'dt', 'posx','posy','velx','vely'],
+                     dtype=(float, float, float, float, float, float)
+                    )
+
+curve1 = plt.plot(data['t'], data['dt'], 'C1-')
+plt.title('Time-step Adaptation')
+plt.xlabel('t')
+plt.ylabel('dt')
+plt.grid(True)
+plt.show()