;+
    ;   GENWEIGTHS : return optimal weigthing coefficients from Powell and Leben 2004<br />
    ;   translated from Matlab genweigths.m program to IDL<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.
    ;   
    ; @author Renaud DUSSURGET, LEGOS/CTOH
    ; @history Created Sep. 2009 from genweights.m (Brian Powell (c) 2004, <br />
    ;   University of Colorado, Boulder)<br />
    ;-

def genweights(p,q,dt,error=None , n=False):
    """
    ;+
    ;   GENWEIGTHS : return optimal weigthing coefficients from Powell and Leben 2004<br />
    ;   translated from Matlab genweigths.m program to IDL<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.
    ;   
    ; @author Renaud DUSSURGET, LEGOS/CTOH
    ; @history Created Sep. 2009 from genweights.m (Brian Powell (c) 2004, <br />
    ;   University of Colorado, Boulder)<br />
    ;-
    """
    
    p = np.abs(p)
    q = np.abs(q)
  
    #check inputs
    if (-p > q) : raise "genweights : P must be lesser than q"
  
    #Build matrices
    N = p + q
    T = N + 1
    A = np.matrix(np.zeros((T,T)))
    A[T-1,:] = np.append(np.repeat(1.0,N),0)
    sn = np.arange(T) - p
    sn = sn.compress(sn != 0)
    for i in np.arange(len(sn)) :
        A[i,:] = np.append(((1./sn)*(-sn[i]/2.)),sn[i]**2.*dt**2./4.) #Eq.11 (PL) 
        A[i,i] = -1.
        
    B = np.zeros(T)
    B[N] = 1.0
  
    #Compute the coefficients
    cn=np.dot(A.I,B)##B
#    cn=cn.transpose()
    cn=np.array([i for i in cn.flat])
    cn = cn[0:N] #Check the indices
  
    #Compute the error
    error = np.sqrt(np.sum(cn.transpose()/(sn*dt))**2. + np.sum( (cn.transpose()/(sn*dt))**2. ) );
  
    return cn, sn if n else cn