PRO beam_cen_plots
;
; Plot the data in beam_cen_results.rdb

; read the results in
results = rdb_read('beam_cen_results.rdb')

; get the names of the gratings and current
; values of period, angle.
xrcf_params, NAM = names
xrcf_params, ANG = angles
xrcf_params, PER = periods
xrcf_params, ROW = rowland
xrcf_params, HC = hc

pre_print_portrait
!p.multi = [0,1,3]

; Loop over the grating types and
;   - calculate the XRCF/AL-K period
;   - make plots
for ig = 0, n_elements(names)-1 do begin
  ; calculate the XRCF weighted average period and error
  xrcf_grat = where(results.grating EQ names(ig) AND $
			STRPOS(results.phase,'XRCF') GE 0  AND $
			STRPOS(results.source,'Al-K') GE 0, nxrcf)
  if nxrcf GE 0 then begin
    xps = results(xrcf_grat).p
    xperrs = results(xrcf_grat).perr
    ; Weighted mean and error
    ; from Bevington 2nd, (4.17) and (4.19)
    xavep = TOTAL(xps/(xperrs)^2)/TOTAL(1./(xperrs)^2)
    xaveperr = SQRT(1./TOTAL(1./(xperrs)^2))
    print, xavep, xaveperr
  end

  ; find this type of measurement
  this_grat = where(results.grating EQ names(ig), nfound)
  if nfound GE 1 then begin
    ; Periods and period errors
    ps = results(this_grat).p
    avep = TOTAL(ps)/n_elements(ps)
    perrs = results(this_grat).perr
    xs = 0.5+indgen(n_elements(ps))
    ; Set the Y range to be +/- 1 part in 1500 for
    ; LEG and MEG, +/- 1 part in 1500 for HEG
    partin = 1500.
    if names(ig) EQ 'HEG' then partin = 1500.
    plot, xs, ps, PSYM=4, $
	YRANGE=avep*[partin/(partin+1.), (partin+1.)/partin], $
	XTITLE = 'Measurement set', YTITLE='Period (A)', $
	TITLE = names(ig)+' Period = '+ $
		STRING(xavep,FORMAT='(F8.2)')+' +/-'+ $
		STRING(xaveperr,FORMAT='(F5.2)') + $
		' A  (XRCF/Al-K data)', $
	CHARSIZE=2.5
    ; add error bars
    plot_errors, xs, xs, xs, ps-perrs, ps+perrs
    ; plot the XRCF-AL-K average period and error range
    oplot, [2.,n_elements(ps)], xavep*[1.,1.], LINESTYLE=2
    oplot, [2.,n_elements(ps)], (xavep-xaveperr)*[1.,1.], LINESTYLE=1
    oplot, [2.,n_elements(ps)], (xavep+xaveperr)*[1.,1.], LINESTYLE=1
    ; label the points
    for il=0,n_elements(ps)-1 do begin
      ytext = avep*partin/(partin+0.8)
      xyouts, il+0.5, ytext, results(this_grat(il)).dataset, $
	ALIGNMENT = 0.5
      ytext = avep*partin/(partin+0.95)
      xyouts, il+0.5, ytext, results(this_grat(il)).source, $
	ALIGNMENT = 0.5
    end
  end
end

device, /close
SPAWN, 'cp idl.ps beam_cen_periods.ps'
SPAWN, 'rm idl.ps'

; OK, now do the ANGLES
pre_print_portrait
!p.multi = [0,1,3]

; Loop over the grating types and
;   - calculate HETG average angle and opening angle
;   - make plots

ang_range = 0.1
mean_offset = -0.25
open_offset = 9.9

; LETG angle plot
xrcf_grat = where(results.grating EQ 'LEG' AND $
			STRPOS(results.phase,'XRCF') GE 0, nxrcf)
