;+
; NAME:
;	smei_star_lsqfit
; PURPOSE:
;	Performs multi-dimensional least-squares fit between
;	observed star and standard star
; CATEGORY:
;	camera/idl/star
; CALLING SEQUENCE:
	FUNCTION smei_star_lsqfit, $
		rr				, $
		standard_stars	, $
		original_star	, $
		in_fullpsf		, $
		in_fitpsf		, $
		in_bckgnd		, $
		back			, $
		bright			, $
		sigma			= sigma			, $
		only			= only			, $
		noplane			= noplane		, $
		noyzero			= noyzero		, $
		bkgnd_model	  	= bkgnd_model	, $
		one_star_model	= one_star_model, $
		use_weights 	= use_weights	, $
		original_weights= original_weights,$
		fixed_bright	= fixed_bright
; INPUTS:
;	rr					array[2,n]; type: float
;							position vectors for all bins in the box around
;							the star, relative to the star centroid
;							(in units of the bin size of the skymap from
;							which the star was extracted).
;	standard_stars		array[n,m]; type: float
;							intensities in equivalent locations in standard
;							star for all bins around star
;							The 2nd dimension indicates the m stars fitted
;							simultaneously (absent if m=1)
;	original_star		array[nx,ny] (with nx*ny=n); type: float
;							intensities in box around star as extracted
;							from the original skymap.
;	in_fullpsf			array; type: long
;							list of array indices into rr, standard_stars and
;							original_star that defineds the full area of the
;							PSF to be used to define the star model
;	in_fitpsf			array; type: long
;							list of array indices into rr, standard_stars and
;							original_star that should be included in the fit
;							for the star PSF.
;							in_fitpsf is a subset of in_fullpsf
;							(see href=smei_star_stdmaps=).
;	in_bckgnd			array; type: long
;							list of array indices to be used to fit the background.
;							If this is a valid list of indices (i.e. if
;							in_bckgnd[0] NE -1 then these pixels are used to fit
;							a background separately from the fit to the PSF.
;							The fit is to a plane or, if /noplane, is SET the
;							mean of the background pixels is used.
; OPTIONAL INPUTS:
;	sigma=sigma 		scalar; type: float
;							if set then a second fit is made after discarding
;							bins that are more than 'sigma' standard deviations
;							in the residual after the first fit.
;	only=only			array; type: integer
;							set of indices into the 'ifit' array
;							Only the subset of bins in 'only' are used in the fit
;	/noplane			if set, then fit constant background (i.e. do
;						not fit a sloped background
;	/noyzero 			passed to href=lsqLinearFit= for fitting of PSF only if
;						background and PSF are fitted separately (i.e. if
;						in_bckgnd contains a valid list of background bins)
;	/use_weights		if set, the fit is done using the inverse of the skymap
;							brightness as weight
;	original_weights	array[nx,ny] (with nx*ny=n); type: float;
;							default: 1.0/(50.0+original_star) (See PROCEDURE)
;							Only used if /use_weights SET.
;							weights to be used for each skybin.
; OUTPUTS:
;	Result				array[nx,ny]; type: float
;							the resulting lsq fit to the PSFs of all stars fitted
;							simultaneously.
;							star_model[nx,ny] = (Sum(m))bright[m]*standard_stars[n,m]
;	back				array[3]; type: float
;							constants for fit of background
;							back[0] intensity at the centroid of the star
;								(i.e. at the origin of the rr position vectors)
;							back[1:2] horizontal and vertical slope of the background
;	bright				array[m]; type: float
;							brightness fit for all stars
; OPTIONAL OUTPUTS:
;	bkgnd_model=bkgnd_model
;				 		array[nx,ny]; type: float
;							the resulting lsq fit for the planar background
;	one_star_model=one_star_model
;						array[nx,ny,n]; type: float
;							the resulting lsq fit to the brightest of the stars
;							one_star_model[nx,ny,i] = bright[i]*standard_stars[n,i]
;							This is equal to the return value if only one star is fitted
; INCLUDE:
	@compile_opt.pro		; On error, return to caller
; CALLS:
;	InitVar, lsqLinearFit, IsType
; PROCEDURE:
;	Note if /use_weights is set then the default weight used
;	is 1.0/(50.0+original_star).
;	The 50.0 is an estimate for the brightness of the darkest sky.
;	The reason for using it is to avoid extremely large weights for
;	background skybins when original_star already has the sidereal
;	background and the zodiacal light removed (resulting in a
;	background that is essentially zero).
; MODIFICATION HISTORY:
;	DEC-2006, Paul Hick (UCSD/CASS)
;	AUG-2006, Paul Hick (UCSD/CASS)
;		Added in_bckgnd argument to be able to fit background and
;		star separately
;	APR-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu)
;		Added keyword original_weights.
;		Changed default weights from 1/original_star
;		to 1/(50+original_star).
;-

InitVar, only	, -1
InitVar, noplane, /key
InitVar, noyzero, 0
InitVar, use_weights, /key

splitfit = in_bckgnd[0] NE -1

ndim = size(original_star,/dim) 		; Nr of positions in standard star

CASE size(standard_stars,/n_dim) OF
1: nsub = 1
2: nsub = (size(standard_stars,/dim))[1]; Nr of simultaneously fitted stars
ENDCASE

CASE only[0] NE -1 OF
0: in_fitpsf_only = in_fitpsf
1: in_fitpsf_only = in_fitpsf[only]		; Limit fit to bins specified in 'only'
ENDCASE

IF use_weights THEN BEGIN
	InitVar, original_weights, $
		1.0/(50.0+(original_star-min(original_star[in_fitpsf_only]))), $
		set=weights
	weights = weights[in_fitpsf_only]
ENDIF

; Fit the background and star brightness.
; Note that when multiple stars are fitted the array 'standard_stars'
; contains brightnesses for multiple stars.


; Check for stars where the brightnesses are fixed.

CASE IsType(fixed_bright,/defined) OF
0: n_fixed = 0
1: fixed = where(finite(fixed_bright),n_fixed,complement=not_fixed,ncomplement=n_not_fixed)
ENDCASE

IF n_fixed NE 0 THEN BEGIN

view, in=magnifyarray(original_star < 150,5), /stretch, /bl

	; Save. Will be restored after the fitting has been done

	original_save = original_star
	standard_save = standard_stars
	nsub_save     = nsub

	; Remove the stars with fixed brightnesses from 'original_star'
	; Only the non-fixed, to-be-fitted stars remain

	FOR i=0,n_fixed-1 DO	$
		original_star -= fixed_bright[fixed[i]]*standard_stars[*,fixed[i]]

	; Only retain the standard_stars for the non-fixed,
	; to-be-fitted stars.

	CASE n_not_fixed NE 0 OF
	0: standard_stars = -1				; No stars left to fit (only background)
	1: standard_stars = standard_stars[*,not_fixed]
	ENDCASE

	nsub = n_not_fixed

view, in=magnifyarray(original_star < 150,5), /stretch, /br

ENDIF


CASE noplane OF

0: BEGIN								; Fit a sloped background

	CASE splitfit OF

	0: BEGIN							; Fit background and star simultaneously

		CASE nsub NE 0 OF

		0: BEGIN						; No stars to fit; only background
			back = lsqLinearFit(original_star[in_fitpsf_only]	, $
				rr[*,in_fitpsf_only], $
				weights = weights)		; Constant background, hor & vert slope

			destroyvar, bright			; No star brightnesses
		END

		1: BEGIN						; Fit nsub stars and background
			tmp = lsqLinearFit(original_star[in_fitpsf_only]	, $
				[rr[*,in_fitpsf_only],transpose(standard_stars[in_fitpsf_only,*])], $
				weights = weights)

			back   = tmp[0:2]			; Constant background, hor & vert slope
			bright = tmp[2+1:2+nsub]	; Star brightnesses
		END

		ENDCASE

	END

	1: BEGIN							; Fit background first, then star separately

		; Fit a plane through the background pixels

		back = lsqLinearFit(original_star[in_bckgnd],rr[*,in_bckgnd])

		; Subtract the planar background from the pixels inside the PSF

		CASE nsub NE 0 OF

		0: destroyvar, bright

		1: BEGIN
	
			tmp = (original_star-(back[0]+back[1]*rr[0,*]+back[2]*rr[1,*]))[in_fitpsf_only]

			; The residual intensity is used to fit the star.
			; The return array 'bright' will have 1+nsub elements.
			; bright[0] is a constant background.
			; If the PSF defined by the standard star is perfect then bright[0] would be zero.

			bright = lsqLinearFit(tmp		, $
					transpose(standard_stars[in_fitpsf_only,*]), $
					weights = weights 		, $
					yzero	= 1-noyzero)

			; Add constant from star fit to constant for background
			; Note that if noyzero is NOT set then bright[0] will be zero

			back[0] += bright[0]
			bright   = bright[1:nsub]	; Fit to star brightness

		END

		ENDCASE

	END
	ENDCASE

END

1: BEGIN								; Fit a constant background

	CASE splitfit OF

	0: BEGIN							; Fit background and star simultaneously

		CASE nsub NE 0 OF

		0: BEGIN						; No stars to fit; only background
			back = mean(original_star[in_fitpsf_only])
			;back = lsqLinearFit(original_star[in_fitpsf_only]	, $
			;	weights = weights)

			back = [back,0.0,0.0]		; Contant background, hor & vert slope zero
			destroyvar, bright			; Star brightnesses
		END

		1: BEGIN						; Fit nsub stars
			tmp = lsqLinearFit(original_star[in_fitpsf_only]	, $
				transpose(standard_stars[in_fitpsf_only,*]) 	, $
				weights = weights)

			back   = [tmp[0],0.0,0.0]	; Contant background, hor & vert slope zero
			bright = tmp[1:nsub]		; Star brightnesses
		END

		ENDCASE

	END

	1: BEGIN							; Fit background first, then star separately

		; The constant background is just the mean over all background pixels
		; Explicitly add slopes of zero.

		back = [mean(original_star[in_bckgnd]),0.0,0.0]

		; Subtract the constant background from the pixels inside the PSF

		CASE nsub NE 0 OF

		0: destroyvar, bright

		1: BEGIN
			tmp = original_star[in_fitpsf_only]-back[0]

			; The residual intensity is used to fit the star.
			; The return array 'bright' will have 1+nsub elements.
			; bright[0] is a constant background.
			; If the PSF defined by the standard star is perfect then bright[0] would be zero.

			bright = lsqLinearFit(tmp	, $
				transpose(standard_stars[in_fitpsf_only,*]), $
				weights = weights		, $
				yzero	= 1-noyzero)

			; Add constant from star fit to constant for background
			; Note that if ynozero is NOT set then bright[0] will be zero

			back[0] += bright[0]
			bright   = bright[1:nsub]		; Fit to star brightness
		END

		ENDCASE

	END

	ENDCASE

END

ENDCASE

; If some of the star brightnesses were fixed, the 'bright'
; array needs to be updated.

IF n_fixed NE 0 THEN BEGIN

	IF n_not_fixed NE 0 THEN tmp = bright	; Fitted star brightnesses
	bright = fixed_bright
	IF n_not_fixed NE 0 THEN bright[not_fixed] = tmp

	; Restore variables saved prior to fit.

	original_star  = original_save
	standard_stars = standard_save
	nsub		   = nsub_save

ENDIF

bkgnd_model = reform(back[0]+back[1]*rr[0,*]+back[2]*rr[1,*],ndim)

; Construct the 'model star' by adding together all brightness fits
; If multiple stars are fit then one_star_model is the model for
; the individual stars

star_model		= fltarr(ndim)
one_star_model	= fltarr([ndim[0]*ndim[1],nsub])
FOR isub=0,nsub-1 DO BEGIN
	tmp = bright[isub]*standard_stars[in_fullpsf,isub]
	star_model	  [in_fullpsf	  ] += tmp
	one_star_model[in_fullpsf,isub]  = tmp
ENDFOR
;star_model = total(one_star_model,1)

IF IsType(sigma,/defined) THEN BEGIN

	; Calculate the residuals (difference between original star and fitted star
	; Find pixels that fit better than the threshold, and use those to refine the fit

	delta_star = (original_star-(bkgnd_model+star_model))[in_fitpsf]
	std  = stddev(delta_star,/nan)
	only = where(abs(delta_star) LE sigma*std)			; Subset of in_fitpsf
	IF only[0] NE -1 THEN	$							; Recursive call WITHOUT sigma keyword
		star_model = smei_star_lsqfit(	$
			rr				, $
			standard_stars	, $
			original_star	, $
			in_fullpsf		, $
			in_fitpsf		, $
			in_bckgnd		, $
			back			, $
			bright			, $
			only			= only		 , $
			noplane			= noplane	 , $
			noyzero			= noyzero	 , $
			bkgnd_model		= bkgnd_model, $
			one_star_model	= one_star_model,$
			use_weights 	= use_weights, $
			original_weights= original_weights)

ENDIF

RETURN, star_model  &  END