function IRIS_ANALYZE_WOBBLE, t0, t1, $
	hdr = hdr, data = data, result = result, thresh = thresh, $
	orbid = orbid, rollang = rollang, wait = wait, aia = aia, sref = sref, $
	save = save, load = load, png = png, runtime = runtime, opath = opath, $
	vertical = vertical, horizontal = horizontal, despike = despike, fuv = fuv, $
	xx1 = xx1, xx2 = xx2, yy1 = yy1, yy2 = yy2, slit = slit, pickorb = pickorb, $
	nocheck = nocheck, win = win, pstop = pstop, up_phase = up_phase, nrt = nrt, $
	chan = chan

;+
;
; Interactively analyze IRIS wobble data and create a wobble correction file
;
; INPUTS:
;	t0	-	Start time
;	t1	-	End time
;
; RETURNS:
;	result	-	a structure with the following fields
;		t0 and t1: Input time parameters
;		tphase: orbital phase (N-element array)
;		offset: measured X and Y pixel shifts from reference image (2xN array)
;		pztoffset: needed A, B and C PZT offsets (3xN array)
;		hdrwobble: X and Y pixel offset applied to the images during
;			data collection (2xN array)
;		
; KEYWORDS:
;	hdr	-	(OPTIONAL INPUT/OUTPUT) a variable holding the structure array of
;			image headers used in the correlation
;	data -	(OPTIONAL INPUT/OUTPUT) a variable holding the array of (processed) 
;			image data
;	result-	(OPTIONAL INPUT) a variable (returned by a previous run of this
;			routine); if passed in, then nothing is calculated, it just replots it	
;	orbid	-	(OPTIONAL INPUT) set to an integer for the ID of the wobble
;				table to compare the result with in the final plot
;	/pickorb-	(SWITCH) set to open a dialog window to select the reference
;				wobble table
;	rollang	-	(OPTIONAL INPUT) set to a float to rotate the images before
;				doing the cross-correlation. Useful if you're looking at or 
;				near a limb. If set, you will have a chance to adjust it iteratively
;	/vertical-	(SWITCH) set if you are trying to align on the E or W limb and 
;				want a vertical reference line
;	/horizontal-	(SWITCH) set if you are trying to align on the N or S limb
;					and want a horizontal reference line
;	/despike -	(SWITCH) set if you want to despike the data during IRIS_PREP
;	/nocheck -	(SWITCH) set if you want to look at a non-standard OBS ID; 
;				If not set, the program returns if it can't find data taken
;				with a wobble table OBS ID
;	win	-	(OPTIONAL INPUT) integer for plot window index ID
;	wait-	(OPTIONAL INPUT) delay (in sec) between frames; defaults to 0.01
;	/save	-	(SWITCH) set to save a genx file with the results of the measurement
;	/png -	(SWITCH) Save png files of output plots
;	opath-	(OPTIONAL INPUT) string specifying output path for png and genx
;	/load	-	(SWITCH) if set, load previously saved results and plot them
;	/pstop	-	(SWITCH) Stop for debugging
;	/slit	-	(SWITCH) (NOT IMPLEMENTED) Apply adjustment for slit wobble
;	/aia	-	(SWITCH) (NOT IMPLEMENTED) Compare self-correlation measurement
;				with a measurement derived from cross-correlation with SDO-AIA
;	/up_phase-	(SWITCH) set if you want to look up the orbital phase from the 
;				online side database and update the value in the headers
;	/nrt	-	(SWITCH) set if you want to use NRT data (because the data were 
;				taken recently and definitive data are not available)
;	xx1/xx2/yy1/yy2	-	Coordinates of the corners of the subregion used for
;						correlation
;	/runtime-	(SWITCH) if set, print the time it took to complete
;
;-

; Set up defaults
tt0 = SYSTIME(/sec)
pper=5856.
answer = ''
if N_ELEMENTS(rollang) eq 0 then rollang = 0.
irollang = rollang
if KEYWORD_SET(nrt) then begin
	ds = 'iris.lev1_nrt' 
	up_phase = 1
