#ifndef COUPLER_BRAID_HPP_INCLUDED
#define COUPLER_BRAID_HPP_INCLUDED
#include <vector>
#include <vortex.hpp>
#include <complex>

template <typename value_t, typename drum_t, typename params_t>
class CouplerBraid: public Coupler<value_t, drum_t>{
public:
  //constructors
  CouplerBraid(const std::vector<Vortex<value_t> >& vortices): vortices_(vortices) {}
  ~CouplerBraid() = default;

  void precompute(const value_t t_end, const value_t dt, const std::vector<drum_t>& drum_vec) noexcept final override
  {
    //we need the drum parameters later on, so extract it
    params_.reserve(drum_vec.size());
    for(auto d: drum_vec){
      params_.push_back(d.get_parameters());
    }

    //we also need the neighbour vectors nicely indexed
    for(auto d: drum_vec){
      if(d.get_parameters().sublattice == 'A'){
        n_vecs_.push_back(d.get_parameters().s0);
        n_vecs_.push_back(d.get_parameters().s1);
        n_vecs_.push_back(d.get_parameters().s2);
        break;
      }
    }
  }
  void step(value_t dt) noexcept final
  {
    //move vortices or alpha here
  }
  value_t operator()(const size_t drum_index, const size_t neighbour_index) const noexcept final override
  {
    Vec2<value_t> v = params_[drum_index].position;
    Vec2<value_t> s = n_vecs_[neighbour_index];

    //the formula is only valid on the 'A' sublattice
    if(params_[drum_index].sublattice == 'B'){
      v = v - s;
    }

    std::complex<value_t> ret_val = 1.; //vortices are multiplicative
    for(auto tex: vortices_){
      ret_val *= tex.distortion(v);
    }
    //add kekule
    ret_val *= vortices_.begin()->kekule(v, s);

    return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
  }

private:
  std::vector<Vortex<value_t> > vortices_;
  std::vector<params_t> params_;
  std::vector<Vec2<value_t> > n_vecs_;
};

#endif