Added rk4 stepper

include/rk4_stepper.hpp

+#ifndef RK4_STEPPER_HPP_INCLUDED
+#define RK4_STEPPER_HPP_INCLUDED
+
+template <typename value_t, template<typename T> class params_t, template<typename R> class vars_t, template<typename Q> class buffer_t, 
+         template <value_t, params_t, vars_t, buffer_t> force_t>
+class Rk4Stepper{
+        public:
+                Rk4Stepper() = default;
+                ~Rk4Stepper() = default;
+
+                //time signifies the state of the system, i.e. the time when the complete step was started.
+                //operations performed at t
+                void step_1(const force_t& force, 
+                                std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums, 
+                                const std::vector<std::vector<int> >& adj_vec,
+                                const value_t dt, 
+                                const value_t time) const noexcept
+                {
+                        //step 1: calculate k1 and l1 for all drums
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+                                drums[i].get_sbuffer().k1 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+                                drums[i].get_sbuffer().l1 = dt*xdot;
+                        }
+                        //update x_temp, xdot_temp
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k1/2.;
+                                drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l1/2.;
+                        }
+                        //now it's time for the owner to update the coupler and driver to t+dt/2.
+
+
+                }
+
+
+                //operations performed at t+dt/2
+                void step_2(const force_t<value_t, params_t, vars_t, buffer_t>& force, 
+                                std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums, 
+                                const std::vector<std::vector<int> >& adj_vec,
+                                const value_t dt, 
+                                const value_t time) const noexcept
+                {
+                        //step 2: calculate k2 and l2 for all drums
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+                                drums[i].get_sbuffer().k2 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+                                drums[i].get_sbuffer().l2 = dt*xdot_temp;
+                        }
+                        //update x_temp, xdot_temp
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k2/2.;
+                                drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l2/2.;
+                        }
+
+
+
+                        //step 3: calculate k3 and l3 for all drums
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+                                drums[i].get_sbuffer().k3 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+                                drums[i].get_sbuffer().l3 = dt*xdot_temp;
+                        }
+                        //update x_temp, xdot_temp
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k3;
+                                drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l3;
+                        }
+                        //now it's time for the owner to update the coupler and driver to t+dt/2.
+                }
+
+
+                //operations performed at t+dt
+                void step_3(const force_t<value_t, params_t, vars_t, buffer_t>& force, 
+                                std::vector<Drum<value_t, params_t, vars_t, buffer_t> >& drums, 
+                                const std::vector<std::vector<int> >& adj_vec,
+                                const value_t dt, 
+                                const value_t time) const noexcept
+                {
+                        //step 4: calculate k4 and l4 for all drums
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                //TODO: Empty neighbours are now signified by adj_vec[i][j] == -1. We could also push a 0-drum to the end and point to that.
+                                drums[i].get_sbuffer().k4 = dt*(force(drums[i], drums[adj_vec[i][0]], drums[adj_vec[i][1]], drums[adj_vec[i][2]], time));
+                                drums[i].get_sbuffer().l4 = dt*xdot_temp;
+                        }
+                        //update x, xdot, x_temp, xdot_temp to t+dt
+                        for(size_t i = 0; i < drums.size()-1; ++i){
+                                drums[i].get_variables().xdot_temp = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+                                drums[i].get_variables().x_temp = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+                                drums[i].get_variables().xdot = drums[i].get_variables().xdot + drums[i].get_sbuffer().k4;
+                                drums[i].get_variables().x = drums[i].get_variables().x + drums[i].get_sbuffer().l4;
+                        }
+                }
+};
+
+#endif