endif else begin
	ds = 'iris.lev1'
	nrt = 0
endelse
if N_ELEMENTS(chan) eq 0 then if KEYWORD_SET(fuv) then chan = 'SJI_1400' else chan = 'SJI_2832'
if not KEYWORD_SET(win) then win = 10
if not KEYWORD_SET(sref) then sref = 0.16635
if not KEYWORD_SET(wait) then wait = 0.01
if not KEYWORD_SET(despike) then despike = 0
if not KEYWORD_SET(opath) then opath = '/Volumes/disk2/data/iris/ops/wobble/'
if KEYWORD_SET(xx1) then insub = 1 else insub = 0
vertical = KEYWORD_SET(vertical)
horizontal = KEYWORD_SET(horizontal)
TVLCT, rr, gg, bb, /get
oldpmulti = !p.multi
PB_SET_LINE_COLOR

; If the /load keyword is set, then pick results from a previously-saved run
if KEYWORD_SET(load) then begin
	savefile = DIALOG_PICKFILE(path = opath, filter = '*.genx')
	RESTGEN, file = savefile, str = result
endif

; Skip to the plotting part if you already have a result structure
if N_TAGS(result) gt 0 then begin
	ttag = TIME2FILE(/sec, result.t0) + '-' + TIME2FILE(/sec, result.t1) + '-' + $
		TIME2FILE(/sec, RELTIME(/now))
	goto, plotit
endif

; Load an existing ORB table into structure for comparison
orbdir = '/irisa/ops/iris_user_tables/ORB/'
if KEYWORD_SET(pickorb) then begin
	orbfile = DIALOG_PICKFILE(path=orbdir)
	orbid = FIX(STRMID(orbfile, STRLEN(orbfile)-6, 3))
endif else begin
	if N_ELEMENTS(orbid) eq 0 then orbid = 22
	orbfile = CONCAT_DIR(orbdir, 'ORB-' + STRING(orbid, form = '(i05)') + '.xml')
endelse
orbdat = SSW_IRIS_XML2STRUCT(orbfile, /add)

; If necessary, get and prep the data
if N_ELEMENTS(data) eq 0 then begin
	files = IRIS_TIME2FILES(t0, t1, drms, /jsoc2, url = remote, nrt = nrt, key='IMG_PATH,ISQOLTID')
;	SSW_JSOC_TIME2DATA, t0, t1, ds = ds, drms, files, /files_only, /jsoc2
	; Only use images taken using the wobble table OBS ID
	if KEYWORD_SET(nocheck) then begin
		wobims = WHERE(drms.img_path eq chan, numwob) 
	endif else begin
		wobims = WHERE(drms.img_path eq chan and drms.isqoltid ge 4202100002ll and $
			drms.isqoltid le 4202100004ll and FILE_EXIST(files), numwob)
	endelse
	if numwob eq 0 then begin
		PRINT, 'The time range specified does not include a wobble sequence!'
		RETURN, 0
	endif
	READ_IRIS, files[wobims], hdr, dat
	IRIS_PREP, hdr, dat, ohdr, odat, despike = despike, /nowarp
	dat = 0
	data = odat[0:1099,*,*]
	odat = 0
endif else begin
	t0 = hdr[0].t_obs
	t1 = hdr[-1].t_obs
endelse

nimg = N_ELEMENTS(hdr)
tphase = (hdr.ophase + FLOOR(MIN(hdr.ophase))) mod 1
if KEYWORD_SET(up_phase) then begin
	mag = IRIS_TIME2HK(hdr, /mag)
	tphase = mag.ophase
endif
ttag = TIME2FILE(/sec, t0) + '-' + TIME2FILE(/sec, t1) + '-' + TIME2FILE(/sec, RELTIME(/now))

