PRO	NGAUSS,X,A,F,PDER
;+
; NAME:
;	NGAUSS
;
; PURPOSE:
;	Evaluate the sum of a Gaussian and a 2nd-order polynomial
;	and optionally return the value of its partial derivatives.
;	Normally, this function is used by CURVEFIT to fit the
;	sum of a line and a varying background to actual data.
;
; CALLING SEQUENCE:
;	NGAUSS, X, A, F [, Pder]
;
; INPUTS:
;	X:	The values of the independent variable.
;	A:	The parameters of the equation described in PROCEDURE below.
;
; OUTPUTS:
;	F:	The value of the function at each X(i).
;
; OPTIONAL OUTPUT PARAMETERS:
;	Pder:	An array of the size (N_ELEMENTS(X),6) that contains the
;		partial derivatives.  Pder(i,j) represents the derivative
;		at the i'th point with respect to j'th parameter.
;
; COMMON BLOCKS:
;        None.
;
; SIDE EFFECTS:
;	None.
;
; RESTRICTIONS:
;	None.
;
; PROCEDURE:
;
;       The function is the sum of N gaussians plus a constant, linear, 
;       or quadratic background.
;
;       Terms 0, 1, and 2 are the amplitude, center, and FWHM of
;       the first gaussian.
;
;       G1 = A(0)*EXP(-C*(X-A(1))^2 / A(2)^2)
;          Where C = 4*ALOG(2.) makes A(2) a FWHM
;
;       Subsequent gaussians are specified with A(3:5), A(6:8), etc.
;
;       The last 1, 2, or 3 terms define a background.  If one
;       term is left, the background is a constant A(3N+0).
;       If two terms left, it is linear = A(3N+0) + A(3N+1)*X
;       If 3 left, quadratic = A(3N+0) + A(3N+1)*X + A(3N+2)*X^2
;
;-
ON_ERROR,2                        ;Return to caller if an error occurs

npts=n_elements(x)
nparms=n_elements(a)
ngauss=(nparms-1)/3
nbkg=(nparms-3*ngauss)

c = 4.*alog(2.)

; Evaluate the function, starting with the exponential parts

ez_a = fltarr(npts, ngauss)
z_a = fltarr(npts,ngauss)

f=replicate(0.,npts)
for g=0,ngauss-1 do begin

    if (a(3*g+2) ne 0.) then $
	z = (x-a(3*g+1)) / a(3*g+2) $
        else z = replicate(10.,npts)
    z_a(*,g)=z

    ; Ignore small terms
    ez = exp(-c*(z^2)) * (abs(z) le 7.)
    ez_a(*,g)=ez

    f = f + a(3*g+0)*ez

    endfor

;       Evaluate the background

bkg = a(3*ngauss+0)
if nbkg gt 1 then bkg = bkg + a(3*ngauss+1)*x
if nbkg gt 2 then bkg = bkg + a(3*ngauss+2)*x^2

f = f+bkg

; Done with the function.  Return if partial not needed.

if n_params(0) le 3 then return 

; Compute array of partial derivatives

pder = fltarr(npts,nparms) 

for g=0,ngauss-1 do begin

    ez=reform(ez_a(*,g))
    z=reform(z_a(*,g))

    pder(0,3*g+0) = ez

    if a(3*g+2) ne 0. then $
    pder(0,3*g+1) = a(3*g+0) * ez * (2.*c/a(3*g+2))*z

    pder(0,3*g+2) = pder(*,3*g+1) * z

    endfor

pder(*,3*ngauss+0) = 1.
if nbkg gt 1 then pder(0,3*ngauss+1) = x
if nbkg gt 2 then pder(0,3*ngauss+2) = x^2

return
end