pro IRIS_PREP_TREND_FID, $
	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_FID (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_FID) 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 1)
;	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 = 100
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-dimension sine wave fits on the fiducial
; positions into a structure array (obsdata)
;-------------------------------------------------
RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_filter_fid_save.genx'), fidfits
numobs = N_ELEMENTS(fidfits)
all_obsdata = IRIS_READ_META()
for i = 0, numobs - 1 do begin
	metamatch = WHERE(all_obsdata.path eq fidfits[i].filename, nummatch)
	if nummatch gt 0 then begin
		this_obsdata = CREATE_STRUCT(fidfits[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 
; fiducial finding into a structure array (imgdata)
;-------------------------------------------------
RESTGEN, file = CONCAT_DIR(outpath, 'iris_prep_all_fid_save.genx'), imgdata
imgtags = TAG_NAMES(imgdata)
xtag = WHERE(imgtags eq 'X1')

if KEYWORD_SET(loud) then PRINT, 'IRIS_PREP_TREND_FID 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 fiducial positions
;-------------------------------------------------
pdb = IRIS_MK_POINTDB()
nomy1 = (pdb.cpx2_1p5 - 1 - 270)
nomy2 = (pdb.cpx2_1p5 - 1 + 270)
nomx = (pdb.cpx1_1p5 - 1)

;+++++++++++++++++++++++++++++++++++++++++++++++++
; Generate a best-overall sine wave by looking at
; only the best-fit OBS and locally smoothing them
;-------------------------------------------------
swaves = STRTRIM([1330, 1400, 2796, 2832], 2)
chans = ['SJI_' + swaves, ['NUV', 'FUVS', 'FUVL']]
axes = ['X', 'Y']
fidnames = obsdata[0].fidnames
params = obsdata[0].paramnames
cols = [255, 2, 6, 5]

for i = 0, N_ELEMENTS(chans) - 1 do begin			;	Loop through img_path

	; First just plot the raw (per-image) fit results for this channel
	WDEF, (win+3) mod 32, 1000, 800
	thistype = WHERE(imgdata.type eq chans[i], numtype)
	if STRMID(chans[i], 0, 3) eq 'SJI' then !p.multi = [0,1,2] else !p.multi = 0
	UTPLOT, imgdata[thistype].t_obs, imgdata[thistype].sumspat * imgdata[thistype].(xtag+1) - nomy1, $
		psym = 3, chars = 1.5, /xstyle, ytitle = 'Pixel [Y] from Nominal', $
		title = chans[i] + ' Y Fiducials', yrange = [-5, 5], /ystyle
		timerange = [imgdata[thistype[0]].t_obs, imgdata[thistype[-1]].t_obs]
	OUTPLOT, imgdata[thistype].t_obs, imgdata[thistype].sumspat * imgdata[thistype].(xtag+3) - nomy2, $
		psym = 3, col = 2
	LEGEND, pos = 10, psym = 3, col = [255, 2], ['Y1', 'Y2'], char = 2
	if STRMID(chans[i], 0, 3) eq 'SJI' then begin
		UTPLOT, imgdata[thistype].t_obs, imgdata[thistype].sumspat * imgdata[thistype].(xtag) - nomx, $
			psym = 3, chars = 1.5, /xstyle, ytitle = 'Pixel [X]', title = chans[i] + ' X Fiducials', $
			timerange = [imgdata[thistype[0]].t_obs, imgdata[thistype[-1]].t_obs], $
			yrange = [-5, 5], /ystyle
		OUTPLOT, imgdata[thistype].t_obs, imgdata[thistype].sumspat * imgdata[thistype].(xtag+2) - nomx, $
			psym = 3, col = 2
		LEGEND, pos = 10, psym = 3, col = [255, 2], ['X1', 'X2'], char = 2				
	endif
	if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_fid_data.png')

	; Generate a structure array (fiddat) containing the results of the good fits
	; on the statistically significant OBS for this img_path
	for j = 0, N_ELEMENTS(fidnames) - 1 do begin		;	Loop through fiducial dimension
		numtag = WHERE(obstag eq chans[i] + '_' + fidnames[j] + '_NUM')
		chisqtag = WHERE(obstag eq chans[i] + '_' + fidnames[j] + '_CHISQ')
		partag = WHERE(obstag eq chans[i] + '_' + fidnames[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	; Loop through usable OBS
			thisdata = CREATE_STRUCT('date_obs', obsdata[goodobs[k]].date_obs, $
				'imgpath', chans[i], $
				'name', fidnames[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(fiddat) eq 0 then fiddat = thisdata else fiddat = [fiddat, thisdata]
		endfor else goto, skipchan
	endfor
	numpoints = N_ELEMENTS(fiddat)
	
	; 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
		for k = 0, 1 do begin
			useobs = WHERE(ABS(ANYTIM2TAI(fiddat.date_obs) - tais[j]) le tairange $
				and STRMID(fiddat.name, 0, 1) eq axes[k], numuse)
			if numuse gt 1 then usespec = fiddat[useobs] else begin
				usespec = CREATE_STRUCT('amplitude', [0d], 'phase', [0d], $
					'offset', [0d])
			endelse
			if k eq 0 then begin
				thissmooth = CREATE_STRUCT('tai', tais[j], 'date', days[j], $
					'amplitudex', MEDIAN(usespec.amplitude), $
					'phasex', MEDIAN(usespec.phase), $
					'offsetx', MEDIAN(usespec.offset) )
			endif else begin
				thissmooth = CREATE_STRUCT(thissmooth, $
					'amplitudey', MEDIAN(usespec.amplitude), $
					'phasey', MEDIAN(usespec.phase), $
					'offsety', MEDIAN(usespec.offset) )
			endelse
		endfor
		if N_TAGS(smoothstr) eq 0 then smoothstr = thissmooth $
			else smoothstr = [smoothstr, thissmooth]
	endfor
		
	;+++++++++++++++++++++++++++++++++++++++++++++++++
	; Generate plots
	;-------------------------------------------------		
	if KEYWORD_SET(loud) then begin

		; Plot amplitude
		WDEF, win, 900, 1000
		!p.multi = [0,1,2]
		for j = 0, 1 do begin
			ax1 = axes[j] + '1'
			rdat = fiddat[WHERE(fiddat.name eq ax1, numr)]
			UTPLOT, rdat.date_obs, rdat.amplitude, psym = 4, chars = 1.5, $
				ytitle = 'Amplitude [pixels]', title = chans[i] + ' Amplitude' + axes[j], $
				yrange = [0, 2]
			ax2 = axes[j] + '2'
			rdat = fiddat[WHERE(fiddat.name eq ax2, numr)]
			OUTPLOT, rdat.date_obs, rdat.amplitude, psym = 4, col = 2
			OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.(2+3*j), thick = 2, col = 6, line = 2, psym = -3
			LEGEND, pos = 10, chars = 1.5, axes[j] + ['1', '2'], col = [255,2], psym = 4
		endfor
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_fid_amplitude.png')

		; Plot offset
		WDEF, win+1, 900, 1000
		for j = 0, 1 do begin
			ax1 = axes[j] + '1'
			rdat = fiddat[WHERE(fiddat.name eq ax1, numr)]
			UTPLOT, rdat.date_obs, rdat.offset, psym = 4, chars = 1.5, $
				ytitle = 'Offset [pixels]', title = chans[i] + ' Offset' + axes[j], $
				yrange = [-8,8]
			ax2 = axes[j] + '2'
			rdat = fiddat[WHERE(fiddat.name eq ax2, numr)]
			OUTPLOT, rdat.date_obs, rdat.offset, psym = 4, col = 2
			OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.(4+3*j), thick = 2, col = 6, line = 2, psym = -3
			LEGEND, pos = 10, chars = 1.5, axes[j] + ['1', '2'], col = [255,2], psym = 4
		endfor
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_fid_offset.png')

		; Plot phase
		WDEF, win+2, 900, 1000
		!p.multi = [0, 1, 2]
		for j = 0, 1 do begin
			ax1 = axes[j] + '1'
			rdat = fiddat[WHERE(fiddat.name eq ax1, numr)]
			UTPLOT, rdat.date_obs, rdat.phase, psym = 4, chars = 1.5, $
				ytitle = 'Phase [radians]', title = chans[i] + ' Phase' + axes[j], $
				yrange = !pi * [0, 2]
			ax2 = axes[j] + '2'
			rdat = fiddat[WHERE(fiddat.name eq ax2, numr)]
			OUTPLOT, rdat.date_obs, rdat.phase, psym = 4, col = 2
			OUTPLOT, FILE2TIME(smoothstr.date), smoothstr.(3+3*j), thick = 2, col = 6, line = 2
			LEGEND, pos = 10, chars = 1.5, axes[j] + ['1', '2'], col = [255,2], psym = 4
		endfor
		if KEYWORD_SET(png) then PB_WIN2PNG, CONCAT_DIR(outpath, chans[i] + '_fid_phase.png')

	endif

	if i eq 0 then begin
		result = CREATE_STRUCT(chans[i], fiddat)
		smooths = CREATE_STRUCT(chans[i], smoothstr)
	endif else begin
		result = CREATE_STRUCT(result, chans[i], fiddat)
		smooths = CREATE_STRUCT(smooths, chans[i], smoothstr)
	endelse
		
	fiddat = 0
	smoothstr = 0
	skipchan:
endfor

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

if KEYWORD_SET(makepngs) then begin
	WDEF, win, 1000, 900
	!p.multi = [0,1,2]
	imtai = ANYTIM2TAI(imgdata.t_obs)

	CD, inpath, current = old_dir
	ff = FILE_SEARCH('*/*/*/*/fid/*.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(chans)-1 do begin		;		Loop through IMG_PATH
			usefiles = WHERE(imtai gt tais[i] and imtai lt tais[i+1] and $
				imgdata.type eq chans[j], numuse)
			if numuse gt 1 then begin
				this_smooth = smooths.(j)
				; Generate a plot for this OBS and wavelength channel			
				plotname = years[i] + '/' + months[i] + '/' + days[i] + '/' + $
					dates[i] + '/fid/' + dates[i] + '_' + chans[j] + '.png'
				usedat = imgdata[usefiles]
				for k = 0, 1 do begin		;	Loop through X/Y
					if k eq 0 then axlab = 'X' else axlab = 'Y'
					reset_flag = 0
					for m = 0, 1 do begin	;	Loop through fiducial 1/2
						goodfid = WHERE(usedat.(xtag + 4 + k + m*2) le chisqlim, numfin)
						if numfin gt 0 then begin
							if reset_flag eq 0 then begin
								flatt = usedat[goodfid].t_obs
								flatfid = usedat[goodfid].(xtag + k + m*2)
							endif else begin
								flatt = [flatt, usedat[goodfid].t_obs]
								flatfid = [flatfid, usedat[goodfid].(xtag + k + m*2)]
							endelse
							reset_flag = 1
							if k eq 1 and m eq 0 then flatfid = flatfid + 541
						endif
					endfor
					if N_ELEMENTS(flatt) lt 2 then begin
						UTPLOT, usedat.t_obs, FLTARR(numuse), ytitle = 'Pixel location [' + axlab + ']', $
							chars = 1.5, psym = 4, title = plotname, yrange = [-5, 5], /ystyle
						goto, skipobs
					endif 
					UTPLOT, flatt, flatfid, /ynoz, ytitle = 'Pixel location [' + axlab + ']', $
						chars = 1.5, psym = 4, title = plotname, $
						yrange = MEDIAN(flatfid) + [-5, 5], /ystyle
					for m = 0, 1 do OUTPLOT, usedat.t_obs, usedat.(xtag+m*2), psym = 4, col = cols[m]
				
					; Compute the smoothed sine wave for this wavelength and time interval
					fparams = this_smooth[WHERE(this_smooth.date eq TIME2FILE(/date, flatt[0]))]
					if k eq 0 then begin
						thisphase = fparams.phasex
						thisoffset = fparams.offsetx
						thisamplitude = fparams.amplitudex
					endif else begin
						thisphase = fparams.phasey
						thisoffset = fparams.offsety
						thisamplitude = fparams.amplitudey
					endelse
					
					; Look up the per-range fits and note the average pixel deviation
					fitind = (WHERE(obsdata.filename eq dates[i]))[0]
					chisqstrings = STRARR(2)
					for m = 0, 1 do begin
						thischisq = obsdata[fitind].(WHERE(obstag eq chans[j] + '_X1_CHISQ') + (k + m*2) * 3)
						if thischisq eq 0 then chisqstrings[m] = 'NaN' else $
							chisqstrings[m] = STRING(thischisq, '(f6.2)')
						chisqstrings[m] = fidnames[k+m*2] + ' : ' + chisqstrings[m]
					endfor

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

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

end