pro pkfit,f,scale,x,y,sky,radius,ronois,phpadu,gauss,psf, $
                  errmag,chi,sharp,niter, DEBUG= debug
;+
; NAME:
;	PKFIT
; PURPOSE:
;	Subroutine of  GETPSF to perform a one-star least-squares fit 
; EXPLANATION:
;	Part of the DAOPHOT PSF photometry sequence
;
; CALLING SEQUENCE:
;	PKFIT, f, scale, x, y, sky, radius, ronois, phpadu, gauss, psf, 
;				errmag, chi, sharp, Niter, /DEBUG 
; INPUTS:
;	F      - NX by NY array containing actual picture data.           
;	X, Y   - the initial estimates of the centroid of the star relative
;		to the corner (0,0) of the subarray.  Upon return, the
;		final computed values of X and Y will be passed back to the
;		calling routine.
;	SKY  -   the local sky brightness value, as obtained from APER
;	RADIUS-  the fitting radius-- only pixels within RADIUS of the
;		instantaneous estimate of the star's centroid will be
;		included in the fit, scalar
;	RONOIS - readout noise per pixel, scalar
;	PHPADU - photons per analog digital unit, scalar
;	GAUSS -  vector containing the values of the five parameters defining
;		the analytic Gaussian which approximates the core of the PSF.
;	PSF   -  an NPSF by NPSF look-up table containing corrections from
;		the Gaussian approximation of the PSF to the true PSF.
;
; INPUT-OUTPUT:
;	SCALE  - the initial estimate of the brightness of the star,
;		expressed as a fraction of the brightness of the PSF.
;		Upon return, the final computed value of SCALE will be
;		passed back to the calling routine.
; OUTPUTS:
;	ERRMAG - the estimated standard error of the value of SCALE
;		returned by this routine.
;	CHI    - the estimated goodness-of-fit statistic:  the ratio
;		of the observed pixel-to-pixel mean absolute deviation from
;		the profile fit, to the value expected on the basis of the
;		noise as determined from Poisson statistics and the 
;		readout noise.
;	SHARP  - a goodness-of-fit statistic describing how much broader  
;		the actual profile of the object appears than the
;		profile of the PSF.
;	NITER -  the number of iterations the solution required to achieve
;		convergence.  If NITER = 25, the solution did not converge.
;		If for some reason a singular matrix occurs during the least-
;		squares solution, this will be flagged by setting NITER = -1.
;
; RESTRICTIONS:
;	No parameter checking is performed
; REVISON HISTORY:
;	Adapted from the official DAO version of 1985 January 25
;	Version 2.0 W. Landsman STX             November 1988
;	Converted to IDL V5.0   W. Landsman   September 1997
;-
 s = size(f)	;Get array dimensions
 nx = s[1] & ny = s[2]
;               ;Initialize a few things for the solution
 redo = 0B		
 pkerr = 0.027/(gauss[3]*gauss[4])^2
 clamp = fltarr(3) + 1.
 dtold = fltarr(3) 
 niter = 0                                
 chiold = 1.

 if keyword_set(DEBUG) then $
     print,'PKFIT: ITER  X      Y      SCALE    ERRMAG   CHI     SHARP'

