Index: /trunk/UTIL/PYTHON/mcd/examples/dustdevil.py
===================================================================
--- /trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 1196)
+++ /trunk/UTIL/PYTHON/mcd/examples/dustdevil.py	(revision 1197)
@@ -9,4 +9,7 @@
 
 from ppplot import writeascii
+
+
+from scipy import stats
 
 
@@ -53,8 +56,9 @@
 ## constant_heights
 #for heights in [20.,100.,500.,1000.,2000.,5000.,7000.,10000.,20000.]:
-for heights in [10.]:
+#for heights in [10.]:
 
 ## multiple_obsheight
 #for heights in [0.1,0.2,0.5,0.75,1.0,1.5,2.0,10.]:
+for heights in [1.0]:
 
  query = mcd()
@@ -72,8 +76,7 @@
 
     ## constant_heights
-    query.xz = heights
-
+    #query.xz = heights
     ## multiple_obsheight
-    #query.xz = ddheight[i-1]*heights
+    query.xz = ddheight[i-1]*heights
 
     query.update() ; uu = query.zonwind ; vv = query.merwind
@@ -87,8 +90,10 @@
        if angle < 0.: angle = angle + 360.
 
-       #query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
-       #query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
-       query.xzs = max(query.xz - 0.1*query.xz,10.)
-       query.xze = query.xz + 0.1*query.xz
+       ## multiple_obsheight
+       query.xzs = max(ddheight[i-1] - 2.*ddheighterror[i-1],10.)
+       query.xze = ddheight[i-1] + 2.*ddheighterror[i-1]
+       ## constant_heights
+       #query.xzs = max(query.xz - 0.1*query.xz,10.)
+       #query.xze = query.xz + 0.1*query.xz
 
        query.profile()
@@ -130,4 +135,14 @@
  writeascii(field=mcdwind,absc=ddwind,name=str(heights)+'.txt')
 
+ #ind = np.isnan(ddwind)
+ #y = ddwind
+ #y[ind] = -9999.
+ #xi = mcdwind[y > -9999.]
+ #y = ddwind[y > -9999.]
+ #slope, intercept, r_value, p_value, std_err = stats.linregress(xi,y)
+ #linear = slope*xi+intercept
+ #mpl.plot(xi,linear)
+ #ecart = y - linear
+ #print np.std(ecart)
 
  mpl.errorbar(mcdwind, ddwind, yerr=[ddwinderror,ddwinderror], xerr=[mcdwindle,mcdwindhe], fmt='bo')
Index: /trunk/UTIL/PYTHON/powerlaw/ddhist.py
===================================================================
--- /trunk/UTIL/PYTHON/powerlaw/ddhist.py	(revision 1196)
+++ /trunk/UTIL/PYTHON/powerlaw/ddhist.py	(revision 1197)
@@ -4,52 +4,4 @@
 import numpy as np
 import scipy.optimize as sciopt
-
-### SAVED
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = False ; typefit=1
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm2_1.txt" ; drop = True ; typefit=2
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = False ; typefit=1
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc_f18.ncm1_1.txt" ; drop = True ; typefit=1
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt" ; drop = False ; typefit=1 #bof. tous les 10 est mieux.
-#
-###########################################################
-##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_1.txt"
-##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_1.txt"
-##namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_1.txt"
-##namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/press_ustm_exomars.nc1.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/press_ustm_exomars.ncm2_1.txt"
-##namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_1.txt"
-
-
-
-
-case = "188324p"
-case = "191798"
-case = "160564p"
-case = "156487"
-case = "2007p"
-case = "13526p"
-case = "172097"
-namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_1.txt"
-
-
-
-drop = False ; typefit=1
-#drop = True ; typefit=1
-#drop = True ; typefit=2
-##########################################################
-## define bins (we expect fit do not depend too much on this -- to be checked)
-nbins = 7 ; www = 3.
-#nbins = 10 ; www = 2.5
-#nbins = 15 ; www = 2. 
-##nbins = 20 ; www = 1.75
-##nbins = 30 ; www = 1.5
-##nbins = 50 ; www = 1.2
-##########################################################
-limrest = 4. # restrict limit (multiple of dx)
-#limrest = 2.
-#limrest = 3.
-#limrest = 0.
-##########################################################
 
 # -- functions for fitting
@@ -62,99 +14,151 @@
 def fitfunc2(x,a,b): return a*np.exp(-b*x)
 
