;+
; NAME:
;  cubic
; PURPOSE: (one line)
;  evaluate a cubic given values and derivatives at two points
; DESCRIPTION:
;  evaluate a cubic given values and derivatives at two points
; CATEGORY:
;  Math
; CALLING SEQUENCE:
;  y = cubic(x0, x1, y0, y1, dy0, dy1, x, dy, dx0, dx1, a)
; INPUTS:
;  x0, x1 - the fixed points
;  y0, y1 - values at the fixed points
;  dy0, dy1 - derivatives at the fixed points, same
;             units as (y1-y0)/(x1-x0)
;  x - new point at which to evaluate the cubic
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  none
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  returns y = value of the cubic at x
; OPTIONAL OUTPUT PARAMETERS:
;  dy = derivative of the cubic at x
;  dx0 = (x-x0)/(x1-x0) and dx1 = (x-x1)/(x1-x0) - for debug only
;  a = coefficients for y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
;     (for debug only)
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
;  Given a cubic function with the following endpoints 
;  y(x0) = y0
;  y(x1) = y1
;  dy(x0)/dx = dy0
;  dy(x1)/dx = dy0
;  calculate the coefficients so that
;  y = sum(ai dx0^(3-i) dx1^(i), i = 0..3)
;  dx0 = (x-x0)/(x1-x0)
;  dx1 = (x-x1)/(x1-x0)
; MODIFICATION HISTORY:
;  Written 2000 October, by Leslie Young, SwRI
;  Nov 2000.  Allow non-scalar argument, LAY
;  11 Mar 2006.  LAY.  Moved to $idl/layoung/math, check nparams
;    before calculating dy
;     
;-

function cubic, x0, x1, y0, y1, dy0, dy1, x, dy, dx0, dx1, a

    del = x1-x0
    dx0 = (x-x0)/del
    dx1 = (x-x1)/del
    a = fltarr(4)
    a[0] = y1
    a[3] = -y0
    a[1] = dy1 * del - 3 * a[0]
    a[2] = dy0 * del - 3 * a[3]
    
    y = 0
    for i = 0,3 do begin
        y = y + a[i]*dx0^(3-i)*dx1^(i)
    end
    if n_params() gt 7 then begin
        dy = (3*a[0]*dx0^2 + a[1]*dx0*(2*dx1+dx0) + $
              a[2]*dx1*(2*dx0+dx1) + 3*a[3]*dx1^2) / del
    end
    return, y
end

pro testcubic
; y = (x+1)^3   dy/dx = 3 (x+1)^2  
x0 = 0. & x1 = 2.
y0 = (x0+1.)^3. & y1 = (x1+1.)^3.
dy0 = 3.*(x0+1)^2. & dy1 = 3.*(x1+1.)^2

x = 0.5

y = cubic (x0, x1, y0, y1, dy0, dy1, x, dy, dx0, dx1, a)
print, [y, (x+1.)^3.]
print, [dy, 3*(x+1)^2.]

end