function gaussian, xi, parms, pderiv, DOUBLE=double;+; NAME:;       GAUSSIAN; PURPOSE:;       Compute the 1-d Gaussian function and optionally the derivative; EXPLANATION:;       Compute the 1-D Gaussian function and optionally the derivative ;       at an array of points.;; CALLING SEQUENCE:;       y = gaussian( xi, parms,[ pderiv ]);; INPUTS:;       xi = array, independent variable of Gaussian function.;;       parms = parameters of Gaussian, 2, 3 or 4 element array:;               parms[0] = maximum value (factor) of Gaussian,;               parms[1] = mean value (center) of Gaussian,;               parms[2] = standard deviation (sigma) of Gaussian.;               (if parms has only 2 elements then sigma taken from previous;               call to gaussian(), which is stored in a  common block).;               parms[3] = optional, constant offset added to Gaussian.; OUTPUT:;       y -  Function returns array of Gaussian evaluated at xi.    Values will;            be floating pt. (even if xi is double) unless the /DOUBLE keyword;            is set.;; OPTIONAL INPUT:;       /DOUBLE - set this keyword to return double precision for both;             the function values and (optionally) the partial derivatives.; OPTIONAL OUTPUT:;       pderiv = [N,3] or [N,4] output array of partial derivatives,;               computed only if parameter is present in call.;;               pderiv[*,i] = partial derivative at all xi absisca values;               with respect to parms[i], i=0,1,2,[3].;;; EXAMPLE:;       Evaulate a Gaussian centered at x=0, with sigma=1, and a peak value;       of 10 at the points 0.5 and 1.5.   Also compute the derivative;;       IDL> f = gaussian( [0.5,1.5], [10,0,1], DERIV );       ==> f= [8.825,3.25].   DERIV will be a 2 x 3 array containing the;       numerical derivative at the two points with respect to the 3 parameters.; ; COMMON BLOCKS:;       None; HISTORY:;       Written, Frank Varosi NASA/GSFC 1992.;       Converted to IDL V5.0   W. Landsman   September 1997;       Use machar() for machine precision, added /DOUBLE keyword,;       add optional constant 4th parameter    W. Landsman   November 2001;-  On_error,2  common gaussian, sigma  if N_params() LT 2 then begin        print,'Syntax - y = GAUSSIAN( xi, parms,[ pderiv, /DOUBLE ])'        print,'         parms[0] = maximum value (factor) of Gaussian'        print,'         parms[1] = mean value (center) of Gaussian'        print,'         parms[2] = standard deviation (sigma) of Gaussian'        print,'         parms[3] = optional constant to be added to Gaussian'        return, -1  endif  common gaussian, sigma        Nparmg = N_elements( parms )        npts = N_elements(xi)         ptype = size(parms,/type)        if (ptype LE 3) or (ptype GE 12) then parms = float(parms)        if (Nparmg GE 3) then sigma = parms[2]        double = keyword_set(DOUBLE)        if double then $       ;Double precision?            gauss = dblarr( npts ) else $            gauss = fltarr( npts )         z = ( xi - parms[1] )/sigma        zz = z*z; Get smallest value expressible on computer.   Set lower values to 0 to avoid; floating underflow        minexp = alog((machar(DOUBLE=double)).xmin)              w = where( zz LT -2*minexp, nw )        if (nw GT 0) then gauss[w] = exp( -zz[w] / 2 )        if N_params() GE 3 then begin                if double then $                 pderiv = dblarr( npts, Nparmg ) else $                pderiv = fltarr( npts, Nparmg )                fsig = parms[0] / sigma                pderiv[0,0] = gauss                pderiv[0,1] = gauss * z * fsig                if (Nparmg GE 3) then  pderiv[0,2] = gauss * zz * fsig                if (Nparmg GE 4) then  pderiv[0,3] = replicate(1, npts)           endif if Nparmg LT 4 then return, parms[0] * gauss else $                     return, parms[0] * gauss + parms[3] end