PRO xrcf_beam_cen_werr, file_name, LPR = lpr
; This procedure reads in beam cen information from a 
; file and i) displays the data, ii) fits the data to
; get grating angles, and iii) determines changes needed in
; period and/or Rowland diameter.
;
; The current parameters are gotten from xrcf_params.pro
;
; 12/20/96 dd
; 7/2/97 dd Modified to use errors for Y and Z and output
;  angle and period errors too.
;  Input file is rdb format also.

; read the data file
bc_data = rdb_read(file_name)
beam_grat = bc_data.grating
beam_order = bc_data.order
beam_Y = DOUBLE(bc_data.Y)
beam_Z = DOUBLE(bc_data.Z)
beam_E = DOUBLE(bc_data.E)
beam_Yerr = bc_data.Yerr
beam_Zerr = bc_data.Zerr
nbeam = n_elements(bc_data)

; open an output file
out_file_name = 'beam_cen.txt'
OPENW, out_unit, out_file_name, /GET_LUN

; print out the input data
printf, out_unit, ''
printf, out_unit, '     Output from  xrcf_beam_cen_werr.pro ,   '+SYSTIME()
printf, out_unit, ''
printf, out_unit, ' Input data file : '+file_name
printf, out_unit, ''
printf, out_unit, 'Grating     Order      Prime Y' + $
		'          Z' + $
		'         Energy       Y_err(um)    Z_err(um)'
printf,out_unit,'. . . . . . . . . . . . . . . . . . . . . . . . . .'+$
	' . . . . . . . . . . . . . . . .'
for ib=0,nbeam-1 do printf, out_unit, beam_grat(ib)+'  '+ $
	STRING(beam_order(ib)) + STRING(beam_Y(ib)) + $
	STRING(beam_Z(ib)) + $
	STRING(beam_E(ib)) + $
	STRING(beam_Yerr(ib)) + STRING(beam_Zerr(ib))
printf,out_unit,'. . . . . . . . . . . . . . . . . . . . . . . . . .'+$
	' . . . . . . . . . . . . . . . .'

; if keyword LPR is set then create postscript plot
if KEYWORD_SET(LPR) then begin
  ; set plot device to postscript
  pre_print_landscape
end

; plot the raw data
plot, beam_Y, beam_Z, PSYM = 4, $
	XRANGE = [MIN(beam_Y),MAX(beam_Y)], $
	YRANGE = [MIN(beam_Z),MAX(beam_Z)], $
	TITLE = 'Beam_cen data Plot, from '+file_name, $
	XTITLE = 'Prime Y (microns)', YTITLE = 'Z (microns)'

; Get information and parameters...
xrcf_params, NAM = names
xrcf_params, ANG = angles
xrcf_params, PER = periods
xrcf_params, ROW = rowland
xrcf_params, HC = hc

printf, out_unit, ''
printf, out_unit, ' Calculate Angles and Periods to best-fit the data:'
printf, out_unit, '  (keeping the Rowland distance at'+$
	STRING(rowland)+' )'
printf, out_unit, ''

; Arrays for results
calc_angles = DBLARR(n_elements(names)) - 100.0
calc_angles_err = DBLARR(n_elements(names))
calc_rowland = DBLARR(n_elements(names))
calc_period = DOUBLE(periods)  ; initial guess
calc_period_err = DOUBLE(periods) 

; go through the grating types
for ig=0,n_elements(names)-1 do begin
  these_grats = where(STRPOS(beam_grat, names(ig)) GE 0, nthese)
  if nthese GE 2 then begin
    ; calculate the angles
    ; find the best fit line for Y and Z
    ; Use LINFIT routine
    result = LINFIT( beam_Y(these_grats), $
		beam_Z(these_grats), $
		SDEV=beam_Zerr(these_grats), SIGMA = sigmaA)
    Z_fit = result(0) + result(1)*beam_Y(these_grats)

    oplot, beam_Y(these_grats), Z_fit, LINESTYLE=1
    calc_angles(ig) = ATAN(result(1))/!DTOR
    ; Angle error
    calc_angles_err(ig) = ATAN(result(1)+sigmaA(1))/!DTOR - calc_angles(ig)
    ; print the measured and current angles
    ang_fmt = '(F8.4)'
    printf, out_unit, '  '+names(ig)+' angle = ' + $
		STRING(calc_angles(ig),FORMAT=ang_fmt)+ $
		' +/-'+STRING(calc_angles_err(ig),FORMAT=ang_fmt) + $
		'  [ currently = '+STRING(angles(ig),FORMAT=ang_fmt)+']'

    ; Adjust the period to have the Rowland value be the value in 'params
    change = 1.10
    while ABS(change-1.0) GE 1.E-8 do begin 
      ; calculate the Rowland distance if the period is accurate
      ; Rowland is the slope between big_mess and Yprime
      big_mess = 0.5*SIN(2.0*ASIN( $
  	beam_order(these_grats)*(hc/beam_E(these_grats)) / $
							calc_period(ig)) )
      big_mess = big_mess/SQRT(1.0+(TAN(!DTOR*calc_angles(ig)))^2)

      ; Use LINFIT routine
      result = LINFIT( big_mess, $
		1.E-3*beam_Y(these_grats), $
		SDEV=1.E-3*beam_Yerr(these_grats), SIGMA = sigmaA)
      Z_fit = result(0) + result(1)*big_mess

      calc_rowland(ig) = result(1)

      ; modify the period to get Rowland back to set value
      change = rowland/calc_rowland(ig)   ; approximation
      this_period = calc_period(ig)
      calc_period(ig) = change * this_period

      ; error estimate for period
      err_rowland = result(1) + sigmaA(1)
      err_change = rowland/err_rowland  ; approximation
      calc_period_err(ig) = calc_period(ig) - err_change * this_period
    end

