function IRIS_PREP_WAVECORR_L2, files, $
	pstop = pstop, showfit = showfit, nosmooth = nosmooth, nosine = nosine

;+
;
; Given some IRIS L2 raster files, measures the position of neutral lines in the 
; NUV and FUV and produces recommended wavelength corrections as a function of time 
; Note that the NUV correction is based on the Ni I line at 2799.474, and the FUV
; correction is based on the O I line at 1355.6, even though other lines are 
; measured and reported.
;
; INPUTS:
;	files	-	string giving path and filename to IRIS L2 raster files. If
;				multiple files are specified, they must all come from the same OBS
;
; RETURNS:
;	result	-	a structure containing some identifying information along with 
;				the results of the fit
;				Use result.corr_tai, result.corr_nuv and result.corr_fuv for
;				correction unless you have reason to do otherwise.
;
; KEYWORDS:
;	/pstop	-	(SWITCH) stop at the end for debugging. Set to 2 to stop in the loop
;	/showfit -	(SWITCH) set to plot the results at the end
;	/nosmooth - (SWITCH) set to avoid smoothing the measurements for the recommended
;				correction. Seems like a bad idea, but may be necessary in some cases.
;				The smoothing first removes outliers from the fitting.
;	/nosine	-	(SWITCH) set to avoid fitting a sine wave for the recommended correction
;				If this is set, then the recommendation is just interpolated through
;				the smoothed version of the data 
;
;-

nfile = N_ELEMENTS(files)
if N_ELEMENTS(pstop) eq 0 then pstop = 0

; Define the lines that will be fit
lname 	=	['Ni I', 	'Mn I', 	'Fe I', 	'O I', 		'Fe II']
wave0s	=	[2799.474,	2801.902,	2805.346,	1355.60,	1392.82]
wmins	=	[2799.3, 	2801.6, 	2805.1, 	1355.4, 	1392.6]
wmaxs	=	[2799.8, 	2802.4, 	2805.7, 	1355.9, 	1393.1]
lmins	=	[5, 		5, 			5, 			0.5, 		0.5]
lsign 	=	[-1, 		-1, 		-1, 		1,			1]
nline = N_ELEMENTS(lname)

; Loop through the files 
for ii = 0, nfile-1 do begin
	l2f = files[ii]
	d = OBJ_NEW('iris_data')
	d -> READ, l2f
	index = d->GETHDR(/struct)

	solarx	=	d->GETXPOS()
	xfac 	= 	[MEAN(solarx), MEAN(DERIV(solarx))]
	nx 		=	N_ELEMENTS(solarx)
	solary	=	d->GETYPOS()
	yfac 	=	[MEAN(solary), MEAN(DERIV(solary))]
	ny 		=	N_ELEMENTS(solary)
	t		=	d->TI2UTC()
	tai		=	ANYTIM2TAI(t)
	nt		=	N_ELEMENTS(t)

	; Initialize some stuff on the first pass
	if ii eq 0 then begin 
		coord = FLTARR(nfile, nx, ny, 2)	;	X and Y coordinates for each image
		corrs = DBLARR(nfile, nt, nline) - !values.d_nan
		chisq = DBLARR(nfile, nt, nline) 
		times = STRARR(nfile, nt)
		tais  = DBLARR(nfile, nt)
		wins  = INTARR(nline) - 1
		for jj = 0, nline-1 do wins[jj] = d -> GETWINDX(wave0s[jj])
	endif

	oldwin = -1
	; Loop through the neutral(-ish) lines that are in this raster
	for jj = 0, nline - 1 do if wins[jj] ge 0 then begin
		winidx= wins[jj]

		; Read the spectral window if needed
		if winidx ne oldwin then begin
			corrlam	=	d->GETLAM(winidx)				;	wavelength grid
			corrwin	=	d->GETVAR(winidx, /load)		;	data array
			corrwin = 	TRANSPOSE(corrwin, [0, 2, 1])	;	flipped time/vertical				
			badpix	= WHERE(FINITE(corrwin) eq 0 or corrwin eq -200, bpc)
			if bpc gt 0 then corrwin[badpix] = 0.
			oldwin = winidx
		endif
		sub_line	=	WHERE(corrlam ge wmins[jj] and corrlam le wmaxs[jj], n_sub)
		if n_sub lt 5 then sub_line = INDGEN(5) + sub_line[0]
		if n_sub lt 7 then nterms = 4 else nterms = 5
		scorrlam = DOUBLE(corrlam[sub_line])
		scorrwin = DOUBLE(corrwin[sub_line, *, *])
		nslitpos = (SIZE(scorrwin))[2]	;	Can't use nx because of roll=90 cases
		if pstop gt 1 then STOP

		; Loop through raster positions and fit line at each one
		for pp = 0, nslitpos-1 do begin
			corrprof = TOTAL(REFORM(scorrwin[*,pp,*]>0d),2) / FLOAT(ny)	;	Mean spectrum along slit 
			if MEAN(corrprof) gt lmins[jj] then begin						;	not too dim
