;+
    ;
    ;   POWEL_LEBEN_FILTER_KM : Compute geostrophic speeds from a sea level dataset <br />
    ;   using a Powell et Leben (2004) along-track filtering to maintain measurement<br />
    ;    noise constant.<br /><br />
    ;
    ;   Reference : Powell, B. S., et R. R. Leben (2004), An Optimal Filter for <br />
    ;   Geostrophic Mesoscale Currents from Along-Track Satellite Altimetry, <br />
    ;   Journal of Atmospheric and Oceanic Technology, 21(10), 1633-1642.
    ;
    ; @param lon {in}{optional}{type:NUMERIC} longitude in degrees
    ; @param lat {in}{optional}{type:NUMERIC} latitude in degrees
    ; @param z {in}{required}{type:NUMERIC} sea level surface. Can either be absolute<br />
    ;          values (SSH) or relative values (SLA). This MUST be given in METERS.
    ; @param p {in}{optional}{type:NUMERIC}{default:12} Filter half-width BEFORE<br />
    ;          center of filtering window in KILOMETERS.
    ; @param q {in}{optional}{type:NUMERIC}{default:12} Filter half-width AFTER <br />
    ;          center of filtering window in KILOMETERS.
    ;          
    ; @returns Geostrophic velocity component, positive eastward
    ;
    ;
    ;
    ; @author Renaud DUSSURGET, LEGOS/CTOH
    ; @history Created Sep. 2009 from genweights.m (Brian Powell (c) 2004, <br />
    ;   University of Colorado, Boulder)<br />
    ;   Modified May 2010 to be compliant with 20Hz datasets (p & n can vary).<br />
    ;     Warining may also be issued for data with holes within the width of the<br />
    ;     window.<br />
    ;   Modified June 2010 to include filtering window width in KM instead of nb. of<br />
    ;     points (Equivalent btw. 1Hz and 20Hz data).<br />
    ;
    ; @uses CALCUL_DISTANCE, EXIST, GENWEIGTHS, SETINTERSECTION, SETUNION, <br />
    ;   OPTIMAL_SLOPE, GRAVITY, CORIOLIS, TRACK_ORIENT
    ;
    ; @example dummy1=powell_leben_filter_KM(lon,lat,sla,p=11,q=11) :<br />
    ;   Return along-track velocity anomalies using a 11km by 11km filter window.
    ;
    ;-

