;+
; NAME:
;          JJGAUSSFIT
;
;
; PURPOSE:
;         Fits a gaussian to a 1D distribution of points, This is a wrapper
;         for Craig Markwardt's infinitely useful fitting package MPFIT.
;
;         Unlike IDL's GAUSSFIT, this function optionally accounts
;         for measurement errors 
;
;         Fits the function:
;
;         f(x) = a0 * exp(-z^2/2) + a3
;
;         where
;
;         z = (x - a1)/a2
;
; CATEGORY:
;
;         Getting the job done (TCB).
;
; CALLING SEQUENCE:
;
;         fit = jjgaussfit(x, y [, param])
;
; INPUTS:
;         x:  independant variable
;         y:  independant variable
;
; KEYWORDS:
;
;         FIXAMP=, FIXCEN=,  FIXWID=, FIXBG=: Scalar value at which the 
;         centroid, amplitude, width or background (offset),
;         respectively, should be held fixed during the fit.
;
;         LIMITAMP=, LIMITCEN=, LIMITWID=, LIMITBG=: 2-element array
;         specifying lower and upper limits of the centroid, amplitude,
;         width or background (offset), respectively. 
;
;         MOVIE: Play a movie during the fit procedure. 
;
;         GUESS: 4-element array containing the parameter guesses. If
;         guesses are not set manually, then parameter guesses are
;         determined automatically. 
;    
;         GUESS[0] = Amplitude
;         GUESS[1] = Centroid 
;         GUESS[2] = Sigma width
;         GUESS[3] = Background offset
;
;         YERR: 1 sigma errors associated with each Y
;
; OUTPUTS:
;
;         FIT: The best-fit Gaussian evaluated at each X
;         PARAM: 4 element vector containing the fit parameters
;                param = [a0, a1, a2, a3]
;
; EXAMPLE:
;
; x = findgen(121) * 0.25 - 15
; y = 1d2 * exp(-x^2) + 2
; yn = 1d2 - (y + randomn(seed, 121)*sqrt(y))
; p0 = [-110, 0.05, 1.1, 90]
; g = jjgaussfit(x, yn, a, yerr=sqrt(y), guess=p0, /movie)
;         
; MODIFICATION HISTORY:
; 12 Dec 2003 Written by JohnJohn
; 05 Dec 2004 JJ - fixed bug in FWHM guess
; 22 Feb 2005 JJ - fixed to work with negative amplitudes. Stole
;                  parameter guesses from GAUSSFIT.PRO
; 06 Mar 2005 JJ - added MOVIE keyword.
; 15 Mar 2005 JJ - modified parameter guess logic to w ork  with
;                  continuum-normalized absorption lines 
; 18 Mar 2005 JJ - added GUESS keyword
;-

function jjgaussfit_func, p, x=x, y=y, err=err, fit=fit, movie=movie
z = (x - p[1])/p[2]
fit = p[0]*exp(-z^2/2d) + p[3]
if n_elements(err) eq 0 then err = 1.
if keyword_set(movie) then begin
    yr = minmax(y)
    if abs(max(y)) gt abs(min(y)) then yr = yr * [0.95, 1.05] else $
      yr = yr * [1.05, 0.95]
    plot, x, y, yr=yr, /ys, ps=4
    oplot, x, fit
    wait,.5
endif
return, (y - fit)/err
end

pro jjgaussfit_plot, fcn, par, iter, fnorm, functargs=fa $
                  , parinfo=parifno, quiet=quiet, dof=dof
wset,30
resid = jjgaussfit_func(par, x=fa.x, y=fa.y, /movie)

mpfit_defiter, fcn, par, iter, fnorm, FUNCTARGS=fa, $
                   quiet=quiet, parinfo=parinfo, dof=dof

wset,0
device,copy = [0,0,!d.x_size,!d.y_size,0,0,30]
end


function jjgaussfit, x, y, p $
                     , yerr=yerr $
                     , fixamp=fixamp , limitamp=limitamp $
                     , fixcen=fixcen , limitcen=limitcen $
                     , fixsig=fixsig, limitsig=limitsig $
                     , fixbg=fixbg, limitbg=limitbg $
                     , loud=loud, movie=movie, guess=guess
ON_ERROR,2
if keyword_set(yerr) then begin
    fa = {x:double(x), y:double(y), err:double(yerr)}
endif else begin
    fa = {x:double(x), y:double(y)}
endelse
;;;Set up parameter guesses, stolen from GAUSSFIT.PRO and modified to
;;;work with the case of a continuum normalized absorption line
if 1-keyword_set(guess) and n_elements(guess) lt 4 then begin
    n = n_elements(y)
    ymed = median(y)
    ytest = y - ymed
    ymax = max(ytest, imax)
    ymin = min(ytest, imin)
    if abs(ymax) gt abs(ymin) then begin ;emiss or absorp?
        i0 = imax 
        amp = ymax - ymin
    endif else begin
        i0 = imin
        amp = ymin - ymax
    endelse
    bg = ymed

    i0 = i0 > 1 < (n-2)         ;never take edges
    del = amp/exp(1.)           ;1/e value
    i=0
    while ((i0+i+1) lt n) and $ ;guess at 1/2 width.
      ((i0-i) gt 0) and $
      (abs(ytest[i0+i]) gt abs(del)) and $
      (abs(ytest[i0-i]) gt abs(del)) do i=i+1

    p0 = [amp, x[i0], abs(x[i0]-x[i0+i]), bg]
endif else p0 = guess

npar = n_elements(p0)
parinfo = replicate( { fixed: 0b, $
                       limited: [0b,0b], $
                       limits: dblarr(2) $
                     }, npar)

if n_elements(fixamp) gt 0 then begin
    parinfo[0].fixed = 1b
    p0[0] = fixamp
endif

if n_elements(fixcen) gt 0 then begin
    parinfo[1].fixed = 1b
    p0[1] = fixcen
endif

if keyword_set(fixsig) then begin
    parinfo[2].fixed = 1b
    p0[2] = fixsig
endif

if n_elements(fixbg) gt 0 then begin
    parinfo[3].fixed = 1b
    p0[3] = fixbg
endif

if keyword_set(limitamp) then begin
    parinfo[0].limited = 1b
    parinfo[0].limits = limitamp
    p0[0] = mean(limitamp)
endif

if keyword_set(limitcen) then begin
    parinfo[1].limited = 1b
    parinfo[1].limits = limitcen
    p0[1] = mean(limitcen)
endif

if keyword_set(limitsig) then begin
    parinfo[2].limited = 1b
    parinfo[2].limits = limitsig
    p0[2] = mean(limitsig)
endif

if keyword_set(limitbg) then begin
    parinfo[3].limited = 1b
    parinfo[3].limits = limitbg
    p0[3] = mean(limitbg)
endif

if keyword_set(movie) then begin
    window, 0
    window, 30, /pixmap
    ip = 'jjgaussfit_plot'
endif
p = mpfit('jjgaussfit_func', p0, parinfo=parinfo, funct=fa, maxiter=20 $
          , quiet=1-keyword_set(loud), iterproc=ip)

z = (x - p[1])/p[2]
fit = p[0]*exp(-z^2/2d) + p[3]
return, fit
end