function GetStarShape, Frame, bScan, Star, Box, Mask, BoxBase, BoxDev, BoxMed,	$
	nocentroid=NoCentroid,file=File, nofit=NoFit, nearbystar=NearbyStar

@compile_opt.pro		; On error, return to caller

;+
; NAME:
;	GetStarShape
; PURPOSE:
;	Get a set of shapes in the form of a 2D rectangular template to be used
;	to detect stars in an image
; CATEGORY:
; CALLING SEQUENCE:
; INPUTS:
;	Frame		2D array, any type
;					Image containing stars
;	Star		scale, type integer
;					Location of suspected star in Frame (usually a local
;					maximum in Frame)
; OPTIONAL INPUT PARAMETERS:
;	/nocentroid	if set then no attempt is made to find the centroid of the
;				suspected star (the centroid is supposed to be at the location
;				of the maximum at position Star).
;	file=File	character, default: starshape.txt in directory $SMEI/sys
;					Name of file containing the basic star template
; 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 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 area in the basic star mask
;					defined as 'inside the star', i.e. the part of the 2D
;					template actually used to match the peak in Frame to a star.
;	BoxBase		scalar, type float
;	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
;
; OPTIONAL OUTPUT PARAMETERS:
; CALLS:
;	ArrayLocation, SuperArray, flt_read
; COMMON BLOCKS:
;	common SaveStarShape, Shape, X, Y, InStar
; PROCEDURE:
; >	The basic template, a 2D array of numbers is stored in an ascii file and
;	read by flt_read. The template is separated in an 'inside star' and
;	'outside star' area. The 'outside' area is used to define a base level,
;	which is subtracted from the 'inside' area. The remaining difference is
;	the basic template used for fitting.
; > If the /nocentroid keyword is not set, then a set of templates is calculated
;	corresponding to set of slightly different centroids for the suspected star.
; MODIFICATION HISTORY:
;	OCT-1998, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;-

common SaveStarShape, Shape, X, Y, InStar

sz = size(Frame)
S  = ArrayLocation(Star[0],size=sz)			; X,Y position of maximum in Frame

NoCentroid	= keyword_set(NoCentroid)

if n_elements(Shape) eq 0 then begin
	if n_elements(File) eq 0 then									$
		File = filepath(root=getenv('SMEI'),subdir='sys','starshape.txt')	$
	else		$
		File = File[0]

	status = flt_read(File, Shape)
	Shape = float(Shape)
	W = size(Shape,/dim)						; Dimension of shape (assumed to be the same in X and Y)
	X = lindgen(W[0])# replicate(1L,W[1])		; X-coordinate
	Y = replicate(1L,W[0])# lindgen(W[1])		; Y-coordinate

	InStar	= Shape gt 131						; Define 'inside star'

	FI = Shape[where(1B-InStar)]				; 'Outside star area
	ShapeMed = median(FI)
	ShapeAve = mean  (FI)
	ShapeDev = stddev(FI)						; Standard deviation outside star

	ShapeDev = stddev(FI[where(FI le ShapeMed)])
	ShapeBase = ShapeMed
;	ShapeBase = ShapeAve

	Shape	= Shape-ShapeBase					; Subtract base from 'Shape'

endif

W = (size(Shape))[1:2]
R = (W-1)/2

XS	= S[0]-R[0]+X								; Index range in Frame around maximum in S
YS	= S[1]-R[1]+Y

Mask = where(0 le XS and XS le sz[1]-1 and 0 le YS and YS le sz[2]-1)

Box	 = replicate(ValuesNaN(sz),W[0],W[1])
Box[Mask] = Frame[XS[Mask],YS[Mask]]

bBox = bytarr(W[0],W[1])
bBox[Mask] = bScan[XS[Mask],YS[Mask]]

Mask = where(finite(Box) and 1B-InStar)			; Finite values outside star

BoxMed = median(Box[Mask])
BoxAve = mean  (Box[Mask])						; Output: Base of Box from pixels outside star
BoxDev = stddev(Box[Mask])						; Output: standard dev in base

