diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/DefiniteIntegral.cpp b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/DefiniteIntegral.cpp new file mode 100644 index 00000000..a70be131 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/DefiniteIntegral.cpp @@ -0,0 +1,44 @@ +#include +#include +#include "SAIntegrator.hh" + +struct DefiniteIntegral { + double lower_limit; + double upper_limit; + double result; + double Coeff[4]; + SA::Integrator* integrator; + double evaluate (); +}; +double DefiniteIntegral::evaluate () { + result = 0.0; + integrator->setIndyVar(lower_limit); + while(integrator->getIndyVar() < upper_limit){ + integrator->integrate(); + } + return result; +} +void Order3Polynomial( double x, double state[], double derivs[], void* udata) { + DefiniteIntegral& DI = *(DefiniteIntegral*)udata; + double x2 = x*x; + double x3 = x*x2; + derivs[0] = DI.Coeff[3] * x3 + + DI.Coeff[2] * x2 + + DI.Coeff[1] * x + + DI.Coeff[0]; +} +int main (int argc, char* argv[]) { + unsigned int NSteps = 20; + DefiniteIntegral DI; + DI.Coeff[0] = 7.0; + DI.Coeff[1] = -5.0; + DI.Coeff[2] = 1.0; + DI.Coeff[3] = 2.0; + DI.lower_limit = 0.0; + DI.upper_limit = 2.0; + double dx = (DI.upper_limit - DI.lower_limit) / NSteps; + double* state[1] = { &DI.result }; + DI.integrator = new SA::RK4Integrator(dx, 1, state, Order3Polynomial, &DI); + double result = DI.evaluate(); + printf("Integral = %g.\n", result); +} diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/README.md b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/README.md new file mode 100644 index 00000000..2e850e14 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/README.md @@ -0,0 +1,21 @@ +# DefiniteIntegral + +This is a code example of using SAIntegrator to compute a definite integral. + +This program will use the SA::RK4Integrator compute an integral of the form : + +![Equation](images/Eq1.png) + +It's specifically parameterized to compute : + +![Equation](images/Eq2.png) + +The correct answer is 14+2/3. + +``` +$ make +$ ./DefiniteIntegral +Integral = 14.6667. +``` + + diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq1.png b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq1.png new file mode 100644 index 00000000..e923aabd Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq1.png differ diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq2.png b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq2.png new file mode 100644 index 00000000..4eb8941a Binary files /dev/null and b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/images/Eq2.png differ diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/makefile b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/makefile new file mode 100644 index 00000000..c136ed78 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/makefile @@ -0,0 +1,19 @@ + +RM = rm -rf +CC = cc +CPP = c++ + +CXXFLAGS = -g -Wall +INCLUDE_DIRS = -I../../include +LIBDIR = ../../lib + +all: DefiniteIntegral + +DefiniteIntegral: + $(CPP) $(CXXFLAGS) DefiniteIntegral.cpp ${INCLUDE_DIRS} -L${LIBDIR} -lSAInteg -o DefiniteIntegral + +clean: + ${RM} DefiniteIntegral.dSYM + +spotless: clean + ${RM} DefiniteIntegral diff --git a/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/plot_trajectory.py b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/plot_trajectory.py new file mode 100644 index 00000000..942d1fe8 --- /dev/null +++ b/trick_source/trick_utils/SAIntegrator/examples/DefiniteIntegral/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()