FUNCTION LocalExtrema, Array, minima=minima, maxima=maxima, 	$
	floatx=floatx, fnc=fnc, count=count, strict=strict

;+
; NAME:
;	LocalExtrema
; PURPOSE:
;	Return location of local maxima and minima in an array
; CATEGORY:
;	Math
; CALLING SEQUENCE:
;	R = LocalExtrema(Array, [/minima, /float])
; INPUTS:
;	Array		array[*]; type: 1-dim array
;					input array to be searched for extrema
; OPTIONAL INPUT PARAMETERS:
;	By default the locations of an alternating sequence of
;	maxima and minima are returned.
;
;	/maxima		if /maxima is set then the locations of maxima are returned
;	/minima		if /minima is set then the locations of minima are returned
;
;	/float		by default integer indices are returned for the location of the extrema.
;				If /float is set then the location is refined by fitting a quadratic to
;				the minimum element and its two neighbours. The extremum of the
;				quadratic is returned as a floating point number.
; OUTPUTS:
;	R			scalar or array; type: integer or float
;					indices of minima or maxima
; OPTIONAL OUTPUT PARAMETERS:
;	fnc=fnc		array[*]; same size as R, same type as 'Array'
;					function values in locations R. If the /float keyword is set it
;					returns the function values in the extremum of the quadratic fit.
;	count=count scalar; type: integer
;					# extrema
; INCLUDE:
	@compile_opt.pro			; On error, return to caller
; CALLS:
;	InitVar, BadValue, LocalExtrema
; RESTRICTIONS:
; >	Only genuine maxima and minima are detected,
; >	The edges of the array are excluded, i.e. only the locations of
;	'internal' extrema (with neighbour on both sides) are returned
; PROCEDURE:
;	Let 'array' represent a function f(x) in points x=[0,1,.]
;	Let f have an extremum in grid point j, i.e. f(j) > f(j-1) and f(j) > f(j+1)
;	Define dx = x-j and g(dx) = f(j+dx)-f(j), i.e. g(0) = 0 is an extremum of g
;	The quadratic for g through 0 and its two neighbours -1 and +1 is:
;		y(dx) = 0.5*dx*{ g(-1)*(dx-1)+g(+1)*(dx+1) }
;	which has a maximum at
;		dx = 0.5*(g(-1)-g(1))/(g(-1)+g(1))
; MODIFICATION HISTORY:
;	MAR-2000, Paul Hick (UCSD/CASS)
;	FEB-2004, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Introduced /maxima keyword to return maxima. This used to be the default.
;		The default now is to return an alternating sequence of minima and maxima.
;-

InitVar, minima, /key
InitVar, maxima, /key
InitVar, floatx, /key
InitVar, strict, /key

x   = -1
fnc = BadValue(Array)
count = 0

IF minima+maxima EQ 0 THEN BEGIN

	x1 = LocalExtrema(Array, floatx=floatx, fnc=fnc1, count=count1, /maxima)
	x2 = LocalExtrema(Array, floatx=floatx, fnc=fnc2, count=count2, /minima)

	count = count1+count2

	IF count NE 0 THEN BEGIN

		CASE 0 OF
		count1: BEGIN
			x   = x2
			fnc = fnc2
		END
		count2: BEGIN
			x   = x1
			fnc = fnc1
		END
		ELSE: BEGIN
			CASE 1 OF
			count1 GT count2: BEGIN
				x   = [x1  [0],reform(transpose([[x2  ],[x1  [1:*]]]), count-1)]
				fnc = [fnc1[0],reform(transpose([[fnc2],[fnc1[1:*]]]), count-1)]
			END
			count1 LT count2: BEGIN
				x   = [x2  [0],reform(transpose([[x1  ],[x2  [1:*]]]), count-1)]
				fnc = [fnc2[0],reform(transpose([[fnc1],[fnc2[1:*]]]), count-1)]
			END
			count1 EQ count2: BEGIN
				x   = reform(transpose([[x1  ],[x2  ]]), count)
				fnc = reform(transpose([[fnc1],[fnc2]]), count)
			END
			ENDCASE
		END
		ENDCASE
	ENDIF

	RETURN, x

ENDIF

;sz = size(Array)
;if sz[0] ge 1 and sz[sz[0]+2] ge 3^sz[0] then begin ; ??????


CASE strict OF

0: BEGIN

	n = n_elements(Array)
	p = where( Array[0:n-2] NE Array[1:n-1] ) ; Unequal neighbours

	IF n_elements(p) GE 2 THEN BEGIN

		;p = [0,p+1]						; Points with different neighbour to the left
		p = [p,n-1] 						; Points with different neighbour to the right
											; (for horizontal stretches use the rightmost point only)
		A = Array[p]						; There are at least 3 elements in A

		CASE minima OF
		0: x = A GT shift(A,-1) AND A GT shift(A,1)
		1: x = A LT shift(A,-1) AND A LT shift(A,1)
		endcase

		x = where(x[1:n_elements(x)-2])		; Discard the edges and check what's left

		IF x[0] NE -1 THEN BEGIN

			x = 1+x							; Index into A
			x = p[x]						; Index into Array

			CASE floatx OF
			0: fnc = Array[x]
			1: begin						; Quadratic fit with two neighbours
				gleft  = Array[x-1]-Array[x]
				gright = Array[x+1]-Array[x]
				dx = (gleft-gright)/(gleft+gright)*0.5
				fnc = Array[x]+0.5*dx*(gleft*(dx-1)+gright*(dx+1))
				x = x+dx
			END
			ENDCASE

			count = n_elements(x)

		ENDIF

	ENDIF

END

1: BEGIN

	IF n_elements(Array) GE 3 THEN BEGIN

		CASE minima OF
		0: x = Array GT shift(Array,-1) AND Array GT shift(Array,1)
		1: x = Array LT shift(Array,-1) AND Array LT shift(Array,1)
		endcase

		x = where(x[1:n_elements(x)-2])		; Discard the edges and check what's left

		IF x[0] NE -1 THEN BEGIN

			x = 1+x							; Add one to compensate for the left discarded edge

			CASE floatx OF
			0: fnc = Array[x]
			1: BEGIN						; Quadratic fit with two neighbours
				gleft  = Array[x-1]-Array[x]
				gright = Array[x+1]-Array[x]
				dx = (gleft-gright)/(gleft+gright)*0.5
				fnc = Array[x]+0.5*dx*(gleft*(dx-1)+gright*(dx+1))
				x = x+dx
			END
			ENDCASE

			count = n_elements(x)

		ENDIF

	ENDIF

END

ENDCASE

RETURN, x  &  END