Implementations for coupling unit test

include/coupler_constant.hpp

+#ifndef COUPLER_CONSTANT_HPP_INCLUDED
+#define COUPLER_CONSTANT_HPP_INCLUDED
+#include <coupler.hpp>
+#include <vector>
+#include <cmath>
+
+template <typename value_t, typename drum_t>
+class CouplerConstant: public Coupler<value_t, drum_t>{
+        public:
+                CouplerConstant(const value_t voltage): voltage_(voltage), time_(0.) {}
+                ~CouplerConstant() = default;
+
+                void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override { }
+
+                void step(value_t dt) noexcept final override
+                {
+                        time_ += dt;
+                }
+
+                value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
+                {
+                        return voltage_;
+                }
+        
+        private:
+                value_t voltage_;
+                value_t time_;
+};
+
+#endif

include/driver_single.hpp

+#ifndef DRIVER_SINGLE_HPP_INCLUDED
+#define DRIVER_SINGLE_HPP_INCLUDED
+#include <driver.hpp>
+#include <vector>
+#include <cmath>
+#include <iostream>
+
+template<typename value_t, typename drum_t>
+class DriverSingle: public Driver<value_t, drum_t>{
+        public:
+                DriverSingle(const value_t amplitude, const value_t frequency, const value_t holdoff_time, const size_t index)
+                        : amplitude_(amplitude), frequency_(frequency), holdoff_time_(holdoff_time), index_(index), time_(0.) {}
+                ~DriverSingle() = default;
+
+                void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
+                {
+                        //precompute all values
+                        precomp_.reserve(static_cast<int>(2.*t_end/dt)+2.);
+                        for(value_t t_pc = 0.; t_pc <= t_end*1.1; t_pc += dt/2.){
+                                precomp_.push_back(std::vector<value_t> ());
+                                precomp_.back().reserve(drum_vec.size());
+                                for(size_t i = 0; i < drum_vec.size(); ++i){
+                                        if(i == index_ && t_pc >= holdoff_time_) [[unlikely]]
+                                                precomp_.back().push_back(amplitude_*(std::sin(frequency_*2.*M_PI*(t_pc-holdoff_time_) - M_PI/2.) + 1.));
+                                        else
+                                                precomp_.back().push_back(0);
+                                }
+                        }
+                        current_ = 0;
+                }
+
+                void step(value_t dt) noexcept final override{
+                        time_ += dt;
+                        ++current_;
+                }
+
+                value_t operator()(const size_t drum_index) const noexcept final override{
+                        return precomp_[current_][drum_index];
+                }
+
+        private:
+                size_t index_;
+                value_t amplitude_;
+                value_t frequency_; //frequency of drive
+                value_t holdoff_time_; //initial waiting time
+                value_t time_; //current simulation time
+
+                //precomputation variables
+                std::vector<std::vector<value_t>> precomp_;
+                size_t current_;
+};
+
+#endif

include/line_generator.hpp

+#ifndef LINE_GENERATOR_HPP_INCLUDED
+#define LINE_GENERATOR_HPP_INCLUDED
+#include <lattice_generator.hpp>
+#include <utility>
+#include <vector>
+#include <vec2.hpp>
+
+template <typename value_t, typename params_t, typename vars_t, typename sbuffer_t>
+class LineGenerator: public LatticeGenerator<value_t, params_t, vars_t, sbuffer_t> {
+        public:
+                LineGenerator(size_t length) noexcept : length_(length){}
+                LineGenerator() = delete;
+                ~LineGenerator() = default;
+
+                //the params_t argument shall contain the position of the first drum that is placed.
+                //the last drum (return.first.back()) signifies inexistent neighbours, it is not part of the lattice.
+                //It shall be of sublattice 'A'.
+                std::pair<std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> >, std::vector<std::vector<int> > > 
+                        operator()(const params_t& drum_params) noexcept final override
+                {
+                        using drum_t = Drum<value_t, params_t, vars_t, sbuffer_t>;
+
+                        //we need drum parameters for A and B sublattices
+                        params_t params_A(drum_params);
+                        params_t params_B(params_A.a, params_A.c, params_A.f, params_A.m, params_A.Q, params_A.gap, params_A.coulomb_prefactor, params_A.position, 'B'); //generate sublattice 'B' type of parameters
+
+                        //connecting vectors
+                        Vec2<value_t> vec_right_up = drum_params.s2;
+                        Vec2<value_t> vec_right_dn = -1. * drum_params.s1;
+
+                        //prepare memory
+                        drums_.reserve(length_);
+                        adjacency_vector_.reserve(length_);
+
+                        //push first drum
+                        drums_.push_back(drum_t(params_A));
+                        adjacency_vector_.push_back(std::vector<int>(3,-1));
+
+                        //initialize next translation vector and next drum parameters
+                        Vec2<value_t> next_tv = vec_right_up;
+                        params_t next_params = params_B;
+
+                        //whoever constructs this object with length=0 is an idiot themselves
+                        for(size_t i = 0; i < length_-1; ++i){
+                                next_params.position = drums_.back().get_parameters().position + next_tv;
+                                drums_.push_back(drum_t(next_params));
+                                adjacency_vector_.push_back(std::vector<int>(3,-1));
+                                if(next_params.sublattice == 'B' || next_params.sublattice == 'b'){
+                                        //update adjacency
+                                        auto it = adjacency_vector_.end();
+                                        --it;
+                                        (*it)[2] = std::distance(adjacency_vector_.begin(), it)-1;
+                                        --it;
+                                        (*it)[2] = std::distance(adjacency_vector_.begin(), it)+1;
+
+                                        //update next parameter and translation vector
+                                        next_tv = vec_right_dn;
+                                        next_params = params_A;
+                                }
+                                //if this else doesn't work, it's the user's fault
+                                else{
+                                        //update adjacency
+                                        auto it = adjacency_vector_.end();
+                                        --it;
+                                        (*it)[1] = std::distance(adjacency_vector_.begin(), it)-1;
+                                        --it;
+                                        (*it)[1] = std::distance(adjacency_vector_.begin(), it)+1;
+
+                                        //update next parameter and translation vector
+                                        next_tv = vec_right_up;
+                                        next_params = params_B;
+                                }
+                        }
+
+                        //push a drum to the end that has no neighbours and signifies the inexistent neighbour
+                        drums_.push_back(drum_t(drum_params));
+                        adjacency_vector_.push_back(std::vector<int>(3,-1));
+                        //make no-neighbours point here
+                        for(size_t i = 0; i < adjacency_vector_.size()-1; ++i){
+                                for(size_t j = 0; j < 3; ++j){
+                                        if(adjacency_vector_[i][j] == -1){
+                                                adjacency_vector_[i][j] = drums_.size()-1;
+                                        }
+                                }
+                        }
+
+			return std::pair<std::vector<drum_t>, std::vector<std::vector<int> > > (drums_, adjacency_vector_);
+                }
+
+        private:
+                size_t length_; //Number of equivalent hexagon rows in x and y direction
+                std::vector<Drum<value_t, params_t, vars_t, sbuffer_t> > drums_;
+                std::vector<std::vector<int> > adjacency_vector_;
+};
+
+#endif