function FitStar, Frame,bScan,Star,Box,Mask,BoxBase,BoxDev,BoxMed,BoxM,	$
	template=Template, DevRatio=DevRatio, nofit=NoFit,nearbystar=NearbyStar

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	FitStar
; PURPOSE:
;	Fit star based on a 2D template of intensities
; CATEGORY:
; CALLING SEQUENCE:
; INPUTS:
;	Frame		2D array, any type
;					image with suspected star at location Star
;	Star		scalar, type integer
;					location of maximum
; OPTIONAL INPUT PARAMETERS:
;	/template	if set then argument Mask returns the exact mask used to fit the star
;				template. By default, Mask retuns a mask covering the central peak in box.
; OUTPUTS:
;	Box			2D array, type float, same size as the star template used to fit the star
;					the suspected star (pixel Star in Frame) is in the center of Box.
;					Box contains the square from Frame around the suspected star.
;					Values !VALUES.F_NAN near the edges of Box indicate pixels outside
;					the range of the Frame array.
;	Mask		1D array, type long integer
;					mask of pixels covering the central peak in Box
;	BoxBase
;	BoxDev		scalar, type float
;					measure for the fluctuations of the values in the 'outside star' pixels
;					in box (~standard deviation)
;	BoxMed		scalar, type float
;					median of the 'outside star' pixels in Box
;	Residue		2D array, type float, same size as Box
;					residual intensities left after subtraction of the Star template
; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
;	FindPeaks, GetStarShape
; PROCEDURE:
; >	A set of shapes for matching the star is collected from the function
;	GetStarShape. The template shape best fitting the star is subtracted.
; >	A single parameter, the total intensity in the star, is fitted.
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

Template = keyword_set(Template)
NoFit	 = keyword_set(NoFit)

Shape	= GetStarShape(Frame,bScan,Star,Box,Mask,BoxBase,BoxDev,BoxMed,nofit=NoFit,nearbystar=NearbyStar)

sz		= size(Shape)
nz		= sz[sz[0]+2]/(sz[1]*sz[2])					; # shapes
Shape	= reform(Shape,sz[1],sz[2],nz,/overwrite)	; Reform to 3D array

BoxM	= Box[Mask]-BoxBase						; Subtract Base from star
FullAve	= mean  (BoxM)
FullDev	= stddev(BoxM)

nM		= n_elements(Mask)
OneM	= replicate(1L,nM)
OneZ	= replicate(1L,nz)
i		= Mask#OneZ+OneM#(sz[1]*sz[2]*lindgen(nz))
i		= reform(i,nM*nz)
Shape	= reform(Shape[i],nM,nz)
i		= total(Shape,1)/nM
Shape	= Shape/(OneM#i)
BoxM	= BoxM#OneZ
BoxM	= BoxM-FullAve*Shape

i = total(BoxM,1)
j = total(BoxM*BoxM,1)
ResidueDev = sqrt((nM*j-i*i)/(nM*(nM-1)))				; Standard deviation: Sum[(Xi-Xmean)^2]/(N-1)

if sz[0] ge 3 then ResidueDev = reform(ResidueDev,sz[3:sz[0]] ,/overwrite)

j = min(ResidueDev,i)

BoxM = BoxM[*,i]
ResidueDev = ResidueDev[i]
DevRatio = ResidueDev/FullDev

; At this point Mask covers the pixels defined by the size of the star template.
; For bright stars this may not be big enough, for faint stars it may be too big.
; Instead set up a mask based on the size of the peak at the center of Box.
; Loop over all maxima in Mask above a certain level and collect all pixels
; assigned to these maxima by FindPeaks.

if not Template then begin
	Peaks = (sz[1]*sz[2]-1)/2							; Peak at the center of Box

	nz = where(finite(Box))
	i  = where(Box[nz] gt BoxMed+2.5*BoxDev)
	if i[0] ne -1 then begin							; Map finite values above threshold only
		nz = FindPeaks(Box,mask=nz[i],map=Map)			; Find local maxima above threshold
		if n_elements(nz) gt 1 then begin				; More then one maximum (1st always is center of box)
			nz = nz[1:*]
			i  = where(Box[nz] gt BoxMed+3.*BoxDev)		; Keep only maxima that are high enough
			if i[0] ne -1 then Peaks = [Peaks,nz[i]]
		endif

		nz = where(Map eq Peaks[0])							; Pixels around central peak in Box
		for i=1,n_elements(Peaks)-1 do	$
			if (where(Peaks[i] eq Mask))[0] ne -1 then nz = [nz,where(Map eq Peaks[I])]

		Mask = nz[uniq(nz,sort(nz))]						; Only keep unique entries

	endif else	$
		Mask = Peaks

endif

; The fraction of Shape subtracted from Box has the same total intensity (and hence
; the same mean) as the peak in Box. So the residue will always have a mean very
; close to zero. The standard deviation is used to decide whether the fit was any good.

return, 1.7*ResidueDev lt FullDev  &  end