; If necessary, measure slit wobble so that it can be corrected in final wobble table
if KEYWORD_SET(slit) then begin
	pdb = IRIS_MK_POINTDB(version = pdb_ver)
	pdb_tags = TAG_NAMES(pdb)
	sji_wv_nm = STRMID(chan, 4, 3)
	nomx = pdb.(WHERE(pdb_tags eq 'CPX1_' + sji_wv_nm))
	nomy1 = (pdb.(WHERE(pdb_tags eq 'CPX2_' + sji_wv_nm)) - 270)
	nomy2 = (pdb.(WHERE(pdb_tags eq 'CPX2_' + sji_wv_nm)) + 270)
	startval = [nomx, nomy1, nomx, nomy2]
	PRINT, 'Measuring slit wobble...'
	slits = IRIS_FID_FINDER(hdr, data, startval = startval)
	!p.multi = [0, 1, 2]
	PLOT, tphase, slits.x1 - MEDIAN(slits.x1), psym = 4, yrange = [-1,1], chars = 1.5, $
		xtit = 'Orbital Phase', ytit = 'Fiducial offset from Median [IRIS pix]', $
		tit = 'Slit X'
	OPLOT, tphase, slits.x1 - MEDIAN(slits.x1), psym = 4, col = 2
	OPLOT, tphase, slits.x2 - MEDIAN(slits.x2), psym = 4, col = 6
	LEGEND, chars = 1.5, pos = 10, ['Bottom', 'Top'], col = [2, 6], psym = 4
	PLOT, tphase, slits.y1 - MEDIAN(slits.y1), psym = 4, yrange = [-1,1], chars = 1.5, $
		xtit = 'Orbital Phase', ytit = 'Fiducial offset from Median [IRIS pix]', $
		tit = 'Slit Y'
	OPLOT, tphase, slits.y1 - MEDIAN(slits.y1), psym = 4, col = 2
	OPLOT, tphase, slits.y2 - MEDIAN(slits.y2), psym = 4, col = 6
	!p.multi = 0
	PRINT, 'How should slit wobble be corrected?'
	PRINT, '1 - No correction'
	PRINT, '2 - Average both fiducials'
	PRINT, '3 - Use only bottom fiducial'
	PRINT, '4 - Use only top fiducial'
	READ, answer, prompt = 'Enter [1-4]: '
	case FIX(answer) of
		2	:	begin
			slitx = slits.x1 - MEDIAN(slits.x1) + slits.y1 - MEDIAN(slits.y2)
			slity = slits.y1 - MEDIAN(slits.y1) + slits.y1 - MEDIAN(slits.y2)
		end
		3	:	begin
			slitx = slits.x1 - MEDIAN(slits.x1)
			slity = slits.y1 - MEDIAN(slits.y1)
		end
		4	:	begin
			slitx = slits.x2 - MEDIAN(slits.x2)
			slity = slits.y2 - MEDIAN(slits.y2)
		end
		else:	begin
			slitx = FLTARR(nimg)
			slity = FLTARR(nimg)
		end
	endcase
	; Set up the parameter structure for using MPFITFUN to fit a sine wave
	numparams = 4
	params = ['Amplitude', 'Frequency', 'Phase', 'Offset']
	values = [-0.3d, (2 * !pi) / 5856d, 1d, 100d]
	fixed = [0, 1, 0, 0]
	limited = [0, 0, 0, 0]
	limits = [1, 1, 1, 1]
	parinfo = REPLICATE({name:params[0], value:values[0], fixed:fixed[0], $
		limited:[1,1]*limited[0], limits:[-1d,1d]*limits[0]}, numparams)
	for i = 1, numparams-1 do begin
		parinfo[i].name 	= params[i]
		parinfo[i].value 	= values[i]
		parinfo[i].fixed 	= fixed[i]
		parinfo[i].limited 	= limited[i] * [1,1]
		parinfo[i].limits 	= limits[i] * [-1d, 1d]
	endfor
	fsineparams = MPFITFUN('IRIS_MP_SINECHISQ', tphase, slitx, FLTARR(nimg)+0.5, $
		parinfo = parinfo, yfit = sinefit, quiet=1 )	;	quiet = (1-loud)	
	PLOT, tphase, slitx, psym = 4, yrange = [-1,1], chars = 1.5, $
		xtit = 'Orbital Phase', ytit = 'Fiducial offset from Median [IRIS pix]', $
		tit = 'Slit X'
	OPLOT, tphase, sinefit, col = 2
	slitx = sinefit
	fsineparams = MPFITFUN('IRIS_MP_SINECHISQ', tphase, slity, FLTARR(nimg)+0.5, $
		parinfo = parinfo, yfit = sinefit, quiet=1 )	;	quiet = (1-loud)	
	PLOT, tphase, slity, psym = 4, yrange = [-1,1], chars = 1.5, $
		xtit = 'Orbital Phase', ytit = 'Fiducial offset from Median [IRIS pix]', $
		tit = 'Slit Y'
	OPLOT, tphase, sinefit, col = 2
	slity = sinefit
	if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-slit.png')
