function AIA_MP_RESPTABLE_CHISQ, x, params
;+
; Helper routine for AIA_BP_RESPONSE_TABLE
; Chi-squared calculator to be used with MPFIT to find the best-fit
; slope for a given time interval.
;-
yout = params[0] + params[1]*x + params[2]*x^2 + params[3]*x^3 
RETURN, yout
end



function AIA_BP_RESPONSE_TABLE, $
	output_data, outdata_all=outdata_all, $
	wavelnth = wavelnth, t_start = t_start, t_stop = t_stop, $
	poly = poly, p_stop = p_stop, silent = silent, loud=loud, nopng=nopng, nofig=nofig, $
	all = all, thicktoo = thicktoo, $
	nosmooth = nosmooth, noevebake = noevebake, no304bake = no304bake, $
	printtab = printtab, ver_num = ver_num, veve = veve, figpath = figpath, scratchpath=scratchpath, $
	fix_prebake = fix_prebake  $
	, nodisc = nodisc, rocket_only = rocket_only  $
	, aia2eve_ratio_path=aia2eve_ratio_path, aia2see_ratio_path=aia2see_ratio_path, aia2model_ratio_path=aia2model_ratio_path  $
    , smooth_width=smooth_width, smooth_threshold=smooth_threshold, post_eve_heavy_smooth=post_eve_heavy_smooth

;+
;
; Generates a response table containing the keywords to be used in the JSOC series aia.response
;
; RETURNS:
;	response	-	a structure that constitutes a single entry (row) in the AIA response table, if no /all set;
;					otherwise, if /all, then it returns a structure that constitutes all entries (multiple rows)
;
; KEYWORD PARAMETERS:
;	wavelnth=	-	(INPUT, long) the wavelength channel to calculate the entry for
;	t_start=	-	(INPUT) the starting time for the epoch described by this entry.
;					Generally will correspond to a change in the response, e.g. a 
;					CCD bakeout. Can be in any time format recognized by ANYTIM
;	t_stop=		-	(INPUT) the ending time for the epoch described by this entry.
;					If not specified, defaults to 1 year after t_start. 
;	poly=		-	The degree of the polynomial used to fit the AIA to EVE ratio during
;					the specified time range. Defaults to 0 (only a normalization). 
;					If set to a number<0, then it is used as a hard-coded scale factor 
;					(the absolute value of poly is multiplied by the nominal effective 
;					area; setting to -1 means it uses the nominal BOL value as-is)
;	/silent		-	if set, then no print info
;	/p_stop		-	if set, then pause in the routine for debugging
;	/nosmooth	-	if not set, then a median filter is used to unspike the data. 
;	/noevebake	-	if not set, then the data from the EVE bakeouts +~2 weeks is left 
;					out of the fit. Set this keyword if you want to include those data
;					(not sure why you would, though)
;	/fix_prebake-	if set, then the data from prior to the first EVE bakeout 
;					(Sept 2010) is left out of the fit (useful for EVE v3, at least for
;					the 304 channel)
;	/no304bake	-	if not set, then the data from 2011-05-19 (during the AIA 304 bakeout)
;					is left out of the fit. Set this keyword if you want to include those
;					data (not sure why you would, though)
;	/thicktoo	-	if set, then generate a response for both thin and thick filter
;					versions of the specified channel (the THICK version will use the same 
;					degradation as the thin version, scaled by the ratio of thin to thick
;					effective area in the initial response functions). Thick filter versions are
;					included for the setkeys output file (to write to DRMS), but not for the
;					text file distributed via SSW
;
;	/all		-	if set, then run this routine recursively on all the EUV wavelengths
;	/printtab	-	if set, then the table is printed out to a text file in the scratchpath
;					'/net/crom/Volumes/disk2/data/aia/test/figs/response_table/'
;					(which can then be copied to somewhere more important if it looks ok)
;					(Note that plots are saved in that path whether or not /printtab is
;					set, unless /nopng is set)
;	ver_num		-	a pass-through to the output table (unless /all is set, in which case
;					it is used to define the intervals of the table)
;	veve		-	set to the integer of the EVE calibration version to use; defaults
;					to 4
;	figpath		-	set to the path where figures should be saved
;
;	--- new since 2017/11/28 -----------
;	scratchpath	-	(added as keyword 2017/11/28) path where to write the output table files
;
;	/nofig		-	if set, then no plot is generated
;	/nopng		-	if set, then the plots are not saved as PNGs
;	/loud		-	set to print more info (much more than w/o /silent), could be messy, use with caution 
;
;	aia2eve_ratio_path, default='/sanhome/weiliu/data/eve/aiafiles/ratio_aia-eve/'
;	aia2fism_ratio_path, default='/sanhome/weiliu/data/eve/aiafiles/ratio_aia-fism/'
;	aia2see_ratio_path, default='/sanhome/weiliu/data/see/ratio_aia-see/'
;
;	smooth_width: default =[9 (from Paul), 27]; two values, before and after 2014-May-25
;   smooth_threshold: default=[0.1, 0.1] (from Paul), two values, used before and after 2014-May-25
;	post_eve_heavy_smooth: set to replace original useratios with median smoothed version entirely, not just outliers as in the default case;
;		this is only for post-eve (after 2014-May-25) epoch times
;
;
; OUTPUT PARAMETER:
;	output_data	-	(OPTIONAL) Structure containing the data used in the fit. 
;					If the /all keyword is set, this is simply set to an array of
;					residuals from all the fits (used for AIA_CALPLOT_RESIDUALS)
;
;	outdata_all: optional output, for all channels, (added 2017/12/09)
;	 ={outdata_94:outdata_94, outdata_131:outdata_131, outdata_171:outdata_171, outdata_193:outdata_193, outdata_211:outdata_211, outdata_304:outdata_304, outdata_335:outdata_335, outdata_1600:outdata_1600, outdata_1700:outdata_1700}
;	 for reach channel, outdata_xxx='{wave: wvls[thiswave], usetime_utc_tai : usetime_all, useratios : useratios_all, yfit : yfit_all}'
;	 Note: can converted usetime_tai from TAI to UTC: tai2utc(usetime_tai, /ccsds) for utplot, etc. Don't convert here for output, b/c (1) otherwise double .genx file size, (2) for consistency here.
;
; History:
;  2017/11/28, Tue: Wei Liu, take over Paul's original version,
;   workspace: ~/sdo/docs/my-job/calib/eve/wkspc/run01_20170926_FISM-V2_PhilC/v1a_batch/cmpr_real-eve_v2_scale-model/response-table/v0_test-Paul 
;   1) changed figpath default from '/Users/boerner/data/aia/trending/response_table/' to '/sanhome/weiliu/data/aia/trending/response_table/', but save *.png to ./png/ subdirectory
;   2) added scratchpath as keyword, default at figpath
;   3) added seefile for eve to see ratio, now seefile= '/sanhome/weiliu/data/aia/cal/20141024_231847_eve_see_ratios.genx'
;   4) hand-added V7 from the screen-shot by Paul's note of AIA-Throughput-trending.html
;      tried hard, but couldn't find Paul's version of this code with V7 epoch dates.
;   5) finally check into SSW:
;      cd /sanhome/weiliu/data/aia/trending/response_table/20171129_195430
; 	5.1) online aia_bp_response_table.pro $SSW/sdo/aia/idl/response/
;      >>> Backing up current version via...:
;      mv -f $SSW/sdo/aia/idl/response//aia_bp_response_table.pro $SSW/offline/ssw_backup/aia_bp_response_table.pro.171129.235133
;   5.2) table files:
;      online aia_V7_20171129_195626_response_table.txt $SSW/sdo/aia/response
;      online aia_V7_20171129_195626_response_table.set_keys  $SSW/sdo/aia/response
;  2017/11/29, Wed: Wei Liu, 
;   workspace: ~/sdo/docs/my-job/calib/eve/wkspc/run01_20170926_FISM-V2_PhilC/v1a_batch/cmpr_real-eve_v2_scale-model/response-table/v1_new-V8
;   1) upgrade to V8, to include 
;     a) latest aia/eve ratio using Phil Chambline's daily average FISM v2 model spectra, for all EUV channels.
;        this requires switching for 304 and 335 channels from using aia/see ratio as Paul added after 2014-May-25 EVE/MEGS-A being out of commission.
;     b) latest aia/see ratio for UV (1600, 1700) channels.
;   2) for epoch times, combine V6 and V7, then revise to add new dates to break up long stretch, roughly once/year in May.
;   3) linear scaling of aia/model ratio to approximate aia/eve ratio after 2014-may-25
;  2017/12/09, Sat (documentation here added/revised 2017/12/11, Mon): 
;   v0: 1) removed the requirement for the last epoch interval being flat (which was added by Paul after EVE/MEGS-A failure, but can under-estimate degradation), 
;     since now we will continuously monitor degradation.
;    2) added plotting and saving all {usetime, useratio, yfit} to outdata_xxx thru new keyword outdata_all, and also into the .genx file
;   v1: revised epoch times, to avoid unnecessary jumps; as baseline before the change, refer plots at:
;	   ~/sdo/docs/my-job/calib/eve/wkspc/run01_20170926_FISM-V2_PhilC/v1a_batch/cmpr_real-eve_v2_scale-model/response-table/v1_new-V8/test-runs/20171210_013000_usetime-TAI
;     General rule to follow, going forward: keep tstarts accumulated to date unchanged, but revise startstr, to minimize changes needed for startstr in the future.
;   v2: heavy smooth - replaced original useratios with median smoothed version entirely, not just outliers as in the default case;
;		this is only for post-eve (after 2014-May-25) epoch times; the result is not much change from before, since polyfit already takes care of noisy data (therefore, made change in v3 below)
;   v3: added option of post_eve_heavy_smooth, set to do heavy smoothing as in v2, default is not to.
;   v4: made 193 and 211 on big epoch interval post EVE/after 2014-may-25, b/c large fluctuations/uncertainties do not justify any additional breaks for piecewise linear fits
;	 => settle with v4 for release to SSW on 2017/12/11, Mon.
;  2020/06/25: PBoerner
;    workspace: /sanhome/boerner/data/aia/v9
;	 Upgrade to v9 of the AIA response table
;-

