function IRIS_PREP_WAVECORR, hdr, dat, msg, $
	angstroms_per_pixel = angstroms_per_pixel, $
	logdir = logdir, showfit = showfit, fuvs = fuvs, fuvl = fuvl, $
	pstop = pstop

;+
;          PURPOSE : Corrections for orbital variation of the spectral line positions
;
; CALLING SEQUENCE : result = iris_prep_wavecorr(hdr, dat) 
; 
;           INPUTS : hdr, dat	-	FITS header and data array for an IRIS NUV SG image
; 
;         KEYWORDS : /indiv - Fit individual profiles and then average the line positions. 
;                    Fit the profile averaged over the selected portion of the slit if not set. 
;					angstroms_per_pixel - spectral scale, defaults to 0.02546
;					logdir	- If set, then the output is added to a log file
;					/fuvs	-	If set (and an FUV image passed in), fit the O I line
;					/fuvl	-	If set (and an FUV image passed in), fit the Fe II line
;
;			OUTPUT : msg	-	a STRING to put in the FITS header
;
;		RETURNS: result	-	a structure containing some identifying information from the
;							header, along with some of the results of the fit:
;                  
;          HISTORY : Sept 7, 2013: iris_orbitvar_corr written by HUI TIAN at CfA
;					Sept 26, 2013: 	P. Boerner
;									iris_prep_wavecorr spun off from iris_orbitvarr_cor
;									and included in iris_prep pipeline
;                    Apr 24, 2017: JPW - changed FUV wrange from 41 to 41-(1-lvl1_in)*5
;
;-

lvl1_in = 1
itags = TAG_NAMES(hdr)
lvlnum_tag = WHERE(itags eq 'LVL_NUM', haslvl)
if haslvl then if hdr.lvl_num ne 1 then lvl1_in = 0

case hdr.img_path of
	'NUV'	:	begin
		; Fit the Ni I 2799.474 line
		wave0 = 2799.474
		; Three spatial pixel ranges
		ymins = FLOOR([280., 500., 720.] / hdr.sumspat)
		ymaxs = CEIL(ymins + 60. / hdr.sumspat)
		; Line is slightly shifted in the three ranges
		wrange = 21. ;;; - (1-lvl1_in)*3. ;  ??
		wmins = FLOOR( (650. + lvl1_in*[-5., 0., 5.]) / hdr.sumsptrl)
		wmaxs = CEIL(wmins + (wrange / hdr.sumsptrl > 4))
		ptag = 'CDLT1_NU'
		imgpath = 'NUV'
		linesign = 1.
	end
	'FUV'	:	begin
		linesign = -1.
		; Spatially, divide into 3 big chunks...
		ymins = FLOOR([20, 380., 740.] / hdr.sumspat)
		ymaxs = CEIL(ymins + 350. / hdr.sumspat)
		if KEYWORD_SET(fuvl) then begin
			; Fit the Fe II 1393 line; pretty hopeless, actually
			wave0 = 1392.82
			wrange = 41.
			wmins = FLOOR( (3007. + lvl1_in*[6., 0., -8.]) / hdr.sumsptrl)
			wmaxs = CEIL(wmins + (wrange / hdr.sumsptrl > 4))
			imgpath = 'FUVL'
			ptag = 'CDLT1_F2'		
		endif else begin
			; Fit the O I line 1356
			wave0 = 1355.60
			wrange = 41.-(1-lvl1_in)*5. ; avoid clipped area in certain cases
			wmins = FLOOR( (1822. + lvl1_in*[6., 0., -8.]) / hdr.sumsptrl)
			wmaxs = CEIL(wmins + (wrange / hdr.sumsptrl > 4))
			imgpath = 'FUVS'
			ptag = 'CDLT1_F1'		
		endelse
	end
	else	:	PRINT, 'IRIS_PREP_WAVECORR : Weird img_path...'
endcase

write = KEYWORD_SET(logdir)

result = {t_obs : hdr.t_obs, fsn : hdr.fsn, img_path : imgpath, roll : hdr.sat_rot, $
	exptime : hdr.exptime, obs_vr : hdr.obs_vr, sumsptrl : hdr.sumsptrl, $
	scvelo_pix : 0d, line_bot : 0d, line_mid : 0d, line_top : 0d, $
	chisq_bot : 0d, chisq_mid : 0d, chisq_top : 0d}