endif else begin
	slits = 0
	slitx = FLTARR(nimg)
	slity = FLTARR(nimg)
endelse

; Check that the data did not have a wobble table running
hdrwobble = FLTARR(2, nimg)
for i = 0, nimg - 1 do hdrwobble[*,i] = IRIS_PZT2XY(hdr[i].IWBPZTA, $
	hdr[i].IWBPZTB, hdr[i].IWBPZTC)/sref
PLOT, tphase, hdrwobble[0,*], psym = 4, xtit = 'Orbital Phase', $
	ytit = 'PZT Wobble Correction [IRIS pix]', yrange = [-20, 20], $
	title = 'Wobble correction from headers'
OPLOT, tphase, hdrwobble[1,*], psym = 6
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-hdrwob.png')
if MAX(ABS(hdrwobble)) gt 1. then begin
	PRINT, 'Non-zero wobble table correction active during data collection...'
endif else begin
	PRINT, "Data collected without wobble table running. That's good..."
endelse

; If a non-zero rollang is specified (or /vert or /horiz is set),
; roll the data and interactively prompt the user to accept or iterate the
; roll angle
roll_it:
	if rollang ne 0. or vertical or horizontal then begin
		LOADCT, 0, /silent
		WDEF, win, 1200, 800
		!p.multi = [0, 2, 1]
		rolled = 0
		refpix = 550
		while not rolled do begin
			PLOT_IMAGE, data[*,*,0], min = 0, max = 400, chars = 1.5, $
				title = 'Original'
			PLOT_IMAGE, ROT(REFORM(data[*,*,0]), rollang, cubic = -0.5), $
				min = 0, max = 400, chars = 1.5, title = 'Rolled'
			if vertical then OPLOT, refpix * [1,1], [0,1095], thick = 2
			if horizontal then OPLOT, [0,1099], refpix * [1,1], thick = 2
			PRINT, 'How does the roll angle look?'
			PRINT, '1 - Accept'
			PRINT, '2 - Adjust roll'
			PRINT, '3 - Adjust reference line position'
			READ, answer, prompt = 'Enter [1-3]: '
			case FIX(answer) of
				2		:	begin
					READ, answer, prompt = 'Enter new roll angle [Deg]: '
					rollang = FLOAT(answer)
				end
				3		:	begin
					READ, answer, prompt = 'Enter new reference line position [pix]: '
					refpix = FIX(answer)
				end
				else	:	begin		;	Anything other than 2 or 3 is treated like 1
					PRINT, 'Accepting current roll...'
					rolled = 1
				end
			endcase
		endwhile
		rotd=data
		for i=0, nimg-1 do begin
			rotd[*,*,i] = ROT(REFORM(data[*,*,i]),rollang,cubic=-0.5)
		endfor
		data = rotd
		rotd = 0
		!p.multi = 0
	endif

