function inside_outside_stats,f,w4

; INPUTS:
;   f=float(nx,nx) = FTransform map (or dirty map)
;   w4=source "where" array defining pixel addresses of source region
;   (from interior_source.pro)
;
; OUTPUTS:
;   stats = 4 x 2 matrix giving the mean, SD, minmax inside and outside the
;   source

  f_n=f
  mean_inside=mean(f(w4))
  sd_inside=stdev(f(w4))
  mnmx_inside=minmax(f(w4))
  f_n(w4)=0  ; zero out the source region   
  n_outside=n_elements(f)-n_elements(w4)
  sum_outside=total(f_n) & mean_outside=sum_outside/n_outside
  sum2outside=total(f_n^2) & sd_outside=sqrt(sum2outside/n_outside-mean_outside^2)
  mnmx_outside=minmax(f_n)

return,[[mean_inside,sd_inside,mnmx_inside[0],mnmx_inside[1]],$
         [mean_outside,sd_outside,mnmx_outside[0],mnmx_outside[1]] ]

end


;+
function hsi_phasor_sum,cbe,det,nx,dx,xyshift,BKGD=bkgd

; PURPOSE:
;   Mapping program independent of modulation patterns
;   Similar to Fourier Transform, but uses calibrated eventlist and
;   phase_map_ctr to create phasor sum for each pixel of a map.

; INPUTS:
;   cbe = POINTER to calibrated eventlist
;   det = SCALAR detector number (0-8)
;   nx = SCALAR size of map in x and y directions (square map)
;   dx = SCALAR pixel size in arc sec
;   xyshift = pixels to shift from map center -- FLOAT(2)
;   optional background input: bkgd=SCALAR

; OUTPUTS:
;   Returns nx x nx array of phasor sums centered around the xyoffset
;   that was used to produce the calibrated eventlist.
;
; NOTE: The inspiration for this program comes from Hurford, Megeath,
; Palmer & Prince, Topical Meeting on Quantum-Limited Imaging,
; Honolulu, Mar 31-Apr 2, 1986.  Optical Society of America, Technical 
; Digest. Should work for single photons.

; 
; EJS, GSFC Code 682, January, 2001 -- mostly unoptimized
; EJS February 26, 2001 -- added correction for shadowing, livetime
; EJS June, 2001 -- changed from complex to real sum,
;                   incorporated bproj normalization
; EJS August, 2001 -- removed the unecessary ideal-grid correction
;-

phase_map_ctr=(*cbe(det,0)).phase_map_ctr
roll_angle=(*cbe(det,0)).roll_angle
count=(*cbe(det,0)).count
tau=(*cbe(det,0)).gridtran * (*cbe(det,0)).livetime
mamp=(*cbe(det,0)).modamp
nx=long(nx)

if n_elements(bkgd) EQ 0 then bkgd=0 
vcount=count-bkgd

; Equivalent ideal-grid count rate:
; ASSUMES: vcount=f0 * tau *(1 + mamp * cos(psi))
; AND:     0=mean(cos(psi))=mean(ct/(tau*mamp))/f0-mean(1/mamp)
; If the background is in error, the maps will be incorrect.
; f0=mean(vcount/(tau*mamp))/mean(1/mamp)
; cospsi=(vcount-f0*tau)/(f0*tau*mamp) ; =vcount/(f0*tau*mamp)-1/mamp
; vcount=f0*(1+(8/!pi^2)*cospsi)  ; fundamental of photon rate


; square arrays centered around shifted map center
x_m=(findgen(nx)+0.5-nx/2-xyshift[0])#replicate(1,nx) ; middle pxls: -0.5,0.5
y_m=-replicate(1,nx)#(findgen(nx)+0.5-nx/2-xyshift[1]); middle pxls: -0.5,0.5
xm=reform(x_m,nx^2)*dx
ym=reform(y_m,nx^2)*dx

N=n_elements(phase_map_ctr)

psi=xm*0 & P=fltarr(nx^2) 
pitch=((hsi_grid_parameters()).pitch)(det)
grid_angle=((hsi_grid_parameters()).orient)(det)

phi=roll_angle+grid_angle ; the + sign makes all the maps the same orientation
k=2*!pi/pitch
kcosphi=k*cos(phi)
ksinphi=k*sin(phi)
P_bar=fltarr(nx^2) & var=P_bar

!x.style=1 & !y.style=1

;  t=systime(1)
; Note: the time waster in this loop is the calculation of P0, which
;       takes 16 times as long as phz_m