;++++++++++++++++++++++++++++++++++++++
; Set up defaults
;--------------------------------------
if not KEYWORD_SET(ver_num) then ver_num = 1
if not KEYWORD_SET(veve) then veve = 6
thicktoo = KEYWORD_SET(thicktoo)
if not KEYWORD_SET(silent) then silent = 0
if not KEYWORD_SET(fix_prebake) then fix_prebake = 0
if not KEYWORD_SET(nopng) then nopng = 0
if not KEYWORD_SET(nofig) then nofig = 0
if not KEYWORD_SET(rocket_only) then rocket_only = 0
if not KEYWORD_SET(p_stop) then p_stop = 0		


	;if N_ELEMENTS(figpath) eq 0 then figpath = '/Users/boerner/data/aia/trending/response_table/'
	;scratchpath = '/Users/boerner/data/aia/trending/response_table/'
checkvar, figpath, '/sanhome/boerner/data/aia/v9/response_table/'
checkvar, scratchpath, figpath
thistime = AIA_BP_DATE_STRING(/time)
scratchpath = scratchpath + thistime + '/'
	;if not silent and KEYWORD_SET(all) then begin
if KEYWORD_SET(all) then begin	; 2017/11/28, commented (not silent) condition, otherwise no result written
  SPAWN, 'mkdir ' + scratchpath, spawnresult, spawnerror
  spawn, 'mkdir ' + scratchpath + '/png'
endif

wvls =	[94,	131,	171,	193,	211,	304,	335,	1600,	1700,	4500]
wvls_used= wvls[0:8]	; skip 4500
ewvls =	[93.9,	131.2,	171.1,	195.1,	211.3,	303.8,	335.4,	1600.0,	1700.0,	4500.0]
if N_ELEMENTS(poly) eq 0 then poly = 0
if poly gt 3 then begin
	BOX_MESSAGE, ['AIA_BP_RESPONSE_TABLE:', $
				'Set poly to a number between 0 and 3! Using 3...', $
				'Or set poly<0 to use scaled EA value without correction' ]
	poly = 3
endif

;--- Wei's addition, starting 2017/11/28 ---------------
;--- PB 2020/06/22: Leave these defined for V7/V8, add a new path for V9 -------
checkvar, aia2eve_ratio_path, '/sanhome/weiliu/data/eve/aiafiles/ratio_aia-eve/'
checkvar, aia2model_ratio_path, '/sanhome/weiliu/data/eve/aiafiles/ratio_aia-fism/'
checkvar, aia2see_ratio_path, '/sanhome/weiliu/data/see/ratio_aia-see/'
eve_end_date='25-may-2014 12:00'
eve_end_tai=ANYTIM2TAI(eve_end_date)

	;@~/idl_lib/bin_wei/reset_bw	; need to do this once
	 !p.color=255     ;0
	 !p.background=0  ;255
PB_SET_LINE_COLOR


checkvar, smooth_width, [9, 27]			; two values, before and after 2014-May-25
checkvar, smooth_threshold, [0.1, 0.1]	; two values, used before and after 2014-May-25	;[0.1, 0.03]	

;++++++++++++++++++++++++++++++++++++++
; Run this routine recursively for a
; predefined set of intervals
;--------------------------------------
if KEYWORD_SET(all) then begin

	window, 9, xs=1000, ys=800
	window, 10, xs=1500, ys=1000

	; Define the intervals depending on the ver_num
	case ver_num of
		1 : begin
				tstarts = [	'24-mar-2010 00:00', '1-may-2030 00:00']
				wavelist 	=	INDGEN(9)
				tstartlist	= 	INTARR(9)
				tstoplist 	=	INTARR(9) + 1
				polylist	=	INTARR(9) - 1
				fname 		=	'aia_preflight_'
			end
		2 : begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $
							'1-may-2030 00:00']
				wavelist 	=	[0, 0, 0, 0, 0, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 5, 5, 5, 6, 6, 7, 7, 8, 8]
				tstartlist	= 	[0, 1, 3, 4, 5, 0, 2, 0, 1, 0, 1, 0, 1, 0, 1, 3, 4, 5, 0, 2, 0, 1, 0, 1]
				tstoplist 	=	[1, 3, 4, 5, 6, 2, 6, 1, 6, 1, 6, 1, 6, 1, 3, 4, 5, 6, 2, 6, 1, 6, 1, 6]
				polylist	=	[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1]
				fname		=	'aia_V2_'
			end
		3 :	begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '12-apr-2012 12:00', '01-may-2030 00:00']		;	6	7	8
				polyperwave = [0, 0, 0, 1, 1, 1, 1, 1, 0, 0]
				startstr = INTARR(50, 9) - 1
				startstr[0, 0] =	[0, 1, 3, 4, 5, 6, 7, 8]	;	94
				startstr[0, 1] =	[0, 2, 6, 7, 8]				;	131
				startstr[0, 2] =	[0, 1, 6, 8]				;	171
				startstr[0, 3] =	[0, 1, 6, 7, 8]				;	193
				startstr[0, 4] =	[0, 1, 6, 7, 8]				;	211
				startstr[0, 5] =	[0, 1, 3, 4, 5, 6, 7, 8]	;	304
				startstr[0, 6] =	[0, 2, 6, 7, 8]				;	335
				startstr[0, 7] =	[0, 1, 6, 7, 8]				;	1600
				startstr[0, 8] =	[0, 1, 6, 7, 8]				;	1700
				for i = 0, 8 do begin
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist	= polyperwave[i]
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
						endelse
					endfor
				endfor
