function hsi_pt_spread_fun,psfsize,subcoll,pixel_size,bessel=bessel

; INPUTS:
;       psfsize=odd integer, size of psf array (square)
;         if psfsize is even, psfsize+1 will be used
;       subcoll=integer 1-9 (rhessi grid number)
;       pixel_size=integer pixel size in arcsec
; OUTPUT:
;       Point spread function for RHESSI subcollimator
; CALLING EXAMPLE:
;       PSF= hsi_pt_spread_fun(127,3,2,/bessel)
; EJS 10/20/05  Fixed error in gau expression
;               Added checks for undersampling & spillover
;     10/21/05  Added bessel function psf

n=fix(psfsize)

if 2*(n/2) eq n then n=n+1
r=shift(dist(n),n/2,n/2)*pixel_size

fwhm=2.26*(3.^(.5*(subcoll-1))) ; fwhm in arcsec
; check for spillover:
if fwhm gt 0.125*psfsize*pixel_size then message,'mapwidth < FWHM*8'
; check for undersampling:
if fwhm lt 2*pixel_size then message,'pixel size > FWHM/2

if keyword_set(bessel) then begin
  p=(hsi_grid_parameters()).pitch
  k=2*!pi/p(subcoll-1)
  psf=beselj(k*r<2.405,0)  ; select innermost lobe to first zero
  psf=psf-min(psf)
endif else begin
 psf=2.^(-(r/(fwhm/2))^2)
endelse

return,psf/total(psf)


end