pro IRIS_MAKE_WOBBLE, vert_file, horiz_file, $
	ref_orbid = ref_orbid, pickorb = pickorb, opath = opath, slit = slit, $
	npoints = npoints, out_orbid = out_orbid, description = description, $
	png = png, win = win, pstop = pstop, save = save, sref = sref, smooth = smooth, $
	phaseshift = phaseshift, vert_angle = vert_angle, write = write, check = check, $
	_extra = extra

;+
;
; Loads the genx files from previously-analyzed wobble runs and generates a
; new wobble.xml file
;
; INPUTS:
;	vert_file-	The path and name of a genx file containing the results of a 
;				previous run of iris_analyze_wobble, for the "east limb" portion
;				of a wobble measurement (where the limb appears vertical in the image)
;	horiz_file-	The path and name of a genx file containing the results of a 
;				previous run of iris_analyze_wobble, for the "north limb" portion
;				of a wobble measurement (where the limb appears horizontal in the image)
;	
; KEYWORDS:
;	opath-	(OPTIONAL INPUT) string specifying path for genx files and xml wobble table
;	ref_orbid-	(OPTIONAL INPUT) can set to an integer specifying the ID
;				of the wobble table that the new one should be compared with
;	npoints-	(OPTIONAL INPUT) can set to an integer specifying how many points
;				are desired in the resulting wobble table file
;	out_orbid-	(OPTIONAL INPUT) can set to an integer specifying the ID 
;				to use for the output wobble table (if left off, no file is generated)
;	description-(OPTIONAL INPUT) can set to a string to be used as a description
;				for the wobble table
;	phaseshift-	(OPTIONAL INPUT) can set to a number between 0 and 1 indicating 
;				the phase shift to apply to the output wobble data
;	vert_angle-	(OPTIONAL INPUT; CURRENTLY HAS NO EFFECT AND SHOULDN'T BE SET) 
;				The angle of the solar limb with respect to the NS axis of the sun
;				in the vert_file. This should be the same as the sat_rot used
;				for the measurement 
;	smooth-		(OPTIONAL INPUT) Defaults to 20; set to 0 if you want no smoothin
;				Set to the number of points over which the raw data should be 
;				smoothed before being binned and averaged. 20-25 seems to help.
;	/check	-	(SWITCH) If set, then the user is given an opportunity to trim
;				the starting and ending points in the dataset
;
;-

; Set up defaults
ttag = TIME2FILE(/sec, RELTIME(/now))
if not KEYWORD_SET(opath) then opath = '/Volumes/disk2/data/iris/ops/wobble/'
if not KEYWORD_SET(win) then win = 10			; 	Window index
if not KEYWORD_SET(sref) then sref = 0.16635	;	Platescale (arcsec per pixel)
if not KEYWORD_SET(pper) then pper = 5856.		;	Orbital period in sec
if N_ELEMENTS(vert_angle) eq 0 then vert_angle = 0
if N_ELEMENTS(npoints) eq 0 then npoints = 30
if N_ELEMENTS(smooth) eq 0 then smooth = 20
if N_ELEMENTS(out_orbid) eq 0 then out_orbid = 0
if N_ELEMENTS(description) eq 0 then $
	description = 'Generated by IRIS_MAKE_WOBBLE on ' + ttag

; Set up for plotting
TVLCT, rr, gg, bb, /get
PB_SET_LINE_COLOR
oldpmulti = !p.multi
!p.multi = 0
WDEF, win, 1200, 800

; Trimming of the input data
h1 = 0
h2 = -1
v1 = 0
v2 = -1
	
; 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)
	ref_orbid = FIX(STRMID(orbfile, STRLEN(orbfile)-6, 3))
endif else begin
	if N_ELEMENTS(ref_orbid) eq 0 then ref_orbid = 22
	orbfile = CONCAT_DIR(orbdir, 'ORB-000' + STRTRIM(ref_orbid, 2) + '.xml')