;				wavelist 	=	[0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 7, 7, 7, 7, 8, 8, 8, 8]
;				tstartlist	= 	[0, 1, 3, 4, 5, 6, 0, 2, 6, 0, 1, 6, 0, 1, 6, 7, 0, 1, 6, 7, 0, 1, 3, 4, 5, 6, 7, 0, 2, 6, 7, 0, 1, 6, 7, 0, 1, 6, 7]
;				tstoplist 	=	[1, 3, 4, 5, 6, 8, 2, 6, 8, 1, 6, 8, 1, 6, 7, 8, 1, 6, 7, 8, 1, 3, 4, 5, 6, 7, 8, 2, 6, 7, 8, 1, 6, 7, 8, 1, 6, 7, 8]
;				polylist	=	[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1]
				fname		=	'aia_V3_'
			end
		4 :	begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2030 00:00']												;	9
				polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 10) - 1
				startstr[0, 0] =	[0, 1, 3, 4, 5, 6, 7, 8, 9]	;	94
				startstr[0, 1] =	[0, 2, 6, 7, 8, 9]			;	131
				startstr[0, 2] =	[0, 1, 6, 8, 9]				;	171
				startstr[0, 3] =	[0, 1, 6, 7, 8, 9]			;	193
				startstr[0, 4] =	[0, 1, 6, 7, 8, 9]			;	211
				startstr[0, 5] =	[0, 1, 3, 4, 5, 6, 7, 8, 9]	;	304
				startstr[0, 6] =	[0, 2, 6, 7, 8, 9]			;	335
				startstr[0, 7] =	[0, 1, 6, 7, 8, 9]			;	1600
				startstr[0, 8] =	[0, 1, 6, 7, 8, 9]			;	1700
				for i = 0, 8 do begin
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist	= polyperwave[i]
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
						endelse
					endfor
				endfor
				fname		=	'aia_V4_'
			end
		5 :	begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2013 12:00', '01-may-2030 00:00']							;	9	10
				polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 11) - 1
				startstr[0, 0] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 10]	;	94
				startstr[0, 1] =	[0, 2, 6, 7, 8, 9, 10]			;	131
				startstr[0, 2] =	[0, 1, 6, 8, 9, 10]				;	171
				startstr[0, 3] =	[0, 1, 6, 7, 8, 9, 10]			;	193
				startstr[0, 4] =	[0, 1, 6, 7, 8, 9, 10]			;	211
				startstr[0, 5] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 10]	;	304
				startstr[0, 6] =	[0, 2, 6, 7, 8, 9, 10]			;	335
				startstr[0, 7] =	[0, 1, 6, 7, 8, 9, 10]			;	1600
				startstr[0, 8] =	[0, 1, 6, 7, 8, 9, 10]			;	1700
				for i = 0, 9 do begin
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist	= polyperwave[i]
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
						endelse
					endfor
				endfor
				fname		=	'aia_V5_'
			end
		6 :	begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2013 12:00', '01-oct-2013 12:00', '01-may-2014 12:00', $	;	9	10	11
							'25-may-2014 12:00', '01-oct-2014 12:00', '01-may-2030 00:00']		;	12	13	14
				polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 10) - 1
				startstr[0, 0] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14]	;	94
				startstr[0, 1] =	[0, 2, 6, 7, 8, 9, 10, 11, 14]			;	131
				startstr[0, 2] =	[0, 1, 6, 8, 9, 10, 11, 14]				;	171
				startstr[0, 3] =	[0, 1, 6, 7, 8, 9, 10, 11, 14]			;	193
				startstr[0, 4] =	[0, 1, 6, 7, 8, 9, 10, 11, 14]			;	211
				startstr[0, 5] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 10, 12, 13, 14]	;	304
				startstr[0, 6] =	[0, 2, 6, 7, 8, 9, 10, 12, 13, 14]			;	335
				startstr[0, 7] =	[0, 1, 6, 7, 8, 9, 10, 12, 14]			;	1600
				startstr[0, 8] =	[0, 1, 6, 7, 8, 9, 10, 12, 14]			;	1700
				for i = 0, 9 do begin
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist = polyperwave[i]
							if (j eq numgood - 2) then polylist = 0 ;	Don't extrapolate last point
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
							if (j eq numgood - 2) then polylist[-1] = 0 ;	Don't extrapolate last point							
						endelse
					endfor
				endfor
				fname		=	'aia_V6_'
			end

		7 :	begin
				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2013 12:00', '01-aug-2013 12:00', '01-oct-2013 12:00', $	;	9	10	11
							'01-may-2014 12:00', '01-may-2015 12:00', '01-may-2030 00:00']		;	12	13	14

		;--- change V7 from V6, Wei's notes 2017/11/29 ---------
		;removed: '25-may-2014 12:00' last EVE/MEGS-A
		;		, '01-oct-2014 12:00'
		;added:	'01-aug-2013 12:00', => looks like for the 1600, 1700 flatfield update, but if so, should really be 16-aug-2013 00:00, as recently noted with John S. and Meng J.
		;		'01-may-2015 12:00'
		;NOTE: 1) hand-added V7 from the screen-shot by Paul's note of AIA-Throughput-trending.html
		;  2) Result is consistent with Paul's V7 run at: /sanhome/boerner/data/aia/trending/response_table/20150821_185751
		;  3) for consistency, online this version of the code to $SSW/sdo/aia/idl/response and the resulting V7 table to $SSW/sdo/aia/response
		;  4) Caution on removing epoch time '25-may-2014 12:00' from V6: 
		;   4.1) while this is good for EUV channels other than 304 & 335, i.e., use the flat average from '01-may-2014 12:00' to '25-may-2014 12:00' for the remaining time up to 2030, when degradation is unknown.
		;   4.2) for 304, 335: do not recommend this version, since makes the next start time of '01-may-2015 12:00'
		;   while only the short interval from '01-may-2014 12:00' to '25-may-2014 12:00' was used to determine the linear fits that is used for the entire 1-year interval (01-may-2014 to 01-may-2015 as in the table file):
		;   2014-05-01T12:00:00.000   2015-05-01T12:00:00.000          7   304_THIN 
		;  5) the *.png plots starting at '01-may-2014 12:00' still end at '25-may-2014 12:00' b/c the aia/eve ratio file ends then, as given by the if block below starting with:
		;   if wavelnth lt 1000 and tai_start lt ANYTIM2TAI('25-may-2014 12:00') then begin)
		;  6) To-Do: the next version should combine V6 and V7 times, recover '25-may-2014 12:00'

				polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 10) - 1
				startstr[0, 0] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 11, 12, 14]	;	94		; changed V6: 10, 11 to 11, 12
				startstr[0, 1] =	[0, 2, 6, 7, 8, 9, 11, 12, 14]			;	131		; changed V6: 10, 11 to 11, 12
				startstr[0, 2] =	[0, 1, 6, 8, 9, 11, 12, 14]				;	171		; changed V6: 10, 11 to 11, 12
				startstr[0, 3] =	[0, 1, 6, 7, 8, 9, 11, 12, 14]			;	193		; changed V6: 10, 11 to 11, 12
				startstr[0, 4] =	[0, 1, 6, 7, 8, 9, 11, 12, 14]			;	211		; changed V6: 10, 11 to 11, 12
				startstr[0, 5] =	[0, 1, 3, 4, 5, 6, 7, 8, 9, 11, 12, 13, 14]	;	304	; changed V6: 10, to 11
				startstr[0, 6] =	[0, 2, 6, 7, 8, 9, 11, 12, 13, 14]			;	335	; changed V6: 10, to 11
				startstr[0, 7] =	[0, 1, 6, 7, 8, 9, 10, 11, 12, 13, 14]		;	1600	; from V6, added 11, 13
				startstr[0, 8] =	[0, 1, 6, 7, 8, 9, 10, 11, 12, 13, 14]		;	1700	; from V6, added 11, 13
				for i = 0, 9 do begin
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist = polyperwave[i]
							if (j eq numgood - 2) then polylist = 0 ;	Don't extrapolate last point
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
							if (j eq numgood - 2) then polylist[-1] = 0 ;	Don't extrapolate last point							
						endelse
					endfor
				endfor
				fname		=	'aia_V7_'
			end

		8 :	begin		; keep latest version at the end of this block to override earlier versions

				;--- 2017/11/29 version ---------------------------------------------------------------------
				;--- V8 change from V7, Wei's notes 2017/11/29 ---------
				;Summary: same as V7 up to #12 (except for #10 revised slightly)
				; change of # only: old #13 -> new #14, old #14 -> new #17
				;removed: N/A
				;added:   '25-may-2014 12:00' (new #13) last EVE/MEGS-A, added back from V6
				;		  '01-may-2016 12:00' (new #15) to break up long stretches where one-piece linear fits do not apply
				;		  '01-nov-2016 12:00' (new #16) for 304 last flat portion
				;		  '01-jan-2017 12:00' (new #17) last data point of FISM-v2 as in /calib/eve/wkspc/run01_20170926_FISM-V2_PhilC/v1a_batch/dat/aia_eve_mission_long_edit_v1_remove-outlier.dat
				;changed: '01-aug-2013 12:00' (old #10) => '16-aug-2013 00:00' (new #10) for 1600, 1700 flatfield update
				;
				;tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
				;			'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
				;			'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
				;			'01-may-2013 12:00', '16-aug-2013 00:00', '01-oct-2013 12:00', $	;	9	10	11
				;			'01-may-2014 12:00', '25-may-2014 12:00', '01-may-2015 12:00', $	;	12	13	14
				;			'01-may-2016 12:00', '01-nov-2016 12:00', '01-jan-2017 12:00', $	;	15	16	17
				;			'01-may-2030 00:00']	; 18
				;polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				;startstr = INTARR(50, 10) - 1
				;startstr[0, 0] =	[0,1,3,4,5,6,7,8,9,   11,12,13,14,15,      18]	;	94		; change from V7: add 13,15,18, use 15-17 fit for 15-18 at constant/flat
				;startstr[0, 1] =	[0,2,      6,7,8,9,   11,12,13,14,15,      18]	;	131		; change from V7: follow 94
				;startstr[0, 2] =	[0,1,      6,  8,9,   11,12,13,14,15,      18]	;	171		; change from V7: follow 94
				;startstr[0, 3] =	[0,1,      6,7,8,9,   11,12,13,14,15,      18]	;	193		; change from V7: follow 94
				;startstr[0, 4] =	[0,1,      6,7,8,9,   11,12,13,14,15,      18]	;	211		; change from V7: follow 94
				;startstr[0, 5] =	[0,1,3,4,5,6,7,8,9,   11,12,13,14,15,16,   18]	;	304		; change from V7: follow 94, except adding 16
				;startstr[0, 6] =	[0,2,      6,7,8,9,   11,12,13,14,15,      18]	;	335		; change from V7: follow 94, 
				;startstr[0, 7] =	[0,1,      6,7,8,9,10,11,12,   14,15,   17,18]	;	1600	; change from V7: add 15,17,18
				;startstr[0, 8] =	[0,1,      6,7,8,9,10,11,12,   14,15,   17,18]	;	1700	; change from V7: follow 1600
				;
				;--- 2017/12/03 version, done 2017/12/09 ---------------------------------------------------------------------
				;Notes: understanding epoch times so far:
				; 1) 1600, 1700 flatfield updates: '01-jan-2012 12:00' '01-aug-2013 12:00'	;  6, 10
				; 2) total 6 bakeouts (see doc AIA03821 on 304 throughput loss, Table 1 on p.5): 
				;  confirmed with search for bakeout at http://aia.lmsal.com/public/jsocobs_info2011.html
				;                      '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	1	2
				; '13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
				;                      '10-apr-2012 12:00'	; 7
				; as noted below:
				; "Leave out the first 48 hours and last 24 hours of the interval from the fit", b/c bakeouts last usually 1-2 days.
				; 3) looks like all '01-may-201* 12:00' are hand-inserted breaks, kind of arbitrary, but they are among flat update dates (usually Feb/May/Aug/Nov-01)
				;
				;Notes on changes:
				; 1) remove #12 '01-may-2014 12:00' for EUV channels, but keep '25-may-2014 12:00'; do this by revising startstr, but keeping tstarts unchanged
				; 2) added '01-sep-2015 12:00' as new #15, all dates from old #15 get +1
				; 3) not yet identified (searched http://aia.lmsal.com/public/jsocobs_info2011.html, and ops daily logs):
				; '15-feb-2013 12:00'

				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2013 12:00', '16-aug-2013 00:00', '01-oct-2013 12:00', $	;	9	10	11
							'01-may-2014 12:00', '25-may-2014 12:00', '01-may-2015 12:00', $	;	12	13	14
	   '01-sep-2015 12:00', '01-may-2016 12:00', '01-nov-2016 12:00', '01-jan-2017 12:00', $	;15	16	17  18
							'01-may-2030 00:00']	; 19
				polyperwave = [0, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 10) - 1
				startstr[0, 0] =	[0,1,3,4,5,6,7,8,9,   11,  13,  15,  19]	;	94		; change from V8 2017/11/29: cut old 12 & 14, add new 15
				startstr[0, 1] =	[0,2,      6,7,8,9,   11,  13,  15,  19]	;	131		; change from V8 2017/11/29: follow 94
				startstr[0, 2] =	[0,1,      6,  8,9,   11,  13,  15,  19]	;	171		; change from V8 2017/11/29: follow 94
				startstr[0, 3] =	[0,1,      6,7,8,9,   11,  13,       19]	;	193		; change from V8 2017/11/29: follow 94, plus cut all after 13 up to the latest 19
				startstr[0, 4] =	[0,1,      6,7,8,9,   11,  13,       19]	;	211		; change from V8 2017/11/29: follow 193
				startstr[0, 5] =	[0,1,3,4,5,6,7,8,9,   11,  13,14, 16,19]	;	304		; change from V8 2017/11/29: cut old 12 & 16
				startstr[0, 6] =	[0,2,      6,7,8,9,   11,  13,14, 16,19]	;	335		; change from V8 2017/11/29: cut old 12, now same as 304
				startstr[0, 7] =	[0,1,      6,7,8,9,10,11,12,      16,19]	;	1600	; change from V8 2017/11/29: cut old 14 & 17
				startstr[0, 8] =	[0,1,      6,7,8,9,10,11,12,      16,19]	;	1700	; change from V8 2017/11/29: follow 1600

				for i = 0, 9 do begin		;--- loop over channels ----
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin		;--- loop over start times for each channel ----
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist = polyperwave[i]
								;if (j eq numgood - 2) then polylist = 0 ;	Don't extrapolate last point	; 2017/12/03, commented, no longer needed
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
								;if (j eq numgood - 2) then polylist[-1] = 0 ;	Don't extrapolate last point	; 2017/12/03, commented, no longer needed						
						endelse
					endfor
				endfor

				fname		=	'aia_V8_'
			end

		9 :	begin		; keep latest version at the end of this block to override earlier versions

				; Notes PB 2020-06-25
				; Added a time interval on 2018-03-05 to enable an alternate slope.

				tstarts = [	'24-mar-2010 00:00', '27-jan-2011 15:00', '24-feb-2011 19:00', $	;	0	1	2
							'13-apr-2011 18:00', '20-may-2011 18:00', '06-oct-2011 12:00', $	;	3	4	5
							'01-jan-2012 12:00', '10-apr-2012 12:00', '15-feb-2013 12:00', $	;	6	7	8
							'01-may-2013 12:00', '16-aug-2013 00:00', '01-oct-2013 12:00', $	;	9	10	11
							'01-may-2014 12:00', '25-may-2014 12:00', '01-may-2015 12:00', $	;	12	13	14
	  						'01-jan-2016 12:00', '01-may-2016 12:00', '01-may-2017 12:00', $	;	15	16	17 
							'05-mar-2018 12:00', '01-oct-2019 12:00', '01-jun-2020 12:00', $	; 	18	19	20 
;							'01-jan-2017 12:00', '05-mar-2018 12:00', '01-jun-2020 12:00', $	; 	18	19	20 
							'01-may-2030 00:00']												;	21
				use_nodisc = [0, 0, 0, $
							0, 0, 0, $
							0, 0, 0, $
							0, 0, 1, $
							1, 1, 1, $
							1, 1, 1, $
							1, 1, 1, $
							1]
				polyperwave = [1, 1, 1, 1, 1, 1, 1, 1, 1, 0]
				startstr = INTARR(50, 10) - 1
				startstr[0, 0] =	[0,1,3,4,5,6,7,8,9,   11,           16,         20,21]	;	94		
				startstr[0, 1] =	[0,2,      6,7,8,9,   11,                 18,   20,21]	;	131		
				startstr[0, 2] =	[0,1,      6,  8,9,   11,  13,   15,   17,	    20,21]	;	171		
				startstr[0, 3] =	[0,1,      6,7,8,9,   11,  13,            18,19,   21]	;	193		
				startstr[0, 4] =	[0,1,      6,7,8,9,   11,  13,         17,   19,   21]	;	211		
				startstr[0, 5] =	[0,1,3,4,5,6,7,  9,   11,  13,14,   16,   18,   20,21]	;	304		
				startstr[0, 6] =	[0,2,      6,7,8,9,   11,  13,      16,   18,   20,21]	;	335		
				startstr[0, 7] =	[0,1,      6,7,8,9,10,   12,        16,   18,   20,21]	;	1600	
				startstr[0, 8] =	[0,1,      6,7,8,9,10,   12,        16,   18,   20,21]	;	1700	