def powell_leben_filter_km(*args,**kwargs):
    """
    ;+
    ;
    ;   POWEL_LEBEN_FILTER_KM : Compute geostrophic speeds from a sea level dataset <br />
    ;   using a Powell et Leben (2004) along-track filtering to maintain measurement<br />
    ;    noise constant.<br /><br />
    ;
    ;   Reference : Powell, B. S., et R. R. Leben (2004), An Optimal Filter for <br />
    ;   Geostrophic Mesoscale Currents from Along-Track Satellite Altimetry, <br />
    ;   Journal of Atmospheric and Oceanic Technology, 21(10), 1633-1642.
    ;
    ; @param lon {in}{optional}{type:NUMERIC} longitude in degrees
    ; @param lat {in}{optional}{type:NUMERIC} latitude in degrees
    ; @param z {in}{required}{type:NUMERIC} sea level surface. Can either be absolute<br />
    ;          values (SSH) or relative values (SLA). This MUST be given in METERS.
    ; @param p {in}{optional}{type:NUMERIC}{default:12} Filter half-width BEFORE<br />
    ;          center of filtering window in KILOMETERS.
    ; @param q {in}{optional}{type:NUMERIC}{default:12} Filter half-width AFTER <br />
    ;          center of filtering window in KILOMETERS.
    ;          
    ; @returns Geostrophic velocity component, positive eastward
    ;
    ;
    ;
    ; @author Renaud DUSSURGET, LEGOS/CTOH
    ; @history Created Sep. 2009 from genweights.m (Brian Powell (c) 2004, <br />
    ;   University of Colorado, Boulder)<br />
    ;   Modified May 2010 to be compliant with 20Hz datasets (p & n can vary).<br />
    ;     Warining may also be issued for data with holes within the width of the<br />
    ;     window.<br />
    ;   Modified June 2010 to include filtering window width in KM instead of nb. of<br />
    ;     points (Equivalent btw. 1Hz and 20Hz data).<br />
    ;
    ; @uses CALCUL_DISTANCE, EXIST, GENWEIGTHS, SETINTERSECTION, SETUNION, <br />
    ;   OPTIMAL_SLOPE, GRAVITY, CORIOLIS, TRACK_ORIENT
    ;
    ; @example dummy1=powell_leben_filter_KM(lon,lat,sla,p=11,q=11) :<br />
    ;   Return along-track velocity anomalies using a 11km by 11km filter window.
    ;
    ;-
    """
    
    
    #From Powell and Leben, 2004. JAOT.
    #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    lon = args[0]
    lat = args[1]
    nu = args[2]
    dst = calcul_distance(lat,lon) * 1e3 #distance in meters
    
    n=nu.size

    #We use a p+q filter (default 2 1Hz points on each side ~ 12km)
    #-> Set p & q distances to defaults if they don't exist    
    p = kwargs.pop('p',12.)
    q = kwargs.pop('q',12.)
    verbose = kwargs.pop('verbose',False)
        
    dt = np.median(dst[1:] - dst[:-1])   #Median point-to-point distance in meters
      
    #Get the corresponding number of points rounded to the nearest integer
    pn=np.round((p*1e3)/dt).astype(int)
    qn=np.round((q*1e3)/dt).astype(int)
    
    if pn + qn > n : raise 'Filtering window is too large wrt array length'
     
    #Compute weights
    #;;;;;;;;;;;;;;;
    
    #Initialise all possible configurations of filters
    #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    cn = []#=np.array([])   #array of weigth arrays
    nout= [] #np.array([]) #array of point indices in filter arrays
    pout=np.ndarray(pn+qn+1) #array containg number of points before any filter
    qout=np.ndarray(pn+qn+1) #array containg number of points after any filter
      
      
    #Compute each filter configuration
    #1)Start with increasing number of points before point
    #2)Compute full resolution filter (p & q points before and after filter)
    #3)Finish with a deacreasing number of points after filter
    
    #1)Weights before
    for i in np.arange(pn):
        w,n=genweights(i,qn,dt,n=True)
        cn.append(w)
        nout.append(n)
        pout[i]=i
        qout[i]=qn
        
    #2)Centered weights
    w,n=genweights(pn,qn,dt,n=True)
    cn.append(w)
    nout.append(n)
    pout[pn]=pn
    qout[pn]=qn
        
    #3)Weights after
    for i in np.arange(pn+1, pn+1+qn):
        w,n=genweights(pn,pn+qn-i,dt,n=True)
        cn.append(w)
        nout.append(n)
        pout[i]=pn
        qout[i]=pn+qn-i
        
    #find track edges and gaps
    #;;;;;;;;;;;;;;;;;;;;;;;;;
    empty = np.where(np.isnan(nu))[0]
    ok = np.where(~np.isnan(nu))[0]
    emptycnt = empty.size
    okcnt = ok.size
    
    st = np.min(ok) #start point
    en = np.max(ok) #end of track
        
    #Initiate dh for slope calculation, Eq. (5)
    dh = np.repeat(np.nan,nu.size)
      
    #Loop over valid points in track
    #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    for i in np.arange(okcnt):
        
        #find valid points within filter width around filter center
        #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        cur = np.arange(pn+qn+1)+ok[i]-pn               #Points around filter center +-p/q points
        cur = np.sort(list(set(cur).intersection(set(ok))))  #intersection of filter window with valid points
        curcnt = cur.size
        
        #check points currently selected
        #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        if (curcnt > 1) :
            #Be careful -> This is intended to select points following directly central point regularly (without gaps within window)
            #We'll see that for a hi-rate case, we'll be able to accept gaps within filter window
            #FIX:this might be buggy...
            a = np.where(np.append(0,cur[1:] - cur[:-1]) == 1)[0]  #select following valid points
            b = np.where(np.append(cur[:-1] - cur[1:],0) ==  -1)[0] #select preceeding valid points
            aftcnt = a.size
            befcnt = b.size
            
#          b = WHERE([cur[:-1] - cur[1:],0] EQ -1,befcnt) ; select preceeding valid points
            #If there are points at least before OR after central point, keep going, otherwise switch to next point
            if (aftcnt != 0) & (befcnt != 0):
                cur = cur[np.sort(list(set(a).union(b)))] #Get the union of valid points before and after central point
                curcnt = cur.size
            else : curcnt =0
        else : curcnt =0
        
        #Adapt filter length to current case
        #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
        #(We know how many points are following directly before and after)
        
        if (curcnt != 0) :
            
            #Get number of points before and after central point
            #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
            nbef = ok[i] - cur[0]
            naft = cur[curcnt - 1] - ok[i]
          
            #If we have less than p or q on each side, we'll set full filter width to
            #the greater size, and reduce the smaller size
            #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
            if (naft > nbef) & (naft+nbef > qn) : naft=pn
            if (naft < nbef) & (naft+nbef > pn) : nbef=qn
          
            #Select adapted filter in filter list
            #(given by pout and qout)
            #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
            which=(np.where((nbef == pout) & (naft == qout)))[0]
            nwhich = which.size
          
            #If we have found a filter, we use it, otherwise we issue a warning message and flag the point
            if (nwhich > 0) : dh[ok[i]] = optimal_slope(cn[which],nout[which],dt,nu,ok[i])
            else :
                if (verbose) : print "[WARNING] No available filter for point no.{0}".format(ok[i])
                dh[ok[i]] = np.NaN #If no suitable filter, flag the data
       
    #Compute currents with given slope
    #;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    #Compute gravitational & coriolis forces
    g = gravity(lat)
    f = coriolis(lat)

    #Get geostrophic velocities
    try :
        ug = - (g*dh) / (f)
    except ValueError :
        print 'error'
  
    #Inverse sign of ug for descending tracks as Coriolis is oriented to the right
    #northward
    if (not track_orient(lon,lat)) : #descending tracks
        ug *=-1
    
    return ug