;				errprof = SQRT(corrprof / MAX(corrprof) * 400.)		;	Error is 5% of peak
;				ggg = MPFITPEAK(scorrlam, corrprof, bbb, nterms = nterms, /double, error = errprof, chisq = thisq)
				ggg = SILENT_GAUSSFIT(scorrlam, corrprof, bbb, nterms = nterms, chisq = thisq)
				abscor = wave0s[jj] - bbb[1]				;	Difference in line position
				; Make sure the fit is in range and the line goes the right way
				if ABS(abscor) le (wmaxs[jj]-wmins[jj])/4. and (bbb[0]*lsign[jj]) gt 0. then begin
					corrs[ii, pp, jj] = abscor
					chisq[ii, pp, jj] = thisq
				endif
			endif
		endfor
	endif
	
	; Populate the ancillary arrays with time and solar X/Y
	times[ii,*]	= t
	tais[ii,*]	= tai
	for pp = 0, ny-1 do coord[ii, *, pp, 0]	= solarx
	for pp = 0, nx-1 do coord[ii, pp, *, 1]	= solary
	
endfor

; Smooth the results of the Ni I and O I fits to obtain recommended 
; corrections for the NUV and FUV as a function of time (TAI)
corr_tai = REFORM(tais, nfile * nt)
tsort = SORT(corr_tai)
corr_tai = corr_tai[tsort]
corr_nuv = (REFORM(corrs[*, *, 0], nfile * nt))[tsort]
corr_fuv = (REFORM(corrs[*, *, 3], nfile * nt))[tsort]
if not KEYWORD_SET(nosmooth) then begin
	; First remove outliers
	nuv_outlier = WHERE( (corr_nuv - MEDIAN(corr_nuv)) gt 0.08, $
		numnout, comp = ngood, ncomp = nngood)
	if numnout gt 0 then corr_nuv[nuv_outlier] = !values.f_nan
	fuv_outlier = WHERE( (corr_fuv - MEDIAN(corr_fuv)) gt 0.05, $
		numfout, comp = fgood, ncomp = nfgood)
	if numfout gt 0 then corr_fuv[fuv_outlier] = !values.f_nan
	; Smooth over a 5 minute sliding window to remove oscillations
	swidth = FLOOR(5 * 60./MEAN(DERIV(corr_tai)))
	if swidth lt nngood then corr_nuv = SMOOTH(corr_nuv, swidth, /edge_truncate, /nan)
	if swidth lt nfgood then corr_fuv = SMOOTH(corr_fuv, swidth, /edge_truncate, /nan)
	if not KEYWORD_SET(nosine) then begin
		; Sine fit the results
		trange = MAX(corr_tai) - MIN(corr_tai)
		orbits = FLOOR(trange / 5856)
		numparams = 4 + (orbits < 3)
		params = ['Amplitude', 'Frequency', 'Phase', 'Offset', 'Linear', 'Squared', 'Cubed', 'Quartic', 'Quintic']
		values = [0.01d, (2 * !pi) / 5856d, 1d, 0.1d, 1d, 0.1d, 0.01d, 0.001d, 0.0001d]
		fixed = [0, 1, 0, 0, 0, 0, 0, 0, 0, 0]
		parinfo = REPLICATE({name:params[0], value:values[0], fixed:fixed[0]}, numparams)
		for i = 1, numparams-1 do begin
			parinfo[i].name 	= params[i]
			parinfo[i].value 	= values[i]
			parinfo[i].fixed 	= fixed[i]
		endfor
		initpars = parinfo
		nsineparams = MPFITFUN('IRIS_MP_SINECHISQ', corr_tai[ngood] - corr_tai[0], $
			corr_nuv[ngood], corr_nuv[ngood] * 0. + 1, /quiet, parinfo = parinfo)
		corr_nuv = IRIS_MP_SINECHISQ(corr_tai - corr_tai[0], nsineparams)
		fsineparams = MPFITFUN('IRIS_MP_SINECHISQ', corr_tai[fgood] - corr_tai[0], $
			corr_fuv[fgood], corr_fuv[fgood] * 0. + 1, /quiet, parinfo = parinfo)
		corr_fuv = IRIS_MP_SINECHISQ(corr_tai - corr_tai[0], fsineparams)
	endif else begin
		; Spline interpolate through the smoothed results
		corr_nuv = INTERPOL(corr_nuv[ngood], corr_tai[ngood], corr_tai, /nan, /spline)
		corr_fuv = INTERPOL(corr_fuv[fgood], corr_tai[fgood], corr_tai, /nan, /spline)
	endelse