;				startstr[0, 0] =	[0,1,3,4,5,6,7,8,9,   11,  13,  15,   17,19,20,21]	;	94		
;				startstr[0, 1] =	[0,2,      6,7,8,9,   11,     14,           20,21]	;	131		
;				startstr[0, 2] =	[0,1,      6,  8,9,   11,12,  14,        19,20,21]	;	171		
;				startstr[0, 3] =	[0,1,      6,7,8,9,   11,     14,        19,20,21]	;	193		
;				startstr[0, 4] =	[0,1,      6,7,8,9,   11,  13,14,        19,20,21]	;	211		
;				startstr[0, 5] =	[0,1,3,4,5,6,7,8,9,   11,  13,14,  16,   19,20,21]	;	304		
;				startstr[0, 6] =	[0,2,      6,7,8,9,   11,  13,14,  16,   19,20,21]	;	335		
;				startstr[0, 7] =	[0,1,      6,7,8,9,10,   12,       16,   19,20,21]	;	1600	
;				startstr[0, 8] =	[0,1,      6,7,8,9,10,   12,       16,   19,20,21]	;	1700	

				for i = 0, 9 do begin		;--- loop over channels ----
					goodstarts = WHERE(startstr[*,i] ne -1, numgood)
					for j = 0, numgood - 2 do begin		;--- loop over start times for each channel ----
						if N_ELEMENTS(wavelist) eq 0 then begin
							wavelist 	= i
							tstartlist	= startstr[j, i]
							tstoplist	= startstr[j+1, i]
							polylist = polyperwave[i]
							if (j eq numgood - 2) then polylist = 0 ;	Don't extrapolate last point	; 2017/12/03, commented, no longer needed ; 2020/06/22 Uncommented, probably still prudent
						endif else begin
							wavelist 	= [wavelist, i]
							tstartlist	= [tstartlist, startstr[j, i]]
							tstoplist	= [tstoplist, startstr[j+1, i]]
							polylist	= [polylist, polyperwave[i]]
							if (j eq numgood - 2) then polylist[-1] = 0 ;	Don't extrapolate last point	; 2017/12/03, commented, no longer needed ; 2020/06/22 Uncommented, probably still prudent
						endelse
					endfor
				endfor

				fname		=	'aia_V9_'
			end

		else : 	begin
					BOX_MESSAGE, ['AIA_BP_RESPONSE_TABLE: /all only works with ver_num=1-9']
					STOP
				end
	endcase

	;--- Main loop: Call this routine (itself) for each interval of each wavelength -----------------------------------
	numentries 	=	N_ELEMENTS(wavelist)

	thiswave= -1	;initialize the current wave counter index: 0-8
	struct_def='={wave: wvls[thiswave], usetime_tai : usetime_all, useratios : useratios_all, yfit : yfit_all}'

	if KEYWORD_SET(p_stop) then STOP		
	for i = 0, numentries-1 do begin
		wave_done=0	; intialize

		;++++++++++++++++++++++++++++++++++++++
		; Version 9: Fix the last interval
		; to be flat at the end of the previous
		; interval
		;--------------------------------------
		if (ver_num eq 9) and (wavelist[(i+1) mod numentries] ne wavelist[i]) then begin
