FUNCTION FOCUS3,IM,DARK=DARK,FLAT=FLAT,PEAK=PEAK,SHARP=SHARP, $
	TV=TV, normalize=norm, quadratic=quad
;+
;NAME:
;       FOCUS3
;PURPOSE:
;       To find 3 positions of best focus by fitting 3 sharpness parameters
;SAMPLE CALLING SEQUENCES:
;       focus=focus3(images,dark=dk,flat=ff,peak=pk,sharp=sh)
;INPUT:
;       IM	array(NX,NY,NF) of image intensities at uniformly spaced
;		focus positions 0 to NF-1
;OUTPUT:
;	FOCUS3	floating array(3), best focus positions based on maxima of:
;		  X-sharpness, mean absolute value of dI/dx
;		  Y-sharpness, mean absolute value of dI/dY
;		  R-sharpness, mean of sqrt(dI/dx^2 + dI/dy^2)
;OPTIONAL KEYWORD INPUT:
;	DARK	dark current array to subtract from IM & FLAT
;	FLAT	way out-of-focus image array, to use for flat field correction
;	TV	if set, then displays blowups of FLAT and each corrected image
;	NORMALIZE  if set, normalize sharpness to average intensity
;	QUADRATIC  if set, uses mean-squared derivatives, not absolute values
;OPTIONAL OUTPUT:
;	PEAK	floating array(3), values of maximum sharpness, at each of the
;		3 best focus positions
;	SHARP	floating array(3,NF), containing each of the sharpness
;		parameters for each input image
;HISTORY:
;       Written Jan 1996 by T. Tarbell

NX=N_ELEMENTS(IM(*,0,0))
NY=N_ELEMENTS(IM(0,*,0))
NF=N_ELEMENTS(IM(0,0,*))
FAC = 1./((NX-1.)*(NY-1.))
IF (KEYWORD_SET(FLAT)) THEN BEGIN
	FF = FLOAT(FLAT)
	IF (KEYWORD_SET(DARK)) THEN FF=FF-DARK
	GG = TOTAL(FF)/N_ELEMENTS(FF)
	G  = GG/(FF > 1.)
	IF (KEYWORD_SET(TV)) THEN tvscl,rebin(ff,3*nx,3*ny),0
ENDIF

SH=FLTARR(3,NF)
FOR I=0,NF-1 DO BEGIN
	II = FLOAT(IM(*,*,I))
	IF (KEYWORD_SET(DARK)) THEN II=II-DARK
	IF (KEYWORD_SET(FLAT)) THEN II=II*G
	IF (KEYWORD_SET(TV)) THEN tvscl,rebin(ii,3*nx,3*ny),i+1
	IF (KEYWORD_SET(NORM)) THEN $ 
	   FAC = 1./TOTAL(II(0:NX-2,0:NY-2))
	DX = II(1:*,0:NY-2)-II(0:NX-2,0:NY-2)
	DY = II(0:NX-2,1:*)-II(0:NX-2,0:NY-2)
	IF (KEYWORD_SET(QUAD)) THEN BEGIN
	  IF (KEYWORD_SET(NORM)) THEN $ 
	   FAC = (NX-1.)*(NY-1.)*FAC^2
	  SH(0,I) = FAC*TOTAL(DX^2)
	  SH(1,I) = FAC*TOTAL(DY^2)
	  SH(2,I) = SH(0,I)+SH(1,I)
	END ELSE BEGIN
	  SH(0,I) = FAC*TOTAL(ABS(DX))
	  SH(1,I) = FAC*TOTAL(ABS(DY))
	  SH(2,I) = FAC*TOTAL(SQRT(DX^2+DY^2))
	END
ENDFOR

IF (KEYWORD_SET(SHARP)) THEN SHARP=SH

FOC = FLTARR(3)
PK = FLTARR(3)
FOR J=0,2 DO BEGIN
  SM = MAX(SH(J,*),II)
  IF ((II EQ 0) OR (II EQ NF-1)) THEN BEGIN
	FOC(J) = II
	PK(J) = SM
  ENDIF ELSE BEGIN
	LIM = II+INDGEN(4)-2
	IF (SH(J,II-1) LT SH(J,II+1)) THEN LIM=LIM+1
	IF (MIN(LIM) LT 0) THEN LIM=LIM-MIN(LIM)
	IF (MAX(LIM) GT NF-1) THEN LIM=LIM-(MAX(LIM)-(NF-1))
	C = POLY_FIT(LIM,SH(J,LIM),2)
	FOC(J) = -C(1)/(2.*C(2))
	PK(J) = POLY(FOC(J),C)
  ENDELSE
ENDFOR
IF (KEYWORD_SET(PEAK)) THEN PEAK=PK
RETURN,FOC+1.
END