;+
;FUNCTION:  MGAUSS
;PURPOSE:
;  Evaluates multiple gaussian function with background.
;  This function may be used with the "fit" curve fitting procedure.
;
;KEYWORDS:
;  PARAMETERS: a structure that contain the parameters that define the gaussians
;     If this parameter is not a structure then it will be created.
;
;USAGE:
;  par = mgauss( [numg =4] ) ; get parameter structure with 4 gaussian peaks
;  x = findgen(100)/10.
;  y = mgauss(x,param = par)
;  plot,x,y
;RETURNS:
;  result
;
;Written by: Davin Larson
;
; $LastChangedBy: davin-win $
; $LastChangedDate: 2011-03-22 12:29:06 -0700 (Tue, 22 Mar 2011) $
; $LastChangedRevision: 8533 $
; $URL: svn+ssh://thmsvn@ambrosia.ssl.berkeley.edu/repos/ssl_general/trunk/tools/fitting/mgauss.pro $
;-

function mgauss, x,a, parameters=p,numg=numg

if not keyword_set(p) then begin
   if not keyword_set(numg) then numg=1
   g1 = {name:'',a:1.d, x0:0.d, s:1.0d}
   g = replicate(g1,numg)
   g.x0 = findgen(numg)
   p = {func:"mgauss",binsize:0.,quantize:1,shift:1,bkg:0d,g:g}
endif

if n_params() eq 0 then return,p
if n_params() eq 2 then begin                    ;; kludge to make routine work with MPFIT
  xfer_parameters,p,'',a,/array_to_struct
endif

n = n_elements(p.g)
xx = x

dx = p.binsize
if dx eq 0 then begin
  if size(/type, p.bkg) eq 8 then f = func(xx,param=p.bkg) else  f= p.bkg
  for i=0,n-1 do begin
    pg = p.g[i]
    z = (x - pg.x0)/pg.s
    e = exp(- z^2/2 )
    height = pg.a / pg.s / sqrt(2*!dpi)
    f +=  height * e
  endfor
endif else begin
  if p.quantize ne 0 then begin
     if p.shift then xx = p.binsize * (floor(xx/p.binsize)+.5d) $
     else xx = p.binsize * round(xx/p.binsize)
  endif
  if size(/type, p.bkg) eq 8 then f = func(xx,param=p.bkg) else  f= p.bkg
  x1= xx - dx/2d
  x2= xx + dx/2d
  s2 = sqrt(2)
  for i=0,n-1 do begin
    pg = p.g[i]
    z1 = (x1-pg.x0)/pg.s/s2
    z2 = (x2-pg.x0)/pg.s/s2
    f += pg.a/2 * (erfc(z1) - erfc(z2)) / dx
  endfor
endelse
;dprint,dlevel=5,[p.bkg,p.g.a,p.g.x0,p.g.s]
return,f
end



;
;function random_mgauss,param=p,seed
;
;numg = n_elements(p.g)
;for i=0,numg-1 do begin
;  r = randomn(seed,p.g[i].a)*p.g[i].s +p.g[i].x0
;  if i eq 0 then rr = r else rr= [rr,r]
;endfor
;return,rr
;end