function wmean, x, xerr, meanerr, exterr=exterr

;+
; NAME:
;  wmean
; PURPOSE: (one line)
;  weighted mean.
; DESCRIPTION:
;  Return weighted mean and the error of the mean.
; CATEGORY:
;  Math
; CALLING SEQUENCE:
;  meanx = wmean(x, xerr, meanerr)
; INPUTS:
;  x - array of values
;  xerr - array of errors for x
; OPTIONAL INPUT PARAMETERS:
;  none
; KEYWORD INPUT PARAMETERS:
;  exterr - calculate the "external error" using scatter
; KEYWORD OUTPUT PARAMETERS:
;  none
; OUTPUTS:
;  return weighted mean.
;  manxerr - error in the mean
; COMMON BLOCKS:
;  None
; SIDE EFFECTS:
; RESTRICTIONS:
;  None
; PROCEDURE:
; MODIFICATION HISTORY:
;  Added to layoung/math 2004 Feb, by Leslie Young, SwRI
;  fixed typo 2005 Aug 4 LAY
;  2009 Oct 20 LAY - mult. residual by w (weight), to wtot
;-

  ; select on positive errors to avoid division by 0
    eps = (machar()).xmin        ; smallest floating point
    rooteps = sqrt(eps)
    gd = where(xerr gt rooteps, ngd)
  
  ; if all points have funny errors, 
  ; return straight mean and calculate error from scatter
    if ngd eq 0 then begin
        xm = mean(x)
        n = n_elements(x)
        if n lt 2 then begin
            meanerr = reform(xerr[0])
        endif else begin
            meanerr = sqrt( total((x-xm)^2.)/(n*(n-1)) )
        endelse
        return, xm
    end
    w = 1./xerr[gd]^2
    wtot = total(w)
    xm = total(w*x[gd])/wtot
    meanerr = 1/sqrt(wtot)

    if keyword_set(exterr) then begin
;        chisq = total( (x[gd]-xm)^2 * wtot[gd] )
        chisq = total( (x[gd]-xm)^2 * w[gd] )
        meanerr = meanerr * sqrt(chisq/(ngd-1))
    end

    return, xm
end