;			tmax = ANYTIM2TAI(thisout.t_stop) - ANYTIM2TAI(thisout.t_start)
;			ea_end = thisout.eff_area + thisout.effa_p1*tmax + $
;					 thisout.effa_p2*tmax^2. + thisout.effa_p3*tmax^3.
;			fit0 = 0 - (ea_end[0] / thisoutdat.ea_0)
;			STOP
			; fix to the scale factor at the end of the previous interval
			fit0 = 0-thisoutdat.yfit[N_ELEMENTS(thisoutdat.yfit)-1]
			thisout = AIA_BP_RESPONSE_TABLE(thisoutdat, wavelnth = wvls[wavelist[i]], $
				t_start = tstarts[tstartlist[i]], figpath = scratchpath, /nosmooth, $
				t_stop = tstarts[tstoplist[i]], poly = fit0, thicktoo = thicktoo, $
				ver_num = ver_num, silent = silent, fix_prebake = fix_prebake, $
				rocket_only=rocket_only, nopng=nopng, nofig=nofig, p_stop=p_stop gt 1)
		endif else begin
			if use_nodisc[tstartlist[i]] then begin
				nodisc = thisoutdat.yfit[N_ELEMENTS(thisoutdat.yfit)-1]
			endif else begin
				nodisc = 0
			endelse
			thisout = AIA_BP_RESPONSE_TABLE(thisoutdat, wavelnth = wvls[wavelist[i]], $
				t_start = tstarts[tstartlist[i]], figpath = scratchpath, nodisc = nodisc, $
				t_stop = tstarts[tstoplist[i]], poly = polylist[i], thicktoo = thicktoo, $
				ver_num = ver_num, silent = silent, fix_prebake = fix_prebake, $
				rocket_only=rocket_only, nopng=nopng, nofig=nofig, p_stop=p_stop gt 1)
		endelse
		if i eq 0 then begin
			outstr = thisout
			output_data = (thisoutdat.useratios - thisoutdat.yfit)
		endif else begin
			outstr = [outstr, thisout]
			output_data = [output_data, (thisoutdat.useratios - thisoutdat.yfit)]
		endelse

		;--- 2017/12/05, added to get all original output_data, not just the difference -----------------
		;note definition => thisoutdata = {usetime : usetime, useratios : useratios, yfit : yfit}

		if thiswave ne wavelist[i] then begin	; initialize for a new channel
			thiswave= wavelist[i]
				;lastwave= thiswave-1
			nextwave= thiswave+1

			usetime_all=thisoutdat.usetime
			useratios_all=thisoutdat.useratios
			yfit_all= thisoutdat.yfit
		endif else begin						; concatenate for the same channel
			usetime_all=[usetime_all, thisoutdat.usetime]
			useratios_all=[useratios_all, thisoutdat.useratios]
			yfit_all= [yfit_all, thisoutdat.yfit]
		endelse

		;--- output and plot each finished channnel --------------------------

		;--- identify finished channel ---------------
		if (i eq numentries-1) then begin	; the last channel
			wave_done=1
		endif else begin					; previous channels
			if (thiswave eq wavelist[i]) and (nextwave eq wavelist[i+1]) then $		; identify channel boundary
			wave_done=1
		endelse

		if wave_done then begin
			output_cmd= 'outdata_' + trim(wvls[thiswave]) + struct_def
				;help, output_cmd & pause	;debug
			tmp= execute( output_cmd )
			;if this cannot be resolved, just do a case block for all waves
			;% Attempt to call undefined procedure/function: 'OUTPUT_DATA_94__DEFINE'.


			if not KEYWORD_SET(nofig) then begin
				wset, 10
				instlabel_this='EVE'
				ylims = LIMITS([yfit_all, useratios_all])
				UTPLOT, TAI2UTC(usetime_all), useratios_all, yrange = ylims, chars = 1.5, title = 'AIA Degradation: Channel '+ trim(wvls[thiswave]), psym = -4, $
					ytitle = 'AIA / ' + instlabel_this + ' ratio', /xstyle
				OUTPLOT, TAI2UTC(usetime_all), yfit_all, col = 2
					;lstring = 'EA = EA_0 * (' + STRING(fits[0], format = '(f7.4)')
					;if N_ELEMENTS(fits) gt 1 then for i = 1, N_ELEMENTS(fits)-1 do $
					;	lstring = lstring + ' ' + STRING(fits[i], format = '(f9.5)') + ' dt^' + STRTRIM(i,2)
				rmse_this = MEAN(SQRT( (yfit_all - useratios_all)^2. ) )
				lstring= ' RMS of Date- Yfit = ' + trim(rmse_this)
				legend = ['Polynomial fit:', lstring + ')']
				PB_LEGEND, /dev, legend, 200, 200, size = 1.5

								;fname = STRING(usetime[0], format = '(i10)') + '-' + STRING(usetime[N_ELEMENTS(usetime)-1], format = '(i10)') + $
				trange_str= tai2utc(usetime_all[[0, N_ELEMENTS(usetime_all)-1]], /ccsds)
				fname_trange= time2file(trange_str, /sec)
				fname_this = fname_trange[0] + '-' + fname_trange[1] + $
					'_' + trim(wvls[thiswave]) + '_v' + STRTRIM(ver_num, 2) + '.png'
				if not KEYWORD_SET(nopng) then PB_WIN2PNG, scratchpath + '/png/' + '0_aia_bp_response_table_' + fname_this
			endif
		endif
		;--- end of output and plot each finished channnel --------------------------

	endfor
	;--- End of main loop -----------------------------------------------------

	numentries = N_ELEMENTS(outstr)
	if KEYWORD_SET(p_stop) then STOP
	
	;--- compile result and save output genx file, 2017/12/05, Tue -----------
	timetag = TIME2FILE(/sec, RELTIME(/now))
	struct_names= 'outdata_' + trim(wvls_used)
	all_struct_names= strjoin(struct_names+ ':' + struct_names+', ')
	strlength= strlen(all_struct_names)
	strput, all_struct_names, '}', strlength-2		; replace the trailing ',' with '}'
	def_cmd= 'outdata_all={'+ all_struct_names
	tmp=execute(def_cmd)
	file_genx= scratchpath + fname + timetag + '_response_table_out-structure.genx'
	savegen, str={outdata_all:outdata_all, outstr:outstr, output_data:output_data}, file=file_genx	; 2017/12/09

		;all_struct_names= '{'+ strjoin(struct_names+ ':' + struct_names+', ')
		;save_cmd= 'savegen, str={' + all_struct_names + 'outstr:outstr, output_data:output_data}, file=file_genx' 		;2017/11/29, Wei added, to save fit result etc.
		;tmp=execute(save_cmd)

	;--- If requested, generate output files ----------
	if KEYWORD_SET(printtab) then begin
			;timetag = TIME2FILE(/sec, RELTIME(/now))
		OPENW, txtlun, scratchpath + fname + timetag + '_response_table.txt', /get_lun
		OPENW, setlun, scratchpath + fname + timetag + '_response_table.set_keys', /get_lun
		otags = TAG_NAMES(outstr)
		PRINTF, txtlun, otags, format = '(3a26, 11a11)'
		for i = 0, numentries - 1 do begin
			if STRPOS(outstr[i].wave_str, 'THICK') lt 0 then begin
				; Don't write thick filter entries to the text file; only to the setkeys file
				PRINTF, txtlun, outstr[i], format = '(3a26, i11, a11, i11, f11.2, 2f11.5, f11.2, 4f11.5)'
			endif
			setstring = "set_keys -c ds='aia.response' "
			for j = 0, N_ELEMENTS(otags) - 2 do begin		;	Don't include RMSE
				thisval = STRTRIM(outstr[i].(j), 2)
				if otags[j] eq 'WAVE_STR' then thisval = "'" + thisval + "'" 
				setstring = setstring + otags[j] + '=' + thisval + ' '
			endfor
			PRINTF, setlun, setstring
		endfor
		FREE_LUN, txtlun
		FREE_LUN, setlun
	endif
	RETURN, outstr		; structure contain all entries/rows
