function AIA_BP_ADD_ERR, inwave, inresp, $
	normout, noiseout, shiftout, $ 
	norm = norm, noise = noise, wshift = wshift, noisesmooth = noisesmooth, $
    normskew = normskew, sigmas = sigmas, loud = loud
    

;+
;
; $Id: aia_bp_add_err.pro,v 1.2 2011/10/11 20:30:55 boerner Exp $
;
;
; Takes an efficiency curve for one of the components of an EUV/SXR telescope
; and adds pseudorandom numbers in a variety of ways in order to simulate
; the effects of calibration error.
;
; INPUTS: 
;	inwave=	-	the wavelength grid of the input efficiency curve (an N-element vector)
;
;	inresp=	-	the input efficiency (generally assumed to be reflectivity or
;				transmittivity, as a fraction (0.1 rather than 10%))
;
; RETURNS:
;	outresp	-	the output efficiency, on the same wavelength grid, with noise
;				effects added
;
; OUTPUTS:
;	normout -	a scalar giving the overall normalization constant applied to
; 				the result (ie the random number chosen from the distribution
;				specified with the norm= keyword)
;
;	noiseout -	an N-element array giving the noise scaling factors applied to
;				each point in the result (ie the pseudo-random array generated
; 				based on the noise and noisesmooth keywords)
;
;	shiftout -	a scalar giving the overall wavelength shift applied to the 
;				result (ie the random number chosen from the distribution
;				specified with the shift= keyword)
;
; KEYWORDS: 
;	norm=	-	the overall normalization error level; specifically, the 
;				standard deviation of the distribution that the normalization
;				randomizer scalar is chosen from
;
;	noise=	-	the point-to-point error level; specifically, the standard 
;				deviation of the distribution that the wavelength-dependent
;				randomizer vector elements are chosen from
;
;	noisesmooth= -	the wavelength scale of the point-to-point error; specifically,
;					the number of Angstroms between independent	noise knots in the
;					noiseout vector
;
;	wshift=	-	the level of the wavelength error; specifically, the standard
;				deviation of the distribution that the wavelength randomizer
;				shiftout is drawn from
;
;   /normskew -	if set, the distribution of normalization factors is skewed towards higher
;				numbers, while maintaining the mode of the distribution at 1. This is
;				sometimes useful if using worst-case numbers (eg for filter oxidization),
;				rather than best-guess numbers.
;
;   sigmas=	-	can be set to the number of sigmas of confidence on the incoming
;				error estimates. Defaults to 1 (ie these are 1-sigma errors, so that
;				the real answer lies within the specified error 68% of the time). Set
;				to 1.65 to give ~90% confidence.
;
;	/loud	-	if set, a plot or two is generated.
;
; $Log: aia_bp_add_err.pro,v $
; Revision 1.2  2011/10/11 20:30:55  boerner
; Removed extra line breaks
;
; Revision 1.1  2010/05/28 22:31:25  boerner
; Initial version
;
;
;-


if not KEYWORD_SET(norm) then norm = 0
if not KEYWORD_SET(noise) then noise = 0
if not KEYWORD_SET(wshift) then wshift = 0
if not KEYWORD_SET(noisesmooth) then noisesmooth = 1

numpoints = N_ELEMENTS(inwave)

if not KEYWORD_SET(sigmas) then sigmas = 1.

; Calculate normout, the overall normalization factor
if norm ne 0 then begin
	normout = 10.^(RANDOMN(seed1) * ALOG10(1 + norm) / sigmas)
	if KEYWORD_SET(normskew) then begin
	    doskew = RANDOMU(seed2)
    	if normout / norm le (-2) * doskew then normout = 0 - normout
	endif
endif else normout = 1.0

; Calculate noiseout, the wavelength-dependent noise factor
if noise ne 0 then begin
	waverange = MAX(inwave) - MIN(inwave)
	numknots = waverange / noisesmooth
	noiseknots = 10.^(RANDOMN(seed3, numknots) * ALOG10(1 + noise) / sigmas)
	noiseknotwaves = MIN(inwave) + (FINDGEN(numknots) * noisesmooth)
	noisewaveshift = RANDOMU(seed4) * noisesmooth / 10.	;	shift the knots in order
	noiseknotwaves = noiseknotwaves + noisewaveshift	;	to avoid funniness
	noiseout = INTERPOL(noiseknotwaves, noiseknots, inwave)
endif else noiseout = FLTARR(numpoints) + 1.0


; Calculate wshift, the wavelength shift scalar
; Also, create outresp, the result vector
if wshift ne 0 then begin
	shiftwave = inwave
	shiftout = RANDOMN(seed5) * (wshift / sigmas)
	shiftwave = inwave + shiftout
	outresp = INTERPOL(inresp, shiftwave, inwave)
endif else begin
	outresp = inresp
endelse

outresp = outresp * noiseout * normout

if KEYWORD_SET(loud) then begin
    PLOT, inwave, outresp / inresp, /ylog, xtitle = 'Wavelength [' + STRING("305B) + ']', ytitle = 'Ratio (Noisy / Original)'
    PRINT, 'Overall offset: ', normout
    if noise ne 0 then PRINT, 'Wavelength offset: ', noisewaveshift
    PRINT, 'Mean ratio: ', MEDIAN(outresp / inresp)
endif

RETURN, outresp

end