; Now prompt the user to select a subregion for cross-correlation
sub_select:
	LOADCT, 0, /silent
	WDEF, win, 800, 800
	; If the user passed in subregion coordinates, apply them
	if insub then begin
		subdat = data[xx1:xx2, yy1:yy2, *]
		subbed = 1
	endif else begin
		subbed = 0
	endelse
	while not subbed do begin
		PLOT_IMAGE, data[*,*,0], min = 0, max = 400, chars = 1.5
		READ, answer, prompt = 'Enter X of lower left corner: '
		xx1 = FLOAT(answer)
		READ, answer, prompt = 'Enter Y of lower left corner: '
		yy1 = FLOAT(answer)
		READ, answer, prompt = 'Enter X of upper right corner: '
		xx2 = FLOAT(answer)
		READ, answer, prompt = 'Enter Y of upper right corner: '
		yy2 = FLOAT(answer)
		; Make the subregion and display it to the user
		subdat = data[xx1:xx2, yy1:yy2, *]
		PLOT_IMAGE, subdat[*,*,0], min = 0, max = 400, chars = 1.5
		for i = 1, nimg-1 do begin
			PLOT_IMAGE, subdat[*,*,i], /noerase, min = 0, max = 400, chars = 1.5
			WAIT, wait
		endfor
		READ, answer, prompt = 'Re-select subregion? '
		if STRUPCASE(STRMID(answer, 0, 1)) ne 'Y' then subbed = 1
	endwhile

; Now run TR_GET_DISP iteratively on the subregion
corred = 0
if KEYWORD_SET(thresh) then cordat = (subdat gt thresh) * 300 $
	else cordat = subdat
WDEF, win, 3*(xx2-xx1), 3*(yy2-yy1)
WDEF, win+1, 1200, 800
while not corred do begin
	if answer ne '2' then begin
		
		; Run cross-correlation
		this_offset = TR_GET_DISP(cordat)

		; Filter out bad points (that deviate substantially from the smoothed average)
		smooth0 = MEDIAN(REFORM(this_offset[0,*]), 10)
		bad0 = WHERE( ABS(this_offset[0,*] - smooth0) gt 5, numbad)
		if numbad gt 0 then begin
			this_offset[0, bad0] = smooth0[bad0]
			PRINT, 'Filtering ', numbad, ' bad images...'
			PRINT, bad0
		endif
		smooth1 = MEDIAN(REFORM(this_offset[1,*]), 10)
		bad1 = WHERE( ABS(this_offset[1,*] - smooth1) gt 5, numbad)
		if numbad gt 0 then begin
			this_offset[1, bad1] = smooth1[bad1]
			PRINT, 'Filtering ', numbad, ' bad images...'
			PRINT, bad1
		endif
		if N_ELEMENTS(all_offset) eq 0 then all_offset = this_offset $
			else all_offset = all_offset + this_offset

		PB_SET_LINE_COLOR
		PLOT, tphase, all_offset[0,*], psym = 4, yrange = [MIN(all_offset), MAX(all_offset)], $
			chars = 1.5, xtitle = 'Orbital Phase', ytitle = 'Offset [IRIS Pixels]'
		OPLOT, tphase, all_offset[1,*], psym = 6
		LEGEND, pos = 6, chars = 1.5, psym = [4, 6], ['X Shift', 'Y Shift']
		cordat = SHIFT_IMG(cordat, this_offset)
		WSET, win
		LOADCT, 0, /silent
		PLOT_IMAGE, cordat[*,*,0], min = 0, max = 400, chars = 1.5
		for i = 0, nimg-1 do begin
			PLOT_IMAGE, cordat[*,*,i], /noerase, min = 0, max = 400, chars = 1.5
			WAIT, wait
		endfor
		WSET, win + 1
	endif
	PRINT, 'How does that look?'
	PRINT, '1 - Accept'
	PRINT, '2 - Re-play movie'
	PRINT, '3 - Iterate cross-correlation'
	PRINT, '4 - Re-start cross-correlation'
	PRINT, '5 - Go back and re-do selection'
	READ, answer, prompt = 'Enter [1-5]: '
	case FIX(answer) of
		2		:	begin
			WSET, win		
			for i = 0, nimg-1 do begin
				PLOT_IMAGE, cordat[*,*,i], /noerase, min = 0, max = 400, chars = 1.5
				WAIT, wait
			endfor
			WSET, win + 1
		end
		3		:	begin
		end
		4		:	begin
			if KEYWORD_SET(thresh) then cordat = (subdat gt thresh) * 300 $
				else cordat = subdat
		end
		5		:	begin
			rollang = irollang
			goto, roll_it
		end
		else	:	begin		;	Anything other than 2-5 is treated like 1
			PRINT, 'Accepting current shifts...'
			corred = 1
		end
	endcase
