Add periodic impulse to MassSpringDamper sim. Improve README. Fix plot script.

trick_source/trick_utils/SAIntegrator/examples/MassSpringDamper/MassSpringDamper.cpp

@@ -1,20 +1,32 @@
 #include <math.h>
 #include <stdio.h>
 #include "SAIntegrator.hh"
-#define N_DIMENSIONS 1
+
 struct MassSpringDamper {
     double pos;
     double vel;
     double k;     // Spring constant
     double c;     // Damping constant
     double mass;
+    double extForce;
+    MassSpringDamper(double k, double c, double m);
 };
-void init_state ( MassSpringDamper& msd ) {
-    msd.pos = 2.0;
-    msd.vel = 0.0;
-    msd.k = 1.0;
-    msd.c = 0.2;
-    msd.mass = 5.0;
+MassSpringDamper::MassSpringDamper( double k, double c, double m)
+: pos(0.0), vel(0.0), k(k), c(c), mass(m), extForce(0.0) {}
+
+struct ImpulseGenerator {
+    double period;
+    double duration;
+    double force;
+    ImpulseGenerator(double period, double duration, double force);
+    double impulse(double time);
+};
+ImpulseGenerator::ImpulseGenerator(double period, double duration, double force)
+: period(period), duration(duration), force(force) {}
+double ImpulseGenerator::impulse(double time) {
+    double tt = fmod(time, period);
+    if ((tt >= 0.0) && (tt < duration)) return force;
+    return 0.0;
 }
 void print_header() {
     printf ("time, msd.pos, msd.vel\n");
@@ -25,19 +37,26 @@ void print_state( double t, MassSpringDamper& msd ) {
 void G( double t, double x[], double v[], double g_out[], void* udata) {
     MassSpringDamper* msd = (MassSpringDamper*)udata;
     g_out[0] = -(msd->k/msd->mass) * x[0]
-               -(msd->c/msd->mass) * v[0];
+               -(msd->c/msd->mass) * v[0]
+               + msd->extForce ;
 }
+#define N_DIMENSIONS 1
 int main ( int argc, char* argv[]) {
-    MassSpringDamper msd;
+    const double mass = 1.0; // kg
+    const double frequency = 2; // Hz
+    const double dampingConstant = 5.0; // kg/s = Ns/m
+    double springConstant = 4.0 * M_PI * M_PI * frequency * frequency * mass; // kg/s^2 = N/m
+    MassSpringDamper msd(springConstant, dampingConstant, mass);
+    ImpulseGenerator ImpGen(5.0, 0.005, 500.0); // 500 Newtons for 0.005 seconds every 5 seconds.
     double time = 0.0;
     double* x_p[N_DIMENSIONS] = { &msd.pos };
     double* v_p[N_DIMENSIONS] = { &msd.vel };
     double dt = 0.001;
-    init_state(msd);
     print_header();
     print_state(time, msd);
     SA::EulerCromerIntegrator integ(dt, N_DIMENSIONS, x_p, v_p, G, &msd);
-    while (time < 50) {
+    while (time < 10) {
+        msd.extForce = ImpGen.impulse(time);
         integ.integrate();
         time = integ.getIndyVar();
         print_state(time, msd);

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

@@ -1,23 +1,45 @@
-# MassSpringDamper
+## MassSpringDamper
 
-This program uses the SA::EulerCromerIntegrator  to simulate a simple mass-spring-damper.
+The MassSpringDamper program uses the SA::EulerCromerIntegrator class to
+simulate a forced mass-spring-damper system.
 
-Generate the results as follows:
+### Initial conditions:
+
+* mass = 1 kg.
+* spring constant = 157.913670 N/m, gives a frequency of 2 hz.
+* damping constant = 5 Ns/m.
+
+### Run-time:
+
+Every 5 seconds, the system is given an impulse of 500 Newtons for 0.005 seconds.
+
+After each numerical integration time-step, the simulation program prints
+
+1. time (s)
+2. position (m)
+3. velocity (m/s)
+
+to ```stdout```, in Comma Separated Values (CSV) format.
+
+### Building & Running the Simulation Program
+To compile MassSpringDamper:
 
 ```
 $ make
+```
+To run it, and generate a CSV file:
+
+```
 $ ./MassSpringDamper > msd.csv
 ```
 
-Plot the results as follows:
+### Plotting the Results
+The Python script, ```plot_position.py``` is provided to plot the results in ```msd.csv``` using (Python) matplotlib.
+
+Plot position vs. time as follows:
 
 ```
 $ python plot_position.py
 ```
 
-and
-
-```
-$ python plot_velocity.py
-```
 ![MSD](images/MSD.png)

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

@@ -12,7 +12,7 @@ data = np.genfromtxt('msd.csv',
 
 curve1 = plt.plot(data['t'], data['pos'], 'C1-')
 plt.title('MSD Position')
-plt.xlabel('position')
-plt.ylabel('time')
+plt.xlabel('time')
+plt.ylabel('position')
 plt.grid(True)
 plt.show()