; script to do linearity correction and then do
; stats and diagnostics
;
; to use this, need to have two variables defined:
;    lin_nb   bias-subtracted linearity frames
;    tint     array of integration times
;
; it makes the plots of the median corrections if 
; any obvious (and populous) set of bad pixels is set
; to zero first, e.g. for NIRC,
;       lin_nb(0:3,*.*) = 0.0
;
; 2/23/96 MCL

;----------------------------------------
; please adjust these:

run = 'CTIO Nov. 97'    ; label for plots
set = 'set 4'           ; label for plots
out = 'ctio-nov97-lin4' ; prefix for output files
out = 'test'
;----------------------------------------


if (!d.name eq 'X') then win, 0, sz = 450
LOWER = 0


print, 'run = "',run,'", set = "',set,'", out = "',out, '"'
if not(getyn('are these ok?')) then stop


print, '# linearity correction stats for linearity data from ', run,', ', set 
print, '#'
spawn, 'date', date
print, '# ', date
print, '#'

; calculate the linearity correction
lincalc, lin_nb, tint, lincoeff, maxmap, cfitmed, record, cfit, $
	tt = run+', '+set , lower = LOWER
print

; run diagnostics
wset,0
lindiag, lincoeff, medout, fracout, title = run+', '+set, $
  outfile = out+'_corr.ps', lower = LOWER

if getyn('do you want to save results & plots?') then begin

    ; make PS file of the pixel by pixel results
    ps_open, out+'_pixels', /port, thick = 3
    fitplot, lin_nb, tint, cfit
    ps_close

    ; make another PS file
    ps_open, out+'_fit', thick = 3
    fitplot, lin_nb, tint, cfit, /med, title = run+', '+set
    ps_close

    writefits, out+'_lincoeff.fits', lincoeff
    writefits, out+'_maxmap.fits', maxmap

endif


; output stats
if (size(lincoeff))(0) eq 2 then begin
    print, '# avg lin coeff stats (from ITERSTAT):'
    print, form = '($,"#")'
    iterstat, lincoeff
endif
print, '# avg of maxmap = ', strc(median(maxmap))

end