-# load data
-data = np.loadtxt(namefile,delimiter=";")
-t = data[:,0] ; s = data[:,1] ; d = data[:,2]
-i = data[:,3] ; j = data[:,4] # int?
+################################################################################
+### HISTODD
+### namefile: text file to be analyzed
+### drop: look at size (False) or pressure drop (True)
+### typefit: use fitfunc above 1, 2, or 3
+### nbins: bins (we expect fit do not depend too much on this -- to be checked)
+### limrest: restrict limit (multiple of dx) -- greater equal 
+### addtitle: add a name for the examined case
+################################################################################
+def histodd(namefile,drop=False,typefit=1,nbins=12,limrest=4,addtitle=""):
 
-# a way to guess resolution
-dx = np.min(s)/2. ; print "resolution: ", dx
+    # width adapted to number of bins
+    widthbin = {7:3.,10:2.5,12:2.2,15:2.0,20:1.75,30:1.5,50:1.2}
+    www = widthbin[nbins]
+    
+    # load data
+    data = np.loadtxt(namefile,delimiter=";")
+    t = data[:,0] ; s = data[:,1] ; d = data[:,2]
+    i = data[:,3] ; j = data[:,4] # int?
+    
+    # a way to guess resolution
+    dx = np.min(s)/2.
+    
+    # choose a variable
+    if drop: var = d
+    else: var = s
+    
+    # restrictions
+    restrict = (s > 0) # initialization (True everywhere)
+    restrict = restrict*(s >= limrest*dx) # remove lowest sizes (detection limit) 
+    restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height) pour ~round shape
+    #restrict = restrict*(d > 0.9) # limit on drop for size 
+                                   # (casse la power law? seulement si keep smaller devils)
+    #restrict = restrict*(d > 0.5)
+    #restrict = restrict*(d > 0.3)
+    poum=False
+    if poum:
+      out = var[np.logical_not(restrict)]
+      outi = i[np.logical_not(restrict)]
+      outj = j[np.logical_not(restrict)]
+      print out, outi, outj
+      print out.size
+    var = var[restrict]
+    
+    ## define bins
+    zebins = [np.min(var)]
+    for i in range(0,nbins):  zebins.append(zebins[i]*(www**0.5))
+    zebins = np.array(zebins)
+    middle = 0.5*(zebins[1:] + zebins[:-1])
+    binwidth = zebins[1:] - zebins[:-1]
+    total = var.shape[0]
+    
+    # plot histogram
+    yeah = mpl.hist(var,log=True,bins=zebins,normed=True,color='white') ; mpl.xscale('log')
+    #yeah = mpl.hist(var,bins=zebins,normed=True,color='white')
+    
+    # add error bars
+    incertitude_nvortex = 10.
+    if not drop:
+      err = incertitude_nvortex/(yeah[0]*binwidth*total+0.0001)
+      ind = err < 1000. ; err = err[ind] ; y = yeah[0][ind] ; x = middle[ind]
+      err = y*err
+      mpl.errorbar(x,y,yerr=[err,err],fmt='k.')
+    
+    ## print info
+    #print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
+    #for iii in range(len(zebins)-1):
+    #    this = zebins[iii] ; next = zebins[iii+1]
+    #    print "%5.2e in [%5.0f %5.0f] %5.0f" % (yeah[0][iii],this,next,np.abs(next-this))
+    
+    # fitting function to superimpose
+    if typefit == 1: xx = sciopt.curve_fit(fitfunc, middle, yeah[0])
+    elif typefit == 2: xx = sciopt.curve_fit(fitfunc2, middle, yeah[0])
+    elif typefit == 3: xx = sciopt.curve_fit(fitfunc3, middle, yeah[0])
+    print "exponent",xx[0][1],"variance %",100.*xx[1][1][1]/xx[0][1]
+    
+    # label
+    lablab = r"$\alpha=$%4.1f"%(xx[0][1])
+    titi = addtitle+"LES "+str(int(dx))+"m. "+str(total)+" detected vortices"
+    
+    # plot obtained fit along with actual points
+    if typefit == 1: func = fitfunc(middle,xx[0][0],xx[0][1])
+    elif typefit == 2: func = fitfunc2(middle,xx[0][0],xx[0][1])
+    elif typefit == 3: func = fitfunc3(middle,xx[0][0],xx[0][1],xx[0][2]) 
+    func[func<0]=np.nan
+    ind = yeah[0]>0
+    mpl.plot(middle[ind],yeah[0][ind],'k.')
+    mpl.plot(middle[ind],func[ind],'r-',label=lablab)
+    mpl.plot(middle[ind],func[ind],'r.')
+    
+    # display bin limits
+    tata = '[bins:'
+    yorgl = ind[np.where(ind)].size + 1 # to add the last limit
+    for el in zebins[0:yorgl]:
+       tata = tata + "%0.f "%(el)
+    tata=tata+']'
+    
+    # print fit vs. actual populations
+    fit = func*binwidth*total
+    pop = yeah[0]*binwidth*total
+    for iii in range(len(fit)):
+        if fit[iii] > 0.99 or pop[iii] != 0:
+            yorgl = 100.*np.abs(fit[iii]-pop[iii])/fit[iii]
+            yargl = 100.*incertitude_nvortex/fit[iii]
+            print "fit %4.0f real %4.0f pc %3.0f ipc %3.0f" % (fit[iii],pop[iii],yorgl,yargl)
+    
+    # plot settings
+    ax = mpl.gca() ; ax.set_xbound(lower=np.min(zebins),upper=np.max(zebins))
+    if drop: 
+        mpl.xlabel("Pressure drop (Pa)")
+        ax.set_xbound(lower=1.e-1,upper=10.)
+    else: 
+        mpl.xlabel("Vortex size (m)")
+        #ax.set_xbound(lower=30.,upper=3000.)
+        ax.set_xbound(lower=30.,upper=1000.)
+        ax.set_ybound(lower=5.e-6,upper=5.e-2)
+    mpl.ylabel('Population density $n / N w_{bin}$')
+    #mpl.title(titi+'\n '+tata)
+    mpl.title(titi)
+    mpl.legend(loc="upper right",fancybox=True)
+        
+    # show plot and end
+    mpl.show()
 