endelse
table_label = 'Orb table ' + STRCOMPRESS([out_orbid, ref_orbid])
nref = N_ELEMENTS(ref_orbid)
for i = 0, nref - 1 do begin
	orbdat = SSW_IRIS_XML2STRUCT(orbfile[i], /add)
	if KEYWORD_SET(phaseshift) then begin
		orbdat.data.time = (pper * 1000. + orbdat.data.time + (phaseshift[i] * pper * 1000.)) mod (pper * 1000.)
		shift_t0 = WHERE(orbdat.data.time eq MIN(orbdat.data.time))
		orbdat.data.pztaoffset = STRING(FIX(orbdat.data.pztaoffset) - FIX(orbdat.data[shift_t0].pztaoffset))
		orbdat.data.pztboffset = STRING(FIX(orbdat.data.pztboffset) - FIX(orbdat.data[shift_t0].pztboffset))
		orbdat.data.pztcoffset = STRING(FIX(orbdat.data.pztcoffset) - FIX(orbdat.data[shift_t0].pztcoffset))
		table_label[i+1] = table_label[i+1] + ' + phase ' + STRING(phaseshift[i], form='(f6.3)')
	endif
	; Generate the X and Y offsets from the reference wobble data
	onpoints = N_ELEMENTS(orbdat.data.time)
	oxy = FLTARR(2,npoints)
	for j = 0, onpoints-1 do oxy[*,j] = IRIS_PZT2XY(orbdat.data[j].pztaoffset, $
		orbdat.data[j].pztboffset, orbdat.data[j].pztcoffset, /no)/sref
	this_orbdat = {time : orbdat.data.time, pztaoffset : orbdat.data.pztaoffset, $
		pztboffset : orbdat.data.pztboffset,  pztcoffset : orbdat.data.pztcoffset, $
	 	onpoints : onpoints, oxy : oxy}
	if i eq 0 then begin
		ref_orbdat = CREATE_STRUCT('orb' + STRTRIM(ref_orbid[i], 2), this_orbdat)
		ref_pztoffs = FIX([this_orbdat.pztaoffset, this_orbdat.pztboffset, this_orbdat.pztcoffset])
		ref_xyoffs = this_orbdat.oxy
	endif else begin
		ref_orbdat = CREATE_STRUCT(ref_orbdat, 'orb' + STRTRIM(ref_orbid[i], 2), this_orbdat)
		ref_pztoffs = [ref_pztoffs, FIX([this_orbdat.pztaoffset, this_orbdat.pztboffset, this_orbdat.pztcoffset])]
		ref_xyoffs = [ref_xyoffs, this_orbdat.oxy]
	endelse
endfor

; Load the results of previous IRIS_ANALYZE_WOBBLE runs, and plot them to allow
; the user to filter out bad points
CD, opath, current = old_dir
if N_ELEMENTS(vert_file) eq 0 then begin
	vert_prompt = 'Select the vertical-limb file (E limb for roll 0)'
	vert_file = DIALOG_PICKFILE(path = opath, filter = '*.genx', title = vert_prompt)
endif
RESTGEN, file = vert_file, str = vert_result
vert_off = REFORM(vert_result.offset[0,*])
vert_phase = vert_result.tphase
numv = N_ELEMENTS(vert_phase)
if KEYWORD_SET(check) then begin
	vert_ok = 'N'
	while vert_ok ne '' do begin
		PLOT, vert_phase[v1:v2], vert_off[v1:v2], psym = 4, chars = 1.5, $
			title = 'IRIS_MAKE_WOBBLE   Vertical Check', xtitle = 'Phase'
		XYOUTS, /data, vert_phase[v1], vert_off[v1], chars = 2, v1, col = 2, charthick = 2
		XYOUTS, /data, vert_phase[v2], vert_off[v2], chars = 2, v2, col = 2, charthick = 2
		WSHOW
		PRINT, 'Total # of vertical points = ', numv, form = '(a40,i10)'
		READ, vert_ok, prompt = 'Trim points? Enter as   first, last   or just enter to proceed:  '
		if vert_ok ne '' then begin
			nums = FIX(STRSPLIT(/extract, vert_ok, ','))
			v1 = nums[0] > 0
			v2 = nums[1] < (numv - 1)
		endif
	endwhile
	vert_phase = vert_phase[v1:v2]
	vert_off = vert_off[v1:v2]
