PRO focus_trends, DO_PS = DO_PS, $
	FT=FT, SHORT_LABEL=short_label, CAPELLA=capella, $
	NOM_AIM=nom_aim
;+
;  PRO focus_trends, DO_PS = DO_PS, $
;	FT=FT, SHORT_LABEL=short_label, CAPELLA=capella, $
;	NOM_AIM=nom_aim
;
; Plot the streak zeroth-order core FWHM fits vs observation date
;
; Can follow this with ft_chipx_plot:
;  IDL> focus_trends,/NOM_AIM, ft=ft
;  IDL> @ft_chipx_plot
;-

; Output plot into ps file
if KEYWORD_SET(DO_PS) then pre_print_landscape
!p.multi = 0

; The data is in an rdb format file, focus_trends.rdb
if STRPOS(!DDLOCATION, 'HAK') LT 0 then begin 
        ; dd's home location of the file:
        ft = rdb_read(!DDIDL+'/marx/focus_trends.rdb')
end else begin
        ; HAK Stand-alone location of the file
        ;;ft = rdb_read(!DDHAKCODE+'/../hak_data/focus_trends.rdb')
        ; or take it from local directory:
        ft = rdb_read('focus_trends.rdb')
end


; Remove any that have negative values of strkFWHM:
posfwhm = where(ft.strkFWHM GT 0.0)
ft = ft(posfwhm)

; Keep only Capella ones if selected...
if KEYWORD_SET(CAPELLA) then begin
  caps = where(strpos(ft.Object,'apella') GT 0)
  ft = ft(caps)
end

; Keep only ones near the nominal +/- 20" aim points...
if KEYWORD_SET(NOM_AIM) then begin
  offsetaim = ABS(ft.chipx - 246.0)
  noms = where( offsetaim GT 25.0 AND offsetaim LT 60.0 )
  ft = ft(noms)
end

; Convert the 90% confidence errors into 1-sigma:
ft.strkFWHMerr = ft.strkFWHMerr/1.65

; Setup the Julian dates along xaxis
juldayval = LONG(indgen(n_elements(ft))+1.0)
for id=0, n_elements(ft)-1 do begin
  juldayval(id) = JULDAY(ft(id).month, ft(id).day, ft(id).year)
end

; Plot the FWHM values vs Date

jul_dates = 0

if jul_dates EQ 1 then begin
; - - -
dummy = LABEL_DATE(DATE_FORMAT = '%N/%D/%Y')
plot, juldayval, ft.strkFWHM, PSYM=4, $
	XRANGE=[min(juldayval)-20.0,max(juldayval)+20.0], XSTYLE=1, $
	XTICKFORMAT = 'LABEL_DATE', $
	XTITLE='Observation Date ( ' + $
	STRCOMPRESS(FIX(ft(0).month),/REMOVE) + '/' + $
		STRCOMPRESS(FIX(ft(0).year),/REMOVE) + ' to ' + $
	STRCOMPRESS(FIX(ft(n_elements(ft)-1).month),/REMOVE) + '/' + $
		STRCOMPRESS(FIX(ft(n_elements(ft)-1).year),/REMOVE) + $
	' )', $
	YRANGE=[0.,55.0], YSTYLE=1, $
	YTITLE='FWHM from streak events (um)', $
	TITLE='Zeroth-order 1D FWHM vs Time'
plot_errors, juldayval, juldayval, juldayval, $
	ft.strkFWHM - ft.strkFWHMerr, ft.strkFWHM + ft.strkFWHMerr
; - - -
end else begin
  julrefdate = JULDAY(1,1,1999)
  juldayval = DOUBLE(juldayval - julrefdate)
  yrval = 1999.0 + juldayval/365.25

  plot, yrval, ft.strkFWHM, PSYM=4, $
	XRANGE=[1999.5, 2007.0], XSTYLE=1, $
	XTITLE='Observation Date ( ' + $
	STRCOMPRESS(FIX(ft(0).month),/REMOVE) + '/' + $
		STRCOMPRESS(FIX(ft(0).year),/REMOVE) + ' to ' + $
	STRCOMPRESS(FIX(ft(n_elements(ft)-1).month),/REMOVE) + '/' + $
		STRCOMPRESS(FIX(ft(n_elements(ft)-1).year),/REMOVE) + $
	' )', $
	YRANGE=[0.,55.0], YSTYLE=1, $
	YTITLE='FWHM from streak events (um)', $
	TITLE='Zeroth-order 1D FWHM vs Time'
  plot_errors, yrval, yrval, yrval, $
	ft.strkFWHM - ft.strkFWHMerr, ft.strkFWHM + ft.strkFWHMerr
