PRO hsi_alk_misaligned
;
; Analyze the HSI - Al-K shell-quadrant MEG +3 data sets
; to fit the mis-aligned MEG parameters precisely.
;
@df_common

; Point to where the Phase 1 data are:
phase_1_data_dir = '/nfs/spectra/d6/'

; Data sets to analyze...
dates1 = ['970104','970104','970104','970104']
runids1 = ['107463','107464','107465','107466']
agauss1 = [140., 400., 170., -100.]
bgauss1 = [0., -100., 100., 0.]

dates2 = ['970109','970109','970109','970109', $
		'970109','970109','970109','970109']
runids2 = ['108026','108027','108028','108029', $
		'108030','108031','108032','108033']
agauss2 = [-100., -100., 110., -40., $
		100., 400., 100., 0.]
bgauss2 = [0., 0., 0., -130., $
		0., -150., 0., 0.]

dates3 = ['970205','970205','970205','970205', $
		'970208','970208','970208','970208']
runids3 = ['111316','111317','111318','111319', $
		'111488','111489','111490','111491']
agauss3 = [0., -50., 130., -70., $
		70., 370., 100., 0.]
bgauss3 = [0., 0., 0., 0., $
		0., -200., 0., 0.]

; Include the ALL,ALL MEG+3's
; Al-K
dates4 = ['970109','970109']
runids4 = ['108034','108035']
agauss4 = [100., -170.]
bgauss4 = [-140., 80.]

; Mg-K : change iteration to i2 !
dates5 = ['970116']
runids5 = ['109002']
agauss5 = [120.]
bgauss5 = [-160.]

dates =  [dates1, dates2, dates3]
runids =  [runids1, runids2, runids3]
agauss = [agauss1, agauss2, agauss3]
bgauss = [bgauss1, bgauss2, bgauss3]

bin_size = 10.0 ; um
reg_size = 500.0

pre_print_portrait
!p.multi = [0,2,4]

for isq=0,n_elements(dates)-1 do begin
  hsi = xrcf_hsi_read(phase_1_data_dir+dates(isq)+ $
		'/hsi'+runids(isq)+'i0.fits',/FAC,/PBLU)

  ; rotate to grating axes
  ; MEG is 4.74 degrees
  rot_angle = -1.0*4.74*!DTOR
  grat_y = hsi.y * COS(rot_angle) - hsi.z * SIN(rot_angle)
  grat_z = hsi.y * SIN(rot_angle) + hsi.z * COS(rot_angle)

;  plot, grat_y, grat_z, PSYM=3, $
;	XRANGE=reg_size*[-1.,1.], XSTYLE=1, $
;	YRANGE=reg_size*[-1.,1.], YSTYLE=1, $
;	TITLE = dates(isq)+'/'+runids(isq)

  ; make the histogram along Z
  lin_hist, grat_z, bin_size, hist_axis, hist_counts
  sel = where(ABS(hist_axis) LE reg_size)
  hist_axis = hist_axis(sel)
  hist_counts = hist_counts(sel)

  plot_io, hist_axis, hist_counts, PSYM=10, $
	XRANGE=reg_size*[-1.,1.], XSTYLE=1, $
	YRANGE=[0.5, 6000.0], YSTYLE=1, $
	TITLE = dates(isq)+'/'+runids(isq)

  ; try fitting it


x = FLOAT(hist_axis)
y = FLOAT(hist_counts > 1)
; initial guess for parameters
; a is three params per gaussian plus the continuum level
; the gaussian params are: height, center, sigma
;
ngauss = 1
if agauss(isq) NE 0.0 then ngauss = 2
if bgauss(isq) NE 0.0 then ngauss = 3

CASE ngauss of
1: a = [MAX(y), 0.0, 15.0, 1.0]
2: a = [MAX(y), 0.0, 15.0, $
	0.15*MAX(y), agauss(isq), 15.0, 1.0]
3: a = [MAX(y), 0.0, 15.0, $
	0.15*MAX(y), agauss(isq), 15.0, $
		0.2*MAX(y), bgauss(isq), 15.0, 1.0]
ENDCASE

; rms error for each bin
rms = SQRT(y)  ; counting statistics
; constant continuum
df_continuum = 1.0 + 0.0*x  ; tricky way to get all 1.'s
df_title = 'Z Projection Fit'
; Fit it
df_fit,ngauss,x,y,a,rms, sig,yfit
; Show the fit result
df_showfit,x,y,a,rms, sig,yfit


end

device,/close
SPAWN, 'cp idl.ps ham.ps'
SPAWN, 'rm idl.ps'
set_plot, 'X'
!p.multi = 0

RETURN
END