; Copyright (c) 1984-2009, ITT Visual Information Solutions. All
    ;       rights reserved. Unauthorized reproduction is prohibited.
    ;
    
    ;+
    ; NAME:
    ;    DERIV
    ;
    ; PURPOSE:
    ;    Perform numerical differentiation using 3-point, Lagrangian 
    ;    interpolation.
    ;
    ; CATEGORY:
    ;    Numerical analysis.
    ;
    ; CALLING SEQUENCE:
    ;    Dy = Deriv(Y)         ;Dy(i)/di, point spacing = 1.
    ;    Dy = Deriv(X, Y)    ;Dy/Dx, unequal point spacing.
    ;
    ; INPUTS:
    ;    Y:  Variable to be differentiated.
    ;    X:  Variable to differentiate with respect to.  If omitted, unit 
    ;        spacing for Y (i.e., X(i) = i) is assumed.
    ;
    ; OPTIONAL INPUT PARAMETERS:
    ;    As above.
    ;
    ; OUTPUTS:
    ;    Returns the derivative.
    ;
    ; COMMON BLOCKS:
    ;    None.
    ;
    ; SIDE EFFECTS:
    ;    None.
    ;
    ; RESTRICTIONS:
    ;    None.
    ;
    ; PROCEDURE:
    ;    See Hildebrand, Introduction to Numerical Analysis, Mc Graw
    ;    Hill, 1956.  Page 82.
    ;
    ; MODIFICATION HISTORY:
    ;    Written, DMS, Aug, 1984.
    ;    Corrected formula for points with unequal spacing.  DMS, Nov, 1999.
    ;-
    ;
    ; on_error,2              ;Return to caller if an error occurs

def deriv(*args):
    """
    ; Copyright (c) 1984-2009, ITT Visual Information Solutions. All
    ;       rights reserved. Unauthorized reproduction is prohibited.
    ;
    
    ;+
    ; NAME:
    ;    DERIV
    ;
    ; PURPOSE:
    ;    Perform numerical differentiation using 3-point, Lagrangian 
    ;    interpolation.
    ;
    ; CATEGORY:
    ;    Numerical analysis.
    ;
    ; CALLING SEQUENCE:
    ;    Dy = Deriv(Y)         ;Dy(i)/di, point spacing = 1.
    ;    Dy = Deriv(X, Y)    ;Dy/Dx, unequal point spacing.
    ;
    ; INPUTS:
    ;    Y:  Variable to be differentiated.
    ;    X:  Variable to differentiate with respect to.  If omitted, unit 
    ;        spacing for Y (i.e., X(i) = i) is assumed.
    ;
    ; OPTIONAL INPUT PARAMETERS:
    ;    As above.
    ;
    ; OUTPUTS:
    ;    Returns the derivative.
    ;
    ; COMMON BLOCKS:
    ;    None.
    ;
    ; SIDE EFFECTS:
    ;    None.
    ;
    ; RESTRICTIONS:
    ;    None.
    ;
    ; PROCEDURE:
    ;    See Hildebrand, Introduction to Numerical Analysis, Mc Graw
    ;    Hill, 1956.  Page 82.
    ;
    ; MODIFICATION HISTORY:
    ;    Written, DMS, Aug, 1984.
    ;    Corrected formula for points with unequal spacing.  DMS, Nov, 1999.
    ;-
    ;
    ; on_error,2              ;Return to caller if an error occurs
    """
    x = args[0]
    n = x.size
    if n < 3 : raise Exception('Parameters must have at least 3 points')


    if (len(args) == 2) :
        y=args[1]
        if n != y.size : raise 'Vectors must have same size'
    
        #;df/dx = y0*(2x-x1-x2)/(x01*x02)+y1*(2x-x0-x2)/(x10*x12)+y2*(2x-x0-x1)/(x20*x21)
        #; Where: x01 = x0-x1, x02 = x0-x2, x12 = x1-x2, etc.
    
        if isinstance(x,np.ma.masked_array) :  x = x.data                # Convert masked arrays to classic arrays
        if not isinstance(x,np.float) : x.astype(np.float) #;If not floating type, ensure floating... 
        
        x12 = x - np.roll(x,-1)                                      #;x1 - x2
        x01 = np.roll(x,1) - x                                       #;x0 - x1
        x02 = np.roll(x,1) - np.roll(x,-1)                           #;x0 - x2
    
        d = np.roll(y,1) * (x12 / (x01*x02)) \
            + y * (1./x12 - 1./x01) \
            - np.roll(y,-1) * (x01 / (x02 * x12))                    #Middle points
        
        
        #Formulae for the first and last points:
        d[0] = y[0] * (x01[1]+x02[1])/(x01[1]*x02[1]) \
            - y[1] * x02[1]/(x01[1]*x12[1]) \
            + y[2] * x01[1]/(x02[1]*x12[1])                          #;First point
        n2 = n-2
        d[n-1] = -y[n-3] * x12[n2]/(x01[n2]*x02[n2]) \
            + y[n-2] * x02[n2]/(x01[n2]*x12[n2]) \
            - y[n-1] * (x02[n2]+x12[n2]) / (x02[n2]*x12[n2])             #;Last point

    #Equally spaced point case
    else :
        d = (np.roll(x,-1) - np.roll(x,1))/2.
        d[0] = (-3.0*x[0] + 4.0*x[1] - x[2])/2.
        d[n-1] = (3.*x[n-1] - 4.*x[n-2] + x[n-3])/2.

    return d