; load CTIO linearity data (taken with reference frames),
; make a plot to examine the data, and then apply the correction to
; the counts in the linearity frames based on the ref frames
;
; resulting arrays of note (in memory) are:
;    lin_nb    bias-subtracted, ref-corrected lin frames
;    tint      array of int. times for lin frames
;    
;    linlist   indicies identifying lin frames
;    reflist   indicies identifying ref frames
;    ut        UT time from FITS header for lin & ref frames
; 
;    medlin    (corrected) median counts of lin frames
;    medref    median counts of ref frames
;    tref      int time of reference frame
;
; 04/02/98 MCL


; start & end frame #s for linearity & dark data
fdir = '../fits-11nov97/'
LIN0 = 4020
LIN0 = 4051
LIN0 = 4082
LIN0 = 4113
LIN1 = LIN0+30
DARK0 = 4144
DARK1 = 4174


; FITS keywords for UT time stamp & int. time
uthead = 'UT'
tinthead = 'INT_S'


; get linearity frames and UT time from FITS header
; we split the ref & linearity frames into separate arrays
;   called 'linimgs' and 'refimgs'
ans1 = 1
if n_elements(lin) gt 0 then $
  ans1 = getyn('linearity array already defined - reload?', 1)
if (ans1) then begin
  loadims, fd=fdir, /ctio, LIN0, LIN1, lin, head0, linfiles
  szlin = size(lin)
  ut = fltarr(szlin(3))
  tint = fltarr(szlin(3))
  for i = 0, szlin(3)-1 do begin
      im = readfits(linfiles(i), head, /silent)
      ut(i) = timeconv(sxpar(head, uthead))
      tint(i) = sxpar(head, tinthead)
      if (i eq 0) then print, 'start UT = ', sxpar(head, uthead)
      if (i eq szlin(3)-1) then print, 'end UT = ', sxpar(head, uthead)
  endfor    
  ut = ut-ut(0)

  ; define list of lin & ref frames
  ; assumes ref frame taken in between each frames, and also
  ; one at the beginning
  nref = 1 + (szlin(3)-1)/2
  nlin = szlin(3) - nref
  reflist = indgen(nref) * 2
  linlist = indgen(nlin) * 2 + 1

  linimgs = lin(*, *, linlist)
  refimgs = lin(*, *, reflist)
  tref = tint(0)
  tint = tint(linlist)

endif




; get biases
ans2 = 1
if n_elements(lindarks) gt 0 then $
  ans2 = getyn('dark array already defined - reload?', 0)
if (ans2) then begin
  loadims, fd=fdir, /ctio, DARK0, DARK1, lindarks, head0, darkfiles
  szdark = size(lindarks)
  utd = fltarr(szdark(3))
  for i = 0, szdark(3)-1 do begin
      im = readfits(darkfiles(i), head, /silent)
      utd(i) = timeconv(sxpar(head, uthead))
  endfor    
  utd = utd-utd(0)
endif
ndarks = (size(lindarks))(3)


; sanity check
if (ndarks ne szlin(3)) then begin
    print, '# of dark & linearity frames do not agree!'
    stop
endif


; do bias subtraction
ans3 = 0
if (ans1 or ans2) then begin
    if getyn('new lin and/or darks: re-do bias sub?', 1) then begin
      lin_nb = linimgs - lindarks(*, *, linlist)
      ref_nb = refimgs - lindarks(*, *, reflist)
      ans3 = 1
  endif
endif


; plot medians of ref & lin frames, normalized to 1 sec int. time
cs = 1.8
medref = meds(ref_nb)
medlin = meds(lin_nb)
plot, ps = 5, /ynoz, ut(reflist), medref/tref, $
  yr = [min([medlin/tint, medref/tref]), max([medlin/tint, medref/tref])], $
  xr = [-10, max(ut)], xs = 1, $
;  yr = [1780, 1880], ys = 1, $
;  xr = [-10, 720], xs = 1, $
  xtit = 'elapsed real time (sec)', ytit = 'count rate', $
  tit = 'lin frames: '+strc(LIN0)+'-'+strc(LIN1), chars = cs
oplot, ps = 1, ut(linlist), medlin/tint


; do robust linear fit to ref counts versus time and overplot
refcc = linfit(ut(reflist), medref)
oplot, line = 5, ut(reflist), refcc(0)+refcc(1)*ut(reflist)


; calculate correction factor for lin data based on ref data drift
lincc = 1.-refcc(1)*ut(linlist)/refcc(0)
if (ans3) then begin
    if getyn('you re-did bias sub. - apply ref drift corrections?') then begin
        for i = 0, nlin-1 do $
          lin_nb(*, *, i) = lin_nb(*, *, i) * lincc(i)
        medlin = medlin * lincc
        print, 'max ref drift correction = ', strc(100.*max(1.-lincc)), ' %'
    endif
endif


; overplot corrected linearity data
oplot, ps = 2, ut(linlist), medlin/tint
legend, ['ref frame ('+string(tref, form = '(f3.1)')+'s) ', $
         'linearity frame ', 'linearity (corr) '], $
  psym = [5, 1, 2], syms = [2, 2, 2], /top, /right, chars = cs


; plot medians for dark frames
if getyn('plot medians of bias frames?', 0) then begin
    plot, ps = 5, /ynoz, utd(reflist), meddarks(reflist), $
      yr = [min(meddarks), max(meddarks)], $
      xr = [-10, 720], xs = 1, $
      xtit = 'elapsed real time (sec)', ytit = 'median bias counts', $
      tit = 'biases for linearity: '+strc(DARK0)+'-'+strc(DARK1), chars = cs
    oplot, ps = 1, utd(linlist), meddarks(linlist)
endif


; write medians vs time to file
if getyn('write medians vs time to file?') then begin
    file = ''
    read, 'enter file name: ', file
    openw, unit0, file, /get_lun
    spawn, 'date', date
    printf, unit0, '# LINREFPREP.IDL: median counts vs time data'
    printf, unit0, '#'
    printf, unit0, '# linearity frames:'
    printf, unit0, '#   file directory = "', fdir, '"'
    printf, unit0, '#   1st file number = ', strc(LIN0)
    printf, unit0, '#   last file number = ', strc(LIN1)
    printf, unit0, '#'
    printf, unit0, '#   1st bias number = ', strc(DARK0)
    printf, unit0, '#   last bias number = ', strc(DARK1)
    printf, unit0, '#'
    printf, unit0, '# median counts have been corrected for ref frame drift'
    printf, unit0, '# (factor is listed in the last column)'
    printf, unit0, '#'
    printf, unit0, '#', date
    printf, unit0, '#'
    printf, unit0, '#   int. time	median counts  med. ct rate  drift corr'
    printf, unit0
    free_lun, unit0
    printarray, [[tint],[medlin],[medlin/tint],[lincc]], outfile = file, /silent
endif

stop

; for CTIO data, we exlude the 1st lin frame
; since it has strange behavior
print, '**  NOTE: excluding 1st lin frame from lin data  **'
lin_nb = lin_nb(*, *, 1:*)
tint = tint(1:*)

end