pro lindiag,coeffs,medout,fracout,LOWER=LOWER, $
            minc=MINC, maxc=MAXC, deltac=DC, $
            title=title, outfile=outfile, silent=silent

;+
;
; diagnostic routine for linearity correction coefficients
; will generate a bunch of useful plots and tables
;
; INPUTS
;	coeffs	coefficients of linearity correction
;		(should be of dimension N for an Nth order fit)
;
; KEYWORD PARAMETERS
;	lower	lower bound for fitting
;       minc    minimum counts for diagnostic plots
;       maxc    max counts for plots
;       dc      step size of counts for plots
;       outfile file for Postscript plotting
;
; OUTPUT
;	medout	median % non-linear versus counts
;	fracout	% of pixels non-linear for a given series of counts
;
; defined so that when at the lower limit the correction = 0.0
; should have two consequtively numbered IDL windows open,
;   with the active window being the lower numbered one when starting
;
; 2/22/96 MCL
;
;-

LINCONST = 1e6  ; coeff array assumed to be multiplied by this number

if not(keyword_set(LOWER)) then LOWER = 0.0
if not(keyword_set(MINC)) then MINC = 0.0
if not(keyword_set(MAXC)) then MAXC = 30000.
if not(keyword_set(DC)) then DC = 3000.
if n_params() lt 1 then begin
	print,'pro lindiag,coeffs,medout,fracout,[LOWER='+strc(lower)+'],'
        print,'            [minc=0.0],[maxc=', strc(maxc), '],[deltac=', strc(dc),'],'
        print,'            [title=],[outfile=],[silent=]'
	retall
endif
;tt =  'Avg lin corr'
tt =  ''
if keyword_set(title) then tt = tt + ' '+title


sz = size(coeffs)
if sz(0) eq 2 then ORDER = 1 else ORDER = sz(3)
npix =  sz(1)*float(sz(2))

cs = 1.2   ; charsize
ss = 1.0   ; symsize

; filled circle plotting symbol
a = findgen(16)*(!pi*2/16.)
usersym,cos(a),sin(a),/fill

message, 'did you remember to set LOWER to correct value??', /info


;
; (1) median % non-linear versus counts
;
nmed = (MAXC - MINC) / DC + 1
medout = fltarr(nmed, 3)
cc = MINC > 1.0
for i = 0, nmed-1 do begin
    if (cc le LOWER) then begin
        medout(i, 0) = cc
        medout(i, 1) = 1.0
        medout(i, 2) = 0.0
    endif else begin
	corr = dblarr(sz(1),sz(2))
        for k = 1, ORDER do $
          corr = corr + coeffs(*, *, k-1)/LINCONST * (cc-LOWER)^(k)
;        corr = - (corr)
        corr = 1D/(1+corr)
        medout(i, 0) = cc
        medout(i, 1) = median(corr)
        iterstat, corr, out, /silent
        medout(i, 2) = out(3) 
;         medout(i, 2) =  stdev(corr(where(corr ne 0))) * 100.
    endelse
  cc = cc + DC
endfor

; print it out
print
print, '# median linearity correction factors'
print, '#'
print, '# LOWER = ', strc(LOWER)
;spawn, 'date', date
;print, '#', date
print, '#'
print, '#     counts	  lin corr   sigma(lin corr)'
print, '#'
printarray,  medout
print
print


; plot it
yr =  [min(medout(*, 1) -medout(*, 2)), max(medout(*, 1)+medout(*, 2))]
plot, medout(*, 0), medout(*, 1), ps = -8, /yno, $
  charsize = cs, symsize = ss, yr = yr, ys = 2, $
  xtitle = 'counts (bias-subtracted)', $
  ytitle = 'median linearity correction factor', $
  title = tt
errplot, medout(*, 0), medout(*, 1)-medout(*, 2), medout(*, 1)+medout(*, 2)


; dump graph to PS file
if keyword_set(outfile) then begin

    ll =  strlen(outfile)
    if (strmid(outfile, ll-3, 3) eq '.ps') then outfile = $
        strmid(outfile, 0, ll-3) 

    ps_open, outfile, thick = 4

    plot, medout(*, 0), medout(*, 1), ps = -8, /yno, $
      charsize = cs, symsize = ss, yr = yr, ys = 2, $
      xtitle = 'counts (bias-subtracted)', $
      ytitle = 'median linearity correction factor', $
      title = tt
    errplot, medout(*, 0), medout(*, 1)-medout(*, 2), medout(*, 1)+medout(*, 2)

    datestamp

    ps_close

endif


;
; (2) fractional non-linearity diagnostic
; outer loop is the counts loop, inner loop is the percentage loop
;
ncf = (MAXC - MINC + DC) / DC + 1
MINP = 1.0                      ; min %age - do not set < 1!
MAXP = 8.0                      ; max %age
STEPP = 2.0                     ; geometrical step for %age
np = round(alog10(MAXP/MINP)/alog10(STEPP) + 1)
fracout = fltarr(np+1, ncf)

; print table heading
pp =  MINP
print, format = '($,A11," | ")','# counts   '
print, format = '(A30)', 'fraction non-linear'
print, format = '($,A11," | ")', '# (no bias)'
for j = 0, np-1 do begin
    print, format = '($,"<",F6.3,"%  ")', strc(pp)
    pp = pp * STEPP
endfor
print, format = '(F6.3,"+%")', strc(pp/STEPP)
print, '# ------------------------------------------------------------'

cc = MINC+DC
for i = 0, ncf-1 do begin

;   calculate linearity correction for all pixels at given count level
    corr = dblarr(sz(1), sz(2))
    if (cc gt LOWER) then begin
        for k = 1, ORDER do $
          corr = corr + coeffs(*, *, k-1)/1e6 * (cc-LOWER)^(k)
        corr = 1D/(1+corr) - 1D
    endif

;   now find correction percentages
    pp =  MINP
    ppold =  0.0
    corr =  abs(corr * 100.)
    w = where(corr eq 0, nbad)
    for j = 0, np-1 do begin
        w =  where((corr lt pp) and (corr ge ppold), nnew)
        fracout(j, i) = (nnew)/npix
        ppold =  pp
        pp =  pp * STEPP
    endfor
    w =  where(corr ge ppold, nnew)
    fracout(j, i) =  nnew/npix

    ff =  '($,A10,"     ",'+strc(fix(np+1))+'(F6.3,"    "))'
    print, format = ff, strc(cc), fracout(*, i)
    print

    cc = cc + DC

endfor
print, '# number of unfittable pixels = ', strc(fix(nbad))


end