BoxDev = stddev((Box[Mask])[where(Box[Mask] le BoxMed)])
BoxBase = BoxMed
;BoxBase = BoxAve

;print, shapebase,shapedev,boxbase,boxdev,boxmed

; 'Mask' always has at least on element (the maximum)
; The difference between finite(Box) and bBox is that bBox also includes points that
; have been modified as the result of fitting a brighter star. These pixels were
; set to a median value during the earlier fit, and distort the shape of the star being
; fitted now.

Mask = where(finite(Box) and InStar)
;Mask = where(bBox and InStar)
NearbyStar = n_elements(where(bBox and InStar)) ne n_elements(where(finite(Box) and InStar))

if NoFit or NoCentroid then return, Shape		; Box, Mask, Base have been set

FI	= Shape[Mask]
dX	= total((X[Mask]-R[0])*FI)/total(FI)		; Template centroid position relative to center
dY	= total((Y[Mask]-R[1])*FI)/total(FI)

;if n_elements(Mask) eq n_elements(where(InStar)) then begin
	FI	= Box[Mask]-BoxBase
	dXS	= total((X[Mask]-R[0])*FI)/total(FI)	; Box centroid position relative to center
	dYS	= total((Y[Mask]-R[1])*FI)/total(FI)
;endif else begin
;	dXS = dX
;	dYS = dY
;endelse

nTable = 9
Step = 0.03
if n_elements(Star) gt 1 then begin
	nTable = 9
	Step = 0.1
endif

One = -(nTable-1)/2.+findgen(nTable)
One = One*Step

dX  = dX-dXS+One								; Centroid difference in middle of dX,dY
dY  = dY-dYS+One

Shapes = SuperArray(Shape,nTable*nTable,/trail)
Shapes = reform(Shapes,W[0],W[1],nTable,nTable,/overwrite)

One = replicate(1.,nTable)
FI	= fltarr(W[0],W[1],nTable)

; First spline in horizontal direction (along rows in Shape)

for j=0,W[1]-1 do begin							; Loop over rows
	V = where(Y eq j and InStar,nV)				; Pixels in star on row j
	if nV ne 0 then begin
		XV = X[V]								; Pixel X-coordinate
		FV = Shape[V]							; Pixel function values
		F2 = spl_init(XV,FV)					; Second derivatives
		XI = XV#One+replicate(1.,nV)#dX			; X-ccordinates for offsets (nTable per pixel)
		FI[XV,j,*] = spl_interp(XV,FV,F2,XI)
	endif
endfor

; Now spline in the vertical direction (along columns)

for i=0,W[0]-1 do begin							; Loop over columns
	V = where(X eq i and InStar,nV)				; Pixels in star on column i
	if nV ne 0 then begin
		XV = Y[V]								; Pixel Y-coordinates
	    for j=0,nTable-1 do begin				; Loop over horizontal offsets
			FV = FI[i,XV,j]
			F2 = spl_init(XV,FV)				; Second derivatives
			XI = XV#One+replicate(1.,nV)#dY
			Shapes[i,XV,j,*] = spl_interp(XV,FV,F2,XI)
		endfor
	endif
endfor

;if n_elements(Star) gt 1 then begin
;	nScale	= 9
;	Scale	= .9^(indgen(nScale))
;
;	Shapex	= Shapes
;	Shapes	= SuperArray(Shapex,nTable*nTable*nScale,/trail)
;	Shapes	= reform(Shapes,W[0],W[1],nTable,nTable,nTable,nTable,nScale,/overwrite)
;
;	for i1=0,nTable-1 do begin
;		for j1=0,nTable-1 do begin
;			for i2=0,nTable-1 do begin
;				for j2=0,nTable-1 do begin
;					for k=0,nScale-1 do begin
;						Shapes[*,*,i1,j1,i2,j2,k] = Shapex[*,*,i1,j1]+Shapex[*,*,i2,j2]*Scale[k]
;					endfor
;				endfor
;			endfor
;		endfor
;	endfor
;endif

return, reform(Shapes)
end