pro IRIS_PREP_TREND_WAVE, $
	win = win, png = png, inpath = inpath, daymedian = daymedian, $
	pstop = pstop, loud = loud, makepngs = makepngs, $
	save = save, outpath = outpath, chisqlim = chisqlim, numlim = numlim
	
;+
;
; Read the output IRIS_PREP_FILTER_WAVE (sine wave fits to each OBS), and 
; then plots the sine parameters to look for a global best-fit sine fit
;
; INPUTS:
;
; OUTPUT:
;
; KEYWORDS:
;	inpath	-	Path to where the individual OBS directories are (defaults to
;		/irisa/data/prep
;	outpath	-	Path to where the trend results (including the output of 
;		IRIS_PREP_FILTER_WAVE) are saved; used for input and output
;	/save	-	If set, then save a genx file with the parameters of the global fit
;		to the outpath
;
;	/pstop	-	standard debug switch
;	/loud	-	If set, then plots the sine parameters from each obs
;	win	-	Optional, starting window index for plots (defaults to 10)
;	/png	-	If set, then the plots are saved
;	daymedian	-	Width of the smoothing to be performed on sine parameters,
;		in days. Defaults to 8.
;	/makepngs - If set, then go back and re-plot each OBS with the global
;		average sine wave on top
;	chisqlim -	The highest allowable chi-squared for the fit to an OBS (defaults
;		to 20)
;	numlim	-	The minimum number of images in an OBS to include in the global
;		fit (defaults to 100)
;
;-

tt0 = SYSTIME(/sec)
if N_ELEMENTS(outpath) eq 0 then outpath = '/net/xema/Volumes/disk2/data/iris/iris_prep/trend/'
if N_ELEMENTS(inpath) eq 0 then inpath = '/irisa/data/prep/'
if N_ELEMENTS(chisqlim) eq 0 then chisqlim = 20
if N_ELEMENTS(numlim) eq 0 then numlim = 100
if N_ELEMENTS(win) eq 0 then win = 10

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Read the results of the per-OBS, per-channel
; and per-spatial range sine wave fits on the
; line positions into a structure array (obsdata)
;-------------------------------------------------
RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_filter_wave_save.genx'), wavefits
numobs = N_ELEMENTS(wavefits)
all_obsdata = IRIS_READ_META()
for i = 0, numobs - 1 do begin
	metamatch = WHERE(all_obsdata.path eq wavefits[i].filename, nummatch)
	if nummatch gt 0 then begin
		this_obsdata = CREATE_STRUCT(wavefits[i], all_obsdata[metamatch])
		if N_ELEMENTS(obsdata) eq 0 then obsdata = this_obsdata $
			else obsdata = [obsdata, this_obsdata]
	endif
endfor
taistart = ANYTIM2TAI(FILE2TIME(obsdata.path))
obstag = TAG_NAMES(obsdata)

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Read the results of the individual image wave
; finding into a structure array (imgdata)
;-------------------------------------------------
RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_all_wave_save.genx'), imgdata
imgtags = TAG_NAMES(imgdata)

if KEYWORD_SET(loud) then PRINT, 'IRIS_PREP_TREND_WAVE loaded data. Elapsed time : ', $
	STRING(SYSTIME(/sec) - tt0, form = '(f10.1)')
if KEYWORD_SET(pstop) then STOP

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Set up plotting
;-------------------------------------------------
if KEYWORD_SET(loud) then begin
	if N_ELEMENTS(win) eq 0 then win = 10
	TVLCT, rr, gg, bb, /get
	PB_SET_LINE_COLOR
	oldpmulti = !p.multi
endif

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Load nominal wavelength positions
;-------------------------------------------------
pdb = IRIS_MK_POINTDB()
rpix = [pdb.cpx1_nuv, pdb.cpx1_fu1, pdb.cpx1_fu2] - 1
rpix[2] = (1392.82 - pdb.cvl1_fu2) / pdb.cdlt1_f2 + pdb.cpx1_fu2 - 1

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Generate a best-overall sine wave by looking at
; only the best-fit OBS and locally smoothing them
;-------------------------------------------------
specs = ['NUV', 'FUVS', 'FUVL']
ranges = obsdata[0].rangenames
params = obsdata[0].paramnames
cols = [255, 2, 6]
for i = 0, N_ELEMENTS(specs) - 1 do begin		;	Loop through img_path (NUV/FUVS/FUVL)

	; First just plot the raw (per-image) fit results for this channel
	WDEF, (win+3) mod 32, 1000, 700
	thistype = WHERE(imgdata.img_path eq specs[i], numtype)
	UTPLOT, imgdata[thistype].t_obs, imgdata[thistype].line_top * imgdata[thistype].sumsptrl, $
		psym = 3, yrange = rpix[i] + [-5, 5] + [-3,3]*i, /xstyle, chars = 1.5, $
		title = specs[i] + ' Line Position', /ysty
	OUTPLOT, imgdata[thistype].t_obs, imgdata[thistype].line_mid * imgdata[thistype].sumsptrl, $
		psym = 3, col = 6
	OUTPLOT, imgdata[thistype].t_obs, imgdata[thistype].line_bot * imgdata[thistype].sumsptrl, $
		psym = 3, col = 2
	OUTPLOT, [imgdata[thistype[0]].t_obs, imgdata[thistype[-1]].t_obs], rpix[i] + [0,0], col = 5, line = 2, thick = 2
	LEGEND, char = 2, pos = 4, col = [255, 6, 2], ['TOP', 'MID', 'BOT']
	if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, specs[i] + '_wave_data.png')

	; Now generate the average fit
	for j = 0, N_ELEMENTS(ranges) - 1 do begin	;	Loop through spatial ranges
		numtag = WHERE(obstag eq specs[i] + '_' + ranges[j] + '_NUM')
		chisqtag = WHERE(obstag eq specs[i] + '_' + ranges[j] + '_CHISQ')
		partag = WHERE(obstag eq specs[i] + '_' + ranges[j] + '_PARAMS')
		goodobs = WHERE(obsdata.(numtag) ge numlim and obsdata.(chisqtag) lt chisqlim $
			and obsdata.(chisqtag) gt 0, numgood)
		if numgood gt 0 then for k = 0, numgood - 1 do begin
			thisdata = CREATE_STRUCT('date_obs', obsdata[goodobs[k]].date_obs, $
				'imgpath', specs[i], $
				'range', ranges[j], $
				'num', obsdata[goodobs[k]].(numtag), $
				'chisq', obsdata[goodobs[k]].(chisqtag))
			for m = 0, N_ELEMENTS(params) - 1 do begin
				thisdata = CREATE_STRUCT(thisdata, params[m], obsdata[goodobs[k]].(partag)[m])
			endfor
			if N_TAGS(specdat) eq 0 then specdat = thisdata else specdat = [specdat, thisdata]
		endfor
	endfor
	numpoints = N_ELEMENTS(specdat)
	
	; Generate running median (smoothed) corrections	
	if not KEYWORD_SET(daymedian) then daymedian = 8
	tairange = daymedian * 86400d
	numdays = (MAX(taistart) - MIN(taistart)) / 86400.
	tais = DINDGEN(numdays)*86400 + taistart[0]
	days = TIME2FILE(TAI2UTC(tais), /date)
	tais = ANYTIM2TAI(FILE2TIME(days))
	for j = 0, numdays - 1 do begin
		useobs = WHERE(ABS(ANYTIM2TAI(specdat.date_obs) - tais[j]) le tairange, numuse)
		if numuse gt 1 then begin
			usespec = specdat[useobs]
			thissmooth = CREATE_STRUCT('tai', tais[j], 'date', days[j], $
				'amplitude', MEDIAN(usespec.amplitude), $
				'phase', MEDIAN(usespec.phase), $
				'offset', MEDIAN(usespec.offset) )
			if N_TAGS(smoothstr) eq 0 then smoothstr = thissmooth $
				else smoothstr = [smoothstr, thissmooth]
		endif
	endfor
	
	if KEYWORD_SET(pstop) then STOP
	;+++++++++++++++++++++++++++++++++++++++++++++++++
	; Generate plots
	;-------------------------------------------------		
	if KEYWORD_SET(loud) then begin

		; Plot amplitude
		WDEF, win, 900, 600
		!p.multi = 0
		UTPLOT, specdat.date_obs, specdat.amplitude, psym = 4, chars = 1.5, $
			ytitle = 'Amplitude [pixels]', title = specs[i] + ' Amplitude', $
			yrange = [0, 3]
