; docformat = 'rst'
;
; NAME:
;   MEANWITHUNC
; PURPOSE:
;   Compute the average of a serie of values, each associated to an uncertainty
;
;+
; :Description:
;   Compute the average of a serie of values, each associated to an uncertainty
;
; :Categories:
;    Maths, Moment
;    
; :Params:
;   VAL: in, required, type=float
;        An array of values for which the average is sought
;   UNC: in, required, type=float
;        An array of the uncertainties for each element in VAL
;
; :Returns: The weighted average of the VALUES, with its uncertainty
;
; :Examples:
;    Compute the average of 5 measurements, weighted by their uncertainties::
;      IDL> measures= [ 1.8,2.0,2.2,1.9,2.1]
;      IDL> errors  = [ 0.1,0.2,0.1,0.3,0.3]
;      IDL> print, meanWithUnc( measures, errors )
;
; :Author:
;    B.Carry (OCA)
;
; :History:
;    Change History::
;      Original Version written in June 2013, B. Carry (IMCCE)
;      2018-11. - B. Carry (OCA) - Set uncertainty to Bevington formulae
;-
function meanWithUnc, val, unc
  compile_opt idl2

 ;--I-- Bullet-proof checks
  nbVAL = n_elements(val)
  nbUNC = n_elements(unc)
  if nbVAL ne nbUNC then return, -1


 ;--II-- Singular array
  if nbVAL eq 1 then return, [val,unc]

 ;--III-- Compute the average value -- Uncertainties are supposed uncorrelated
  avg = total( val/unc^2, /NAN ) / total( 1./unc^2, /NAN )
    tVAL = val/avg
    tUNC = unc/avg

 ;--IV-- Evaluate the uncertainty on the resulting mean value
  mvar= sqrt( nbVal/ total( 1./unc^2, /NAN ) )
  var = mvar*sqrt( total( (val-avg)*(val-avg)/unc^2, /NAN )/(nbVal-1.) )
  if var eq 0 then var = mean( unc, /NAN )

 ;--V-- Return the mean and its standard deviation
  return, [avg, var]

end