if nxrcf GE 1 then begin
  leg_angs = results(xrcf_grat).ang
  leg_angs_err = results(xrcf_grat).angerr
  xs = 0.5+indgen(n_elements(leg_angs))
  plot, xs, leg_angs, PSYM=4, $
	YRANGE = ang_range*[-1.,1.], $
	XTITLE = 'Measurement set', YTITLE='LEG Angle (deg.)', $
	TITLE = 'LETG XRCF Dispersion Angle', $
	CHARSIZE=2.5
    ; add error bars
  plot_errors, xs, xs, xs, leg_angs - leg_angs_err, $
	leg_angs + leg_angs_err
  ; label the points
  for il=0,n_elements(leg_angs)-1 do begin
    ytext = -0.75*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).dataset, $
	ALIGNMENT = 0.5
    ytext = -0.92*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).source, $
	ALIGNMENT = 0.5
  end
end

; HETG mean angle and opening angle plots
; Assume that an HEG measurement follows each MEG one...
xrcf_grat = where(results.grating EQ 'MEG' AND $
			STRPOS(results.phase,'XRCF') GE 0, nxrcf)
if nxrcf GE 1 then begin
  meg_angs = results(xrcf_grat).ang
  meg_angs_err = results(xrcf_grat).angerr
  heg_angs = results(xrcf_grat+1).ang
  heg_angs_err = results(xrcf_grat+1).angerr
  mean_angs = (meg_angs+heg_angs)/2.0
  mean_errs = SQRT(meg_angs_err^2 + heg_angs_err^2)/2.0
  xs = 0.5+indgen(n_elements(meg_angs))
  plot, xs, mean_angs, PSYM=4, $
	YRANGE = mean_offset + ang_range*[-1.,1.], $
	XTITLE = 'Measurement set', YTITLE='Mean Angle (deg.)', $
	TITLE = 'HETG XRCF Mean Dispersion Angle', $
	CHARSIZE=2.5
    ; add error bars
  plot_errors, xs, xs, xs, mean_angs - mean_errs, $
	mean_angs + mean_errs
  ; label the points
  for il=0,n_elements(leg_angs)-1 do begin
    ytext = mean_offset - 0.75*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).dataset, $
	ALIGNMENT = 0.5
    ytext = mean_offset - 0.92*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).source, $
	ALIGNMENT = 0.5
  end

  open_angs = meg_angs - heg_angs
  open_errs = SQRT(meg_angs_err^2 + heg_angs_err^2)

  ; Weighted mean and error
  ; from Bevington 2nd, (4.17) and (4.19)
  ave_open = TOTAL(open_angs/(open_errs)^2)/TOTAL(1./(open_errs)^2)
  ave_open_err = SQRT(1./TOTAL(1./(open_errs)^2))
  print, ave_open, ave_open_err

  plot, xs, open_angs, PSYM=4, $
	YRANGE = open_offset + ang_range*[-0.9,0.9], $
	XTITLE = 'Measurement set', YTITLE='Opening Angle (deg.)', $
	TITLE = 'HETG Opening Angle =' + $
		STRING(ave_open,FORMAT='(F7.4)')+' +/-'+ $
		STRING(ave_open_err,FORMAT='(F7.4)') + $
		' degrees', $
	CHARSIZE=2.5, XTICKS = 4
  ; add error bars
  plot_errors, xs, xs, xs, open_angs - open_errs, $
	open_angs + open_errs

  ; plot the XRCF average opening angle and error range
  oplot, [0.,n_elements(open_angs)], ave_open*[1.,1.], LINESTYLE=2
  oplot, [0.,n_elements(open_angs)], (ave_open - ave_open_err)*[1.,1.], LINESTYLE=1
  oplot, [0.,n_elements(open_angs)], (ave_open + ave_open_err)*[1.,1.], LINESTYLE=1

  ; label the points
  for il=0,n_elements(leg_angs)-1 do begin
    ytext = open_offset - 0.75*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).dataset, $
	ALIGNMENT = 0.5
    ytext = open_offset - 0.92*ang_range
    xyouts, il+0.5, ytext, results(xrcf_grat(il)).source, $
	ALIGNMENT = 0.5
  end

end

device, /close
SPAWN, 'cp idl.ps beam_cen_angles.ps'
SPAWN, 'rm idl.ps'

set_plot, 'X'
!p.multi = 0

RETURN
END