endif

if N_ELEMENTS(horiz_file) eq 0 then begin
	horiz_prompt = 'Select the horizontal-limb file (N limb for roll 0)'
	horiz_file = DIALOG_PICKFILE(path = opath, filter = '*.genx', title = horiz_prompt)
endif
RESTGEN, file = horiz_file, str = horiz_result
horiz_off = REFORM(horiz_result.offset[1,*])
horiz_phase = REFORM(horiz_result.tphase)
numh = N_ELEMENTS(horiz_phase)
if KEYWORD_SET(check) then begin
	horiz_ok = 'N'
	while horiz_ok ne '' do begin
		PLOT, horiz_phase[h1:h2], horiz_off[h1:h2], psym = 4, chars = 1.5, $
			title = 'IRIS_MAKE_WOBBLE   Horizontal Check', xtitle = 'Phase'
		XYOUTS, /data, horiz_phase[h1], horiz_off[h1], chars = 2, h1, col = 2, charthick = 2
		XYOUTS, /data, horiz_phase[h2], horiz_off[h2], chars = 2, h2, col = 2, charthick = 2
		WSHOW
		PRINT, 'Total # of horizontal points = ', numh, form = '(a40,i10)'
		READ, horiz_ok, prompt = 'Trim points? Enter as   first, last   or just enter to proceed:  '
		if horiz_ok ne '' then begin
			nums = FIX(STRSPLIT(/extract, horiz_ok, ','))
			h1 = nums[0] > 0
			h2 = nums[1] < (numh - 1)
		endif
	endwhile
	horiz_phase = horiz_phase[h1:h2]
	horiz_off = horiz_off[h1:h2]
endif

CD, old_dir

if KEYWORD_SET(slit) then begin
	slitx = vert_result.slitx - MEDIAN(vert_result.slitx)
	slity = horiz_result.slity - MEDIAN(horiz_result.slity)
endif else begin
	slitx = vert_phase * 0.
	slity = horiz_phase * 0.
endelse

result = CREATE_STRUCT('vert_file', vert_file, 'horiz_file', horiz_file, $
	'vert_angle', vert_angle, 'description', description, $
	'tphase', FLTARR(npoints), 'xoff', FLTARR(npoints), 'yoff', FLTARR(npoints), $
	'h1', h1, 'h2', h2, 'v1', v1, 'v2', v2 )
		
; First, regrid and average everything onto a common time axis
halfbin = 1./(npoints * 2.)
tstep = ROUND(DOUBLE(pper) / npoints)
tgrid = FINDGEN(npoints) * tstep * 1000d
result.tphase = tgrid / 1000d / pper 
orbdat_out = REPLICATE({time : '', pztaoffset : '', pztboffset : '', pztcoffset : ''}, npoints)
gap = INTARR(npoints)
if KEYWORD_SET(smooth) then begin
	vert_off = SMOOTH(vert_off, smooth)
	horiz_off = SMOOTH(horiz_off, smooth)
endif

if KEYWORD_SET(pstop) then STOP

