pro IRIS_CHECK_WOBBLE, l2dir, opath, $
	disps, subdat, corrdat, $
	nocorr = nocorr, noplot = noplot, pixels = pixels, pstop = pstop, bitrate=bitrate, trim=trim

;+
;
; Examine the wobble in a set of level 2 data, in a number of ways:
;	1) Make low-cadence (~3 minute) SJI movies, which tend to make the
;		residual wobble much easier to identify;
;	2) Measure frame-to-frame cross-correlation of selected SJI channels
;		and plot the results
;
; INPUTS:
;	l2dir	-	Set to the path to the level2 data to look at
;				(e.g. /irisa/data/level2_nrt/2015/05/31/20150531_075615_3820356897/)
;	opath	-	Set to the path where the output files (movies, maybe a PNG
;				of a plot) should be written
;
; OUTPUTS:
;	disps
;
; KEYWORDS:
;	/nocorr	-	If set, then the correlation is not performed
;	/noplot	-	If set, then the output plots are not made (currently doesn't matter)
;	/pstop	-	Stop in the routine for debugging
;   /trim   -   If set, movie is trimmed to include only area that has data in all frames
;	pixels=		The size of the output movie, in un-summed IRIS pixels (i.e. set
;				pixels=400 to play a movie of a 400x400 region of the image)
;-

if n_elements(bitrate) eq 0 then bitrate=2e6

wobble_cadence = 180.

f = FILE_SEARCH(CONCAT_DIR(l2dir, 'iris_l2_*_SJI_*.fits'))
numf = N_ELEMENTS(f)
if f[0] eq '' then begin
	PRINT, 'IRIS_CHECK_WOBBLE: No SJI files found!'
	RETURN
endif
s2832 = STRPOS(f, 'SJI_2832')
if MAX(s2832) lt 0 then PRINT, $
	'IRIS_CHECK_WOBBLE: Try to use something with 2832 SJI'

; Read in the first SJI channel and check it out
READ_IRIS_L2, f[0], hdr, dat
exptime = hdr[0].exptime
nframe = N_ELEMENTS(hdr)
sumspat = hdr[0].sumspat
sumsptrl = hdr[0].sumsptrl
duration = MAX(hdr.time)
pztxs = LIMITS(hdr.pztx)


; Check that the OBS is reasonable for checking wobble
if duration lt 7200d then PRINT, $
	'IRIS_CHECK_WOBBLE: Try to use something with at least 2 hour duration'
if exptime le 2.0 then PRINT, $
	'IRIS_CHECK_WOBBLE: Try to use something with exposure time longer than 2 seconds'
if (pztxs[1] - pztxs[0]) gt 50 then PRINT, $
	'IRIS_CHECK_WOBBLE: Try to use something that doesnt raster so far'
	 
; Set up color tables/backing store
width = 700
height = 700
fps = 20
TVLCT, rr, gg, bb, /get
oldpmulti = !p.multi
WINDOW, 12, retain = 2, xsize = width, ysize = height
!p.multi = 0

; Make a wobble movie
for i = 0, numf - 1 do begin
	if i gt 0 then READ_IRIS_L2, f[i], hdr, dat
	nframe = N_ELEMENTS(hdr)
	sumspat = hdr[0].sumspat
	sumsptrl = hdr[0].sumsptrl
	tai0 = ANYTIM2TAI(hdr[0].date_obs)
	thispath = hdr[0].tdesc1
	filebase = CONCAT_DIR(opath, thispath + '_')
	dustdat = IRIS_DUSTBUSTER(hdr, dat)
	IRIS_LCT, hdr
	
	; Downsample in time to accentuate the wobble
	cadence = MEDIAN(DERIV(hdr.time))
	cadence_sample = ROUND(wobble_cadence / cadence) > 1
	nsub = nframe / cadence_sample
	index = INDGEN(nsub) * cadence_sample

	if keyword_set(trim) then begin
	    ; Trim down to only that part of the movie that contains data in all frames
	    dmin = MIN(dustdat, dim = 3)
	    dmask = dmin gt -200
	    dmx = TOTAL(dmask, 1)
	    dmy = TOTAL(dmask, 2)
	    subx = LIMITS(WHERE(dmx gt (MAX(dmx) * 0.8)))
	    suby = LIMITS(WHERE(dmy gt (MAX(dmy) * 0.8)))
	endif else begin
	    subx = [0,n_elements(dat[0,*,0])-1]
	    suby = [0,n_elements(dat[*,0,0])-1]
	endelse
	subdat = IRIS_INTSCALE(dustdat[suby[0]:suby[1], subx[0]:subx[1],*], hdr)
	sdsize = SIZE(subdat)	
	
	if KEYWORD_SET(pixels) then begin
		; Find the brightest part of the image and subframe it
		profx = SMOOTH(TOTAL(TOTAL(subdat, 3), 2), 20)
		profy = SMOOTH(TOTAL(TOTAL(subdat, 3), 1), 20)
		centx = WHERE(profx eq MAX(profx))
		centy = WHERE(profy eq MAX(profy))
		xmin = 0 > (centx - (pixels / 2. / sumsptrl))
		ymin = 0 > (centx - (pixels / 2. / sumspat))
		xmax = (xmin + (pixels / sumsptrl) ) < (sdsize[1]-1)
		ymax = (ymin + (pixels / sumspat) ) < (sdsize[2]-1)
	endif else begin
		; Just use the whole image
		xmin = 0
		ymin = 0
		xmax = sdsize[1]-1
		ymax = sdsize[2]-1
	endelse

	vfile = CONCAT_DIR(opath, thispath + '_wobble.mp4')
	CD, opath, current = old_dir
	video = IDLffVideoWrite(FILE_BASENAME(vfile))
	vidStream = video.AddVideoStream(width, height, fps, bit_rate=bitrate)
	for j = 0, nsub - 1 do begin
		PLOT_IMAGE, subdat[xmin:xmax,ymin:ymax,index[j]], chars = 1.5, scale = [sumsptrl, sumspat], $
			xtitle = 'Un-summed pixels [X]', ytitle = 'Un-summed pixels [Y]', $
			title = thispath + ' ' + TAI2UTC(tai0 + hdr[index[j]].time, /ccsds)
		thisframe = TVRD(/true)
		dummy = video.Put(vidStream, thisframe)
	endfor
	video.Cleanup
	CD, old_dir
endfor

; Make correlation-tracking plots (using the last SJI channel, 
; which is typically the 2832
if not KEYWORD_SET(nocorr) then begin
	corrdat = subdat
	; De-emphasize the slit so it doesn't confuse the correlation
;	corrdat[WHERE(subdat lt (MEDIAN(subdat)*0.6))] = MEDIAN(subdat)
	disps = TR_GET_DISP(corrdat, /shift)
	; Make correlation-corrected movie
	WINDOW, 14, xsize = width*2, ysize = height, retain = 2
	!p.multi = [0, 2, 1]
	vfile = CONCAT_DIR(opath, thispath + '_correl.mp4')
	CD, opath, current = old_dir
	video = IDLffVideoWrite(FILE_BASENAME(vfile))
	vidStream = video.AddVideoStream(width*2, height, fps, bit_rate=bitrate)
	for j = 0, nsub - 1 do begin
		PLOT_IMAGE, subdat[xmin:xmax,ymin:ymax,index[j]], chars = 1.5, scale = [sumsptrl, sumspat], $
			xtitle = 'Un-summed pixels [X]', ytitle = 'Un-summed pixels [Y]', $
			title = 'Original ' + TAI2UTC(tai0 + hdr[index[j]].time, /ccsds), $
			min = 0, max = 255
		PLOT_IMAGE, corrdat[xmin:xmax,ymin:ymax,index[j]], chars = 1.5, scale = [sumsptrl, sumspat], $
			xtitle = 'Un-summed pixels [X]', ytitle = 'Un-summed pixels [Y]', $
			title = 'Corrected ' + TAI2UTC(tai0 + hdr[index[j]].time, /ccsds), $
			min = 0, max = 255
		thisframe = TVRD(/true)
		dummy = video.Put(vidStream, thisframe)
	endfor
	video.Cleanup
	CD, old_dir
	; Make plot of correlation
	!p.multi = 0
	LOADCT, 0
	WINDOW, 12, xsize = width, ysize = height, retain = 2
	PLOT, disps[0,*] * sumsptrl, disps[1,*] * sumspat, chars = 1.5, psym = -4, $
		xtitle = 'Un-summed Pixels X', ytitle = 'Un-summed Pixels Y', /iso, $
		title = 'Frame-to-frame correlation'
	thisplot = TVRD(/true)
	WRITE_PNG, CONCAT_DIR(opath, 'correl.png'), thisplot
endif

; Make plots of the expected wobble correction
if not KEYWORD_SET(noplot) then begin
;	IRIS_COMPARE_WOBBLE
endif

if KEYWORD_SET(pstop) then STOP
TVLCT, rr, gg, bb
!p.multi = oldpmulti

end