PRO JAN20_511OPTIMIZER, VERBOSE=verbose
;
; Optimizes selection of 511 keV energy/time range for imaging Jan 20, 2005 event.
;
; 2-June-05 gh
;
; /VERBOSE prints some of the details
;
; Load input light curve data for 511 line
; Solar line flux, influx, (from GS 5-31-05 email) is specified in irregular time bins centered at 'intime'
; intime=time in minutes after 06:00
intime = ([419,425,431,437,443,449,455,461,470,482,494,506,518,530,542,554,566,578,590,602,617,637.5,649.5,660.5,673]*10+20000-6*3600)/60
influx = [-49,167,350,905,670,822,665,569,510,388,246,240,261,275,168,253,196,130,126,225,202,188,  190,  191,  178]*.001
; PLOT, intime, influx, XTITLE='minutes after 06:00', YTITLE='gammas/cm2/s', PSYM=6
;
; Interpolate into uniform 60 second bins
; time = output time base, in 60 second bins, starting at 06:43.  Times are at midbin.
; flux = line flux, interpolated to 'time'
ntime=44
time = FINDGEN(ntime)+43.5				; 06:43 - 07:27
flux = INTERPOL(influx, intime, time)
OPLOT, time, flux
; FOR n=0, ntime-1 DO PRINT, time[n], flux[n]
;
; Load line spectral data for use before/after 07:08
kev = FINDGEN(20)+500.5
spectpeak=[0.031,0.037,0.038,0.037,0.038,0.045,0.046,0.061,0.065,0.078,0.082,0.078,0.077,0.065,0.056,0.039,0.036,0.033,0.031,0.026]
spectlate=[0.021,0.025,0.025,0.029,0.024,0.034,0.031,0.056,0.071,0.115,0.111,0.103,0.104,0.087,0.056,0.035,0.022,0.024,0.014,0.012]
gm = 0.048		; gridtran*drm for rear segments
effarea = 40*gm
linecount = FLTARR(20, ntime)
IF KEYWORD_SET(verbose) THEN BEGIN
		PRINT
		PRINT, 'TIME   ENERGIES-->      LINE COUNT MATRIX'
ENDIF
FOR n=0, ntime-1 DO BEGIN
	IF time(n) LT 68 THEN spect = spectpeak ELSE spect=spectlate
	linecount(*,n) = effarea * spect * flux(n) * 60					; *60 to convert gammas/s to gammas in 1 minute
	IF KEYWORD_SET(verbose) THEN PRINT, time[n], linecount(*,n),FORMAT='(F4.1, 20F6.2)'
ENDFOR
;
; Load count data
; csp = count spectra
csp=[606.9, 640.2, 634.4, 744.3, 758.6, 842.1, 967.8,1216.6,1608.2,2137.4,2466.5,2480.6,2188.5,1445.6, 835.2, 535.5, 360.9, 340.4, 353.3, 318.9]
relcsp = csp/TOTAL(csp)
; cpm = counts per minute between 500 and 520 keV
cpm = [522.7, 558.3, 757.3, 843.1, 781.9, 752.8, 744.3, 668.4, 619.5, 585.1, 528.4, 493.3, 527.6, 465.6, 486.8, $
 		468.3, 457.4, 455.4, 433.9, 414.9, 389.2, 458.1, 405.1, 457.1, 426.6, 403.3, 446.8, 408.8, 387.3, 414.3, $
 		413.4, 404.7, 405.1, 436.2, 455.0, 431.4, 437.8, 433.2, 443.7, 463.7, 432.0, 462.6, 501.2, 412.2]
totalcount = FLTARR(20, ntime)
FOR n=0, ntime-1 DO BEGIN
	totalcount[*,n] = cpm[n]*relcsp
	IF KEYWORD_SET(verbose) THEN PRINT, time[n], totalcount[*,n], format='(f4.1, 20F6.2)'
ENDFOR
;
; Generate plots
!P.MULTI=[0,1,1]
PLOT, time, TOTAL(linecount, 1)/MAX(TOTAL(linecount,1)),  PSYM=7, YRANGE=[0,1.2], $
	XTITLE='MINUTES AFTER 06:00', yTITLE='RELATIVE COUNTS', TITLE = '511 keV LIGHT CURVE COMPARISON'
OPLOT, time, TOTAL(totalcount, 1)/MAX(TOTAL(totalcount,1)), PSYM=10
XYOUTS, 60, 1.0, 'X Bkgd-subtracted predictions'
XYOUTS, 60, 1.1, '- Observed counts including bkgd'
PLOT, kev, TOTAL(linecount,2)/MAX(TOTAL(linecount,2)), YRANGE=[0, 1.2], $
	XTITLE='ENERGY (keV)', YTITLE='RELATIVE COUNTS', PSYM=10, TITLE='511 SPECTRA COMPARISON'
OPLOT, kev, spectpeak/MAX(spectpeak), PSYM=7
OPLOT, kev, spectlate/MAX(spectlate), PSYM=6
XYOUTS, 501, 1.14, '- Avg observed spectrum, including bkgd'
XYOUTS, 501, 1.09, 'X Bkgd-subtracted prediction at peak'
XYOUTS, 501, 1.04, 'O Bkgd-subtracted prediction late'
;
; Do optimization
IF KEYWORD_SET(verbose) THEN BEGIN
	PRINT
	PRINT, 'TIME   ENERGIES-->      OBSERVED COUNT MATRIX'
ENDIF
PRINT
estart = [8,13]					; 508-514 keV
tstart = [13,42]				; 06:56-07:26
ratio 	= FLTARR(9)
e		= INTARR(2,9)
t 		= INTARR(2,9)
e[*,0] = estart
t[*,0] = tstart
FOR nit=0, 50 DO BEGIN
	e[*,1] = e[*,0]
	e[*,2] = e[*,0]
	e[*,3] = e[*,0]
	e[*,4] = e[*,0]
	e[*,5] = (e[*,0] +[-1, 0]) >0 <19
	e[*,6] = (e[*,0] +[ 1, 0]) >0 <19
	e[*,7] = (e[*,0] +[ 0,-1]) >0 <19
	e[*,8] = (e[*,0] +[ 0, 1]) >0 <19
	t[*,1] = (t[*,0] +[-1, 0]) >0 <(ntime-1)
	t[*,2] = (t[*,0] +[ 1, 0]) >0 <(ntime-1)
	t[*,3] = (t[*,0] +[ 0,-1]) >0 <(ntime-1)
	t[*,4] = (t[*,0] +[ 0, 1]) >0 <(ntime-1)
	t[*,5] = t[*,0]
	t[*,6] = t[*,0]
	t[*,7] = t[*,0]
	t[*,8] = t[*,0]
	FOR i = 0, 8 DO ratio[i] =   TOTAL(linecount [e[0,i]:e[1,i], t[0,i]:t[1,i]]) $
						  / SQRT(TOTAL(totalcount[e[0,i]:e[1,i], t[0,i]:t[1,i]]))
	best = MAX(ratio, ibest)
	line0 	= TOTAL ( linecount[e[0,0]:e[1,0], t[0,0]:t[1,0]])
	total0 	= TOTAL (totalcount[e[0,0]:e[1,0], t[0,0]:t[1,0]])
	PRINT, nit+1, e[*,0], t[*,0], line0, total0, ratio, FORMAT='(5I4, 2F10.2, 9F9.4)'
	IF ibest EQ 0 THEN STOP
	e[*,0] = e[*,ibest]
	t[*,0] = t[*,ibest]
ENDFOR
STOP
END