for i = 0, npoints - 1 do begin
	thisphase = result.tphase[i]
	vert_inpoints = WHERE(vert_phase gt thisphase-halfbin and $
		vert_phase le thisphase + halfbin, vert_numin)
	if vert_numin eq 1 then vert_inpoints = [vert_inpoints, vert_inpoints]
	if vert_numin gt 0 then begin
		this_voff = vert_off[vert_inpoints]
		this_vstdev = STDEV(this_voff)
		if this_vstdev gt 2 then begin
			PRINT, 'Too much variance in vertical offset...Auto-filtering'
			this_voff = this_voff[WHERE( ABS(this_voff - MEDIAN(this_voff)) lt 5)]
		endif
		this_vmean = MEAN(this_voff)
		result.xoff[i] = 0 - this_vmean
		if KEYWORD_SET(slit) then result.xoff[i] = result.xoff[i] + MEAN(slitx[vert_inpoints])
	endif else this_vstdev = 0
	horiz_inpoints = WHERE(horiz_phase gt thisphase-halfbin and $
		horiz_phase le thisphase + halfbin, horiz_numin)
	if horiz_numin eq 1 then horiz_inpoints = [horiz_inpoints, horiz_inpoints]
	if horiz_numin gt 0 then begin
		this_hoff = horiz_off[horiz_inpoints]
		this_hstdev = STDEV(this_hoff)
		if this_hstdev gt 2 then begin
			PRINT, 'Too much variance in horizontal offset...Auto-filtering'
			this_hoff = this_hoff[WHERE( ABS(this_hoff - MEDIAN(this_hoff)) lt 5)]
		endif
		this_hmean = MEAN(this_hoff)
		result.yoff[i] = 0 - this_hmean
		if KEYWORD_SET(slit) then result.xoff[i] = result.xoff[i] + MEAN(slity[horiz_inpoints])
	endif else this_hstdev = 0
	PRINT, thisphase, this_vstdev, this_hstdev, form = '(f12.3,2f12.2)'
	if (horiz_numin eq 0) or (vert_numin eq 0) then gap[i] = 1
endfor

if KEYWORD_SET(pstop) then STOP

; Now go back and clean up the results (interpolate across any gaps, shift to 0)
if TOTAL(gap) gt 0 then begin
	bads = WHERE(gap eq 1, numbad)
	goods = WHERE(gap eq 0, numgood)
	xgood = [result.xoff[goods], result.xoff[goods], result.xoff[goods]]
	ygood = [result.yoff[goods], result.yoff[goods], result.yoff[goods]]
	tgood = [result.tphase[goods]-1, result.tphase[goods], result.tphase[goods]+1]
	for i = 0, numbad - 1 do begin
		result.xoff[bads[i]] = INTERPOL(xgood, tgood, result.tphase[bads[i]], _extra = extra)
		result.yoff[bads[i]] = INTERPOL(ygood, tgood, result.tphase[bads[i]], _extra = extra)
	endfor
endif

; Shift everything to 0 at t0
xstart = result.xoff[0]
ystart = result.yoff[0]
result.xoff = result.xoff - xstart
result.yoff = result.yoff - ystart

; Apply the roll angle
xdata = 0 - (xstart + REFORM(vert_result.offset[0,*]))		;	Measurements used
ydata = 0 - (ystart + REFORM(horiz_result.offset[1,*]))
;if vert_angle ne 0 then begin
;	xin = result.xoff
;	yin = result.yoff
;	vert_rad = vert_angle * !pi / 180.
;	result.xoff = xin * COS(vert_rad) - yin * SIN(vert_rad)
;	result.yoff = xin * SIN(vert_rad) + yin * COS(vert_rad)
;	; Roll the measurements as well
;	tmp_xdata = xdata * COS(vert_rad) - ydata * SIN(vert_rad)
;	tmp_ydata = xdata * SIN(vert_rad) + ydata * COS(vert_rad)
;	xdata = tmp_xdata
;	ydata = tmp_ydata
;endif

if KEYWORD_SET(pstop) then STOP

; Now detrend the results to ensure that they don't jump on phase wrapping
xhop = INTERPOL(result.xoff, INDGEN(npoints), npoints)
yhop = INTERPOL(result.yoff, INDGEN(npoints), npoints)
result.xoff = result.xoff - xhop * result.tphase
result.yoff = result.yoff - yhop * result.tphase