; . . . . . . 
; this calculation for paranoid individuals! - skip it
if 1 EQ 0 then begin
    ; calculate the Rowland distance if the period is accurate
    ; (double check that the rowland distance is back where
    ; it should be
    ; Rowland is the slope between big_mess and Yprime
    big_mess = 0.5*SIN(2.0*ASIN( $
  	beam_order(these_grats)*(hc/beam_E(these_grats)) / $
							calc_period(ig)) )
    big_mess = big_mess/SQRT(1.0+(TAN(!DTOR*calc_angles(ig)))^2)
    result = LINFIT(big_mess , 1.E-3*beam_Y(these_grats), $
		SDEV=1.E-3*beam_Yerr(these_grats), SIGMA = sigmaA)
    calc_rowland(ig) = result(1)
    ; print the measured and current Rowland distance
    printf, out_unit, '  '+names(ig)+' rowland = '+ $
		STRING(calc_rowland(ig))+ $
		'  [ currently = '+STRING(rowland)+']'
end
; . . . . . . . 

    ; print the measured and current Period
    per_fmt = '(F8.2)'
    printf, out_unit, '  '+names(ig)+' period = '+ $
		STRING(calc_period(ig),FORMAT=per_fmt) + $
		' +/-'+STRING(calc_period_err(ig),FORMAT='(F5.2)') + $
		'  [ currently = '+STRING(periods(ig),FORMAT=perfmt)+']'

    printf, out_unit, ''
    printf, out_unit, '     Errors between measured '+ $
		'and calculated Y,Z (microns):'
    mlamp = beam_order(these_grats)* $
		(hc/beam_E(these_grats))/calc_period(ig)
    xrcf_disp_calc, mlamp, rowland, calc_angles(ig), 0.0,0.0,0.0, $
		Xpred, Ypred, Zpred
    ; convert predictions to um and shift to have same averages
    Ypred = 1000.0*Ypred + $
		TOTAL(beam_Y(these_grats))/n_elements(these_grats)
    Zpred = 1000.0*Zpred + $
		TOTAL(beam_Z(these_grats))/n_elements(these_grats)
    for it = 0,n_elements(these_grats)-1 do begin
      printf, out_unit, '    Y: ',beam_Y(these_grats(it)), '   -' , $
		Ypred(it) , ' = ', $
		STRING(beam_Y(these_grats(it)) - Ypred(it),FORMAT='(F6.1)')
    end
    for it = 0,n_elements(these_grats)-1 do begin
      printf, out_unit, '    Z: ',beam_Z(these_grats(it)), '   -' , $
		Zpred(it) , ' = ', $
		STRING(beam_Z(these_grats(it)) - Zpred(it),FORMAT='(F6.1)')

    end    
    printf, out_unit, '        -   -   -   -   -'
    printf, out_unit, ''
  end
end
printf,out_unit,'. . . . . . . . . . . . . . . . . . . . . . . . . .'+$
	' . . . . . . . . . .'

; determine the common period change

; print the recommended new rowland value and the recommended
; new period values

; print instructions for updating the constants


; close the output file
CLOSE, out_unit
FREE_LUN, out_unit

; if LPR then print the file and plot
if KEYWORD_SET(LPR) then begin
  ; set plot device to postscript
  post_print
  SPAWN, 'enscript -fCourier9 '+out_file_name
end


RETURN
END