-# choose a variable
-if drop: var = d
-else: var = s
+#################################################################################
+#################################################################################
+#################################################################################
+#
+### UNCOMMENT
+###import plfit,plpva,plplot
+#
+##[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)
+#[alpha, xmin, L] = plfit.plfit(var)
+#print alpha,xmin
+#
+##a = plpva.plpva(var,0.75,'xmin',0.75)
+##print a
+#
+#
+#h = plplot.plplot(var,xmin,alpha)
+#
+#ppplot.save(mode="png")
+#
+##mpl.loglog(h[0], h[1], 'k--',linewidth=2)
+#
+##mpl.show()
+#
+#
 
-# restrictions
-restrict = (s > 0) # initialization (True everywhere)
-restrict = restrict*(s >= limrest*dx) # remove lowest sizes (detection limit) 
-#restrict = restrict*(np.abs(i-j) <= 6.*dx) # condition sur i,j (width,height)
-#restrict = restrict*(d > 0.9) # limit on drop for size (casse la power law? seulement si keep smaller devils)
-#restrict = restrict*(d > 0.5)
-var = var[restrict]
-
-## define bins
-zebins = [np.min(var)]
-for i in range(0,nbins):  zebins.append(zebins[i]*(www**0.5))
-zebins = np.array(zebins)
-middle = 0.5*(zebins[1:] + zebins[:-1])
-binwidth = zebins[1:] - zebins[:-1]
-
-# plot histogram
-yeah = mpl.hist(var,log=True,bins=zebins,normed=True,color='white') ; mpl.xscale('log')
-#yeah = mpl.hist(var,bins=zebins)#,normed=True)
-
-# print info
-print "min %5.2e // max %5.2e" % (np.min(var),np.max(var))
-for iii in range(len(zebins)-1):
-    print "%5.2e in [%5.0f %5.0f] %5.0f" % (yeah[0][iii],zebins[iii],zebins[iii+1],np.abs(zebins[iii+1]-zebins[iii]))
-
-# fitting function to superimpose
-if typefit == 1: xx = sciopt.curve_fit(fitfunc, middle, yeah[0])
-elif typefit == 2: xx = sciopt.curve_fit(fitfunc2, middle, yeah[0])
-elif typefit == 3: xx = sciopt.curve_fit(fitfunc3, middle, yeah[0])
-print "exponent",xx[0][1]
-print "variance %",100.*xx[1][1][1]/xx[0][1]
-
-# plot obtained fit along with actual points
-if typefit == 1: func = fitfunc(middle,xx[0][0],xx[0][1])
-elif typefit == 2: func = fitfunc2(middle,xx[0][0],xx[0][1])
-elif typefit == 3: func = fitfunc3(middle,xx[0][0],xx[0][1],xx[0][2]) 
-func[func<0]=np.nan
-ind = yeah[0]>0
-mpl.plot(middle[ind],yeah[0][ind],'k.')
-mpl.plot(middle[ind],func[ind],'r-')
-mpl.plot(middle[ind],func[ind],'r.')
-
-# print fit vs. actual populations
-total = var.shape[0]
-fit = func*binwidth*total
-pop = yeah[0]*binwidth*total
-for iii in range(len(fit)):
-    if fit[iii] > 0.99 or pop[iii] != 0:
-        print "fit %5.0f actual %5.0f" % (fit[iii],pop[iii])
-
-# plot settings
-ax = mpl.gca() ; ax.set_xbound(lower=np.min(zebins),upper=np.max(zebins))
-if drop: mpl.xlabel("Pressure drop (Pa)") ; ax.set_xbound(lower=1.e-1,upper=10.)
-else: mpl.xlabel("Vortex size (m)") ; ax.set_xbound(lower=10.,upper=5000.)
-mpl.ylabel('Population density $n / N w_{bin}$')
-mpl.title("Statistics on "+str(total)+" detected vortices")
-
-# show plot and end
-mpl.show()
-exit()
-
-################################################################################
-################################################################################
-################################################################################
-
-## UNCOMMENT
-##import plfit,plpva,plplot
-
-#[alpha, xmin, L] = plfit.plfit(var,'xmin',0.3)
-[alpha, xmin, L] = plfit.plfit(var)
-print alpha,xmin
-
-#a = plpva.plpva(var,0.75,'xmin',0.75)
-#print a
-
-
-h = plplot.plplot(var,xmin,alpha)
-
-ppplot.save(mode="png")
-
-#mpl.loglog(h[0], h[1], 'k--',linewidth=2)
-
-#mpl.show()
-
-
-
Index: /trunk/UTIL/PYTHON/powerlaw/ddstat.py
===================================================================
--- /trunk/UTIL/PYTHON/powerlaw/ddstat.py	(revision 1196)
+++ /trunk/UTIL/PYTHON/powerlaw/ddstat.py	(revision 1197)
@@ -3,70 +3,43 @@
 import numpy as np
 
