FUNCTION max_pos, f

;+
; NAME:
;	max_pos	
; PURPOSE:
;	Finds the position of maximum in a 2D array
; CALLING SEQUENCE:
;	max = max_pos(f)
; INPUTS:
;	f	array
; OUTPUTS:
;	max	fltarr(2)	2D array containing the position of the maximum
;				max(0) = x-coordinate of maximum
;				max(1) = y-coordinate of maximum
;				(values based on array index)
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS
;	ArrayLocation
; PROCEDURE:
;	The indices of the element with the max. value are obtained from
;	a call to `max'.
;	If the element has a neighbour on either side (in horizontal or 
;	vertical direction) the x and y positions of the maximum are refined
;	by fitting a quadratic polynomial through the maximum array element
;	and its two neighbours.
; MODIFICATION HISTORY:
;	H. Cohl, 20 Nov, 1990 --- Initial implementation. 
;-

dim = size(f,/dim)

fx = max(f,ip)
ip = ArrayLocation(ip, dim=dim)

ix = ip[0]
iy = ip[1]

pp = float(ip)

IF ip[0] NE 0 AND ip[0] NE dim[0]-1 THEN	$
	pp[0] = pp[0]+.5*(f[ix+1,iy]-f[ix-1,iy])/(2*fx-f[ix-1,iy]-f[ix+1,iy])

IF ip[1] NE 0 AND ip[1] NE dim[1]-1 THEN	$
	pp[1] = pp[1]+.5*(f[ix,iy+1]-f[ix,iy-1])/(2*fx-f[ix,iy-1]-f[ix,iy+1])

RETURN, pp  &  END