;+ Fits data in data container specified by index and fits multiple gaussians
;  to regions in this data, according to the values set in !g.gauss
;
; @keyword index {in}{optional}{type=long}{default=0} index of data container in !g whose data is to be fit.
; @keyword modelindex {in}{optional}{type=long}{default=-1} index of
; data container in !g to hold the resulting model.  If not set (the
; default) then no index is filled with the model.  
; @keyword ok {out}{optional}{type=long} result; 1 - good, 0 - bad
; @keyword quiet {in}{optional}{type=boolean}{default=F} When set, the
; report of the fits and initial guesses is not made.
; @param fit {out}{optional}{type=array} results of fit. identical to !g.gauss.fit
; @param fitrms {out}{optional}{type=array} errors of fit. identical to !g.gauss.fitrms
;
; @uses gauss_fits
; @uses report_gauss
;
; @examples
; <pre>
; ; get some continuum data
; cont
; filein,'peaks.fits'
; get,index=1
; ; setup the gauss fit
; gregion,[60,80]
; ngauss,[1]
; gmaxiter,500
; gparamvalues, 0, [400000.,70.,100.]
; ; find and show the fit
; gauss
; gshow
; </pre>
;-
pro gauss, index=index, modelindex=modelindex, ok=ok, quiet=quiet, fit, fitrms
    compile_opt idl2

    ; init output
    fit = -1
    fitrms = -1

    ; argument checks
    ok = 0
    if (n_elements(index) eq 0) then index=0

    maxindex = !g.line ? n_elements(!g.s) : n_elements(!g.c)

    if  (index lt 0 or index gt maxindex) then begin
        message, 'requested index does not exist',/info
        return
    endif

    npts = !g.line ? data_valid(!g.s[index]) : data_valid(!g.c[index])
    if (npts lt 1) then begin
        message, 'no data at requested index', /info
        return
    endif

    if (n_elements(modelindex) eq 0) then modelindex=-1

    if (modelindex gt maxindex) then begin
        message, 'requested model index does not exist', /info
        return
    endif

    if (check_gauss_settings() eq 0) then message, "cannot fit gaussians with guide structure in bad state"
    
    ; look at only those regions specified
    nregions = !g.gauss.nregion
    if (nregions le 0) then return

    regions = !g.gauss.regions[*,0:(nregions-1)]
    ngauss = !g.gauss.ngauss

    if (ngauss le 0) then return
    
    ; prepare for fitting data
    allChans = dindgen(npts)

    data = !g.line ? *!g.s[index].data_ptr : *!g.c[index].data_ptr

    ; convert the parameter values to a simpler array[3,total(ngauss)]
    inits = dblarr(3,ngauss)
    for i = 0, (ngauss-1) do begin
        inits[0,i] = !g.gauss.params[0,i].value
        inits[1,i] = !g.gauss.params[1,i].value
        inits[2,i] = !g.gauss.params[2,i].value
    endfor

    ; dont let users forget to set the maximum number of iterations. use 500 if they have
    maxiter = (!g.gauss.maxiter eq 0) ? 500 : !g.gauss.maxiter
    
    ; fit all gaussians using all regions, each fits model is returned in an appended array. params are not used yet
    yfits = gauss_fits(allChans,data,nregions,regions,inits,ngauss,maxiter,coefficients,errors,quiet=1)
  
    ; add up yfits and write results to !g
    model=make_gauss_data(allChans,coefficients,0.0)
    ; add these results to !g

    !g.gauss.fit[*,0:(ngauss-1)] = coefficients
    !g.gauss.fitrms[*,0:(ngauss-1)] = errors

    ; if specified, then copy over the model data to it's location
    if (modelindex ge 0) then begin
        copy, 0, modelindex
        setdata,model,index=modelindex
    endif

    ; print out the results
    if (not keyword_set(quiet)) then report_gauss

    ; optional output
    fit = coefficients
    fitrms = errors

    ok = 1

end