-##########################################################
-namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm1_2.txt"
-namefile = "/home/aymeric/Big_Data/LES_dd/psfc_f18.ncm2_2.txt"
-##########################################################
-namefile = "/home/aymeric/Big_Data/LES_dd/psfc.LMD_LES_MARS.160564.ncm1_2.txt"
-#namefile = "/home/aymeric/Big_Data/LES_dd/sav/psfc.LMD_LES_MARS.160564.ncm2_2.txt"
-##########################################################
-#namefile = "/home/aymeric/Big_Data/LES_dd/psfc_oldinsight100m.ncm1_2.txt"
-##########################################################
+### STATDD
+def statdd(namefile):
 
-
- 
-
-#case = "188324p"
-#case = "191798"
-#case = "160564p"
-#case = "156487"
-#case = "2007p"
-#case = "13526p"
-#case = "172097"
-
-
-namefile = "/planeto/aslmd/LESdata/"+case+".ncm1_2.txt"
-
-
-
-
-
-# load data
-data = np.loadtxt(namefile,delimiter=";")
-t = data[:,0] ; n = data[:,1] ; s = data[:,2] ; d = data[:,3]
-
-# remove size and drop point when no vortex detected
-d[np.where(n==0)] = np.nan
-s[np.where(n==0)] = np.nan
-
-## PLOTS
-
-number = ppplot.plot1d()
-number.f = n
-number.x = t
-number.linestyle = ''
-number.marker = '.'
-number.color = 'b'
-number.xlabel = "Local time (hour)"
-number.ylabel = "Detected vortices"
-number.makeshow()
-
-drop = ppplot.plot1d()
-drop.f = d
-drop.x = t
-drop.linestyle = ''
-drop.marker = '.'
-drop.color = 'r'
-drop.fmt = "%.1f"
-drop.xlabel = "Local time (hour)"
-drop.ylabel = "Maximum drop of detected vortices (Pa)"
-drop.makeshow()
-
-size = ppplot.plot1d()
-size.f = s
-size.x = t
-size.linestyle = ''
-size.marker = '.'
-size.color = 'g'
-size.xlabel = "Local time (hour)"
-size.ylabel = "Maximum size of detected vortices (m)"
-size.makeshow()
+    # load data
+    data = np.loadtxt(namefile,delimiter=";")
+    t = data[:,0] ; n = data[:,1] ; s = data[:,2] ; d = data[:,3]
+    
+    # remove size and drop point when no vortex detected
+    d[np.where(n==0)] = np.nan ; s[np.where(n==0)] = np.nan
+    
+    ## PLOTS
+    number = ppplot.plot1d()
+    number.f = n
+    number.x = t
+    number.linestyle = ''
+    number.marker = '.'
+    number.color = 'b'
+    number.xlabel = "Local time (hour)"
+    number.ylabel = "Detected vortices"
+    number.makeshow()
+    
+    drop = ppplot.plot1d()
+    drop.f = d
+    drop.x = t
+    drop.linestyle = ''
+    drop.marker = '.'
+    drop.color = 'r'
+    drop.fmt = "%.1f"
+    drop.xlabel = "Local time (hour)"
+    drop.ylabel = "Maximum drop of detected vortices (Pa)"
+    drop.makeshow()
+    
+    size = ppplot.plot1d()
+    size.f = s
+    size.x = t
+    size.linestyle = ''
+    size.marker = '.'
+    size.color = 'g'
+    size.xlabel = "Local time (hour)"
+    size.ylabel = "Maximum size of detected vortices (m)"
+    size.makeshow()
Index: /trunk/UTIL/PYTHON/powerlaw/finddd.py
===================================================================
--- /trunk/UTIL/PYTHON/powerlaw/finddd.py	(revision 1196)
+++ /trunk/UTIL/PYTHON/powerlaw/finddd.py	(revision 1197)
@@ -6,4 +6,9 @@
 import matplotlib.pyplot as mpl
 import ppcompute