endwhile

if KEYWORD_SET(pstop) then STOP

; Generate the output structure
;offset = 0 - (all_offset + hdrwobble)
offset = (all_offset + hdrwobble)
orbphase = orbdat.data.time / pper / 1000.
pztoffset = FLTARR(3, nimg)
secoffx = (offset[0,*] - MEDIAN(offset[0,*])) * sref 
secoffy = (offset[1,*] - MEDIAN(offset[1,*])) * sref
for i = 0, nimg - 1 do pztoffset[*,i] = IRIS_XY2PZT(secoffx[i], secoffy[i], /no)

orbxys = FLTARR(2, N_ELEMENTS(orbphase))
for i = 0, N_ELEMENTS(orbphase) - 1 do $
	orbxys[*,i] = IRIS_PZT2XY(orbdat.data[i].pztaoffset, $
		orbdat.data[i].pztboffset, orbdat.data[i].pztcoffset)/sref

result = CREATE_STRUCT('t0', t0, 't1', t1, 'tphase', tphase, $
	'offset', offset, 'hdrwobble', hdrwobble, 'pztoffset', pztoffset, $
	'rollang', rollang, 'xx1', xx1, 'yy1', yy1, 'xx2', xx2, 'yy2', yy2, $
	'orbid', orbid, 'orbphase', orbphase, 'orbdat', orbdat, 'orbxys', orbxys, $
	'hdr', hdr, 'cordat', BYTSCL(cordat, min = 0, max = 400), $
	'slits', slits, 'slitx', slitx, 'slity', slity)

; Now plot the new data, as well as the existing ORB table 
plotit:
PB_SET_LINE_COLOR
WDEF, win+2, 1000, 800

