function hsi_polar_psf,R_src,AZ_src,roll_angle,R0,det_index,mapHwidth,$
         MAP_AZIMUTH=map_azimuth,PSF_XY=psf_xy,PLOT=plot
;+
; NAME: hsi_polar_psf
;       
; PURPOSE:  Creates a polar map for a unit point source at (R0,AZ)
;               
; METHOD:
;      Compute the convolution of visibility and exp(i*k*r*cos(roll_angle))
;        This is the Fourier Transform in polar coordinates (r, roll_angle)
;
; INPUTS:
;       R_src = radial distance (asec) of point source from Sun center
;       AZ_src = azimuth of point source, measured from equatorial plane
;       det_index=vector of integers taken from [0-8]
;       mapHwidth = scalar width of map to be made (asec)
;       R0=radius of center of map
;       MAP_AZIMUTH = azimuth of map center.  If defined, only 90
;       degrees of angle will appear in the map.
;          (For some cockeyed reason, this azimuth starts on the
;          negative y azis, and goes clockwise. So if the source is 30 
;          degrees below the positive x axis, set map_azimuth=!pi+!pi/6.)
;
;       PSF_XY -- if set, the polar map will be interpolated to the xy plane
;       PLOT   -- if set, plots will be made to the screen
;
; OUTPUTS:
;       psf=point spread function in polar coordinates, normalized to
;       1 at the peak. 
;        Default Dimensions = n_elements(roll_angle) X 2*mapHwidth
;        If map_azimuth set, dimensions=n_elements(roll_angle)/4 X 2*mapHwidth
;       Lots of plots if optional keyword PLOT is set
;       Additional plot of xy projection if PSF_XY is set
;
; EXAMPLE:
;       N=512 & det_index=[4,6]
;       roll_angle=findgen(N)*2*!pi/N 
;       mapHwidth=256. & R0=700
;       R_src=600. & AZ_src=!pi/6
;       psf=hsi_polar_psf(R_src,AZ_src,roll_angle,R0,det_index,mapHwidth,$
;                MAP_AZIMUTH=!pi+!pi/6,/PSF_XY,/PLOT)
;
;
; FUNCTIONS CALLED:
;      hsi_convl.pro
;      hsi_psf_xy
;
; RESTRICTIONS:
;       Visibility must span 360 degrees of roll, even though maps may
;       span only 90 degrees is map_azimuth is set.
;
;
; VERSION HISTORY:
;        VERS. 1 EJS Nov 24, 1999
;        VERS. 1.1 Modified to output 90-degree sector EJS Dec 2, 1999
;        VERS. 1.2 EJS Jan 11, 2000
;              Made top the inner radius to be consistent with hsi_polar_mapper
;
; NOTES:
;-
;function hsi_polar_psf,R_src,AZ_src,R0,roll_angle,det_index,mapHwidth,$
;         MAP_AZIMUTH=map_azimuth,PSF_XY=psf_xy,PLOT=plot


t0=systime(1)
i=complex(0,1)
pitch=((hsi_grid_parameters()).pitch)(det_index)

  r2=r0+mapHwidth 
  r1=r0-mapHwidth 
  Nroll=n_elements(roll_angle)
  r=r1+findgen(r2-r1)   ; =[r1,r1+1,...,r2-1]
;  message,'Doing'+string(r2-r1)+' radii...',/inform
  Rsrc2picos=R_src*2*!pi*cos(roll_angle-AZ_src)
;  r2pi=r*2*!pi
  cos_roll=cos(roll_angle)
  cos_roll_x_r=cos_roll#r
  psf=fltarr(Nroll,r2-r1)

w_az=findgen(Nroll)
if keyword_set(map_azimuth) then begin
  map_az_start=map_azimuth*Nroll/(2*!pi) - Nroll/2
  w_az=(map_az_start + findgen(Nroll/4)) mod Nroll
endif

for jpitch=0,n_elements(pitch)-1 do begin
  vis=exp(i*Rsrc2picos/pitch[jpitch])          
  k=2*!pi/pitch[jpitch]
  v=exp(-i*k*cos_roll_x_r)
  t0=systime(1)
  for j=0,r2-r1-1 do begin
;       psf[*,j]=hsi_convl(vis,v[*,j]);1st arg is azimuth, 2nd is radius
;       psf1=hsi_convl(vis,v[*,j]);1st arg is azimuth, 2nd is radius
       psf1=hsi_convl(vis,v[*,j]);1st arg is azimuth, 2nd is radius

       psf[w_az,j]=psf1(w_az)  ; make top the inner radius
  endfor

;    print,"convl time=",systime(1)-t0
  
  ;t3=systime(1)
; This is an alternative, but slower form:  (source should be in 4th quadrant)
; The FFT form can do 512 angles in the time the SHIFT form does 11.
;  k0=256
;  for jr=0,r2-r1-1 do begin
;      vshift=shift(v[*,jr],0)
;      for k=k0,k0+127 do begin
;          psf(k,jr)=total(vis*vshift)
;          vshift=shift(vshift,1)
;      endfor
;  endfor
;  print,'multi-shift time=',systime(1)-t3

   if jpitch eq 0 then psf_sum=psf else psf_sum=psf_sum+psf
endfor ; jpitch



psf=psf_sum/max(psf_sum)

;print,'Time for psf=',systime(1)-t0,' sec'

if (keyword_set(plot)) then begin
  erase
  n_az=n_elements(psf(*,0))  & n_R=n_elements(psf(0,*)) 
;  window,0,xsize=n_az,ysize=n_R,title='polar_psf.pro: Real Part of PSF'
  erase
  
  tvscl,psf

endif

if keyword_set(psf_xy) then begin
  message,'Calling pmap2xy.pro...' ,/inform
;              psf: 1st index is azimuth (0 to 2*!pi), 2nd is  radius

  psf_xy=pmap2xy(psf,r1,r2) ;,/min_box)  
  message,'r1,r0,r2='+string(r1)+string(r0)+string(r2),/inform
  dxsize=(size(psf_xy))(1) &  dysize=(size(psf_xy))(2)
  if (!d.name NE 'PS') then window,1,xsize=dxsize,ysize=dysize,title='polar_psf.pro'
  tvscl,psf_xy
endif

return,psf
end