function cubic_min,x,y,max=max,yfit=yfit

;+
;
; CUBIC_MIN
; A function to find the minimum of a curve Y(X) that is best fit by a
; cubic polynomial.
;
; Syntax:
; Result = CUBIC_MIN(X, Y [, /MAX] [, YFIT=Variable])
;
; Input:
; X - The independent variable / abscissa.
; Y - The dependent variable / ordinate.
;
; Returns:
; Fits coefficients for y = ax^3 + bx^2 + cx + d, then returns the
; root of dy/dx = 0 = 3ax^2 + 2bx + c that corresponds to the local
; minimum in y.  If no minimum exists, returns a value of 0.0 with a
; warning error message.  If minimum exists but is outside the domain
; specified by x, returns that minimum with a warning error message.
;
; Keyword:
; MAX - Set /MAX to return maximum instead of minimum.
; YFIT - Set equal to variable name to have the best-fit y-values
;        passed back to user.
;
; Katie Peek / Sept 2007
; Dec 2007 / Added YFIT keyword and syntax printing.
; Oct 2008 / Added STATUS keyword in poly_fit to suppress printing to screen.
;
;-

; Print syntax if no variables are provided.
if (n_elements(x) eq 0) then begin
  print,' Cubic_Min syntax: '
  print,'   Minimum = Cubic_Min(x, y [, /max] [, yfit=Variable]) '
  goto,theend
end

; Do polynomial fit.
aa = poly_fit(x,y,3,/double,yfit=yfit,status=status)
a = aa[3]
b = aa[2]
c = aa[1]
d = aa[0]

; Find the roots of the resulting equation.
; dy/dx = 3ax^2 + 2bx + c = 0
radicand = 4*b^2 - 12*a*c
if (radicand lt 0.0) then begin
  message,'CUBIC_MIN: No minimum exists for cubic function.  Returning 0.0.' $
         ,/continue
  goodroot=0.0
  goto,theend
endif
roots = [(-2*b + sqrt(radicand))/(6*a), (-2*b - sqrt(radicand))/(6*a)]

; Which corresponds to a minimum?
; d2y/dx2 = 6ax + 2b <- plug in roots.
concavity = 6*a*roots + 2*b
if keyword_set(max) then begin
  goodroot = roots[where(concavity lt 0.0)]
endif else begin
  goodroot = roots[where(concavity gt 0.0)]
endelse
if (n_elements(goodroot) eq 2) then $
  message,'CUBIC_MIN: Cubic functions seem to be broken.'

; Is the good root within the domain of x?
if ((goodroot lt min(x)) or (goodroot gt max(x))) then $
  message,'CUBIC_MIN: Warning--cubic minimum outside specified domain.' $
         ,/continue

theend: ; Escape to here if something breaks.
return,goodroot[0]
end