;		for j = 0, N_ELEMENTS(ranges) - 1 do begin
		for j = 1, N_ELEMENTS(ranges) - 1 do begin
			rdat = specdat[WHERE(specdat.range eq ranges[j])]
			OUTPLOT, rdat.date_obs, rdat.amplitude, col=cols[j], psym = 4
		endfor
		OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.amplitude, thick = 2, col = 7, line = 2
		if i gt 0 then begin
			OUTPLOT, FILE2TIME(smooths.nuv.date), smooths.nuv.amplitude, $
				thick = 2, col = 6, line = 3
			LEGEND, pos = 12, ['FUV trend', 'NUV trend'], col = [7, 6], chars = 1.5, linesty = [2,3]
		endif
		LEGEND, pos = 10, chars = 1.5, ranges, col = cols, psym = 4
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, specs[i] + '_wave_amplitude.png')
		
		; Plot offset
		WDEF, win+1, 900, 600
		pdb = IRIS_MK_POINTDB(version = pdb_ver)
		case i of
			0	:	xnom = pdb.cpx1_nuv - 1
			1	:	xnom = pdb.cpx1_fu1 - 1
;			2	:	xnom = pdb.cpx1_fu2 - 1
			2	:	xnom = 3023
		endcase
		UTPLOT, specdat.date_obs, specdat.offset, psym = 4, chars = 1.5, $
			ytitle = 'Offset [pixels]', title = specs[i] + ' average line position', $
			yrange = xnom + [-5,5]
		OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.offset, thick = 2, col = 7, line = 2
		if i gt 0 then begin
			OUTPLOT, FILE2TIME(smooths.nuv.date), thick = 2, col = 6, line = 3, $
				0 - smooths.nuv.offset + ( MEAN(smoothstr.offset) + MEAN(smooths.nuv.offset) )
			LEGEND, pos = 6, col = [7, 6], linesty = [2,3], $
				['FUV trend', 'NUV trend (sign flipped)'], chars = 1.5
		endif
		for j = 1, N_ELEMENTS(ranges) - 1 do begin
			rdat = specdat[WHERE(specdat.range eq ranges[j])]
			OUTPLOT, rdat.date_obs, rdat.offset, col=cols[j], psym = 4
		endfor
		OUTPLOT, specdat.date_obs, specdat.offset * 0. + xnom, col = 5, line = 1
		LEGEND, pos = 4, chars = 1.5, ranges, col = cols, psym = 4
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, specs[i] + '_wave_offset.png')

		; Plot phase
		WDEF, win+2, 900, 600
		UTPLOT, specdat.date_obs, (((specdat.phase/!pi)+1) mod 2) - 1, psym = 4, chars = 1.5, $
			ytitle = 'Phase [pi radians]', title = specs[i] + ' Phase', yrange = [-1, 1]
		OUTPLOT, FILE2TIME(smoothstr.date), (((smoothstr.phase/!pi)+1) mod 2) - 1, thick = 2, col = 7, line = 2
		if i gt 0 then begin
			OUTPLOT, FILE2TIME(smooths.nuv.date), (((smooths.nuv.phase/!pi)+1) mod 2) - 1, $
				thick = 2, col = 6, line = 3
			LEGEND, pos = 12, ['FUV trend', 'NUV trend'], col = [7, 6], chars = 1.5, linesty = [2,3]
		endif			
		for j = 1, N_ELEMENTS(ranges) - 1 do begin
			rdat = specdat[WHERE(specdat.range eq ranges[j])]
			OUTPLOT, rdat.date_obs, (((rdat.phase/!pi)+1) mod 2) - 1, col=cols[j], psym = 4
		endfor
		LEGEND, pos = 6, chars = 1.5, ranges, col = cols, psym = 4
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, specs[i] + '_wave_phase.png')
			
	endif

	if i eq 0 then begin
		result = CREATE_STRUCT(specs[i], specdat)
		smooths = CREATE_STRUCT(specs[i], smoothstr)
	endif else begin
		result = CREATE_STRUCT(result, specs[i], specdat)
		smooths = CREATE_STRUCT(smooths, specs[i], smoothstr)
	endelse
		
	specdat = 0
	smoothstr = 0
