function gaussian2d, xi1, xi2, parms, pderiv, DOUBLE=double ;+ ; NAME: ; GAUSSIAN2d ; PURPOSE: ; Compute the 2-d Gaussian function and optionally the derivative ; ; based on Goddard IDL Astro's "Gaussian.pro" ; ; RESTRICTION: Right now only circularly symmetric gaussians ; ; EXPLANATION: ; Compute the 2-D Gaussian function and optionally the derivative ; at an array of points. ; ; CALLING SEQUENCE: ; y = gaussian( x, y, parms,[ pderiv ]) ; ; INPUTS: ; x,y = arrays, 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 in X, ; parms[2] = mean value (center) of Gaussian in Y, ; parms[3] = 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[4] = 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: ; gaussian ; why is this here?? ; HISTORY: ; 2005-May-06 Forked from gaussian.pro by Marshall Perrin ; ;- 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 ) npts1 = N_elements(xi1) npts2 = N_elements(xi2) ptype = size(parms,/type) if (ptype LE 3) or (ptype GE 12) then parms = float(parms) if (Nparmg GE 4) then sigma = parms[3] double = keyword_set(DOUBLE) if double then \$ ;Double precision? gauss = dblarr( npts1,npts2 ) else \$ gauss = fltarr( npts1,npts2 ) xx = rebin(xi1-parms[1],npts1,npts2) yy = rebin(reform(xi2-parms[2],1,npts2),npts1,npts2) xi = sqrt(xx^2+yy^2) z = ( xi )/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 ) ; TODO make the following work in 2D if N_params() GE 4 then begin if double then \$ pderiv = dblarr( npts1, npts2, Nparmg ) else \$ pderiv = fltarr( npts1, npts2, 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 5 then return, parms[0] * gauss else \$ return, parms[0] * gauss + parms[4] end