+
+## method 3
+#from skimage import filter,transform,feature
+#import matplotlib.patches as mpatches 
+###########
 
 ###############################################################################
@@ -71,4 +76,6 @@
     mean = np.mean(psfc,dtype=np.float64)
     std = np.std(psfc,dtype=np.float64)
+    damax = np.max(psfc)
+    damin = np.min(psfc)
     ## some information about inferred limits
     print "**************************************************************"
@@ -94,12 +101,12 @@
 
      ## MAIN ANALYSIS. LOOP ON FAC. OR METHOD.
-     #for fac in faclist:
-     fac = 3.75
-     for method in [2,1]:
+     for fac in faclist:
+     #fac = 3.75
+     #for method in [2,1]:
   
       ## initialize arrays
       tabij = [] ; tabsize = [] ; tabdrop = []
       tabijcenter = [] ; tabijvortex = [] ; tabijnotconv = [] ; tabdim = []
-    
+
       ################ FIND RELEVANT POINTS TO BE ANALYZED
       ## lab is 1 for points to be treated by minimum_position routine
@@ -110,5 +117,5 @@
           # method 1: standard deviation
           lab[np.where(psfc2d < mean-fac*std)] = 1
-      else:
+      elif method == 2:
           # method 2: polynomial fit
           # ... tried smooth but too difficult, not accurate and too expensive
@@ -131,11 +138,66 @@
           limlim = fac*std ## same as method 1
           lab[np.where(anopsfc2d < -limlim)] = 1
