function IRIS_INTSCALE, image, header, $
	img_path = img_path, nobytescale = nobytescale, nonorm = nonorm, $
	maxfrac = maxfrac, minfrac = minfrac, minval = minval, maxval = maxval, $
	limb = limb, lmax = lmax, radlimb = radlimb, xoff = xoff, yoff = yoff, $
	exptime = exptime, pstop = pstop

;+
;
; Function for coarse first-estimate intensity scaling of the images
;
; INPUT: 
;	image	-	IRIS data cube, all with the same img_path (e.g. level2 data)
;	header (OPTIONAL) - FITS header (structure array)
;
; RETURNS: 
;	result	-	Data cube with the intensity adjusted for prettiness
;
; KEYWORDS:
;	/nobytescale-	By default, the result is bytescaled. Set this keyword if
;					you want to keep the same type as the input (though it still
;					does have a top and bottom threshold applied)
;	/nonorm		-	By default, the result is exposure-time normalized. Set this
;					keyword if you don't want to normalize it. (This does save
;					some memory space)
;	exptime=	-	Set to a float array of exposure time values for each frame
;					in the movie (in case you don't pass in a header structure;
;					this keyword is ignored if a header is passed in)
;	img_path=	-	Set to a string indicating the img_path of the input image
;					(in case you don't pass in a header structure; if a header
;					is passed in, this keyword is ignored)
;
; 	maxval		-	Pixel value to be set to maximum intensity (255); if set, ignore
;					maxfrac
;	maxfrac		-	Percentage of the image that should not be saturated; defaults
;					to 0.999. Lower values make things brighter. Don't set to 1.0
;					or more	or the image will be totally posterized.
;	minval		-	Pixel value to be set to minimum intensity (0); if set, ignore
;					minfrac
;	minfrac		-	Percentage of the image that should be clipped to black;
;					defaults to 0.002. Higher values make things darker. Don't set
;					below 0, or bad things will happen
;
;	/limb		-	If set, then it tries to apply a radial filter for limb
;					observations (linear scaling above limb, and whatever the
;					channel-dependent default is on the disk)
;	radlimb=	-	Radius of the limb, in arcsec, if /limb is set; defaults to 972
;	xoff/yoff	-	Offset of disk center from the header value, in arcsec; used
;					for fine-tuning limb alignment
;	lmax		-	Percentage of the portion above the limb that should not be
;					saturated (if /limb is set); defaults to 0.999. 
;
;	/pstop		-	Stop for debugging
;
; Paul Boerner        2014/12/11 spun off from AIA_INTSCALE
;
;-

imsize = SIZE(image)
if not KEYWORD_SET(nobytescale) then nobytescale = 0
if N_ELEMENTS(header) gt 0 then begin
	tags = TAG_NAMES(header[0])
	isspec = header[0].instrume eq 'SPEC'
	tdescind = WHERE(tags eq 'TDESC1', hastdesc)
	if hastdesc then begin
		if isspec then img_path = STRMID(header[0].tdet1, 0, 3) else img_path = header[0].tdesc1
	endif else begin
		img_path = header[0].img_path
	endelse
	exptimesind = WHERE(tags eq 'EXPTIMES', hasexps)
	if hasexps then begin
		if isspec then exptimes = hdr.exptimen > hdr.exptimef else exptimes = header.exptimes
	endif else begin
		exptimes = header.exptime
	endelse
endif

; Do exposure-time normalization unless exptimes unavailable or /nonorm is set
if imsize[0] eq 2 then numframes = 1 else numframes = imsize[3]
if not KEYWORD_SET(nonorm) and N_ELEMENTS(exptimes) eq numframes then begin
	norm_image = FLOAT(image)
	for i = 0, numframes - 1 do norm_image[*,*,i] = image[*,*,i] / exptimes[i]
	imvar = 'norm_image'
	logmin = 1. / MAX(exptimes)
endif else begin
	imvar = 'image'
	logmin = 1.
endelse

; Allowed values of img_path
; img_paths = ['SJI5000W', 'SJI2832', 'SJI2796', 'SJI1600W', 'SJI1400', 'SJI1330', $
;	'FUV', 'NUV', 'SJI_NUV']

; Set up intensity scaling based on image path
scale_log = 0		;	Use log10 scaling
scale_pow = 0		;	Use SQRT or some other power law scaling
dminfrac = 0.002	;	Default histogram level to make black
dmaxfrac = 0.999	;	Default histogram level to make white
case img_path of
	'SJI_1330'	:	begin
		scale_log = 1
	end
	'SJI_1400'	:	begin
		scale_log = 1
	end
	'SJI_2796'	:	begin
		scale_pow = 0.75
		dminfrac = 0.005
	end
	'SJI_2832'	:	begin
		scale_pow = 1.00
	end
	'NUV'	:	begin
		scale_pow = 0.50
	end
	'FUV'	:	begin
		scale_log = 1
	end
	else		: begin
	  PRINT, 'IRIS_INTSCALE: no intensity scaling defined for selected img_path'
	  RETURN, -1
	end