endif 

; Plot the results if requested
if KEYWORD_SET(showfit) then begin
	oldpmulti = !p.multi
	oldpcol = !p.color
	!p.multi = [0, 1, 2]
	TVLCT, rr, gg, bb, /get
	LOADCT, 12, /silent
	bluecol		= 110
	redcol 		= 200
	greencol	= 50
	forecol 	= 255
	!p.color 	= forecol	;	UTPLOT requires this for the axis...yuck
	
	cnuv = corrs[*,*,0:2]
	cnfin = WHERE(FINITE(cnuv), numfin)
	if numfin gt 0 then begin
		yrange = LIMITS(cnuv[cnfin])
		UTPLOT, times, corrs[*,*,0], yrange = yrange, psym = 4, chars = 1.5, $
			ytitle = 'Wavelength error (' + STRING(197B) + ')', $
			title = 'NUV wavelength correction', col = forecol
		OUTPLOT, times, corrs[*,*,1], psym = 4, col = greencol
		OUTPLOT, times, corrs[*,*,2], psym = 4, col = bluecol
		if not KEYWORD_SET(nosmooth) then $
			OUTPLOT, TAI2UTC(corr_tai), corr_nuv, col = redcol, thick = 2, line = 2
		AL_LEGEND, lname[0:2], col = [forecol, greencol, bluecol], psym = [4,4,4], $
			/top, /left, chars = 1.5
	endif

	cfuv = corrs[*,*,3:4]
	cffin = WHERE(FINITE(cfuv), numfin)
	if numfin gt 0 then begin
		yrange = LIMITS(cfuv[cffin])
		UTPLOT, times, corrs[*,*,3], yrange = yrange, psym = 4, chars = 1.5, $
			ytitle = 'Wavelength error (' + STRING(197B) + ')', $
			title = 'FUV wavelength correction', col = forecol
		OUTPLOT, times, corrs[*,*,4], psym = 4, col = greencol
		if not KEYWORD_SET(nosmooth) then $
			OUTPLOT, TAI2UTC(corr_tai), corr_fuv, col = redcol, thick = 2, line = 2
		AL_LEGEND, lname[3:4], col = [forecol, greencol], psym = [4,4], $
			/top, /left, chars = 1.5
	endif
	
	!p.multi = oldpmulti
	!p.color = oldpcol
	TVLCT, rr, gg, bb
endif

; Generate the output array
if KEYWORD_SET(pstop) then STOP
note = 'corrs[file, rasterstep, line] gives target_wave - measured_wave in angstroms'
result = {note : note, files : files, lname : lname, wave0s : wave0s, $
	times : times, tais : tais, nx : nx, ny : ny, coord : coord, $
	corrs : corrs, chisq : chisq, $
	corr_tai : corr_tai, corr_nuv : corr_nuv, corr_fuv : corr_fuv}
RETURN, result

end