delta=-1.0
 for m=0,nx^2-1 do begin
   phz_m=-xm[m]*ksinphi+ym[m]*kcosphi ;= kxm*sin(thetam-phi)
   ;   P0=(1+mamp*cos(phase_map_ctr - phz_m))*tau  ; old way
   P0=mamp*cos(phase_map_ctr-delta - phz_m)           ; new way
   mom=moment(P0,maxmoment=2)
   P_bar[m]=mom[0]     ;mean(P0)
   var[m]=mom[1]*mean(tau)    ; variance(P0)
   if (var[m] LT .0001) then begin
     message,'Singularity at pixel',locate_val(reform(var,64,64),var[m]),/info
     message,'Replacing var[m]',var[m],' with 1.0',/info
     var[m]=1.0
   endif
   psi[m]=total(vcount*(P0-P_bar[m])/var[m])  
endfor

; print,'phasor sum took',systime(1)-t,' sec'


return,reform(psi,nx,nx)

end

;+
FUNCTION HSI_OVERRES_EVAL,CBE,NX,DMAX,ftmap,NODISPLAY=nodisplay,$
                          PIXEL_SIZE=pixel_size,VERBOSE=verbose,  $
                          STATS=stats,TITLE=title
; PURPOSE:
;   Tries to evaluate a HESSI calibrated event list to determine which
;   subcollimators over resolve the source.
;
; METHOD:
;   Starts with detector dmax, makes a Fourier transform map
;   (hsi_phasor_sum), selects the half-power contour of the source,
;   then goes to the finer subcollimators, making FT maps and 
;   determining the ratio of the maximum inside the contour 
;   to the maximum outside the contour.
;   Usually, if this ratio is less than one, the subcollimator
;   overresolves the source, but noise, erronious pointing, or other
;   factors can throw it off.  Therefore the default is to show the
;   single-collimator maps.

; INPUTS:
;    CBE = pointer to a calibrated eventlist
;    NX  = dimensions of map
;    DMAX = coarsest collimator to use
;    ff  = FT maps
;    NODISPLAY -- if set, the maps will not be displayed
;    VERBOSE  -- if set, the inside/outside maxima will be printed
;
; OUTPUTS
;    ratio = FLOAT array for all selected detectors of the
;             inside/outside ratio
; OPTIONAL:: 
;     stats = FLOAT array of inside/outside statistics (minmax,rms,mean)
;     pixelsize =scalar float (arcseconds), defaults to 1.0
;
; RESTRICTIONS:
;     Maps must have  power-of-2 dimensions, no bigger than 128.
;
; VERSION HISTORY:
;    EJ Schmahl, GSFC Code 682 August, 2001
;-
if n_elements(pixel_size) eq 0 then pixel_size=1.0
dx=pixel_size
ftmap=fltarr(128,128,dmax)

i=0 & xyshift=[0.0,0.0]
if (keyword_set(PLOT)) then loadct, file=loc_file(path=hessi_data_paths(), 'colors_hessi.tbl'), 41

x_y=[1,1]
if (NOT keyword_set(NODISPLAY)) then window,xsize=128*dmax,ysize=128,title=title

for det=dmax-1,0,-1 do begin
  print,'DETECTOR',strcompress(det)
  f1map=hsi_phasor_sum(cbe,det,nx,dx,xyshift)  
  ftmap(*,*,det)=rebin(f1map,128,128)
  if det EQ dmax-1 then begin
    fm=ftmap(*,*,dmax-1)
    w4=interior_source(fm, 0.5*max(fm),xycont=x_y)
  endif
    ftmap1=ftmap(*,*,det)
    ftmap1(x_y(*,0),x_y(*,1))=max(ftmap1)
    ; help,x_y
  if (NOT keyword_set(NODISPLAY)) then begin
    tvscl,ftmap1,det
    loadct, file=loc_file(path=hessi_data_paths(), 'colors_hessi.tbl'), 41
  endif
endfor


nw4=n_elements(w4)

stats=fltarr(4,2,dmax)

for j=0,dmax-1 do begin
  fj=ftmap(*,*,j)
  stats(*,*,j)=inside_outside_stats(fj,w4)  ; 4 x 2 array
endfor

if keyword_set(VERBOSE) then begin
  print,'   INSIDE/OUTSIDE RATIO OF MAXIMA'
  ratio=stats(3,0,*)/stats(3,1,*)
  for j=0,dmax-1 do print,j,ratio[j]
endif

RETURN,ratio
END

;MAIN
restore,'flare9_cbe.sav',/verb & id='FLARE'+strcompress(id_number)
nx=64
dmax=6
; window,/free,xsize=dmax*128,ysize=128,title=id+'  Phasor_Sum Maps (0-'+strcompress(dmax-1,/rem)+')'

fmap=fltarr(128,128,dmax)
title=id+' PHASOR SUM MAPS, dets 1-'+strcompress(dmax,/rem)
stats=HSI_OVERRES_EVAL(CBE,NX,DMAX,fmap,/verbose,title=title)


end