end

; Calc and print the weighted mean value
print, ''
print, ' Weighted Mean:'
vals = ft.strkFWHM
errs = ft.strkFWHMerr
; weights are 1/sigma^2
weights = 1.0/(errs)^2
mean = TOTAL(weights*vals)/TOTAL(weights)
mean_err = SQRT(1.0/TOTAL(weights))
print, 'mean and error = ', mean, mean_err
; and Chi-square
chi2n = TOTAL( ((ft.strkFWHM-mean)/ft.strkFWHMerr)^2 ) / $
	FLOAT(n_elements(ft)-1)
print, ' chi-squared / N-1  = '+string(chi2N)
print, ''

if jul_dates EQ 1 then begin
; - - -
; Add the mean to the plot
oplot, [0.0,max(juldayval)+1.0], mean*[1.0,1.0], LINESTYLE=2
oplot, [0.0,max(juldayval)+1.0], (mean-mean_err)*[1.0,1.0], LINESTYLE=1
oplot, [0.0,max(juldayval)+1.0], (mean+mean_err)*[1.0,1.0], LINESTYLE=1
xyouts, juldayval(0),50.0, 'Ave FWHM = '+STRING(mean,FORMAT='(F5.1)') + $
	' +/- ' + STRING(mean_err,FORMAT='(F5.2)') + ' um'
xyouts, juldayval(0),48.0, 'Chi^2 / (N-1) = '+STRING(chi2n,FORMAT='(F5.2)')

; and add observation info too
lab_offset = 3.0
for id=0,n_elements(ft)-1 do begin
  if KEYWORD_SET(SHORT_LABEL) then begin
    label_str = strpad(ft(id).obsid,5,/RIGHT) + ' - '
    ; and a "*" if the events were filtered for low E
    if STRPOS(ft(id).object,'<') GT 0 then begin
      label_str = label_str + '*'
    end
  end else begin
    label_str = strpad(ft(id).obsid,5,/RIGHT) + ' - ' + $
	ft(id).object + ' ('+ft(id).PItype+')'
  end
  xyouts, juldayval(id)-0.5, lab_offset, label_str, $
	orientation=90.0, CHARSIZE=0.8
  if lab_offset LT 10.0 then lab_offset=17.0 else lab_offset=3.0
end
; - - -
end else begin

  ; Add the mean to the plot
  oplot, [1999.0,max(yrval)+1.0], mean*[1.0,1.0], LINESTYLE=2
  oplot, [1999.0,max(yrval)+1.0], (mean-mean_err)*[1.0,1.0], LINESTYLE=1
  oplot, [1999.0,max(yrval)+1.0], (mean+mean_err)*[1.0,1.0], LINESTYLE=1
  xyouts, yrval(0),50.0, 'Ave FWHM = '+STRING(mean,FORMAT='(F5.1)') + $
	' +/- ' + STRING(mean_err,FORMAT='(F5.2)') + ' um'
  xyouts, yrval(0),48.0, 'Chi^2 / (N-1) = '+STRING(chi2n,FORMAT='(F5.2)')

  ; and add observation info too
  lab_offset = 3.0
  for id=0,n_elements(ft)-1 do begin
    if KEYWORD_SET(SHORT_LABEL) then begin
      label_str = strpad(ft(id).obsid,5,/RIGHT) + ' - '
      ; and a "*" if the events were filtered for low E
      if STRPOS(ft(id).object,'<') GT 0 then begin
        label_str = label_str + '*'
      end
    end else begin
      label_str = strpad(ft(id).obsid,5,/RIGHT) + ' - ' + $
	ft(id).object + ' ('+ft(id).PItype+')'
    end
    xyouts, yrval(id)-0.002, lab_offset, label_str, $
	orientation=90.0, CHARSIZE=0.8
    if lab_offset LT 10.0 then lab_offset=17.0 else lab_offset=3.0
  end

end

; label the plot...
plot_creator_label, '~/idl/marx/focus_trends.pro'

if KEYWORD_SET(DO_PS) then begin
  device, /close
  SPAWN, 'cp idl.ps focus_trends.ps'
  set_plot, 'X'
end

RETURN
END