;+
; NAME: hsi_polar_clean.pro
;
;
; PURPOSE:
;       To extract a "clean" map from a "dirty" map
;
; METHOD:
;       Iteratively decomposes the image into single bright pixels,
;       subtracting out the side lobes for each point.
;
;       hsi_polar_clean does not know about collimators, pitches or angles. 
;
; INPUTS:
;       Dirty_map = polar map created from visibilities
;                   --1st index=azimuth,
;                     2nd index=radius (r1 to r2)
;
;       psf = 3D array of polar point-spread functions for points on
;           the ray azimuth=0 (positive x axis) for radii ranging from 
;           r1 to r2
;                1st index = azimuth of psf value (spans 360 degrees)
;                2nd index = radius of psf value  (spans r1 to r2)
;                3rd index = radius of point source (spans r1 to r2)
;
; OPTIONAL KEYWORDS:
;       NITER = number of iterations to perform (default = 100)
;       FRAC = amount of PSF to subtract in each iteration (default = 0.1)
;       FWHM = 2-element fwhm of convolving kernel
;       CHI_sq_crit = max allowed chi_sq (default = 1)
;
; OUTPUTS:
;        A structure containing:
;          clean_map
;          clean_components
;          clean_amplitudes
;          clean_vis 
;          niter
;          frac
;          chi-square statistic
;
; CALLING SEQUENCE
;     EXAMPLE:       
;      
;     CLEAN_MAP=hsi_polar_clean(DMAP, PSF, NITER = 300, fwhm=[2,8]$
;           CHI_SQ_CRIT = 1.)
;
;-

; ALGORITHM:
;      1. Copy dirty_map to map, set N=0
;      2. Find the brightest point (r,az) in map, if negative, quit.
;      3. Look at component list to see if (r,az) was done before
;         If not, then shift the PSF
;      4. Subtract PSF*frac from map, set N=N+1, save (r,az), increm clean_map
;      5. If compute model visibility and chi_sq
;      6. If chi_sq >  chi_sq_crit or N < Niter, go to 2;      
;      7. Add up all the clean components, replace vis with the
;         model vis and return.
;
;
;
; THIS VERSION ASSUMES ALL REQUIRED PSF ARRAYS EXIST--NO CACHE!

;CLEAN PROGRAM

function hsi_polar_clean,DIRTY_MAP,PSF,         $
                     NITER = niter,             $
                     FRAC = frac,               $
                     FWHM = fwhm,               $
                     CHI_SQ_CRIT = chi_sq_crit

        if NOT keyword_set(niter) then niter=100
        if NOT keyword_set(frac) then frac=0.1
        if NOT keyword_set(chi_sq_crit) then chi_sq_crit=1.


; 1. Copy dirty_map to map, initialize stuff

        Nmap=n_elements(dirty_map)
	map=dirty_map
        clean_map=map*0
	iter=0
        ncomp=0
        map_max=1.
        N_az=n_elements(dirty_map(*,0))
        N_r= n_elements(dirty_map(0,*))
        cc=[-1]
        ca=[-1]
        if keyword_set(fwhm) then $
         kernel=psf_gaussian(npixel=6*max(fwhm),fwhm=fwhm,/normalize)

message,'Starting first iteration of '+string(Niter),/inform
loop:
        blank=dirty_map*0   
; 2. Find the brightest point (x,y) in map, if negative, then exit loop

   max_pt:
        map_max=max(map)    & print,'map max=',max(map)
   if (map_max GT 0) then begin
        R_Az=(where(map eq max(map)))(0)
        R=R_Az/N_Az
        Az=R_Az-R*N_Az

; 3. Select the right PSF, then shift it to the current azimuth 
        PSF1=psf(*,*,R) ;  N_r-R-1)  
        PSF1=shift(psf1,Az-N_az/2,0) 

; 4. Subtract PSF*frac*map_max from map, increment iteration count and map
; print,'max(map) at  ',locate_val(map,max(map)),'        AZ,R=',AZ,R
if (where(psf1 eq max(psf1)))(0) NE (where(map eq max(map)))(0) then begin
      print,'psf1 max at ',locate_val(psf1,max(psf1)),$
           ' map max at ',locate_val(map,max(map))
      message,'psf1 is not peaking at location of map maximum!',/info
endif
          map=map-psf1*(frac*map_max) & print,' frac*map_max=', frac*map_max

;         xy=locate_val(map,max(map))
;         map1=map & map1(xy[0],*)=max(map) & map1(*,xy[1])=max(map)
;         wset,0 & tvscl,map1 & wset,1 & tvscl,psf1

          clean_map(R_Az)=clean_map(R_Az)+frac*map_max
          cc=[cc,R_Az]       ; clean component location
          ca=[ca,map_max*frac]   ; clean component amplitude
          iter=iter+1
          ncomp=ncomp+1
     print,'iter=',iter,' R=',R,' Az=',AZ,' max(map)=',max(map)

;     blank(R_Az)=map_max & contour,blank

; 5. If vis exists, compute model_vis and chi_sq

          if (Keyword_set(vis)) then begin
             plot,vis,xsty=1,tit='VIS AND MODEL_VIS)
             model_vis=_hsi_clean_model2vis(clean_map)
             oplot,model_vis,line=1
	     clean_flux=total(ca)
;             chi_sq = chi_square(vis,model_vis)
          endif

; 6. If chi_sq >  chi_sq_crit or N < Niter OR clean_flux < obs_flux, go to 2
          
     if (iter lt Niter) then goto, loop
   endif else begin
            print,"Negative Maximum. Exiting Loop."
   endelse

done:
   print,'iter=',iter
          

; 7. Return the result

      if keyword_set(FWHM) then clean_map=convolve(clean_map,kernel)
      output={clean_map:clean_map, $
              cc: cc[1:*], $
              ca: ca[1:*] } ;, $
;              chi_sq: chi_sq, $
;              model_vis: model_vis }

      return, output

 end