endif


;++++++++++++++++++++++++++++++++++++++
; Get the AIA response functions that 
; were used by the EVE team to generate
; the predicted AIA band emissions
;--------------------------------------
aiaresp = AIA_GET_RESPONSE(/full, /all, version = 1, /phot)
resptags = TAG_NAMES(aiaresp)
fuvresp = AIA_GET_RESPONSE(/uv, /full, version = 1, /phot)
fuvresptags = TAG_NAMES(fuvresp)

;++++++++++++++++++++++++++++++++++++++
; Set up the defaults and result format
;--------------------------------------
tai_start = ANYTIM2TAI(t_start)
if N_ELEMENTS(t_stop) eq 0 then t_stop = TAI2UTC(tai_start + (3600d*24d*365d) )
tai_stop = ANYTIM2TAI(t_stop)

template_str = {date : 'a', t_start : 'a', t_stop : 'a', ver_num : 1, $
	wave_str : 'a', wavelnth : 1, $
	eperdn : 1.0, dnperpht : 1.0, $
	eff_area : 1.0, eff_wvln : 1.0, $
	effa_p1 : 0.0, effa_p2 : 0.0, effa_p3 : 0.0, $
	rmse : 0.0}

tags = TAG_NAMES(template_str)
ntags = N_TAGS(template_str)

