Changeset 629 in lmdz_wrf for trunk

Timestamp:

Sep 16, 2015, 11:28:35 AM (9 years ago)

Author:

lfita

Message:

Working on the radial averages. Now points function works!

File:

• trunk/tools/nc_var_tools.py (modified) (8 diffs)

trunk/tools/nc_var_tools.py

@@ -15201 +15201 @@
         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 +15216 @@
              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 +15222 @@
         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 +15242 @@
     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 +15262 @@
       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 +15301 @@
 # 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 +15314 @@
         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 +15486 @@
             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()