;+
; NAME:
;  roboscale
; PURPOSE: (one line)
;  Robust scaling of a template to match an array.
; DESCRIPTION:
;   Robust fitting of y = scaleFac * x
; CATEGORY:
;  Statistics
; CALLING SEQUENCE:
;  scaleFac = roboscale(x, y, sigma, thresh, bplist, )
; INPUTS:
;  x - the template
;  y - the data array, same length as x
;  sigma - estimated noise for each pixel in y
;  thesh - threshold for outliers 
;  bplist - optional bad pixel list (-1 for all OK)
; OPTIONAL INPUT PARAMETERS:
; INPUT KEYWORD PARAMETERS:
; OUTPUT KEYWORD PARAMETERS:
;  gdindex - named variable that will contain an array of indices used
; OUTPUTS:
;  returns scaleFac such that y ~= scaleFac * x
; COMMON BLOCKS:
;  None.
; SIDE EFFECTS:
;  None.
; RESTRICTIONS:
;  None.
; PROCEDURE:
;   First iteration finds minimum absolute deviation.
;   Outliers are identified after the first iteration.
;   The second iteration minimizes least squares 
;
; MODIFICATION HISTORY:
;  Written July 22, 2002 Eliot Young, SwRI
;  Rewritten by Leslie A. Young, SwRI, 2002 Oct 29.
;  - first iteration minimizes absolute deviation
;  - exchanged order of x and y in argument list
;  - only two iterations
;  - optionally returns the error and the list of good pixel indices
;-

function roboscale, x, y, sigma, thresh, bplist, GDINDEX = gdindex, scaleErr = scaleErr

  n = n_elements(x)
  scaleErr = 0.

; Get the weights
  w = fltarr(n)
  gd = where(sigma ge 0)
  w[gd] = 1/sigma[gd]^2.
  if n_params() gt 4 then begin
    if n_elements(bplist) gt 0 then begin
      if bplist[0] ne -1 then w[bplist] = 0.
    end 
  end

; First iteration: minimize weighted absolute deviation
  scaleFac = madscalew(x, y, w)
  yfit = scaleFac * x
  
; Find good points
  gdindex = where( (abs(y-yfit) lt sigma * thresh) and (w ne 0), ngd)
  if ngd eq 0 then begin
    print, 'ROBOSCALE: All points rejected in first iteration.'
    return, 0
  end
  
; Second iteration: minimize sum of squares for good points
  swxy = total(w[gdindex]*x[gdindex]*y[gdindex])
  swxx = total(w[gdindex]*x[gdindex]^2)
  if swxx eq 0 then begin
    print, 'ROBOSCALE: all xs are zero, scale factor is undefined'
    return, 0 
  end
  scaleFac = swxy / swxx
  scaleErr = sqrt(1./swxx)
  return, scaleFac
end

pro roboscaleTEST

  x = (findgen(21)-10)*0.2
  f = exp(-x^2)
  a = 1000.
  y = a * f
  ysig = sqrt(y)
  yerr = randomn(0L, 21) * ysig
  ydat=y+yerr
  for i = 0, 4 do begin
    ibad = floor(randomu(0L)*20.99999)
    ydat[ibad] = ydat[ibad]+2*a*randomu(0L)
  end
  thresh = 5.
  bplist = -1
  roboa = roboscale(f, ydat, ysig, thresh, bplist, GDINDEX = gdindex, scaleErr=roboe)
  print, a, roboa,"/-",roboe
  plot, x, y, xr=[-2.1,2.1], yr=[-0.1,1.2]*max(ydat), /xs,/ys
  oplot, x, ydat, ps=4
  oplot, x[gdindex], ydat[gdindex], ps=1
  oplot, x, roboa*f, line=2, thick=2
end