Index: trunk/tools/nc_var_tools.py
===================================================================
--- trunk/tools/nc_var_tools.py	(revision 628)
+++ trunk/tools/nc_var_tools.py	(revision 629)
@@ -15201,14 +15201,13 @@
         if maxradii is not None and rad > maxradii: break
         Nrad = int(values[2])
-        couples = None
+        couples = []
         Nradgood = 0
         if Nrad > 1:
            Ncouples = 0
-           couples = []
            for ic in range(Nrad):
                ix = int(values[3 + ic*2])
                iy = int(values[3 + ic*2 + 1])
                if ix <= Npt and iy <= Npt:
-                   couples.append(np.array([ix,iy], dtype=int))
+                   couples.append(np.array([ix,iy]))
                    Nradgood = Nradgood + 1
                elif ix > Npt and iy > Npt: break
@@ -15217,5 +15216,5 @@
              iy = int(values[4])
              if ix <= Npt and iy <= Npt:
-                 couples = np.array([ix,iy], dtype=int)
+                 couples.append(np.array([ix,iy]))
                  Nradgood = Nradgood + 1
              elif ix > Npt and iy > Npt: break
@@ -15223,5 +15222,8 @@
         if Nradgood > 0:
             Nradii[rad] = Nradgood
-            gridradii[rad] = np.array(couples)
+            finalcouples = np.zeros((Nradgood,2), dtype=int)
+            for ic in range(Nradgood):
+                finalcouples[ic,:] = couples[ic]
+            gridradii[rad] = finalcouples
 
 #    for rad in sorted(gridradii.keys()):
@@ -15240,13 +15242,13 @@
     gridcomb = []
     gridcomb.append([x,y])
-    gridcomb.append([-x,y])
-    gridcomb.append([-x,-y])
-    gridcomb.append([x,-y])
+    if x != 0: gridcomb.append([-x,y])
+    if x != 0 and y != 0: gridcomb.append([-x,-y])
+    if y != 0: gridcomb.append([x,-y])
 
     if x != y:
         gridcomb.append([y,x])
-        gridcomb.append([-y,x])
-        gridcomb.append([-y,-x])
-        gridcomb.append([y,-x])
+        if y != 0: gridcomb.append([-y,x])
+        if x != 0 and y != 0: gridcomb.append([-y,-x])
+        if x != 0: gridcomb.append([y,-x])
 
     return gridcomb
@@ -15260,8 +15262,20 @@
       dx= x size of the domain of values
       dy= y size of the domain of values
+      >>> radpos, SQradpos = radii_points(12,12,1,8,25,25)
+      radpos:
+      [[[ 13.          12.        ]]
+       [[ 12.70710678  12.70710678]]
+       [[ 12.          13.        ]]
+       [[ 11.29289322  12.70710678]]
+       [[ 11.          12.        ]]
+       [[ 11.29289322  11.29289322]]
+       [[ 12.          11.        ]]
+       [[ 12.70710678  11.29289322]]]
+      SQradpos:
+      {1.0: [[13, 12], [11, 12], [12, 13], [12, 11]], 1.41421356237: [[13, 13], [11, 13], [11, 11], [13, 11]]}
     """
     fname = 'radii_points'
 
-    radiipos = np.zeros((Nang, Lrad, 2), dtype = int)
+    radiipos = np.zeros((Nang, Lrad, 2), dtype = np.float)
     SQradiipos = {}
 
@@ -15287,6 +15301,6 @@
 # Creation of a matrix with all the grid points at each time-step
             if iipos >= 0 and iipos < dx-1 and jjpos >= 0 and jjpos < dy-1:
-                radiipos[ia,ip,0] = iipos
-                radiipos[ia,ip,1] = jjpos
+                radiipos[ia,ip,0] = posangle[ip,0]
+                radiipos[ia,ip,1] = posangle[ip,1]
 
 # SQ radii values
@@ -15300,16 +15314,9 @@
         Npairs = len(pairs.flatten())/2
         pairsw = []
-        print 'Lluis ',SQradii[ir],'pairs:',pairs, type([pairs])
         for ip in range(Npairs):
             if Npairs == 0: 
-                print '  Lluis NO pair:', SQradii[ir],':',pairs
                 break
-            elif Npairs == 1: 
-                print '  Lluis single pair:',pairs[ip*2],pairs[ip*2+1]
-                allpairs = grid_combinations(pairs[ip*2],pairs[ip*2+1])
             else:
-                print ip,'  Lluis 1 pair:',pairs[ip,0],pairs[ip,1]
-                allpairs = grid_combinations(pairs[ip,0],pairs[ip,1])
-                
+                allpairs = grid_combinations(pairs[ip,0],pairs[ip,1])                
             for iap in range(len(allpairs)):
                 ppos = allpairs[iap]
@@ -15479,7 +15486,11 @@
             trjradii[it,:,:,:], trjSQradii[it] = radii_points(trajvals[it,0],        \
               trajvals[it,1],radii,Nangle,dimx,dimy)
+            print trjradii[it,:,:,:]
+            print '-----------'
+            print trjSQradii[it,:,:,:]
 
             it = it + 1
             iline = iline + 1
+
             quit()
 
