1 files changed, 100 insertions, 0 deletions
diff --git a/torus_symmetries.hpp b/torus_symmetries.hpp
new file mode 100644
index 0000000..17772de
--- /dev/null
+++ b/torus_symmetries.hpp
@@ -0,0 +1,100 @@
+
+#pragma once
+
+#include <array>
+#include <list>
+#include <vector>
+
+#include "matrix.hpp"
+#include "vector.hpp"
+
+template <int D> void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+  if (level > 0) {
+    unsigned new_level = level - 1;
+    unsigned old_length = sequences.size();
+    for (std::array<double, D>& sequence : sequences) {
+      std::array<double, D> new_sequence = sequence;
+      new_sequence[new_level] = -1;
+      sequences.push_front(new_sequence);
+    }
+    one_sequences<D>(sequences, new_level);
+  }
+}
+
+template <class U, int D> std::vector<Vector<U, D>> torus_vecs(U L) {
+  std::vector<Vector<U, D>> vecs;
+  std::array<double, D> ini_sequence;
+  ini_sequence.fill(1);
+  std::list<std::array<double, D>> sequences;
+  sequences.push_back(ini_sequence);
+  one_sequences<D>(sequences, D);
+  sequences.pop_front(); // don't want the identity matrix!
+
+  for (std::array<double, D> sequence : sequences) {
+    Vector<U, D> v;
+    for (unsigned i = 0; i < D; i++) {
+      if (sequence[i] == 1) {
+        v(i) = 0;
+      } else {
+        v(i) = L / 2;
+      }
+    }
+
+    vecs.push_back(v);
+  }
+
+  return vecs;
+}
+
+template <class U, int D> std::vector<Matrix<U, D>> torus_mats() {
+  std::vector<Matrix<U, D>> mats;
+
+  std::array<double, D> ini_sequence;
+  ini_sequence.fill(1);
+  std::list<std::array<double, D>> sequences;
+  sequences.push_back(ini_sequence);
+
+  one_sequences<D>(sequences, D);
+
+  sequences.pop_front(); // don't want the identity matrix!
+
+  for (std::array<double, D> sequence : sequences) {
+    Matrix<U, D> m;
+    for (unsigned i = 0; i < D; i++) {
+      for (unsigned j = 0; j < D; j++) {
+        if (i == j) {
+          m(i, j) = sequence[i];
+        } else {
+          m(i, j) = 0;
+        }
+      }
+    }
+
+    mats.push_back(m);
+  }
+
+  for (unsigned i = 0; i < D; i++) {
+    for (unsigned j = 0; j < D; j++) {
+      if (i != j) {
+        Matrix<U, D> m;
+        for (unsigned k = 0; k < D; k++) {
+          for (unsigned l = 0; l < D; l++) {
+            if ((k == i && l == j) || (k == j && l == i)) {
+              if (i < j) {
+                m(k, l) = 1;
+              } else {
+                m(k, l) = -1;
+              }
+            } else {
+              m(k, l) = 0;
+            }
+          }
+        }
+        mats.push_back(m);
+      }
+    }
+  }
+
+  return mats;
+}
+