; Now generate the PZT offsets
for i = 0, npoints - 1 do begin
	thispzt = IRIS_XY2PZT(result.xoff[i]*sref, result.yoff[i]*sref, /no)
	orbdat_out[i].time = STRING(tgrid[i], form = '(i10)')
	orbdat_out[i].pztaoffset = STRING(thispzt[0], form = '(i4)')
	orbdat_out[i].pztboffset = STRING(thispzt[1], form = '(i4)')
	orbdat_out[i].pztcoffset = STRING(thispzt[2], form = '(i4)')
endfor

; Then generate a file using the common time axis
orbhead_out = CREATE_STRUCT('ID', STRING(out_orbid, form = '(i05)'), $
	'SIZE', STRING(FIX(npoints*4 + 3)), 'NUMENTRIES', STRING(FIX(npoints)))
orbinfo_out = CREATE_STRUCT('Description', description)
orbstr_out = CREATE_STRUCT('info', orbinfo_out, 'header', orbhead_out, $
	'data', orbdat_out, 'xmldata', orbdat.xmldata)

; Write the file out and read it back in to use for plotting (or just plot using
; the data queued for output
if KEYWORD_SET(write) then begin
	out_file = CONCAT_DIR(opath, ttag + '_ORB-' + STRING(out_orbid, form = '(i05)') + '.xml' )
	dummy = SSW_IRIS_STRUCT2XML(orbstr_out, out_file = out_file)
	rresult = SSW_IRIS_XML2STRUCT(out_file)
endif else begin
	rresult = orbstr_out
endelse
nnpoints = N_ELEMENTS(rresult.data.time)
nxy = FLTARR(2,nnpoints)
for i = 0, nnpoints-1 do nxy[*,i] = IRIS_PZT2XY(rresult.data[i].pztaoffset, $
	rresult.data[i].pztboffset, rresult.data[i].pztcoffset, /no)/sref

; Generate some plots to check the results
; First, re-read the output xml file and compare it with the reference xml file
; In PZT space...
pztsym = [4,6,5]
cols = [255, 2, 6]
PLOT, rresult.data.time, rresult.data.pztaoffset, psym = 3, chars = 1.5, $
	xtitle = 'Relative time [msec]', ytitle = 'PZT offset [DN]', $
	yrange = LIMITS(FLOAT([rresult.data.pztaoffset, rresult.data.pztboffset, $
		rresult.data.pztcoffset, ref_pztoffs])), title = description + ' comparison with old wobble table : PZT ABC'
OPLOT, rresult.data.time, rresult.data.pztaoffset, psym = 3, col = 0
for i = 0, 2 do begin
	OPLOT, rresult.data.time, rresult.data.(i+1), psym = pztsym[i], thick = 2
	for j = 0, nref - 1 do begin
		this_orbdat = ref_orbdat.(j)
		OPLOT, this_orbdat.time, this_orbdat.(i+1), psym = pztsym[i], col = cols[j+1]
	endfor
endfor
LEGEND, pos = 10, chars = 1.5, 'PZT ' + ['A', 'B', 'C'], psym = pztsym
LEGEND, pos = 12, chars = 1.5, table_label, col = cols[0:nref], psym = 4+INTARR(nref+1), thick = [2, INTARR(nref)+1]
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-total_pzt_comp.png')
; And in XY space...
WDEF, win+1, 1200, 800
PLOT, rresult.data.time, nxy[0,*], psym = 4, chars = 1.5, $
	xtitle = 'Relative time [msec]', ytitle = 'Image offset [IRIS pixels]', $
	yrange = LIMITS([nxy, ref_xyoffs]), title = description + ' comparison with old wobble table : XY'
for j = 0, N_ELEMENTS(ref_orbid) - 1 do begin
	this_orbdat = ref_orbdat.(j)
	OPLOT, this_orbdat.time, this_orbdat.oxy[0,*], psym = 4, col = cols[j+1]
	OPLOT, this_orbdat.time, this_orbdat.oxy[1,*], psym = 6, col = cols[j+1]
endfor
OPLOT, rresult.data.time, nxy[1,*], psym = 6, thick = 2
OPLOT, rresult.data.time, nxy[0,*], psym = 4, thick = 2
LEGEND, pos = 10, chars = 1.5, ['X', 'Y'], psym = [4, 6]
LEGEND, pos = 12, chars = 1.5, table_label, col = cols[0:nref], psym = 4+INTARR(nref+1), thick = [2, INTARR(nref)+1]
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-total_xy_comp.png')

; Then, compare the output data with the raw analysis results it is based on
WDEF, win+2, 1200, 800
PLOT, result.tphase, result.xoff, psym = 4, xtit = 'Orbital phase', chars = 1.5, $
	ytit = 'Offset [IRIS pix]', title = description + ' comparison with raw data', $
	yrange = LIMITS([result.xoff, result.yoff, xdata, ydata])
OPLOT, vert_result.tphase, xdata, psym = 4, col = 6
OPLOT, horiz_result.tphase, ydata, psym = 6, col = 6
if KEYWORD_SET(check) then begin
	OPLOT, vert_result.tphase[v1:v2], xdata[v1:v2], psym = 4, col = 5
	OPLOT, horiz_result.tphase[h1:h2], ydata[h1:h2], psym = 6, col = 5
	LEGEND, pos = 12, chars = 1.5, ['Data Used', 'Data Excluded'], $
		psym = [6, 6], col = [5, 6], thick = [2, 2]
endif
OPLOT, result.tphase, result.yoff, psym = 6, thick = 2
OPLOT, result.tphase, result.xoff, psym = 4, thick = 2
LEGEND, pos = 10, chars = 1.5, ['Combined X', 'Combined Y', 'Vert X', 'Horiz Y'], $
	psym = [4, 6, 4, 6], col = [255, 255, 6, 6], thick = [2, 2, 1, 1]
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-total_smoothcomp.png')
; And compare across the orbital period boundary...
WDEF, win+3, 1200, 800
PLOT, ((result.tphase+0.5) mod 1) - 0.5, result.xoff, psym = 4, xtit = 'Orbital phase', chars = 1.5, $
	ytit = 'Offset [IRIS pix]', title = description + ' comparison with raw data (check period wrap)', $
	yrange = LIMITS([result.xoff, result.yoff]) + [-5, 5]
OPLOT, ((vert_result.tphase+0.5) mod 1) - 0.5, xdata, psym = 4, col = 6
OPLOT, ((horiz_result.tphase+0.5) mod 1) - 0.5, ydata, psym = 6, col = 6
if KEYWORD_SET(check) then begin
	OPLOT, ((vert_result.tphase[v1:v2]+0.5) mod 1) - 0.5, xdata[v1:v2], psym = 4, col = 5
	OPLOT, ((horiz_result.tphase[h1:h2]+0.5) mod 1) - 0.5, ydata[h1:h2], psym = 6, col = 5
	LEGEND, pos = 12, chars = 1.5, ['Data Used', 'Data Excluded'], $
		psym = [6, 6], col = [5, 6], thick = [2, 2]
endif
OPLOT, ((result.tphase+0.5) mod 1) - 0.5, result.yoff, psym = 6, thick = 2
OPLOT, ((result.tphase+0.5) mod 1) - 0.5, result.xoff, psym = 4, thick = 2
LEGEND, pos = 10, chars = 1.5, ['Combined X', 'Combined Y', 'Vert X', 'Horiz Y'], $
	psym = [4, 6, 4, 6], col = [255, 255, 6, 6], thick = [2, 2, 1, 1]
if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(opath, ttag + '-total_wrap.png')

if KEYWORD_SET(pstop) then STOP

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

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

end