Fixed vortex calculation

projects/braidingTightBinding/include/coupler_braid.hpp

@@ -2,6 +2,7 @@
 #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>{
@@ -40,12 +41,14 @@ class CouplerBraid: public Coupler<value_t, drum_t>{
                         if(params_[drum_index].sublattice == 'B')
                                 v = v - s;
 
-                        value_t ret_val = 1.; //vortices are multiplicative
+                        std::complex<value_t> ret_val = 1.; //vortices are multiplicative
                         for(auto tex: vortices_){
-                                ret_val *= tex(v, s);
+                                ret_val *= tex.distortion(v);
                         }
+                        //add kekule
+                        ret_val *= vortices_.begin()->kekule(v, s);
 
-                        return 1. + ret_val; //return with added offset, as in eliska's code
+                        return 1. + std::real<value_t>(ret_val); //return with added offset, as in eliska's code
                 }
 
         private:

projects/braidingTightBinding/include/vortex.hpp

@@ -2,6 +2,7 @@
 #define VORTEX_HPP_INCLUDED
 #include <vec2.hpp>
 #include <cmath>
+#include <complex>
 
 template <typename value_t>
 class Vortex{
@@ -31,11 +32,12 @@ class Vortex{
 
                 //functional
                 //there is no divide_by_zero check for efficiency reasons. l0_ can not be zero upon call.
-                value_t operator()(const Vec2<value_t>& v, const Vec2<value_t>& s)
+                std::complex<value_t> distortion(const Vec2<value_t>& v) const
                 {
-                        value_t prefactor = delta_*std::tanh(v.r_wrt(position_)/l0_);
-                        value_t cosine = std::cos(alpha_ - v.phi_wrt(position_) + K_*(s + 2.*(v - position_)));
-                        return prefactor*cosine;
+                        return std::polar<value_t>(delta_*std::tanh(v.r_wrt(position_)/l0_), alpha_ - v.phi_wrt(position_));
+                }
+                std::complex<value_t> kekule(const Vec2<value_t>& v, const Vec2<value_t>& s) const{
+                        return std::polar<value_t>(1., K_*(s + 2.*(v - position_)));
                 }
         private:
                 Vec2<value_t> position_;