Added utility for taking magnetization moments from distribution data.HEADmaster

1 files changed, 155 insertions, 0 deletions

diff --git a/process_mags.cpp b/process_mags.cpp
new file mode 100644
index 0000000..ab7112f
--- /dev/null
+++ b/process_mags.cpp
@@ -0,0 +1,155 @@
+
+#include <cinttypes>
+#include <cmath>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include <iostream>
+#include <iomanip>
+
+signed M(unsigned N, unsigned i) {
+  return -2 * ((signed)(N - 1) / 2 - (signed)i);
+}
+
+int main(int argc, char *argv[]) {
+  std::vector<std::string> all_args;
+
+  if (argc > 1) {
+    all_args.assign(argv + 1, argv + argc);
+  }
+
+  for (std::string filename : all_args) {
+    std::ifstream in_file;
+    in_file.open(filename);
+
+    if (!in_file.is_open()) {
+      std::cout << "Failed to open " << filename << ".\n";
+    } else {
+
+      std::vector<uint64_t> dat;
+
+      std::string line;
+      while (std::getline(in_file, line)) {
+        std::istringstream iss(line);
+        uint64_t a;
+        while (iss >> a) {
+          dat.push_back(a);
+        }
+      }
+
+      in_file.close();
+
+      double Mt = 0;
+      double δMt = 0;
+
+      double Ma = 0;
+      double δMa = 0;
+
+      double Me = 0;
+      double δMe = 0;
+
+      double Mm = 0;
+      double δMm = 0;
+
+      double χt = 0;
+      double δχt = 0;
+
+      double χa = 0;
+      double δχa = 0;
+
+      double χe = 0;
+      double δχe = 0;
+
+      double χm = 0;
+      double δχm = 0;
+
+      uint64_t normt = 0;
+      uint64_t norme = 0;
+      uint64_t normm = 0;
+
+      for (unsigned i = 0; i < dat.size(); i++) {
+        int64_t MM = M(dat.size(), i);
+        int64_t ni = dat[i];
+        normt += ni;
+
+        Mt += ni * MM;
+        δMt += ni * pow(MM, 2);
+
+        Ma += ni * std::abs(MM);
+        δMa += ni * pow(MM, 2);
+
+        if (i + 1 >= (dat.size()+1) / 2) {
+          norme += ni;
+          Me += ni * MM;
+          δMe += ni * pow(MM, 2);
+        }
+
+        if (i + 1 <= (dat.size()+1) / 2) {
+          normm += ni;
+          Mm += ni * MM;
+          δMm += ni * pow(MM, 2);
+        }
+      }
+
+      Mt /= normt;
+      Ma /= normt;
+      Me /= norme;
+      Mm /= normm;
+
+      δMt = sqrt(δMt) / normt;
+      δMa = sqrt(δMa) / normt;
+      δMe = sqrt(δMe) / norme;
+      δMm = sqrt(δMm) / normm;
+
+      for (unsigned i = 0; i < dat.size(); i++) {
+        int64_t MM = M(dat.size(), i);
+        int64_t ni = dat[i];
+
+        χt += ni * MM * (MM - Mt);
+        δχt += ni * pow(MM, 2) * (pow(MM - Mt, 2) + ni * pow(δMt, 2));
+
+        χa += ni * std::abs(MM) * (std::abs(MM) - Ma);
+        δχa += ni * pow(MM, 2) * (pow(std::abs(MM) - Ma, 2) + ni * pow(δMa, 2));
+
+        if (i + 1 >= (dat.size()+1) / 2) {
+          χe += ni * MM * (MM - Me);
+          δχe += ni * pow(MM, 2) * (pow(MM - Me, 2) + ni * pow(δMe, 2));
+        }
+
+        if (i + 1 <= (dat.size()+1) / 2) {
+          χm += ni * MM * (MM - Mm);
+          δχm += ni * pow(MM, 2) * (pow(MM - Mm, 2) + ni * pow(δMm, 2));
+        }
+      }
+
+      χt /= normt;
+      χa /= normt;
+      χe /= norme;
+      χm /= normm;
+
+      δχt = sqrt(δχt) / normt;
+      δχa = sqrt(δχa) / normt;
+      δχe = sqrt(δχe) / norme;
+      δχm = sqrt(δχm) / normm;
+
+      std::ofstream out_file;
+      out_file.open(filename + ".processed");
+
+      out_file << Me << " " << δMe << "\n";
+      out_file << Mm << " " << δMm << "\n";
+      out_file << Mt << " " << δMt << "\n";
+      out_file << Ma << " " << δMa << "\n";
+      out_file << χe << " " << δχe << "\n";
+      out_file << χm << " " << δχm << "\n";
+      out_file << χt << " " << δχt << "\n";
+      out_file << χa << " " << δχa << "\n";
+
+      out_file.close();
+
+      std::cout << "Finished processing " << filename << ".\n";
+    }
+  }
+
+  return 0;
+}