BIGLOOP:  		        ;Begin the big least-squares loop
 niter = niter+1

 ixlo = fix(x-radius) > 0	;Choose boundaries of subarray containing
 iylo = fix(y-radius) > 0        ;points inside the fitting radius
 ixhi = fix(x+radius) +1 < (nx-1)
 iyhi = fix(y+radius) +1 < (ny-1)
 ixx  = ixhi-ixlo+1
 iyy  = iyhi-iylo+1
 dy   =  findgen(iyy) + iylo - y    ;X distance vector from stellar centroid
 dysq = dy^2
 dx   = findgen(ixx) + ixlo - x
 dxsq = dx^2
 rsq  = fltarr(ixx,iyy)  ;RSQ - array of squared 

 for J = 0,iyy-1 do rsq[0,j] = (dxsq+dysq[j])/radius^2

 ; The fitting equation is of the form
 ;
 ; Observed brightness =
 ;      SCALE + delta(SCALE)  *  PSF + delta(Xcen)*d(PSF)/d(Xcen) +
 ;                                           delta(Ycen)*d(PSF)/d(Ycen)
 ;
 ; and is solved for the unknowns delta(SCALE) ( = the correction to
 ; the brightness ratio between the program star and the PSF) and
 ; delta(Xcen) and delta(Ycen) ( = corrections to the program star's
 ; centroid).
 ;
 ; The point-spread function is equal to the sum of the integral under
 ; a two-dimensional Gaussian profile plus a value interpolated from            
 ; a look-up table.

 good = where(rsq lt 1.,ngood)
 ngood = ngood > 1

 t = fltarr(ngood,3)
 dx = dx[good mod ixx]
 dy = dy[good/ixx]
 model = dao_value(dx, dy, gauss, psf, dvdx, dvdy)

 if keyword_set(DEBUG) then begin print,'model created ' & stop & end

 t[0,0] = model
 t[0,1] = -scale*dvdx  
 t[0,2] = -scale*dvdy
 fsub = f[ixlo:ixhi,iylo:iyhi]
 fsub = fsub[good]
 rsq = rsq[good]
 df = fsub - scale*model - sky     ;Residual of the brightness from the PSF fit

 ; The expected random error in the pixel is the quadratic sum of
 ; the Poisson statistics, plus the readout noise, plus an estimated
 ; error of 0.75% of the total brightness for the difficulty of flat-
 ; fielding and bias-correcting the chip, plus an estimated error of 
 ; of some fraction of the fourth derivative at the peak of the profile,
 ; to account for the difficulty of accurately interpolating within the
 ; point-spread function.  The fourth derivative of the PSF is
 ; proportional to H/sigma**4 (sigma is the Gaussian width parameter for
 ; the stellar core); using the geometric mean of sigma(x) and sigma(y),
 ; this becomes H/ sigma(x)*sigma(y) **2.  The ratio of the fitting
 ; error to this quantity is estimated from a good-seeing CTIO frame to
 ; be approximately 0.027 (see definition of PKERR above.)

 fpos = (fsub-df) > 0	;Raw data - residual = model predicted intensity
 sigsq = fpos/phpadu + ronois + (0.0075*fpos)^2 + (pkerr*(fpos-sky))^2
 sig = sqrt(sigsq)
 relerr = df/sig

 ; SIG is the anticipated standard error of the intensity
 ; including readout noise, Poisson photon statistics, and an estimate
 ; of the standard error of interpolating within the PSF.  

 rhosq = fltarr(ixx,iyy)

 for j = 0,iyy-1 do rhosq[0,j] = (dxsq/gauss[3]^2+dysq[j]/gauss[4]^2)

 rhosq = rhosq[good]
 if (niter GE 2) then begin    ;Reject any pixel with 10 sigma residual
        badpix = where( ABS(relerr/chiold) GE 10.,nbad ) 
        if nbad GT 0 then begin 
            remove, badpix, fsub, df, sigsq, sig
            remove, badpix, relerr, rsq, rhosq
            ngood = ngood-badpix
        endif
 endif

 wt = 5./(5.+rsq/(1.-rsq))
 lilrho = where(rhosq LE 36.)	;Include only pixels within 6 sigma of centroid
 rhosq[lilrho] = 0.5*rhosq[lilrho]
 dfdsig = exp(-rhosq[lilrho])*(rhosq[lilrho]-1.)
 fpos = ( fsub[lilrho]-sky) >0 + sky

 ; FPOS-SKY = raw data minus sky = estimated value of the stellar
 ; intensity (which presumably is non-negative).

 sig  = fpos/phpadu + ronois + (0.0075*fpos)^2 + (pkerr*(fpos-sky))^2
 numer = total(dfdsig*df/sig)
 denom = total(dfdsig^2/sig)

 ; Derive the weight of this pixel.  First of all, the weight depends
 ; upon the distance of the pixel from the centroid of the star-- it
 ; is determined from a function which is very nearly unity for radii
 ; much smaller than the fitting radius, and which goes to zero for
 ;  radii very near the fitting radius.  

 chi = total(wt*abs(relerr))
 sumwt = total(wt)

 wt = wt/sigsq   ;Scale weight to inverse square of expected mean error
 if niter GE 2 then $	;Reduce weight of a bad pixel
         wt = wt/(1.+(0.4*relerr/chiold)^8)

 v = fltarr(3)       ;Compute vector of residuals and the normal matrix. 
 c = fltarr(3,3)

 for kk = 0,2 do begin   
    v[kk] = TOTAL(df*t[*,kk]*wt)
    for ll = 0,2 do C[kk,ll] = TOTAL(t[*,kk]*t[*,ll]*wt)
 end 

 ; Compute the (robust) goodness-of-fit index CHI.
 ; CHI is pulled toward its expected value of unity before being stored
 ; in CHIOLD to keep the statistics of a small number of pixels from
 ; completely dominating the error analysis.

 if sumwt GT 3.0 then begin  
   chi = 1.2533*chi*sqrt(1./(sumwt*(sumwt-3.)))
   chiold = ((sumwt-3.)*chi+3.)/sumwt
 endif

 C = INVERT(C)	;Invert the normal matrix
 dt = c#v 	;Compute parameter corrections

; In the beginning, the brightness of the star will not be permitted
; to change by more than two magnitudes per iteration (that is to say,
; if the estimate is getting brighter, it may not get brighter by
; more than 525% per iteration, and if it is getting fainter, it may
; not get fainter by more than 84% per iteration).  The x and y
; coordinates of the centroid will be allowed to change by no more
; than one-half pixel per iteration.  Any time that a parameter
; correction changes sign, the maximum permissible change in that
; parameter will be reduced by a factor of 2.

 div = where( dtold*dt LT -1.e-38, nbad )
 if nbad GT 0 then clamp[div] = clamp[div]/2.
 dtold = dt
 adt = abs(dt)

 scale = scale+dt[0]/    $
  (1.+(( dt[0]/(5.25*scale)) > (-1*dt[0]/(0.84*scale)) )/clamp[0])
 x = x + dt[1]/(1.+adt[1]/(0.5*clamp[1]))
 y = y + dt[2]/(1.+adt[2]/(0.5*clamp[2]))
 redo = 0B

; Convergence criteria:  if the most recent computed correction to the
; brightness is larger than 0.1% or than 0.05 * sigma(brightness),
; whichever is larger, OR if the absolute change in X or Y is
; greater than 0.01 pixels, convergence has not been achieved.

 sharp = 2.*gauss[3]*gauss[4]*numer/(gauss[0]*scale*denom)
 errmag = chiold*sqrt(c[0,0])
 if ( adt[0] GT ( 0.05*errmag > 0.001*scale )) then redo = 1b
 if ((adt[1] > adt[2] ) GT 0.01) then redo = 1b

 if keyword_set(DEBUG) then print,format='(1H  ,I9,2F7.2,2F9.3,F8.2,F9.2)', $ 
                       niter,x,y,scale,errmag,chiold,sharp
 if niter LT 3 then goto, BIGLOOP 	;At least 3 iterations required

; If the solution has gone 25 iterations, OR if the standard error of 
; the brightness is greater than 200%, give up.

 if (redo and (errmag LE 1.9995) and (niter LT 25) ) then goto, BIGLOOP
 sharp = sharp>(-99.999)<99.999

 return
 end