function getrms, fluxes, errors

;; Keith Horne's method to estimate intrinsic variability.
;; This is a modified version of getrms initially written by Jon Trump
;; V = sum( (xi-u)^2 gi^2 ) / sum(gi)
;;   V = intrinsic variance  -- Eqn (8) of sample char paper
;;   sig0 = sqrt(V)
;;   xi = flux measurement
;;   u = mean flux with optimal estimator -- Eqn (9) of sample char paper; 
;;     note there was a typo on the error of u in the original Eqn (9) of the paper, missing a square sign
;;   gi = 1 / (1+e^2/V)
;;   e = flux error
;;
;; Iterate to solve for V and u.  Note that e should be a combination of
;; the standard flux error plus the spectrophotometry error.
;; It returns: [u, uerror, sig0, sig0err]

  if stddev(fluxes) eq 0 or n_elements(fluxes) le 1 then return,[0,-1]

  uguess = median(fluxes)
;;   errors = sqrt(errors^2 + (0.04*meanflux)^2)  ; add 0.04% specphoterror

  varguess = (stddev(fluxes)^2 - median(errors)^2)>0.05

  repeat begin
     varold = varguess
     uold = uguess
     gi = 1. / (1+errors^2/varold)
     varnew = total( (fluxes-uold)^2 * gi^2 ) / total(gi)
     unew = total( fluxes*gi ) / total(gi)
     varguess = varnew
     uguess = unew
  endrep until abs(varold-varnew)/varnew le 1e-5 or varnew lt 1e-30

  gi = 1. / (1+errors^2/varnew)
  varerror = varnew / $
     sqrt((total(gi)*total((fluxes-unew)^2*gi^3)/total((fluxes-unew)^2*gi^2) - total(gi^2)/2.))
  sig0 = sqrt(varnew>0)
  sig0error = 0.5*sqrt(varnew>0) / $
     sqrt((total(gi)*total((fluxes-unew)^2*gi^3)/total((fluxes-unew)^2*gi^2) - total(gi^2)/2.))
  uerror = (varnew>0)/total(gi) ; there was a typo in the second-half of Eqn (9) in the sample char paper, i.e., missing a square sign. 
  uerror = sqrt(uerror)

  return, [unew, uerror, sig0, sig0error]

end