;+
; NAME:
;       FORCEFIT_SIGCLIP
; PURPOSE:
;       Force a fit in either slope or zero-point of a line where the
;       points lying outside a certain sigma are not included in the
;       fit. This occurs iteratively until no points are rejected anymore.
; CALLING SEQUENCE:
;       result = FORCEFIT_SIGCLIP([x, y[,a=a][,b=b][,niter=niter][,sigma=sigma][,verb=verb][,/fixed]][/help])
; EXAMPLES:
;       result = FORCEFIT_SIGCLIP(x, y, a=1, sigma=3)
;           fits a line of slope 1 to the input values where the
;           cut-off is at 3 sigma
; KEYWORDS:
;       x = the input x values
;       y = the input y values
;       a = the fixed slope
;       b = the fixed zero-point
;       niter = the number of iterations to reach a stable fit
;       sigma = the cut-off boundaries for points out to which points
;       are included in the fit
;       verb = verbose on when non-zero, omitting means no verbose
;       /fixed = omitting this keyword does a normal LINFIT at the
;       end, including it forced the fit to use the input a or b
;       /help = displays this help
; MODIFICATION HISTORY:
;       before 27/08/04 written by Eduard Westra
;       27/08/04 added help header by Eduard Westra
;       14/12/04 if both a and b are not set, then just return the
;       forcefit result
;-
FUNCTION forcefit_sigclip, x, y, a=a, b=b, niter=i, sigma=sigma, verb=verb, fixed=fixed, help=help
IF (KEYWORD_SET(HELP)) THEN BEGIN & doc_library,'forcefit_sigclip' & RETURN, [0.,0.] & ENDIF
    IF NOT KEYWORD_SET(sigma) THEN sigma = 3.
    IF NOT KEYWORD_SET(verb) THEN verb = 0
    IF NOT KEYWORD_SET(fixed) THEN fixed = 0
    IF (KEYWORD_SET(a) AND KEYWORD_SET(b)) THEN BEGIN
        PRINT, "You already set both parameters. What do you want me to do?!?"
        RETURN, [b,a]
    ENDIF
    IF KEYWORD_SET(a) THEN BEGIN
        ox = x
        oy = y
        nprev = N_ELEMENTS(ox)
        i = 0
        WHILE (1 EQ 1) DO BEGIN
            result = FORCEFIT(ox, oy, a=a)
            a      = result[1]
            b      = result[0]
            c      = oy + ox/a
            XX     = (c-b)/(a+(1/a))
            YY     = a*XX + b
            dist   = SQRT((ox - XX)^2 + (oy - YY)^2)

            idx    = WHERE( (dist - MEAN(dist))/STDDEV(dist) LE sigma )

            IF (verb) THEN BEGIN
                print, i, ": ", N_ELEMENTS(idx), " out of ", N_ELEMENTS(dist), " are within ", TRIM(sigma)," sigma."
                print, i, ": ", nprev/DOUBLE(N_ELEMENTS(x))
            ENDIF

            ix = ox[idx]
            iy = oy[idx]

            IF (fixed) THEN $
              result = FORCEFIT(ix, iy, a=a) $
            ELSE $
              result = FORCEFIT(ix, iy)

            IF (verb) THEN BEGIN
                print, i, ": ", result
            ENDIF

            IF (nprev EQ N_ELEMENTS(idx)) THEN BREAK
            IF (nprev/DOUBLE(N_ELEMENTS(x)) LE 0.10) THEN BEGIN
                print, "Failed to converge to a solution. Less than 10% of the points would be used."
                print, "The result is not worth trusting."
                RETURN, result
            ENDIF
            nprev = N_ELEMENTS(idx)
            ox = ix
            oy = iy
            i = i + 1
        ENDWHILE

        RETURN, result
    ENDIF
    IF KEYWORD_SET(b) THEN BEGIN
        ; implement later
    ENDIF
    IF NOT (KEYWORD_SET(a) AND KEYWORD_SET(b)) THEN BEGIN
        RETURN, FORCEFIT(x, y)
    ENDIF
END