if N_ELEMENTS(wavelnth) eq 1 then begin

	;++++++++++++++++++++++++++++++++++++++
	; Get the AIA and EVE data, select the
	; time range of interest, and do some
	; initial cleanup of bad points, etc.
	;--------------------------------------

	;--- depending on epoch intervals, determine and read appropriate ratios --------------
	case 1 of
		ver_num eq 9 : begin	;	For version 9, we're loading the ratios from a different source
								;	and much of the complexity in cleaning and parsing those
								; 	ratios has been farmed out to the comparison routine
			if wavelnth gt 1000 then begin
				AIA_SOLSTICE_COMPARE, /load, results, /silent, /mask, $
					file = '/sanhome/boerner/crom/disk2/data/solstice/20200620_145659_aia-solstice_ratios.genx' 
				chans = results.chans
				instlabel = "SOLSTICE"
			endif else begin
				results = READ_EUV_RATIOS()
				chans = results.chans
				instlabel = "EVE/FISM"
			endelse
		end
		(wavelnth lt 1000) and (tai_start lt eve_end_tai): begin	;--- for EUV and tai_start < eve_end_date, use aia/eve ratio --------------
				;--- Paul's original, changed below 2017/11/28 -----------
				;ff = FILE_SEARCH('/Users/boerner/data/eve/*EVLshort' + STRTRIM(veve, 2) + '_aia-eve_ratios*')
				;filetime = AIA_BP_DATE_STRING2TAI(STRMID(ff, 24, 15))
			ff = FILE_SEARCH(aia2eve_ratio_path+'/*EVLshort' + STRTRIM(veve, 2) + '_aia-eve_ratios*')
			filetime = AIA_BP_DATE_STRING2TAI(STRMID(file_basename(ff), 0, 15))		; 2017/11/28, Wei Liu, use file_basename() to get rid of the path, which is somewhat user-dependent.
			latestfile = ff[WHERE(filetime eq MAX(filetime))]	;use latest file
			instlabel = 'EVE'
			AIA_EVE_PLOT_RATIOS, /load, /silent, results, file = latestfile, /mask
			chans = results.chans
		end
		(wavelnth lt 1000) and (tai_start ge eve_end_tai):  begin	;--- for EUV and tai_start >= eve_end_date, use aia/model ratio ---------------
			ff = FILE_SEARCH(aia2model_ratio_path + '/*aia-*-ratio_*.genx')
			filetime = AIA_BP_DATE_STRING2TAI(STRMID(file_basename(ff), 0, 15))
			latestfile = ff[WHERE(filetime eq MAX(filetime))]	;use latest file
			instlabel = 'MODEL'
			restgen, str=results, file = latestfile
			chans = results.chans
		end
		(wavelnth ge 1000):  begin									;--- for UV, wavelnth >= 1000 \AA\ , use aia/see ratio ------------------------
				;--- Paul's original, changed below 2017/11/28 -----------
				;ff = FILE_SEARCH('/Users/boerner/data/see/*aia-see_ratios*')
				;filetime = AIA_BP_DATE_STRING2TAI(STRMID(ff, 24, 15))
			ff = FILE_SEARCH(aia2see_ratio_path + '/*aia-see_ratios*')
			filetime = AIA_BP_DATE_STRING2TAI(STRMID(file_basename(ff), 0, 15))
			latestfile = ff[WHERE(filetime eq MAX(filetime))]	;use latest file
			instlabel = 'SEE'
			AIA_SEE_COMPARE, /load, /silent, results, file = latestfile, /mask
			chans = STRMID(results.chans, 1, 4)
		end
	endcase

	;--- now, do initial cleanup of bad points, etc. -----------------
	tai_data = ANYTIM2TAI(results.date)
	intime = WHERE(tai_data ge tai_start and tai_data le tai_stop)
	inband = WHERE(chans eq wavelnth)

	if (wavelnth lt 1000) or (ver_num eq 9) then begin		; for aia/eve or aia/solstice ratio, [nwave, ntime]		;  2017/11/29, commented: and (tai_start lt eve_end_tai)
		bandratios = results.ratios[inband[0], intime]
		if ver_num eq 9 then begin
			goodratios = WHERE(results.ratios[inband[0], intime] gt 0) 	; take only > 0 ratios; doesn't catch negative-negative points
		endif else begin
			goodratios = WHERE(results.eve[inband[0], intime] gt 0) 	; take only > 0 ratios
		endelse
	endif else begin
		bandratios = results.ratios[intime, inband[0]]				; for aia/see ratio, [ntime, nwave]
		goodratios = WHERE(results.ratios[intime, 0] gt 0)
	endelse

	useratios = bandratios[goodratios]

	;++++++++++++++++++++++++++++++++++++++
	; Version 8: hard-coded scaling of
	; FISM to agree with EVE
	;--------------------------------------
	if (ver_num eq 8) and (wavelnth lt 1000) and (tai_start ge eve_end_tai) then begin		; for EUV beyond EVE/MEGS-A failure on eve_end_date='25-may-2014 12:00'

		;--- 2017/11/29, copy/paste from /net/sunflower/disk2/scr0/weiliu/sdo/docs/my-job/calib/eve/wkspc/run01_20170926_FISM-V2_PhilC/v1a_batch/cmpr_real-eve_v2_scale-model/aia_eve_modeled_compare.dat ----
		;Channels, Linear Corr.Coef. & Uncertainty, (AIA/EVE)/(AIA/Modeled): Median, Mean;       Linfit: A,      B;                 Fitexy: A,     B
		model2eve_scale= [  $
		[ 94d ,         0.61584041 ,     0.025932508 ,       1.0555048 ,       1.0556922 ,      0.47111445 ,      0.51337972 ,      0.19038501 ,      0.83314272 ], $
		[ 131 ,         0.84323426 ,     0.025941230 ,       1.0392760 ,       1.0375196 ,      0.23657549 ,      0.73133414 ,      0.13084770 ,      0.86721330 ], $
		[ 171 ,         0.74055215 ,     0.025932508 ,      0.78199676 ,      0.78023902 ,      0.61757295 ,      0.34164870 ,      0.44716241 ,      0.46168375 ], $
		[ 193 ,         0.33682442 ,     0.025932508 ,      0.99951437 ,       1.0069643 ,      0.92267919 ,      0.18695042 ,      0.50721890 ,      0.55463748 ], $
		[ 211 ,         0.58473927 ,     0.025932508 ,       1.0011734 ,       1.0145264 ,      0.49462705 ,      0.48690817 ,      0.16817678 ,      0.83229673 ], $
		[ 304 ,         0.99771736 ,     0.025932508 ,      0.99075332 ,      0.99480978 ,    -0.021722123 ,       1.0309228 ,    -0.022775077 ,       1.0325413 ], $
		[ 335 ,         0.97064818 ,     0.025941230 ,       1.1047379 ,       1.1127064 ,    -0.020537844 ,       1.1445834 ,    -0.045512218 ,       1.1783370 ] ]
		ifit=[8,9]-1	; array index of the above linear fit result to be used

		useratios = model2eve_scale[ifit[0], inband[0]] + useratios * model2eve_scale[ifit[1], inband[0]]		; scale aia/model ratio to be a proxy of aia/eve ratio

			;--- 2017/11/29, changed to aia/fism ratio, use linear correlation instead ----------------
			;	  ;RESTGEN, file = '~/data/aia/cal/20141024_231847_eve_see_ratios.genx', eve_see_ratios	; changed below 2017/11/28 
			;seefile= '/sanhome/weiliu/data/aia/cal/20141024_231847_eve_see_ratios.genx'
			;RESTGEN, file = seefile, eve_see_ratios
				;;IDL> print, eve_see_ratios
				;;      0.19005049     0.035992068      0.28383121      0.87576466      0.75325248      0.76398102       1.7181232
			;useratios = useratios * eve_see_ratios[inband[0]]		; scale aia/see ratio to be a proxy of aia/eve ratio
	endif
		
	;++++++++++++++++++++++++++++++++++++++
	; Smooth the data a little bit
	;--------------------------------------
	if not KEYWORD_SET(nosmooth) then begin
			;rawratios = useratios	;never used
		if (tai_start lt eve_end_tai) then ismooth=0 $
			else ismooth=1
		smoothratios = MEDIAN(useratios, smooth_width[ismooth])		;Median smoothing replaces each point with the median of the neighborhood of a given width.
	
		;--- replace outliers beyond smooth_threshold with the median of the running window of smooth_width ----
		unsmooth = WHERE( ABS(useratios - smoothratios) gt smooth_threshold[ismooth], count)
		if count gt 0 then useratios[unsmooth] = smoothratios[unsmooth]
			;if not KEYWORD_SET(silent) then $
		if KEYWORD_SET(loud) then $
		PRINT, 'Number of points smoothed: ', count, format = '(a45, i10)'

		;--- 2017/12/09, more aggressive smoothing after EVE failure, do the same for UV channels, when TIMED/SEE data gets noisy; Result: not much change, since polyfit already takes care of noisy data -----
		if keyword_set(post_eve_heavy_smooth) and (tai_start ge eve_end_tai) then begin
			useratios = smoothratios
		endif
	endif
	usetime = tai_data[intime[goodratios]]
	
	;++++++++++++++++++++++++++++++++++++++
	; Filter the data to exclude some known
	; periods of funny AIA or EVE data
	;--------------------------------------
	if not KEYWORD_SET(noevebake) then begin
		; Eliminate data taken during the 2 weeks following each of the EVE bakeouts
		rawtimes = usetime
		eb1_start = ANYTIM2TAI('2010-06-18T12:00:00Z')
		eb1_stop = ANYTIM2TAI('2010-07-01T12:00:00Z')
		eb2_start = ANYTIM2TAI('2010-09-24T12:00:00Z')
		eb2_stop = ANYTIM2TAI('2010-10-11T12:00:00Z')
		eb3_start = ANYTIM2TAI('2012-03-11T12:00:00Z')
		eb3_stop = ANYTIM2TAI('2012-04-11T12:00:00Z')
		usetime_ind = WHERE( 	(usetime lt eb1_start) or (usetime gt eb3_stop) or $
							(usetime gt eb1_stop and usetime lt eb2_start) or $
							(usetime gt eb2_stop and usetime lt eb3_start) , count )
		if count gt 0 then begin
			usetime = usetime[usetime_ind]
			useratios = useratios[usetime_ind]
		endif else STOP
	endif

	if not KEYWORD_SET(no304bake) and (wavelnth eq 304 or wavelnth eq 94) then begin
		; Eliminate the bad datapoints taken during the ATA4 bakeouts (on 2011-05-19 or 2011-10-04-6)
		badtimes = ['2011-05-19T19:00:00Z', '2011-10-04T19:00:00Z', '2011-10-05T19:00:00Z']
		for i = 0, N_ELEMENTS(badtimes) - 1 do begin
			usetime_ind = WHERE(ABS(usetime - ANYTIM2TAI(badtimes[i])) ge 86400d, count)
			if count gt 0 then begin
				usetime = usetime[usetime_ind]
				useratios = useratios[usetime_ind]
			endif else STOP
		endfor
	endif
	
	;++++++++++++++++++++++++++++++++++++++
	; Now do the polynomial fit to the data
	; Leave out the first 48 hours and last
	; 24 hours of the interval from the fit
	;--------------------------------------
	fittime = (usetime - tai_start)/86400d
	fitend = MAX(fittime) - 1
	fitstart = MIN(fittime) + 2

	if KEYWORD_SET(fix_prebake) and (wavelnth eq 304) then begin
		; Eliminate the data prior to the first EVE bakeout from the fit (not the plots)
		badtime = (ANYTIM2TAI('2010-09-24T00:00:00Z') - tai_start) / 86400d
		fitstart = badtime > fitstart
	endif

	fitindex = WHERE(fittime ge fitstart and fittime le fitend, numfitind)
	
	; Set up the parameter structure for using MPFITFUN 
	numparams = 4
	params = ['Offset', 'Linear', 'Quadratic', 'Cubic']
	values = [1.0d, 0.01d, 0d, 0d]
	fixed = [0, 0, 1, 1]
	limited = [0, 0, 0, 0]
	limits = [1, 1, 1, 1]
	parinfo = REPLICATE({name:params[0], value:values[0], fixed:fixed[0], $
		limited:[1,1]*limited[0], limits:[-1d,1d]*limits[0]}, numparams)
	for i = 1, numparams-1 do begin
		parinfo[i].name 	= params[i]
		parinfo[i].value 	= values[i]
		parinfo[i].fixed 	= fixed[i]
		parinfo[i].limited 	= limited[i] * [1,1]
		parinfo[i].limits 	= limits[i] * [-1d, 1d]
	endfor

	case 1 of
		poly lt 0 : $
			begin
				fits = ABS(poly) 
			end
		KEYWORD_SET(nodisc) : $
			begin
				if poly eq 0 then begin
					fits = nodisc
				endif else begin
					parinfo[0].value = nodisc
					parinfo[0].fixed = 1
					fits = MPFITFUN('AIA_MP_RESPTABLE_CHISQ', fittime[fitindex], $
									useratios[fitindex], err, $
									parinfo = parinfo, yfit = thisyfit, quiet=1 )