plotted = 0
pixoff = 0
pztoff = 0
phaseoff = 0
while not plotted do begin
	plotphase = ( (result.tphase + phaseoff) - FLOOR(MIN(result.tphase + phaseoff)) ) mod 1
	
	; Generate a plot comparing the frame-to-frame offset as measured for the
	; present wobble table with the frame-to-frame offset that is being compensated
	; for in an existing wobble table
	WSET, win + 2
	PLOT, plotphase, result.offset[0,*]+pixoff, psym = 4, chars = 1.5, $
		xtit = 'Orbital Phase', xrange = [0, 1], ytitle = 'Offset to compensate for [IRIS Pixels]', yrange = [-25, 25], $
		title = 'Comparison with ORB-000' + STRTRIM(result.orbid,2), /ystyle
	OPLOT, plotphase, result.offset[1,*]+pixoff, psym = 6, col = 2
	OPLOT, result.orbphase, 0-result.orbxys[0,*], col = 6, psym = 4, thick = 2
	OPLOT, result.orbphase, 0-result.orbxys[1,*], col = 5, psym = 6, thick = 2
	LEGEND, col = [255, 2, 6, 5], psym = [4, 6, 4, 6], pos = 10, chars = 1.5, $
		['Recommended X', 'New Y', 'Old X', 'Old Y']
	if phaseoff ne 0 or pixoff ne 0 then $
		LEGEND, pos = 12, chars = 1.5, ['Phase offset = ' + STRTRIM(phaseoff,2), $
			'Pixel offset = ' + STRTRIM(pixoff, 2)]
	if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-pixcomp.png')

	WSET, win + 1
	PLOT, plotphase, result.pztoffset[0,*]+pztoff, psym = 4, chars = 1.5, xtit = 'Orbital Phase', $
		ytitle = 'PZT Offset To Apply', title = 'Comparison with ORB-000' + STRTRIM(result.orbid,2), $
		xrange = [0,1], yrange = [MIN(result.pztoffset), MAX(result.pztoffset)]+pztoff+[-1,1]
	OPLOT, plotphase, result.pztoffset[1,*]+pztoff, psym = 6
	OPLOT, plotphase, result.pztoffset[2,*]+pztoff, psym = 2
	OPLOT, result.orbphase, result.orbdat.data.pztaoffset, col = 6, psym = 4, thick = 2
	OPLOT, result.orbphase, result.orbdat.data.pztboffset, col = 6, psym = 6, thick = 2
	OPLOT, result.orbphase, result.orbdat.data.pztcoffset, col = 6, psym = 2, thick = 2
	LEGEND, col = [255, 255, 255, 6, 6, 6], psym = [4, 6, 2, 4, 6, 2], pos = 10, chars = 1.5, $
		['New PZT A', 'New PZT B', 'New PZT C', 'Old PZT A', 'Old PZT B', 'Old PZT C']
	if phaseoff ne 0 or pztoff ne 0 then $
		LEGEND, pos = 12, chars = 1.5, ['Phase offset = ' + STRTRIM(phaseoff,2), $
			'PZT offset = ' + STRTRIM(pztoff, 2)]
	if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-pztcomp.png')

	PRINT, 'How do the plots look?'
	PRINT, '1 - Accept'
	PRINT, '2 - Adjust phase'
	PRINT, '3 - Adjust pixel offset'
	PRINT, '4 - Adjust PZT offset'
	READ, answer, prompt = 'Enter [1-4]: '
	case FIX(answer) of
		2		:	begin
			READ, answer, prompt = 'Enter new phase shift [-1 to +1]: '
			phaseoff = FLOAT(answer)
		end
		3		:	begin
			READ, answer, prompt = 'Enter new pixel offset [pix]: '
			pixoff = FLOAT(answer)
		end
		4		:	begin
			READ, answer, prompt = 'Enter new PZT offset [DN]: '
			pztoff = FLOAT(answer)
		end
		else	:	begin		;	Anything other than 2-4 is treated like 1
			PRINT, 'Keeping current plot...'
			plotted = 1
		end
	endcase
endwhile

WDEF, win+3, 1000, 800
PLOT, result.tphase, result.offset[0,*]-result.hdrwobble[0,*], psym = 4, $
	chars = 1.5, xtit = 'Orbital Phase', xrange = [0, 1], $
	ytitle = 'Offset [IRIS Pixels]', yrange = [MIN(result.offset), MAX(result.offset)]-result.hdrwobble, $
	title = 'Offset on top of current correction (ORB-000' + STRTRIM(result.orbid,2) + ')'
OPLOT, result.tphase, result.offset[1,*]-result.hdrwobble[1,*], psym = 6, col = 2
PLOTS, [0, 1], [0, 0], col = 5, line = 2
LEGEND, col = [255, 2], psym = [4, 6], pos = 10, chars = 1.5, ['X', 'Y']
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-pixdelt.png')

TVLCT, rr, gg, bb
!p.multi = oldpmulti

if KEYWORD_SET(pstop) then STOP

if KEYWORD_SET(save) then begin
	savefile = CONCAT_DIR(opath, ttag + '-iwobble.genx')
	SAVEGEN, file = savefile, str = result
endif

if KEYWORD_SET(runtime) then begin
	deltat = STRING(SYSTIME(/sec) - tt0, form = '(f7.2)')
	PRINT, 'Elapsed time (IRIS_ANALYZE_WOBBLE): ', deltat, ' seconds'
endif

RETURN, result

end