endcase

; Determine top and bottom intensity cutoffs
if N_ELEMENTS(maxval) * N_ELEMENTS(minval) eq 0	then begin		;	Need image histogram
	hist_cmdstr = 'finimg = ' + imvar + '[WHERE(FINITE(image))]'
	hist_status = EXECUTE(hist_cmdstr)
	if MAX(finimg) lt 0 then dhist = [1] else $			;	Avoid crash on pathological data
		dhist = HISTOGRAM(finimg, min = 0, bin = 1)
	chist = TOTAL(dhist, /cum)/DOUBLE(TOTAL(dhist))
endif
if N_ELEMENTS(maxval) eq 0 then begin
	if N_ELEMENTS(maxfrac) eq 0 then maxfrac = dmaxfrac				;	Percentage cutoff
	maxbin = INTERPOL(INDGEN(N_ELEMENTS(chist)), chist, maxfrac); intensity bin of cutoff
endif else begin
	maxbin = maxval
endelse
if N_ELEMENTS(minval) eq 0 then begin
	if N_ELEMENTS(minfrac) eq 0 then minfrac = dminfrac				;	Percentage cutoff
	minbin = INTERPOL(INDGEN(N_ELEMENTS(chist)), chist, minfrac); intensity bin of cutoff
endif else begin
	minbin = minval
endelse

; Scale the image (build a command string to pass to EXECUTE)
if scale_log eq 1 then begin
	arg1 = 'ALOG10( ' + imvar + '>logmin )' 
	argmin = 'ALOG10( minbin>logmin )'
	argmax = 'ALOG10( maxbin )'
endif else begin
	arg1 = '(' + imvar + '>0.)^(' + STRING(scale_pow, form = '(f4.2)') + ')
	argmin = '(minbin>0.)^(' + STRING(scale_pow, form = '(f4.2)') + ')
	argmax = '(maxbin>0.)^(' + STRING(scale_pow, form = '(f4.2)') + ')
endelse
if not KEYWORD_SET(nobytescale) then begin
	cmdstr = 'result = BYTSCL( ' + arg1 + ', min = ' + argmin + ', max = ' + argmax + ' )'
endif else begin
	cmdstr = 'result = ' + arg1 + ' > ' + argmin + ' < ' + argmax
endelse

if KEYWORD_SET(pstop) then STOP

status = EXECUTE(cmdstr)

; Apply linear scaling above the limb (useful, maybe, for looking at spicules in Mg II)
if KEYWORD_SET(limb) then begin
	if N_ELEMENTS(xoff) eq 0 then xoff = 0
	if N_ELEMENTS(yoff) eq 0 then yoff = 0
	xx = REBIN(FINDGEN(imsize[1]), imsize[1], imsize[2])
	yy = REBIN(TRANSPOSE(FINDGEN(imsize[2])), imsize[1], imsize[2])
	xx = (xx - header[0].crpix1 ) * header[0].cdelt1 + header[0].crval1 + xoff
	yy = (yy - header[0].crpix2 ) * header[0].cdelt2 + header[0].crval2 + yoff
	rad = SQRT(xx^2 + yy^2)
	if imsize[0] eq 3 then rad = REBIN(rad, imsize[1], imsize[2], imsize[3])
	if N_ELEMENTS(radlimb) eq 0 then radlimb = 972
	offlimbpix = WHERE(rad ge radlimb, numoff)
	way_offlimbpix = WHERE(rad ge radlimb+5 and FINITE(image), numwayoff)
	if numoff gt 0 then begin
		if numwayoff gt 1 then begin
			lhist = HISTOGRAM(image[way_offlimbpix], min = 0)
			lchist = TOTAL(lhist, /cum)/DOUBLE(TOTAL(lhist))
			if N_ELEMENTS(lmax) eq 0 then lmax = 0.999				;	Percentage cutoff
			limbmax = INTERPOL(INDGEN(N_ELEMENTS(lchist)), lchist, lmax) ; intensity bin of cutoff
		endif else begin
			lhist = HISTOGRAM(image[offlimbpix], min = 0)
			lchist = TOTAL(lhist, /cum)/DOUBLE(TOTAL(lhist))
			if N_ELEMENTS(lmax) eq 0 then lmax = 0.999			;	Percentage cutoff
			limbmax = INTERPOL(INDGEN(N_ELEMENTS(lchist)), lchist, lmax) ; intensity bin of cutoff
		endelse
		limbresult = 0. > image < limbmax	;	Use linear scaling and a low threshold above the limb
		if not KEYWORD_SET(nobytescale) then limbresult = BYTSCL(limbresult, min = 0, max = limbmax)
		result[offlimbpix] = limbresult[offlimbpix]
	endif
endif

RETURN, result

end