Merge branch 'master' of treefish.org:~/src/latlib

[phys/latlib.git] / neigh.cpp

diff --git a/neigh.cpp b/neigh.cpp
index a25d1e70de0f91deb940b9155c83107cfb183eda..8e4084cce461d4f2318d453bd8e0984f94f0bfa8 100644 (file)
--- a/neigh.cpp
+++ b/neigh.cpp
@@ -1,28 +1,44 @@
 #include "neigh.h"
 
-neigh::neigh(const int& dimension, const int& length)
+#include <iostream>
+
+using namespace std;
+
+int neigh::dirstep(int dir)
 {
-  nfield = new int[ipow(dimension,length)*2*dimension];
+  int step=1;
+  for(int idir=0; idir<dir; idir++) step*=len[idir];
+  return step;
+}
+
+neigh::neigh(const int& dimension, ...)
+{
+  int totsize=1;
+  len = new int[dimension];
+  va_list length;
+
+  va_start(length, dimension);
+  for(int idim=0; idim < dimension; idim++) {
+    len[idim] = va_arg(length, int);
+    totsize *= len[idim];
+  }
+  va_end(length);
+
+  nfield = new int[totsize*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);
-           }
-         
-       }
-    }
+  for( int ipos=0; ipos<totsize; ipos++ )
+    for( int idir=0; idir<2*dimension; idir++) {
+      if( idir < dimension ) { // positive direction
+       if ( ( ipos/dirstep(idir%dimension) ) % len[idir%dimension] == len[idir%dimension]-1 )
+         nfield[ipos*2*dimension + idir] = ipos - dirstep(idir%dimension)*(len[idir%dimension]-1);
+       else
+         nfield[ipos*2*dimension + idir] = ipos + dirstep(idir%dimension);
+      }
+      else { // negative direction
+       if ( ( ipos/dirstep(idir%dimension) ) % len[idir%dimension] == 0 )
+         nfield[ipos*2*dimension + idir] = ipos + dirstep(idir%dimension)*(len[idir%dimension]-1);
+       else
+         nfield[ipos*2*dimension + idir] = ipos - dirstep(idir%dimension);
+      }
+    }  
 }