Add double integral example to trick_utils/SAIntegrator #1058

trick_source/trick_utils/SAIntegrator/README.md

@@ -25,6 +25,7 @@ Some examples of using these integrators can be found in the **examples/** direc
 * [CannonBall](examples/CannonBall/README.md) uses the RK2Integrator.
 * [MassSpringDamper](examples/MassSpringDamper/README.md) uses the EulerCromerIntegrator.
 * [Orbit](examples/Orbit/README.md) uses the EulerCromerIntegrator.
+* [DoubleIntegral](examples/DoubleIntegral/README.md) shows an example of a double integral.
 
 <a id=class-Integrator></a>
 ## class Integrator

trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/DoubleIntegral.cpp

@@ -0,0 +1,53 @@
+#include <math.h>
+#include <stdio.h>
+#include "SAIntegrator.hh"
+
+struct IContext {
+    SA::RK2Integrator* integ;
+    double start;
+    double end;
+};
+
+void deriv_y( double y, double state_y[], double derivs_y[], void* udata) {
+    double& x = state_y[0];
+    derivs_y[0] = 0.0;           // (dx/dy) = 0, Hold x constant.
+    derivs_y[1] = cos(x)*cos(y); // area = ∫ f(X,y) dy
+}
+
+void deriv_x( double x, double state_x[], double derivs_x[], void* udata) {
+    IContext* context = (IContext*)udata;
+    double area = 0.0;
+    double* state_y[2] = { &x, &area };
+    SA::RK2Integrator* integ_y = context->integ;
+    integ_y->set_in_vars(state_y);
+    integ_y->set_out_vars(state_y);
+    integ_y->setTime(context->start);
+    while (integ_y->getTime() <= context->end) {
+        integ_y->load();
+        integ_y->step();
+        integ_y->unload();
+    }
+    derivs_x[0] = area; // volume = ∫ area dx
+}
+
+double doubleIntegral( double x_start, double x_end, double y_start, double y_end ) {
+    double volume = 0.0;
+    double dx = 0.01;
+    double dy = 0.01;
+    double* state_x[1] = { &volume };
+    SA::RK2Integrator integ_y (dy, 2, NULL, NULL, deriv_y, NULL);
+    IContext y_integ_context = {&integ_y, y_start, y_end};
+    SA::RK2Integrator integ_x (dx, 1, state_x, state_x, deriv_x, &y_integ_context);
+    integ_x.setTime(x_start);
+    while (integ_x.getTime()<=x_end) { // x-end
+        integ_x.load();
+        integ_x.step();
+        integ_x.unload();
+    }
+    return volume;
+}
+
+int main ( int argc, char* argv[]) {
+    double vol = doubleIntegral( -M_PI/2.0, M_PI/2.0, -M_PI/2.0, M_PI/2.0 );
+    printf("Volume = %g.\n", vol);
+}

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

@@ -0,0 +1,17 @@
+# Double Integral
+
+This is a code example of how one can do a double integral.
+
+This example code computes the following definate integral :
+
+![Equation](images/Equation.png)
+
+using the SA::RK2Integrator, with dx = dy = 0.01
+
+```
+$ make
+$ ./DoubleIntegral
+Volume = 3.99989.
+```
+
+

trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/makefile

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

trick_source/trick_utils/SAIntegrator/examples/DoubleIntegral/plot_trajectory.py

@@ -0,0 +1,18 @@
+#!/usr/bin/env python
+  
+import matplotlib.pyplot as plt
+import numpy as np
+data = np.genfromtxt('cannon.csv',
+                     delimiter=',',
+                     skip_header=1,
+                     skip_footer=1,
+                     names=['t', 'px', 'py', 'vx', 'vx'],
+                     dtype=(float, float, float, float, float)
+                    )
+
+curve1 = plt.plot(data['px'], data['py'], 'C1-')
+plt.title('Cannonball Trajectory')
+plt.xlabel('pos-x')
+plt.ylabel('pos-y')
+plt.grid(True)
+plt.show()

trick_source/trick_utils/SAIntegrator/examples/makefile

@@ -3,13 +3,16 @@ examples:
 	@make -C CannonBall
 	@make -C MassSpringDamper
 	@make -C Orbit
+	@make -C DoubleIntegral
 
 clean:
 	@make -C CannonBall clean
 	@make -C MassSpringDamper clean
 	@make -C Orbit clean
+	@make -C DoubleIntegral clean
 
 spotless:
 	@make -C CannonBall spotless
 	@make -C MassSpringDamper spotless
 	@make -C Orbit spotless
+	@make -C DoubleIntegral spotless

trick_source/trick_utils/SAIntegrator/include/SAIntegrator.hh

@@ -13,6 +13,7 @@ namespace SA {
         virtual void load() = 0;
         virtual void unload() = 0;
         double getTime() {return time;}
+        void setTime(double t) {time = t;}
     };
 
     typedef void (*derivsFunc)( double t, double state[], double derivs[], void* user_data);
@@ -32,10 +33,11 @@ namespace SA {
     public:
         FirstOrderODEIntegrator( double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func, void* user_data);
         ~FirstOrderODEIntegrator();
-        // void init(double dt, int N, double* in_vars[], double* out_vars[], derivsFunc func);
         virtual void undo_step();
         void load();
         void unload();
+        void set_in_vars( double* in_vars[]);
+        void set_out_vars( double* out_vars[]);
     };
 
     class FirstOrderODEVariableStepIntegrator : public FirstOrderODEIntegrator {

trick_source/trick_utils/SAIntegrator/src/SAIntegrator.cpp

@@ -21,15 +21,29 @@ SA::FirstOrderODEIntegrator::~FirstOrderODEIntegrator() {
 }
 void SA::FirstOrderODEIntegrator::load() {
     reset = false;
-    for (int i=0 ; i<state_size; i++ ) {
+    if (orig_vars != NULL) {
-        istate[i] = *(orig_vars[i]);
+        for (int i=0 ; i<state_size; i++ ) {
+            istate[i] = *(orig_vars[i]);
+        }
+    } else {
+        std::cerr << "Error: SA::FirstOrderODEIntegrator::load(). orig_vars is not set." << std::endl;
     }
 }
 void SA::FirstOrderODEIntegrator::unload() {
-    for (int i=0 ; i<state_size; i++ ) {
+    if (dest_vars != NULL) {
-        *(dest_vars[i]) = ostate[i];
+        for (int i=0 ; i<state_size; i++ ) {
+            *(dest_vars[i]) = ostate[i];
+        }
+    } else {
+        std::cerr << "Error: SA::FirstOrderODEIntegrator::load(). dest_vars is not set." << std::endl;
     }
 }
+void SA::FirstOrderODEIntegrator::set_in_vars( double* in_vars[]){
+    orig_vars = in_vars;
+}
+void SA::FirstOrderODEIntegrator::set_out_vars( double* out_vars[]){
+    dest_vars = out_vars;
+}
 void SA::FirstOrderODEIntegrator::undo_step() {
     if (!reset) {                           // If we've not already reset the integrator, then reset it.
         for (int i=0 ; i<state_size; i++ ) {