diff --git a/projects/braidingTightBinding/bin/main_braiding_zerostate.cpp b/projects/braidingTightBinding/bin/main_braiding_zerostate.cpp
index 124af20946d236989e5114e380ac8fb7e355e057..b9e892b274ed10e65e103666e0aae0b92cbe1d64 100644
--- a/projects/braidingTightBinding/bin/main_braiding_zerostate.cpp
+++ b/projects/braidingTightBinding/bin/main_braiding_zerostate.cpp
@@ -25,111 +25,121 @@
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();
+ //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 = CouplerBraid<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 = 493480220.0544679;
- value_t k2 = 2467401100.2723393;
- value_t c = 5.235987755982988;
-
- vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it
- vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it
- vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it
- //std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
- std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
-
- 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);
-
- //change edge potentials to push edge modes around
+ 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 = CouplerBraid<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 = 493480220.0544679;
+ value_t k2 = 2467401100.2723393;
+ value_t c = 5.235987755982988;
+
+ vortex_t vortex0 (vec2_t(10.5,24.), 5., 1., 1., 1.); //as Eliska did it
+ vortex_t vortex1 (vec2_t(31.5,24.), 5., 1., 1., 1.); //as Eliska did it
+ vortex_t vortex2 (vec2_t(52.5,24.), 5., 1., 1., 1.); //as Eliska did it
+ //std::vector<vortex_t> vortices ({vortex0, vortex1, vortex2});
+ std::vector<vortex_t> vortices ({vortex1});
+
+ 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);
+
+ //change edge potentials to push edge modes around
value_t onsite_A = 27388152213.022964;
value_t onsite_B = 21959869792.42382;
- for(size_t i = 0; i < lat_pair.first.size()-1; ++i){
+ for(size_t i = 0; i < lat_pair.first.size()-1; ++i){
if(lat_pair.second[i][0] == lat_pair.first.size() || lat_pair.second[i][1] == lat_pair.first.size() || lat_pair.second[i][2] == lat_pair.first.size()){
if(lat_pair.first[i].get_parameters().sublattice == 'A'){ //'A' site
lat_pair.first[i].get_parameters().k0 = onsite_A;
- for(int j: lat_pair.second[i])
- if(j != lat_pair.first.size())
+ for(int j: lat_pair.second[i]){
+ if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_B;
+ }
+ }
}
else{ //'B' site
lat_pair.first[i].get_parameters().k0 = onsite_B;
- for(int j: lat_pair.second[i])
- if(j != lat_pair.first.size())
+ for(int j: lat_pair.second[i]){
+ if(j != lat_pair.first.size()){
lat_pair.first[j].get_parameters().k0 = onsite_A;
+ }
+ }
}
}
}
-
+
+ //set initial conditions
+ double x_0 = 0.;
+ double v_0 = 0.;
+ for(auto drum: lat_pair.first){
+ drum.get_variables().x = x_0;
+ drum.get_variables().xdot = v_0;
+ }
- sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
- force_t force;
- stepper_t stepper;
- //generate grabber
- 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);
+ sysparams_t system_parameters (coupler_t(vortices), driver_t(), lat_pair.second);
+ force_t force;
+ stepper_t stepper;
- system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
- system.save();
+ //generate grabber
+ 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);
- //set up matrix element calculator
- matelecalc_t mec;
+ system_t system (1., 0.00001, lat_pair.first, stepper, force, system_parameters, grabber);
+ system.save();
- //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);
+ //set up matrix element calculator
+ matelecalc_t mec;
- //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);});
+ //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);
- //find the interesting eigenvector range
- size_t ev_begin = 913 * diag_pair.first.size();
- size_t ev_end = 943 * diag_pair.first.size();
+ //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());
+ //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());
}