endfor

if KEYWORD_SET(save) then begin
	output = CREATE_STRUCT('fitdat', result, 'smoothdat', smooths)
	SAVEGEN, file = CONCAT_DIR(outpath, 'wave_fit.genx'), str = output
endif

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Done generating overall best-fit sine wave
;-------------------------------------------------
if KEYWORD_SET(loud) then PRINT, 'IRIS_PREP_TREND_WAVE generated average sine. Elapsed time : ', $
	STRING(SYSTIME(/sec) - tt0, form = '(f10.1)')
if KEYWORD_SET(pstop) then STOP

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Now go back and plot the global best-fit wave 
; over each obs (relies on the genx file written
; by IRIS_PREP_FILTER_WAVE, /write_logdat)
;-------------------------------------------------
if KEYWORD_SET(makepngs) then begin
	WDEF, win, 900, 700
	!p.multi = 0
	bottag = WHERE(imgtags eq 'LINE_BOT')
	imtai = ANYTIM2TAI(imgdata.t_obs)

	CD, inpath, current = old_dir
	ff = FILE_SEARCH('*/*/*/*/wavecorr/*.txt') 
	numf = N_ELEMENTS(ff)
	dates = STRMID(ff, 11, 15)
	years = STRMID(dates, 0, 4)
	months = STRMID(dates, 4, 2)
	days = STRMID(dates, 6, 2)
	tais = ANYTIM2TAI(FILE2TIME(dates))
	for i = 0, numf - 2 do begin					;		Loop through OBS
		for j = 0, N_ELEMENTS(specs)-1 do begin		;		Loop through IMG_PATH
			usefiles = WHERE(imtai gt tais[i] and imtai lt tais[i+1] and $
				imgdata.img_path eq specs[j], numuse)
			if numuse gt 1 then begin
			
				; Generate a plot for this OBS and wavelength channel			
				plotname = years[i] + '/' + months[i] + '/' + days[i] + '/' + $
					dates[i] + '/wavecorr/' + dates[i] + '_' + specs[j] + '.png'
				usedat = imgdata[usefiles]
				thissum = usedat[0].sumsptrl
				reset_flag = 0
				for k = 0, N_ELEMENTS(ranges) - 1 do begin
