--- /dev/null
+#include "neigh.h"
+
+neigh::neigh(const int& dimension, const int& length)
+{
+ nfield = new int[ipow(dimension,length)*2*dimension];
+
+ for (int ipos = 0; ipos < ipow(length,dimension); ipos++)
+ {
+ for (int idir = 0; idir < 2*dimension; idir++)
+ {
+ if ( idir < dimension ) //pos direction
+ {
+ if ( ( ipos/ipow(length,idir%dimension) ) % length == length-1 )
+ nfield[ipos*2*dimension + idir] = ipos - ipow(length,idir%dimension)*(length-1);
+ else
+ nfield[ipos*2*dimension + idir] = ipos + ipow(length,idir%dimension);
+ }
+ else //negative dir
+ {
+ if ( ( ipos/ipow(length,idir%dimension) ) % length == 0 )
+ nfield[ipos*2*dimension + idir] = ipos + ipow(length,idir%dimension)*(length-1);
+ else
+ nfield[ipos*2*dimension + idir] = ipos - ipow(length,idir%dimension);
+ }
+
+ }
+ }
+}
--- /dev/null
+#ifndef NEIGH_H
+#define NEIGH_H
+
+class neigh
+{
+ private:
+ int *nfield;
+ static int ipow(const int& x, const int& p){
+ if (p == 0) return 1;
+ if (p == 1) return x;
+ return x * ipow(x, p-1);
+ }
+ public:
+ neigh(const int& length, const int& dimension);
+ int& operator[] (int i) {return nfield[i];}
+};
+
+#endif