diff --git a/meson.build b/meson.build
index 961527543bb6db6e22b9bc303e76c78ccfb17863..928852aac311eedb009061edc1b981838a93a929 100644
--- a/meson.build
+++ b/meson.build
@@ -1,5 +1,8 @@
-project('usadelndsoc', 'cpp')
+project('usadelndsoc', 'cpp',
+        default_options : ['cpp_std=c++17'])
 
 cppad_dep = dependency('cppad')
-deps = [cppad_dep]
-executable('main', 'src/main.cpp', dependencies: deps)
+eigen3_dep = dependency('eigen3')
+
+deps = [cppad_dep, eigen3_dep]
+executable('main', 'src/main.cpp', dependencies : deps)
diff --git a/src/action.hpp b/src/action.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..6d9127f8b1f65122a62754c198a998e885c52ded
--- /dev/null
+++ b/src/action.hpp
@@ -0,0 +1,15 @@
+#ifndef ACTION_HPP_
+#define ACTION_HPP_
+
+#include "common.hpp"
+
+
+inline ADComplex S_2(Array<ADComplex,3> Q)
+{
+    size_t i = 2;
+    return Q(0u, i, 2u);
+}
+
+
+
+#endif
diff --git a/src/common.h b/src/common.h
deleted file mode 100644
index de55b0f15cc4a83074279ecb7e6a6b043b2719d1..0000000000000000000000000000000000000000
--- a/src/common.h
+++ /dev/null
@@ -1,9 +0,0 @@
-#ifndef COMMON_H_
-#define COMMON_H_
-
-#include <cppad/cppad.hpp>
-
-typedef CppAD::AD<double> ADDouble;
-typedef std::complex<ADDouble> ADComplex;
-
-#endif
diff --git a/src/common.hpp b/src/common.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..12d7b98fd98e2866e9f3e2807c783a629ade2aa2
--- /dev/null
+++ b/src/common.hpp
@@ -0,0 +1,62 @@
+#ifndef COMMON_H_
+#define COMMON_H_
+
+#include <array>
+#include <tuple>
+#include <iostream>
+
+#include <cppad/cppad.hpp>
+#include <cppad/example/cppad_eigen.hpp>
+
+#include <Eigen/Core>
+
+typedef CppAD::AD<double> ADDouble;
+typedef std::complex<ADDouble> ADComplex;
+
+template <typename Scalar, int NDim>
+class Array
+{
+private:
+    std::vector<Scalar>& data_;
+    std::array<size_t, NDim> shape_;
+    std::array<size_t, NDim> stride_;
+
+protected:
+    template <typename... Idx>
+    size_t index(Idx... idxs) const
+        {
+            const std::array<size_t,NDim> m{idxs...};
+            size_t idx = 0;
+
+            for (size_t i = 0; i < NDim; ++i)
+                idx += stride_[i] * m[i];
+
+            return idx;
+        }
+    
+public:
+    Array(std::vector<Scalar>& data, std::array<size_t, NDim> shape) : data_(data), shape_(shape)
+        {
+            if (NDim > 0) {
+                size_t i = NDim - 1;
+                stride_[i] = 1;
+                while (i-- > 0)
+                    stride_[i] = shape[i+1] * stride_[i+1];
+            }
+        }
+
+    template <typename... Idx>
+    Scalar& operator()(Idx... idxs)
+        {
+            return data_[index(idxs...)];
+        }
+
+    template <typename... Idx>
+    const Scalar&& operator()(Idx... idxs) const
+        {
+            return data_[index(idxs...)];
+        }
+};
+
+
+#endif
diff --git a/src/main.cpp b/src/main.cpp
index 783cb8ed508b5e8aa21bd6a74a022b615f1702ea..7d79f02540a20c72de71579a9277195ee8ba253a 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,8 +1,27 @@
-#include "common.h"
+#include "common.hpp"
+#include "action.hpp"
 
 
 int main()
 {
+    std::vector<double> v(3*3*3);
+    Array<double,3> Q(v,{3,3,3});
+
+    for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; i < 3; ++i) {
+            for (size_t k = 0; i < 3; ++i) {
+                Q(i,j,k) = i+j+k;
+            }
+        }
+    }
+
+    for (size_t i = 0; i < 3; ++i) {
+        for (size_t j = 0; i < 3; ++i) {
+            for (size_t k = 0; i < 3; ++i) {
+                std::cout << Q(i,j,k) << std::endl;
+            }
+        }
+    }
 
     return 0;
 }