Add copying/readme

.gitignore

 /builddir
 /build
+/Makefile
 *~

README.md

+# usadelndsoc
+
+Solver for Usadel equations in 2D with SOC gauge fields.
+
+Uses a discretization with local gauge invariance.
+
+## Requirements
+
+- C++ compiles
+- C++ libraries: Eigen, CppAD
+- Meson (build requirement)
+
+## Compiling
+
+See https://mesonbuild.com/Quick-guide.html#compiling-a-meson-project
+

src/action.hpp

+/* usadelndsoc
+ *
+ * Copyright © 2022 Pauli Virtanen
+ *    @author Pauli Virtanen <pauli.t.virtanen@jyu.fi>
+ *
+ * SPDX-License-Identifier: AGPL-3.0-or-later
+ */
 #ifndef ACTION_HPP_
 #define ACTION_HPP_
 
@@ -5,8 +12,18 @@
 #include "array.hpp"
 
 template <typename Scalar, typename Shape>
-inline Scalar S_2(array::Array<Scalar, Shape> Q)
+inline auto get_Q_matrix(array::Array<Scalar, Shape> Q)
 {
+    return 0;
+}
+
+template <typename Scalar, typename Shape, typename Shape2>
+inline Scalar S_2(array::Array<Scalar, Shape> Q, array::Array<Complex, Shape2> U, double h)
+{
+    static const int N_i[4] = {0, 1, -1, 0};
+    static const int N_j[4] = {1, 0, 0, -1};
+    Scalar S = Scalar(0);
+
     size_t nx = Q.dim(0);
     size_t ny = Q.dim(1);
 
@@ -15,29 +32,32 @@ inline Scalar S_2(array::Array<Scalar, Shape> Q)
             auto Q1 = Q.part(i,j).matrix();
 
             /* Sum over neighbors */
-            for (int di = -1; di <= 1; ++di) {
-                for (int dj = -1; dj <= 1; ++dj) {
-                    if (i == 0 && di < 0)
+            for (size_t k = 0; k < 4; ++k) {
+                size_t di = N_i[k];
+                size_t dj = N_j[k];
+                
+                if (i == 0 && di < 0)
                         continue;
-                    if (i == nx-1 && di > 0)
+                if (i == nx-1 && di > 0)
                         continue;
-                    if (j == 0 && dj < 0)
-                        continue;
-                    if (j == ny-1 && dj > 0)
-                        continue;
-                    if (di == 0 && dj == 0)
+                if (j == 0 && dj < 0)
                         continue;
+                if (j == ny-1 && dj > 0)
+                    continue;
+
+                size_t i2 = i + di;
+                size_t j2 = j + dj;
 
-                    size_t i2 = i + di;
-                    size_t j2 = j + dj;
+                auto Q2 = Q.part(i2,j2).matrix();
+                auto U12 = U.part(i,j,k);
+                auto U21 = U.part(j,i,k);
 
-                    auto Q2 = Q.part(i2,j2).matrix();
-                }
+                S += 2/(4*h) * (Q1 - U12*Q2*U21).trace();
             }
         }
     }
 
-    return Scalar(0);
+    return S;
 }
 
 

src/common.hpp

+/* usadelndsoc
+ *
+ * Copyright © 2022 Pauli Virtanen
+ *    @author Pauli Virtanen <pauli.t.virtanen@jyu.fi>
+ *
+ * SPDX-License-Identifier: AGPL-3.0-or-later
+ */
 #ifndef COMMON_HPP_
 #define COMMON_HPP_
 
@@ -8,6 +15,7 @@
 
 #include <Eigen/Core>
 
+typedef std::complex<double> Complex;
 typedef CppAD::AD<double> ADDouble;
 typedef std::complex<ADDouble> ADComplex;