;***************************************************************************
;+
;*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.
;
;
;*EXAMPLE:
;	Create a 31 x 31 array containing a normalized centered gaussian 
;	with an X FWHM = 4.3 and a Y FWHM = 3.6
;
;	IDL> gausspsf,31,array,fwhm=[4.3,3.6],/normal
;  or
;       IDL> p = fltarr(3,2)
;       IDL> p(0) = [1,15.5,4.3,1,15.5,3.6]
;	IDL> gausspsf,31,array,params=p,/normal
;
;*SUBROUTINES USED:
;       gauss
;
;*NOTES:
;       - Either fwhm or stdev must be specified if params is
;       not included.
;       - Specifying xy_correl only affects profiles with ndim = 2.
;    
;
;*MODIFICATION HISTORY:
;       PSF_GAUSSIAN Written by Frank Varosi NASA/GSFC 1991.
;       11/30/94 RWT converted to procedure gausspsf, and replaced
;                    gaussian function with gauss procedure
;       
;-
;***************************************************************************
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
 parcheck,npar,2,'gausspsf'
;
 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
;		psf = gaussian( x, [1,0,stdev] )
           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
;			psfx = gaussian( x, [ 1, 0, stdev(0) ] )
;			psfy = gaussian( y, [ 1, 0, stdev(1) ] )
			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
;		psfx = gaussian( x, [ 1, 0, stdev(0) ] )
;		psfy = gaussian( y, [ 1, 0, stdev(1) ] )
;		psfz = gaussian( z, [ 1, 0, stdev(2) ] )
		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