-    
+      elif method == 3:
+          # method 3 : find centers of circle features using a Hough transform
+          # initialize the array containing point to be further analyzed
+          lab = np.zeros(psfc2d.shape)
+          ## prepare the field to be analyzed by the Hough transform
+          ## normalize it in an interval [-1,1]
+          field = psfc2d
+          # ... test 1. local max / min used for normalization.
+          mmax = np.max(field) ; mmin = np.min(field) ; dasigma = 2.5
+          ## ... test 2. global max / min used for normalization.
+          #mmax = damax ; mmin = damin ; dasigma = 1.0 #1.5 trop restrictif
+          spec = 2.*((field-mmin)/(mmax-mmin) - 0.5)
+          # perform an edge detection on the field
+          # ... returns an array with True on edges and False outside
+          # http://sciunto.wordpress.com/2013/03/01/detection-de-cercles-par-une-transformation-de-hough-dans-scikit-image/    
+          edges = filter.canny(filter.sobel(spec),sigma=dasigma)
+          # initialize plot for checks
+          if plotplot:
+            fig, ax = mpl.subplots(ncols=1, nrows=1, figsize=(10,10))
+            ax.imshow(field, cmap=mpl.cm.gray)
+          # detect circle with radius 3dx. works well.
+          # use an Hough circle transform
+          radii = np.array([3]) 
+          hough_res = transform.hough_circle(edges, radii)
+          # analyze results of the Hough transform
+          nnn = 0 
+          sigselec = neighbor_fac
+          #sigselec = 3.
+          for radius, h in zip(radii, hough_res):
+            # number of circle features to keep
+            # ... quite large. but we want to be sure not to miss anything.
+            nup = 30 
+            maxima = feature.peak_local_max(h, num_peaks=nup)
+            # loop on detected circle features
+            for maximum in maxima:
+              center_x, center_y = maximum - radii.max()
+              # nup is quite high so there are false positives.
+              # ... but those are easy to detect
+              # ... if pressure drop is unclear (or inexistent)
+              # ... we do not take the point into account for further analysis
+              # ... NB: for inspection give red vs. green color to displayed circles
+              diag = field[center_x,center_y] - (mean-sigselec*std)
+              if diag < 0:  
+                  col = 'green'
+                  nnn = nnn + 1
+                  lab[center_x,center_y] = 1
+              else:
+                  col = 'red'
+              # draw circles
+              if plotplot:
+                circ = mpatches.Circle((center_y, center_x), radius,fill=False, edgecolor=col, linewidth=2)
+                ax.add_patch(circ)
+          if plotplot:
+            mpl.title(str(nnn)+" vortices") ; mpl.show()
+            mpl.close()
+
       ## while there are still points to be analyzed...
       while 1 in lab:
         ## ... get the point with the minimum field values
-        if method == 1:
+        if method == 1 or method == 3:
             ij = minimum_position(psfc2d,labels=lab)
-        else:
+        elif method == 2:
             ij = minimum_position(anopsfc2d,labels=lab)
         ## ... store the indexes of the point in tabij
@@ -153,7 +215,7 @@
       ## --> or even lower as shown by plotting reslab 
       reslab = np.zeros(psfc2d.shape)
-      if method == 1:
+      if method == 1 or method == 3:
           reslab[np.where(psfc2d < mean-neighbor_fac*std)] = 1
-      else:
+      elif method == 2:
           reslab[np.where(anopsfc2d < -neighbor_fac*std)] = 1
      
@@ -227,5 +289,5 @@
             ## OK. this is most likely an actual vortex. we get the drop.
             ## we multiply by mesh area, then square to get approx. size of vortex
-            if method == 1:
+            if method == 1 or method ==3:
                 drop = -psfc2d[i,j]+mean
             else:
@@ -301,9 +363,9 @@
     
       #### PLOT PLOT PLOT PLOT
-      if nvortex>0 and plotplot:
+      if (nvortex>0 and plotplot) or (nvortex>0 and maxsize > 800.):
        mpl.figure(figsize=(12,8))
        myplot = plot2d()
-       myplot.absc = np.array(range(psfc2d.shape[1]))*dx/1000.
-       myplot.ordi = np.array(range(psfc2d.shape[0]))*dx/1000.
+       myplot.x = np.array(range(psfc2d.shape[1]))*dx/1000.
+       myplot.y = np.array(range(psfc2d.shape[0]))*dx/1000.
        myplot.title = str(nvortex)+" vortices found (indicated diameter / pressure drop)"
        myplot.xlabel = "x distance (km)"
@@ -312,5 +374,5 @@
        #if method == 1:
            #myplot.field = ustm 
-           myplot.field = psfc2d
+           myplot.f = psfc2d
            #myplot.field = polypsfc2d
            myplot.vmin = -2.*std 
@@ -324,5 +386,5 @@
        myplot.fmt = "%.1f"
        myplot.div = 20
-       myplot.colorb = "spectral"
+       myplot.colorbar = "spectral"
        myplot.make()
       
Index: /trunk/UTIL/PYTHON/powerlaw/performdd.py
===================================================================
--- /trunk/UTIL/PYTHON/powerlaw/performdd.py	(revision 1196)
+++ /trunk/UTIL/PYTHON/powerlaw/performdd.py	(revision 1197)
@@ -2,6 +2,12 @@
 from finddd import finddd
 
-finddd("/planeto/aslmd/LESdata/case_A_50m_145_145_201_12km.nc",dt_out=100.,lt_start=8.)
 
+finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=3)
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",timelist=[250,385],method=3)
+
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=2)
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",method=1)
+
+#finddd("/home/aymeric/Big_Data/LES_dd/psfc_f18.nc",timelist=range(0,591,20),method=1,save=False,plotplot=True)
 exit()