;; averaged positions of the neutral line at different spatial positions
mean_wpos = DBLARR(N_ELEMENTS(ymins))
chisqs = DBLARR(N_ELEMENTS(ymins))

for k = 0, N_ELEMENTS(ymins) - 1 do begin

	ymin = ymins[k]
	ymax = ymaxs[k]
	wmin = wmins[k]
	wmax = wmaxs[k]
	odat = dat[wmin:wmax,ymin:ymax]

	nw = wmax - wmin + 1 
	xx = FINDGEN(nw) + wmin  

	;++++++++++++++++++++++++++++++++++++
	; Average all line profiles with good
	; data in the selected range of the
	; slit, then do Gaussian fit
	;------------------------------------
	cuts = MEAN(odat, dim=1, /nan)
	goodcut = WHERE(FINITE(cuts), numgood)
	if numgood le 1 then ii = !values.f_nan else $
		ii = MEAN(odat[*,goodcut], dim=2, /nan)
	nans = WHERE(ii ne ii, numnan)	;	Don't try to fit if you don't have data
	if numnan ne 0 then begin
		a = FLTARR(5) - !values.f_nan
		chisq = !values.f_nan
	endif else begin
		ifit=mpfitpeak(xx, ii, a, nterms=5, chisq=chisq)
		;; find bad fit
		if (a[0]*linesign) gt 0. or $
		   finite(a[0]) ne 1 or $
		   a[1] gt ((wmin+wmax)/2.+(wrange/2.)/hdr.sumsptrl) or $
		   a[1] lt ((wmin+wmax)/2.-(wrange/2.)/hdr.sumsptrl) then a[*]=-!values.f_nan
	endelse
	mean_wpos[k] = a[1]
	chisqs[k] = chisq
	if (KEYWORD_SET(showfit)) then begin
		PLOT, xx, ii, xstyle=1, ystyle=1, $
			title='Gaussian Fit', xtitle='Wavelength pixel', ytitle='Intensity'
		OPLOT, xx, ifit, color=120
	endif

endfor
 
; Define the spectral pixel size
if not KEYWORD_SET(angstroms_per_pixel) then begin
	pdb = IRIS_MK_POINTDB()
	ptag_ind = WHERE(TAG_NAMES(pdb) eq ptag, hastag)
	angstroms_per_pixel = pdb.(ptag_ind)
endif
specsize=angstroms_per_pixel * hdr.sumsptrl  ;; spectral pixel size

; convert the spacecraft velocity to wavelength pixel, 
; then subtract spacecraft velocity from the line position
scvelo = hdr.obs_vr / 3.e8	;	v/c
result.scvelo_pix = scvelo * wave0 / specsize	;	Doppler shift in image pixels

; FUV produces a single summed result; stretch it to the same dimension as the NUV
mean_wpos = REBIN(mean_wpos, 3)
chisqs = REBIN(chisqs, 3)

result.line_bot = mean_wpos[0]
result.line_mid = mean_wpos[1]
result.line_top = mean_wpos[2]
result.chisq_bot = chisqs[0]
result.chisq_mid = chisqs[1]
result.chisq_top = chisqs[2]

msg = STRING('Orbital shift: ', result.scvelo_pix, $
	'pixels; measured line position = [', mean_wpos, ']', $
	form = '(a15, f6.2, a38, f7.2, 2f8.2, a1)')

if write then begin
	date_obs = STRMID(TIME2FILE(hdr.t_obs), 0, 8)
	txtpath = CONCAT_DIR(logdir, 'wavecorr')
	if not FILE_EXIST(txtpath) then begin
		BOX_MESSAGE, ['IRIS_PREP_WAVECORR: Output path does not exist!', txtpath]
	endif else begin
		filename = CONCAT_DIR(txtpath, date_obs + '_wave.txt')
		if not FILE_EXIST(filename) then begin
			OPENW, lun, /get, filename
			PRINTF, lun, TAG_NAMES(result), 'WRITTEN', form = '(A25,4A10,9A12,A25)'
		endif else begin
			OPENW, lun, /get, filename, /append
		endelse
		PRINTF, lun, result, RELTIME(/now, out='CCSDS'), form = '(A25,I10,A10,2F10.3,9F12.3,A25)'
		FREE_LUN, lun
	endelse
endif

if KEYWORD_SET(pstop) then STOP

RETURN, result

end