function wmean, x, xerr, meanerr

;+
; 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:
;  none
; 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
;-

  ; select on positive errors to avoid division by 0
  eps = (machar()).eps ; smallest floating point
  rooteps = sqrt(eps)
  gd = where(xerr gt root, ngd)
  
  ; if all points have funny errors, 
  ; return straight mean and calculate error from scatter
  if ng eq 0 then begin
    xm = mean(x)
    n = n_elements(x)
    if n eq 0 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)
  return, xm
end