New files

projects/braidingTightBinding/bin/main_nocoupling.cpp

+#include <coupler.hpp>
+#include <diagonalizer.hpp>
+#include <driver.hpp>
+#include <drum.hpp>
+#include <force.hpp>
+#include <grabber.hpp>
+#include <lattice_generator.hpp>
+#include <rk4_buffer.hpp>
+#include <rk4_stepper.hpp>
+#include <system.hpp>
+#include <system_parameters.hpp>
+#include <vec2.hpp>
+#include <coupler_const.hpp>
+#include <driver_braid.hpp>
+#include <drum_parameters_braid.hpp>
+#include <drum_variables_braid.hpp>
+#include <force_braid.hpp>
+#include <matrix_element_calculator_braid.hpp>
+#include <rbcomb_generator_braid.hpp>
+#include <vortex.hpp>
+#include <string>
+#include <fstream>
+#include <cassert>
+
+template <typename iterator_t>
+void writeToFile(const std::string filename, const std::string header, int line_length, iterator_t begin, iterator_t end) noexcept
+{
+  //no line breaks if line_length == -1
+  if(line_length == -1){
+    line_length = std::distance(begin, end);
+  }
+  std::fstream file (filename, file.out);
+  file << header;
+  file << std::endl;
+  for(auto it = begin; it != end; ++it){
+    if(std::distance(begin, it) % line_length == 0){
+      file << std::endl;
+    }
+    file << *it << " ";
+  }
+  file << std::flush;
+  file.close();
+
+}
+
+int main(int argc, char** argv){
+  assert(argc == 2);
+  //shorthands
+  using value_t = double;
+  using params_t = DrumParametersBraid<value_t>;
+  using vars_t = DrumVariablesBraid<value_t>;
+  using buffer_t = RK4Buffer<value_t>;
+  using drum_t = Drum<value_t, params_t, vars_t, buffer_t>;
+  using coupler_t = CouplerConst<value_t, drum_t, params_t>;
+  using driver_t = DriverBraid<value_t, drum_t>;
+  using sysparams_t = SystemParameters<coupler_t, driver_t>;
+  using force_t = ForceBraid<value_t, params_t, vars_t, buffer_t>;
+  using stepper_t = Rk4Stepper<value_t, params_t, vars_t, buffer_t, force_t>;
+  using generator_t = RbcombGeneratorBraid<value_t, params_t, vars_t, buffer_t>;
+  using grabber_t = Grabber<value_t, drum_t>;
+  using matelecalc_t = MatrixElementCalculatorBraid<value_t, params_t, vars_t, drum_t>;
+  using system_t = System<value_t, drum_t, grabber_t, sysparams_t, force_t, coupler_t, driver_t, stepper_t, matelecalc_t>;
+  using vortex_t = Vortex<value_t>;
+  using vec2_t = Vec2<value_t>;
+
+  value_t k0 = 24674011002.723392;
+  value_t k1 = 0.;
+  value_t k2 = 2467401100.2723393;
+  value_t c = 5.235987755982988;
+
+  params_t drum_parameters (k0, k1, k2, c, vec2_t (0., 0.), 'A');
+  generator_t lattice_generator (37,15);
+  auto lat_pair = lattice_generator(drum_parameters);
+
+  //set initial conditions
+  double x_0 = 1.;
+  double xdot_0 = 0.;
+  for(drum_t& drum: lat_pair.first){
+    drum.get_variables().x = x_0;
+    drum.get_variables().xdot = xdot_0;
+  }
+
+  sysparams_t system_parameters (coupler_t(0.), driver_t(), lat_pair.second);
+  force_t force;
+  stepper_t stepper;
+
+  //generate grabber
+  double period = 0.00004;
+  double points_per_period = 100.;
+  size_t save_every = 10;
+  std::string parameters_file =  std::string(argv[1]) + "PARAMS.TXT";
+  std::string adjacency_file =  std::string(argv[1]) + "ADJACENCY.TXT";
+  std::string dynamics_file = std::string(argv[1]) + "DYNAMICS.TXT";
+  grabber_t grabber (save_every, parameters_file, adjacency_file, dynamics_file);
+
+  system_t system (10000*period, period/points_per_period, lat_pair.first, stepper, force, system_parameters, grabber);
+  system.simulate();
+  system.save();
+
+  //set up matrix element calculator
+  matelecalc_t mec;
+
+  //set up diagonalizer
+  Diagonalizer diagonalizer;
+  //calculate eigenvalues
+  std::vector<value_t> matrix = system.get_matrix(mec);
+  //std::vector<value_t> evals = diagonalizer.ev(system.get_matrix(mec), lat_pair.first.size()-1);
+  std::pair<std::vector<value_t>, std::vector<value_t> > diag_pair = diagonalizer.evv(system.get_matrix(mec), lat_pair.first.size()-1);
+
+  //make eigenvalues to frequencies
+  std::transform(diag_pair.first.begin(), diag_pair.first.end(), diag_pair.first.begin(), [](value_t x) -> value_t {return std::sqrt(x)/(2.*M_PI);});
+
+  //write to file
+  writeToFile(std::string(argv[1]) + "eigenvalues.txt", "", 1, diag_pair.first.begin(), diag_pair.first.end());
+  writeToFile(std::string(argv[1]) + "eigenvectors.txt", "", diag_pair.first.size(), diag_pair.second.begin(), diag_pair.second.end());
+}

projects/braidingTightBinding/bin/plot_states.py

+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+readpath = "../results/data/zerostate/"
+writepath = "../results/plots/zerostate/"
+
+if(len(sys.argv) < 2):
+    print(f"Usage: python3.9 {sys.argv[0]} mode1 mode2 mode3 ...")
+
+else:
+    num_modes = int(len(sys.argv)-1)
+    #eigenvalues
+    ev = np.loadtxt(readpath+"eigenvalues.txt")
+    #eigenvectors
+    data = np.loadtxt(readpath+"eigenvectors.txt")
+    pos_data = np.loadtxt(readpath+"PARAMS.txt", dtype=str) #positions are rows 4 and 5
+    desired_length = len(data[0,:])
+    pos_x = pos_data[0:desired_length,4].astype('float')
+    pos_y = pos_data[:desired_length,5].astype('float')
+
+    plt.figure(figsize=(8,5*num_modes))
+    count = 1
+    for str_i in sys.argv[1:]:
+        i = int(str_i)-1
+        plt.subplot(num_modes, 1, count)
+        cmax = np.amax(data[i,:])
+        cmin = np.amin(data[i,:])
+        if(-cmin > cmax):
+            cmax = -cmin
+        elif(-cmin <= cmax):
+            cmin = -cmax
+        plt.scatter(pos_x, pos_y, c=data[i,:], edgecolors='k', linewidth=0.1, cmap='seismic', s=16, vmin=cmin, vmax=cmax)
+        plt.colorbar()
+        plt.title(f"Mode {i+1}, ({ev[i]} Hz)")
+        plt.xlabel("x")
+        plt.ylabel("y")
+        count += 1
+
+name = "Modes"
+for str_i in sys.argv[1:]:
+    name += f"_{str_i}"
+plt.savefig(writepath + name, dpi=500, bbox_inches='tight', pad_inches=0)