;***************************************************************************
;
;*PROCDURE NAME:
;   Gausspsf
;
;*PURPOSE:
;   Return a point spread function having Gaussian profiles,
;   as either a 1D vector, a 2D image, or 3D volumetric-data.
;
;*CALLING SEQUENCE:
;       GAUSSPSF,NPIX,GARRAY,params=gpar, ndim=ndim, fwhm=fwhm,  $
;        centroid=cntrd, stdev=stdev, xy_correl=xy_corr, normalize=normalize
; or:
;   GAUSSPSF,NPIX,GARRAY,params=gpar
;
;*PARAMETERS:
;   npix   (REQ) (I) (012) (BILF)
;             number of pixels for each dimension, specify as an array,
;         or just one number to make all sizes equal.
;
;       garray  (REQ) (O) (123) (FD)
;             output Gaussian profile as either a 1-, 2- or 3-dimensional
;             array
;
;   params   (OPT) (KEY) (2) (BILF)
;             a 3 x ndim array giving for each dimension:
;             [ maxval, center, stdev ],  overrides other keywords.
;
;   ndim (OPT) (KEY) (0) (BIL)
;             dimension of result: 1 (vector), 2 (image), or 3 (volume),
;         default = 2 (an image result).
;
;   fwhm  (OPT) (KEY) (0) (BIL)
;             the desired Full-Width Half-Max (pixels) in each dimension,
;         specify as an array, or single number to make all the same.
;             (Must be specified if params or stdev is not included.)
;
;   centroid (OPT) (KEY) (0) (BIL)
;              pixels numbers of PSF maximum ( 0.5 is center of a pixel ),
;          default is exact center of requested vector/image/volume.
;
;   stdev (OPT) (KEY) (0) (BIL)
;              optional way to specify width by standard deviation param.
;              Must be specified if params or fwhm is not specified.
;
;   xy_correl (OPT) (KEY) (0) (BIL)
;               scalar between 0 and 1 specifying correlation coefficient
;     Use this keyword, for example, to specify an elliptical
;     gaussian oriented at an angle to the X,Y axis
;
;   normalize (OPT) (KEY) (0) (BIL)
;               causes resulting PSF to be normalized so Total( psf ) = 1.
;
;***************************************************************************
pro gausspsf,npix,garray,params=gpar, ndim=ndim, fwhm=fwhm,  $
    centroid=cntrd, stdev=stdev, xy_correl=xy_corr, normalize=normalize
;
 On_error,2
 npar = n_params()
 if (npar eq 0) then begin
    print,' PRO GAUSSPSF,NPIX,GARRAY,params=gpar, ndim=ndim, '
    print,'  fwhm=fwhm, centroid=cntrd, stdev=stdev, xy_correl=xy_corr, '
    print,'  normalize=normalize'
    retall
 endif

 if keyword_set(gpar) then begin
     sp = size(gpar)
     ndim = sp(0)
     if (ndim eq 2) then ndim = sp(2)
     factor = total( gpar(0,*) )/float( ndim )
     cntrd = gpar(1,*)
     stdev = gpar(2,*)
 endif
;
 if (n_elements(ndim) ne 1) then ndim = 2
 ndim = ndim > 1
;
 if (n_elements(npix) le 0) then begin
    print,'number of pixels not specified'
    retall
 endif else if (n_elements(npix) lt ndim) then npix = replicate(npix(0),ndim)
;
 if (n_elements(cntrd) lt ndim) and (n_elements(cntrd) gt 0) then $
    cntrd = replicate(cntrd(0), ndim)
;
 if (n_elements(cntrd) le 0) then cntrd=(npix-1)/2. else cntrd=cntrd-0.5
 if (n_elements(fwhm) gt 0) then stdev = fwhm/( 2* sqrt( 2* aLog(2) ) )
 if (n_elements(stdev) le 0) then begin
    message,"must specify stdev or fwhm=",/INFO
    retall
 endif
 if (n_elements(stdev) lt ndim) then stdev = replicate( stdev(0), ndim )
;
case ndim of

    1: begin
       x = findgen(npix(0)) - cntrd(0)
                gauss,x,0,stdev(0),1,psf
;
           end

    2: begin
       psf = make_array( DIM=npix(0:ndim-1), /FLOAT )
       x = findgen( npix(0) ) - cntrd(0)
       y = findgen( npix(1) ) - cntrd(1)
                sigfac = 1 / (2. * stdev * stdev )
       if n_elements(xy_corr) eq 1 then begin
         y2 = sigfac(1) * y * y
         x1 = sigfac(0) * x
         yc = y * ( xy_corr/(stdev(0)*stdev(1)) )
         for j=0,npix(1)-1 do begin
          zz = x * (yc(j) + x1) + y2(j)
          w = where( zz LT 86, nw )
          if (nw GT 0) then psf(w,j) = exp( -zz(w) )
               endfor
                endif else begin
                        gauss,x,0,stdev(0),1,psfx
                        gauss,y,0,stdev(1),1,psfy

         for j=0,npix(1)-1 do psf(0,j) = psfx * psfy(j)
                endelse
           end

    3: begin
       psf = make_array( DIM=npix(0:ndim-1), /FLOAT )
       x = findgen( npix(0) ) - cntrd(0)
       y = findgen( npix(1) ) - cntrd(1)
       z = findgen( npix(2) ) - cntrd(2)
                gauss,x,0,stdev(0),1,psfx
                gauss,y,0,stdev(1),1,psfy
                gauss,z,0,stdev(2),1,psfz

       for k=0,npix(2)-1 do begin
           for j=0,npix(1)-1 do psf(0,j,k) = psfx * psfy(j) * psfz(k)
        endfor
           end
 endcase
;
 if keyword_set(normalize) then garray = psf/total(psf)
 if (n_elements(factor) eq 1) then begin
    if (factor ne 1) then garray = factor*psf else garray = psf
 endif else garray = psf
;
 return
 end    ; gausspsf