;					goodrange = WHERE(FINITE(usedat.(bottag + 3 + k)), numfin)
					goodrange = WHERE(usedat.(bottag + 3 + k) le chisqlim, numfin)
					if numfin gt 0 then begin
						if reset_flag eq 0 then begin
							flatt = usedat[goodrange].t_obs
							flatwave = usedat[goodrange].(bottag + k) * thissum
						endif else begin
							flatt = [flatt, usedat[goodrange].t_obs]
							flatwave = [flatwave, usedat[goodrange].(bottag + k) * thissum] 
						endelse
						reset_flag = 1
					endif
				endfor
				if N_ELEMENTS(flatt) lt 2 then goto, skipobs
				
				UTPLOT, flatt, flatwave, /ynoz, ytitle = 'Pixel location', $
					chars = 1.5, psym = 4, title = plotname, $
					yrange = MEDIAN(flatwave) + [-8, 8], /ystyle
				for k = 0, 2 do	OUTPLOT, usedat.t_obs, usedat.(bottag+k), psym = 4, col = cols[k]
				
				; Compute the smoothed sine wave for this wavelength and time interval
				nparams = smooths.nuv[WHERE(smooths.nuv.date eq TIME2FILE(/date, flatt[0]))]
				fparams = smooths.fuvs[WHERE(smooths.fuvs.date eq TIME2FILE(/date, flatt[0]))]
				thisphase = (nparams.phase + imgdata[usefiles[0]].roll / 180. * !pi) mod (2 * !pi)
				case j of
					0	:	begin
						thisoffset = nparams.offset
						thisamp = nparams.amplitude
					end
					1	:	begin
						thisoffset = fparams.offset
						thisamp = 0 - nparams.amplitude
					end
					2	:	begin
						thisoffset = 3020
						thisamp = 0 - nparams.amplitude
					end
				endcase
				
				; Look up the per-range fits and note the average pixel deviation
				fitind = (WHERE(obsdata.filename eq dates[i]))[0]
				chisqstrings = STRARR(3)
				for k = 0, 2 do begin
					thischisq = obsdata[fitind].(WHERE(obstag eq specs[j] + '_BOT_CHISQ') + k*3)
					if thischisq eq 0 then chisqstrings[k] = 'NaN' else $
						chisqstrings[k] = STRING(thischisq, '(f6.2)')
					chisqstrings[k] = ranges[k] + ' : ' + chisqstrings[k]
				endfor

				; include the SINE wave fit
				if numuse gt 10 then begin
					xtai = ANYTIM2TAI(flatt)
					taigrid = MIN(xtai) + DINDGEN(MAX(xtai) - MIN(xtai))
					ysine = thisoffset + thisamp * SIN( ((2 * !pi) / 5856d * taigrid) + thisphase )
					OUTPLOT, TAI2UTC(taigrid), ysine, thick = 2
					sinepred = INTERPOL(ysine, taigrid, ANYTIM2TAI(flatt))
					fitdev = MEAN(ABS(flatwave - sinepred), /nan)
					chisqstrings = ['Fit : ' + STRING(fitdev, form = '(f6.2)'), chisqstrings]
					chisqpsym = [-3,4,4,4]
					ccol = [255, cols]
					clin = [0,1,1,1]
					for k = 0, 2 do begin
						thisparams = obsdata[fitind].(WHERE(obstag eq specs[j] + '_BOT_PARAMS') + k*3)
						ysine = thisparams[0] * SIN( ((2 * !pi) / 5856d * taigrid) + thisparams[2] )
						ysine = ABS(thisparams[3]) + ysine
						OUTPLOT, TAI2UTC(taigrid), ysine, col = cols[k], line = 1
						if KEYWORD_SET(pstop) then STOP
					endfor
				endif else begin
					cpsym = [-4,-4,-4]
					ccol = cols
					clin = [1,1,1]
				endelse
				
				LEGEND, pos = 6, chisqstrings, col = ccol, chars = 1.5, psym = cpsym
				PB_WIN2PNG, plotname
				skipobs : if KEYWORD_SET(pstop) then STOP

			endif
		endfor
	endfor
		
	CD, old_dir
endif

if KEYWORD_SET(loud) then begin
	TVLCT, rr, gg, bb
	!p.multi = oldpmulti
	PRINT, 'IRIS_PREP_TREND_WAVE elapsed time : ', STRING(SYSTIME(/sec) - tt0, form = '(f10.1)')
endif

if KEYWORD_SET(pstop) then STOP

end