diff --git a/include/dsp/ode.hpp b/include/dsp/ode.hpp index c94dfe19..f50e3aea 100644 --- a/include/dsp/ode.hpp +++ b/include/dsp/ode.hpp @@ -3,10 +3,14 @@ namespace rack { -typedef void (*stepCallback)(float x, const float y[], float dydt[]); +/** The callback function `f` in each of these stepping functions must have the signature + void f(float x, const float y[], float dydt[]) +A capturing lambda is ideal for this. +*/ -/** Solve an ODE system using the 1st order Euler method */ -inline void stepEuler(stepCallback f, float x, float dx, float y[], int len) { +/** Solves an ODE system using the 1st order Euler method */ +template +void stepEuler(float x, float dx, float y[], int len, F f) { float k[len]; f(x, y, k); @@ -15,8 +19,28 @@ inline void stepEuler(stepCallback f, float x, float dx, float y[], int len) { } } -/** Solve an ODE system using the 4th order Runge-Kutta method */ -inline void stepRK4(stepCallback f, float x, float dx, float y[], int len) { +/** Solves an ODE system using the 2nd order Runge-Kutta method */ +template +void stepRK2(float x, float dx, float y[], int len, F f) { + float k1[len]; + float k2[len]; + float yi[len]; + + f(x, y, k1); + + for (int i = 0; i < len; i++) { + yi[i] = y[i] + k1[i] * dx / 2.f; + } + f(x + dx / 2.f, yi, k2); + + for (int i = 0; i < len; i++) { + y[i] += dx * k2[i]; + } +} + +/** Solves an ODE system using the 4th order Runge-Kutta method */ +template +void stepRK4(float x, float dx, float y[], int len, F f) { float k1[len]; float k2[len]; float k3[len]; @@ -26,14 +50,14 @@ inline void stepRK4(stepCallback f, float x, float dx, float y[], int len) { f(x, y, k1); for (int i = 0; i < len; i++) { - yi[i] = y[i] + k1[i] * dx / 2.0; + yi[i] = y[i] + k1[i] * dx / 2.f; } - f(x + dx / 2.0, yi, k2); + f(x + dx / 2.f, yi, k2); for (int i = 0; i < len; i++) { - yi[i] = y[i] + k2[i] * dx / 2.0; + yi[i] = y[i] + k2[i] * dx / 2.f; } - f(x + dx / 2.0, yi, k3); + f(x + dx / 2.f, yi, k3); for (int i = 0; i < len; i++) { yi[i] = y[i] + k3[i] * dx; @@ -41,7 +65,7 @@ inline void stepRK4(stepCallback f, float x, float dx, float y[], int len) { f(x + dx, yi, k4); for (int i = 0; i < len; i++) { - y[i] += dx * (k1[i] + 2.0 * k2[i] + 2.0 * k3[i] + k4[i]) / 6.0; + y[i] += dx * (k1[i] + 2.f * k2[i] + 2.f * k3[i] + k4[i]) / 6.f; } } diff --git a/include/util/common.hpp b/include/util/common.hpp index 535259c3..b0d0e990 100644 --- a/include/util/common.hpp +++ b/include/util/common.hpp @@ -116,14 +116,14 @@ Example: fclose(file); }); */ -template +template struct DeferWrapper { F f; DeferWrapper(F f) : f(f) {} ~DeferWrapper() { f(); } }; -template +template DeferWrapper deferWrapper(F f) { return DeferWrapper(f); }