pro sot_focus_scan, index, scandata, best_focus, movie_graphic, $
   stddevx=stddevx, moments=moments, use_moment=use_moment, typey=typey , $
   closest=closest, movie_size=movie_size , tf_scan=tf_scan,$
   indiv_scale=indiv_scale, no_stat=no_stat
;
;+
;   Name: sot_focus_scan
;
;   Purpose: calculate statistics and estimate best focus position for SOT scand
;
;   Input Parameters:
;      index,scandata - assumed SOT focus scan index,data
;   
;   Output Parameters:
;      best_focus - estimated best focus position
;      movie_graphic - optionally return a movie summary of scan/calulation  
;                      (may be used to visually confirm or deny estimated best)
;
;   Keyword Parameters:
;      stddevsx (output) - standard deviations for input scan images (1:1)
;      moments (output) - corresponding moments (only defined if USE_MOMENT is set
;      use_moment (input) - if set, use moments instead of standard deviation
;                           (0,1,2,3) -> Select Mean,Variance(def),Skew,Kurtosis
;                           My guess is that only 1 & 3 will work
;      typey - (output string) - what param was use (Standard Dev, Variance...)
;      closest (output) - index closest in time to optimal
;      movie_size - optional size [nx,ny] if movie_graphic is "requested" (def =nx from scandata)
;      tf_scan - if set, look at TFn paramters...
;
;   Calling Examples:
;      IDL> sot_focus_scan,index,data,best,movieout ; use standard deviation
;      IDL> sot_focus_scan,index,data,best,movieout, /use_moment ; use/plot variance
;      IDL> sot_focus_scan,index,data,best,movieout, use_moment=3 ; use/plot kurtosis
;
;   Calling Context:
;      IDL> sot_focus_scan,index,data,bb,movie,use_moment=3
;      IDL> linecolors ; set color table for 'fancy colors'
;      IDL> xstepper,movie ; show returned focus scan movie+plot 
;   
;
;    Method:
;      Use image statistics comparisons to estimate best focus
;      Default=maximum of standard deviation; optionally Variance & Skewness
;
;   History:
;      26-October-2006 - Auto sotscans -> www applications
;
;-

no_stat=keyword_set(nostat)
use_mom=n_elements(use_moment) gt 0 ; default is standard deviation

tf_scan=keyword_set(tf_scan)
focus=gt_tagval(index,/focus)
range=[min(focus),max(focus)]
if n_elements(all_vals(range)) eq 1 and (1-tf_scan) then begin 
   box_message,'FOCUS is a constant over your data, bailing...'
   return
endif
allfoc=lindgen(deriv_arr(range))+range(0)
nimages=data_chk(scandata,/nimage)

typey='Standard Deviation'
momtypes=str2arr('Mean,Variance,Skewness,Kurtosis')
whichmom=1 ; init
if use_mom then begin 
   moments=fltarr(4,nimages)
   whichmom=use_moment<3  ; {0,1,2,3}
   typey=momtypes(whichmom)
   for i=0,nimages-1 do moments(0,i)=moment(scandata(*,*,i),$
     maxmoment=whichmom+1)
   fy=moments(whichmom,*) ; dependent variable
   moments=fy 
endif else begin 
   fy=fltarr(nimages)
   for i=0,nimages-1 do fy(i)=stddev(scandata(*,*,i))
   stddevx=fy
endelse

; now calculate best focus position using selected statistics
fx=interp1d(focus,fy,allfoc) ; interp to all focus positions in range
ssbest=where(fx eq call_function((['max','min'])(whichmom gt 2),fx))
best_focus=allfoc(ssbest(0))
resid=abs(focus-best_focus)
closest=index( (where(resid eq min(resid)))(0))

if n_params() ge 4 then begin ; user wants movie graphic
   ptemp=!d.name
   mindex=struct2ssw(index,/nopoint)
   mdata=confac(scandata,.5)
   if n_elements(movie_size) eq 0 then movie_size=[data_chk(mdata,/nx),256]
   wdef,xx,movie_size(0),movie_size(1),/zbuffer
   set_plot,'z'
   erase
   linecolors
   tit=([' Focus Scan',' Tunable Filter Scan'])(tf_scan)
   utplot,struct2ssw(mindex,/nopoint),fy,back=11,$
      color=([5,11])(no_stat),title='Hinode/SOT' +tit ,$
      ytitle=([typey,''])(no_stat) + (['',' & TF Pos'])(tf_scan) ,ynozero=(1-tf_scan),psym=2
   if tf_scan then begin 
      cols=[4,5,7,9,4,5,7,9] ; TF1,2,...8
      psyms=[replicate(4,4),replicate(2,4)]
      tftags='TF'+strtrim(indgen(8)+1,2)
      for i=0,7 do begin 
         outplot,mindex,mindex.(tag_index(mindex,tftags(i))), $
            color=cols(i), symsize=.7,psym=psyms(i)
      endfor
      legend,tftags,psym=psyms, colors=cols,textcolors=4,/center,/box,back=11,$
        charsize=.5,symsize=replicate(.7,8),/bottom
      
   endif else begin 
      evt_grid,mindex,label=strtrim(focus,2),$
      tickpos=.5,ticklen=.001,/vertical,labcol=4,labsize=.8, color=11 
      evt_grid,struct2ssw(closest,/nopoint),/arrow, color=9,labcol=9,ticklen=.01,tickpos=.92
      evt_grid,struct2ssw(closest,/nopoint),label='Best? '+strtrim(best_focus,2),labsize=.7,$
         color=9,labpos=.6
  
 
   endelse
   if keyword_set(indiv_scale) then begin 
      sdata=float(mdata)
      for i=0,data_chk(mdata,/nimage)-1 do $
         sdata(0,0,i)=hist_equal(sdata(*,*,i))
   endif else sdata=ssw_sigscale(mindex,float(mdata))
   sdata=bytscl(sdata,top=255-15)+15
   event_movie3,mindex,movie_graphic,data=sdata,outsize=movie_size,tcolor=11
   set_plot,ptemp ; restore plot device
endif

return
end