;						derivratios = DERIV(useratios[fitindex]-useratios[fitindex[0]])
;						bfits = POLY_FIT(fittime[fitindex], derivratios, poly-1, /double)
;						ofits = POLY_FIT(fittime[fitindex], useratios[fitindex], poly, /double)
;						bfits[0] = bfits[0] + (ofits[0]-nodisc)/(MAX(fittime[fitindex])-MIN(fittime[fitindex]))
;						fits = [nodisc, bfits]
				endelse
			end
		else : $
			begin
				fits = MPFITFUN('AIA_MP_RESPTABLE_CHISQ', fittime[fitindex], $
								useratios[fitindex], err, $
								parinfo = parinfo, yfit = thisyfit, quiet=1 )
;				fits = POLY_FIT(fittime[fitindex], useratios[fitindex], poly, /double)
				if poly gt 0 then begin
				;	We don't want any positive slopes...  ; Commented out 2020-06; positive slopes are ok now
		;			if fits[1] gt 0 then begin
		;				fits = POLY_FIT(fittime[fitindex], useratios[fitindex], 0, /double)
		;			endif
				endif
			end
	endcase
	; Note that the time-dependent polynomial coefficients need to be divided by the
	; 0th order coefficient, because we want change in time WRT the EA at the start of
	; the epoch, not WRT the "nominal" EA
	pfits = fits	;	Hold the "raw" coefficients 
	if N_ELEMENTS(fits) gt 1 then fits[1:*] = fits[1:*] / fits[0]

	;++++++++++++++++++++++++++++++++++++++
	; Generate a plot showing the data and
	; the fit over the interval
	;--------------------------------------	
	yfit = useratios * 0.
	for i = 0, N_ELEMENTS(fits)-1 do yfit = yfit + pfits[i] * fittime^i
	rmse = MEAN(SQRT( (yfit[fitindex] - useratios[fitindex])^2. ) )
	
		;if not KEYWORD_SET(silent) then begin
	if not KEYWORD_SET(nofig) then begin
		wset, 9
		ylims = LIMITS([yfit, useratios])
		UTPLOT, TAI2UTC(usetime), useratios, yrange = ylims, chars = 1.5, title = wavelnth, psym = -4, $
			ytitle = 'AIA / ' + instlabel + ' ratio', /xstyle
		OUTPLOT, TAI2UTC(usetime), yfit, col = 2
		lstring = 'EA = EA_0 * (' + STRING(fits[0], format = '(f7.4)')
		if N_ELEMENTS(fits) gt 1 then for i = 1, N_ELEMENTS(fits)-1 do $
			lstring = lstring + ' ' + STRING(fits[i], format = '(f9.5)') + ' dt^' + STRTRIM(i,2)
		legend = ['Polynomial fit:', lstring + ')']
		PB_LEGEND, /dev, legend, 200, 200, size = 1.5

						;fname = STRING(usetime[0], format = '(i10)') + '-' + STRING(usetime[N_ELEMENTS(usetime)-1], format = '(i10)') + $
		trange_str= tai2utc(usetime[[0, N_ELEMENTS(usetime)-1]], /ccsds)
		fname_trange= time2file(trange_str, /sec)
		fname_this = fname_trange[0] + '-' + fname_trange[1] + $
			'_' + STRTRIM(wavelnth, 2) + '_v' + STRTRIM(ver_num, 2) + '.png'
		if not KEYWORD_SET(nopng) then PB_WIN2PNG, figpath + '/png/' + 'aia_bp_response_table_' + fname_this
	endif

	if KEYWORD_SET(p_stop) then STOP

	;++++++++++++++++++++++++++++++++++++++
	; Look up the response structure in 
	; order to determine the CCD gain 
	; parameters and the "nominal" EA at 
	; the effective wavelength
	;--------------------------------------	
	chantag = 'A' + STRTRIM(wavelnth,2) + '_FULL'
	wvlind = WHERE(wvls eq wavelnth, count)
	if count eq 1 then wvlind = wvlind[0] else STOP
	if wavelnth lt 1000 then begin
		chanstrind = WHERE(resptags eq chantag, count)
		if count eq 1 then chanstr = aiaresp.(chanstrind[0]) else STOP
		wavestring = STRTRIM(wavelnth,2) + '_THIN'
	endif else begin
		chanstrind = WHERE(fuvresptags eq chantag, count)
		if count eq 1 then chanstr = fuvresp.(chanstrind[0]) else STOP
		wavestring = STRTRIM(wavelnth,2)
	endelse

	;++++++++++++++++++++++++++++++++++++++
	; Generate the output structure
	;--------------------------------------		
	thisstr = template_str
;	thisstr.date = TAI2UTC(/ccsds, tai_start)
	thisstr.date = ANYTIM(RELTIME(/now), /ccsds)
	thisstr.t_start = TAI2UTC(/ccsds, tai_start)
	thisstr.t_stop = TAI2UTC(/ccsds, tai_stop)
	thisstr.ver_num = ver_num
	
	thisstr.wave_str = wavestring
	thisstr.wavelnth = wavelnth
	thisstr.eff_wvln = ewvls[wvlind]

	thisstr.eperdn = chanstr.elecperdn
	thisstr.dnperpht = 12398. / thisstr.eff_wvln / 3.65 / thisstr.eperdn
	thisstr.rmse = rmse
	
	ea_0 = INTERPOL(chanstr.effarea, chanstr.wave, thisstr.eff_wvln)
	thisstr.eff_area = ea_0 * fits[0]
	if N_ELEMENTS(fits) gt 1 then thisstr.effa_p1 = fits[1]
	if N_ELEMENTS(fits) gt 2 then thisstr.effa_p2 = fits[2]
	if N_ELEMENTS(fits) gt 3 then thisstr.effa_p3 = fits[3]

	;++++++++++++++++++++++++++++++++++++++
	; Also pass out the fit data as an 
	; optional output parameter
	;--------------------------------------		
	output_data = {usetime : usetime, useratios : useratios, yfit : yfit, ea_0 : ea_0}
	
	;++++++++++++++++++++++++++++++++++++++
	; Print a short summary of the fit
	;--------------------------------------			
	if not KEYWORD_SET(silent) then begin
			;if KEYWORD_SET(loud) then begin
		PRINT, 'wave_str, t_start, t_stop: ', thisstr.wave_str, thisstr.t_start, thisstr.t_stop, format='(10(a, 2x))'
		PRINT, 'Effective area at central wavelength [cm2]:', ea_0, format = '(a45, f10.5)'
		PRINT, 'Best-fit polynomial coefficients:', fits, format = '(a45, ' + STRTRIM(N_ELEMENTS(fits), 2) + 'f10.5)'
		if keyword_set(loud) then begin
			PRINT
			HELP, thisstr
		endif ;else begin
		;	help, thisstr.wave_str, thisstr.t_start, thisstr.t_stop
		;endelse
	endif

	;++++++++++++++++++++++++++++++++++++++
	; Generate the output for the thick
	; focal plane filter as well
	;--------------------------------------			
	if thicktoo and wavelnth lt 1000 then begin
		thickstr = thisstr
		thick_chantag = 'A' + STRTRIM(wavelnth,2) + '_THICK_FULL'
		thick_chanstrind = WHERE(resptags eq thick_chantag, count)
		if count eq 1 then thick_chanstr = aiaresp.(thick_chanstrind[0]) else STOP
		thick_ea0 = INTERPOL(thick_chanstr.effarea, thick_chanstr.wave, thickstr.eff_wvln)
		thick_over_thin = thick_ea0 / ea_0
		thickstr.eff_area = thisstr.eff_area * thick_over_thin
		thickstr.wave_str = STRTRIM(wavelnth, 2) + '_THICK'
		thisstr = [thisstr, thickstr]
	endif
	
endif

RETURN, thisstr		; structure contain one entry/row for one time interval of one wavelength

end