;+
; NAME:
;  center
;
; PURPOSE: (one line)
;  return a star center
;
; DESCRIPTION:
;  return a star center
;
; CATEGORY:
;  IMAGES
;
; CALLING SEQUENCE:
;  center, d, xin, yin, xout, yout, method=method, nomax=nomax, $
;            miss=miss, ap=ap, skyin=skyin, skyout=skyout, count=counts, 
;            skyval=skyval, fwhm=fwhm
; INPUTS:
;  d : the two-dimensional image
;  xin, yin: initial guess at the star position
; OPTIONAL INPUT PARAMETERS:
;  method: boxm, centrod, cntrd, or psf_1_gauss
;    [default = centrod]
;  nomax : don't begin by centering up on the local maximum
;  miss : allowable miss distance between xin,yin and the 
;         x, y of maximum flux [default = 10]
;  ap : 
;    boxm - not used
;    centrod - radius used for calculating the center of max
;    cntrd - size of box to use for calculating center
;    psf_1_gauss - radius of pixels to include
;    [default = 3]
; 
; skyin, skyout: inner and outer radii for calculating the 
;    sky background [default = 10,20]
;
; OUTPUTS:
;  xout, yout : calculated x,y position
;     boxm - location of the maximum (or passed xin,yin
;            if nomax=1
;     centrod - "center of mass" using Buie's centrod
;     cntrd - position where X and Y derivatives go to zero using
;         daophot's cntrd
;     psf_1_gauss - fit of an unweighted circular gaussian
;
; OUTPUT KEYWORDS:
;  ckyval : value of the sky from the sky aperture
;    calculated using robomean
;  counts : counts in the star
;     boxm - the value at xout,yout
;     centrod - "mass" of the object in counts from buie's centrod
;     cntrd - aperture photometry using daophot's cntrd
;     psf_1_gauss - fit of an unweighted circular gaussian
;         returns mathematical total under the PSF, not just
;         from the center to ap
;  fwhm :
;     for psf_1_gauss, the fitted FWHM
;
; SIDE EFFECTS:
;
; EXAMPLE:
;
; MODIFICATION HISTORY:
;  Written  14 Jun 2007 Leslie Ann Young SwRI
;  Modified 13 Sep 2007 LAY
;  Modified 07 Oct 2007 LAY set counts, fwhm if off frame
;-

pro center, d, xin, yin, xout, yout, method=method, nomax=nomax, $
            miss=miss, ap=ap, skyin=skyin, skyout=skyout, $
            counts=counts, skyval=skyval, fwhm=fwhm, debug=debug

    if not keyword_set(miss) then miss = 10.
    if not keyword_set(ap) then ap = 4.
    if not keyword_set(skyin) then skyin = 10.
    if not keyword_set(skyout) then skyout = 20.
    if not keyword_set(method) then method = 'centrod'
    debug = keyword_set(debug)
   
    nxnyn, d, nx, ny

    if xin lt 1 or xin ge nx-2 or yin lt 1 or yin ge ny-2 then begin
        xout = xin
        yout = yin
        counts = 0
        fwhm = -1
        return
    end
 
    if debug then print, 'center : xin, yin = ', xin, yin

    ; center up on the maximum, unless nomax = 1
    if not keyword_set(nomax) then begin
        boxm, d, xin, yin, miss, miss, x, y
        x = float(x) & y = float(y)
    endif else begin
        x = xin & y = yin
        if debug then print, 'center : x, y = ', x, y
    endelse

    ; calculate the background
    Getannul, d, x, y,  skyin, skyout, skybuf
    Robomean, skybuf, 3.0, 0.5, sky, avgdev, skysig, var, skew, kurt, nsky
    db = d - sky

    case method of
        'boxm': begin
            counts = db[round(x), round(y)]
            xout = x 
            yout = y
        end

        'centrod': begin
            ; set sky radii to 0 to force sky
            centrod, db, x, y, ap, 0, 0, 0., xout, yout, counts
        end
    
        'cntrd': begin
            cntrd, db, x, y, ap/0.637, xout, yout
            aper, db, xout, yout, mags, errap, sky, skyerr, $
              1., [ap], [0,0], [0,0], /NAN, /EXACT, FLUX=1, /SILENT, $
              SETSKYVAL = sky
            counts = mags[0]
        end

        'psf_1_gauss': begin
            param = psf_1_gauss_initparam(db, x, y, ap)
            if debug then print, 'center : param init x,y = ', param[0:1]
            arg = psf_1_gauss_arg(db, x, y, ap)
            z =  db[arg.indx]
            zfit = curvefit(arg, z, sqrt( 1. > z ), $
              param, paramerr, chisq=chisq_dof, $
              /double, func = 'psf_1_gauss_func')
            if debug then print, 'center : param fit x,y = ', param[0:1]
            p = psf_1_gauss_param2struct(param)
            xout = p.xctr
            yout = p.yctr
            counts = p.total
            fwhm = p.fwhm
        end
    
        else: begin
            print, $
              'center: method must be one of '+$
              'boxm, centrod, cntrd, or psf_1_gauss'
            xout = xin
            yout = yin
        end
    endcase

    if n_elements(counts) eq 0 then stop

    skyval = sky

end