Improve the organization and description of the Orbit example sim for SAIntegrator, in prep for tutorial.

trick_source/trick_utils/SAIntegrator/examples/Orbit/Orbit.cpp

@@ -10,13 +10,14 @@ struct Orbit {
     double pos[2];
     double vel[2];
     double planet_mass;
+    Orbit(double px, double py, double vx, double vy, double m);
 };
-void init_state ( Orbit& orbit ) {
-    orbit.pos[0] = 0.0;
-    orbit.pos[1] = 6578000.0;
-    orbit.vel[0] = 7905.0;
-    orbit.vel[1] = 0.0;
-    orbit.planet_mass = EARTH_MASS;
+Orbit::Orbit(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, orbit.pos[0], orbit.pos[1], orbit.vel[0], orbit.vel[1]\n");
@@ -32,16 +33,25 @@ void G( double t, double pos[], double vel[], double acc_out[], void* udata) {
     acc_out[1] = -pos[1] * GRAVITATIONAL_CONSTANT * orbit->planet_mass / (d*d*d);
 }
 int main ( int argc, char* argv[]) {
-    Orbit orbit;
-    double time = 0.0;
+
+    const double altitude_AGL = 408000.0;            // m
+
+    double mu = GRAVITATIONAL_CONSTANT * EARTH_MASS; // m^3s^2
+    double R_bar_mag = EARTH_RADIUS + altitude_AGL;  // m
+    double V_bar_mag = sqrt(mu / R_bar_mag);         // m/s
+    double Orbital_period = 2 * M_PI * sqrt((R_bar_mag * R_bar_mag * R_bar_mag) / mu); // s
+
+    double dt = 0.1; // s
+    Orbit orbit(0.0, R_bar_mag, V_bar_mag, 0.0, EARTH_MASS);
     double* x_p[2] = { &orbit.pos[0], &orbit.pos[1] };
     double* v_p[2] = { &orbit.vel[0], &orbit.vel[1] };
-    double dt = 0.1;
-    init_state(orbit);
+
+    double time = 0.0; // s
+
     print_header();
     print_state(time, orbit);
     SA::EulerCromerIntegrator integ(dt, 2, x_p, v_p, G, &orbit);
-    while (time < 5600) {
+    while (time < Orbital_period) {
         integ.integrate();
         time = integ.getIndyVar();
         print_state(time, orbit);

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

@@ -1,17 +1,41 @@
 # Orbit
 
-This program uses the *SemiImplicitEuler* integrator to simulate the orbit of an Earth satellite.
+The Orbit program uses the *SA::EulerCromerIntegrator* class to simulate
+the orbit of a satellite around the Earth.
+
+The altitude of the satellite is chosen to be 408000.0 meters, the
+approximate altitude of the International Space Station. The initial
+position, and velocity of the satellite is calculated to maintain a
+circular orbit around the Earth. The orbital period is also calculated,
+to determine how long the simulation should run to complete a full orbit.
+
+For each numerical integration time-step, the simulation program prints:
+
+1. time (s)
+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
 $ ./Orbit > orbit.csv
 ```
+### Plotting the Results
+The Python script, ```plot_position.py``` is provided to plot the results
+in ```orbit.csv``` using (Python) matplotlib.
 
-Plot the results as follows:
+Plot position vector for the orbit duration as follows:
 
 ```
 $ python plot_position.py
 ```
 
 ![Orbit](images/Orbit.png)
+
+### References:
+[https://en.wikipedia.org/wiki/Circular_orbit](https://en.wikipedia.org/wiki/Circular_orbit)