pro sag_ltc_plots, out0, infil_ref, v, tit20, skip_right=skip_right, fmax=fmax, $
	outstr=outstr, progver=progver, info=info
;
;
;HISTORY:
;V1.0	Written 19-Nov-99 by M.Morrison taking sxi_ltc.pro as start
;V1.01	22-Mar-02 (MDM) - Modified to plot the excluded points as dots and to
;			  spread out the printing for the matrix results
;
progver = 'SAG_LTC  Ver 1.01
;
outstr = {gain: 0., avg_gain: 0., dev_gain: 0., range_gain: 0.}
;
out = out0
if (keyword_set(skip_right)) then begin
    nx = n_elements(out(*,0))
    out = out(0:nx/2-2,*,*)
    tit2 = 'Skipping the right half of CCD'
end
if (keyword_set(tit20)) then tit2 = tit20	;copy it to avoid passing it out
;
if (keyword_set(fmax)) then begin
    ss = where(out.avg le fmax, nss)
    if (nss eq 0) then begin
	print, 'No data with signal below: ', fmax
	print, 'Returning...
	return
    end
    out2 = out(ss)
end else begin
    out2 = out
end
;
ss33 = avoid_spikes(out2.noise)
ss44 = where(out2.stat eq 1, nss44)	;MDM 22-Mar-02
if (nss44 ne 0) then ss33 = ss33( rem_elem(ss33, ss44) )
out2 = out2(ss33)
;
if (n_elements(v) eq 0) then v = 5	;1		;255
bits, v, barr
gain = out.avg / out.noise	;not quite right -- camera read noise should be subtracted
gain2 = out2.avg / out2.noise	;not quite right -- camera read noise should be subtracted
;
if (!d.name ne 'PS') then wdef, !d.window>0, 600, 800, /already
;
if (barr(0)) then begin		;standard signal vs noise^2 plot
    !p.multi = 0
    ;
    nnn = n_elements(out(0,*))
    if (nnn ne 1) then begin
	tmp = out(*,(nnn-3)>0:*,*).noise
;;	yr = [0, median(tmp)*1.15]
	ss33 = avoid_spikes(out.noise)
	yr = [0, max(out(ss33).noise)]
    end
    ;
    plot, out.avg, out.noise, psym=2, tit = infil_ref, $
		xtit='Signal (DN)', ytit = 'Noise^2 (DN^2)', yr=yr
    ;
    coeff=poly_fit(out2.avg, out2.noise, 1, yfit)
    oplot, out2.avg, yfit
    g = 1/coeff(1)
    outstr.gain = g
    plottimes, .1, .80, 'Gain = ' + string(g, format='(f10.3)') + ' e/DN'
    plottimes, .1, .77, 'Y-Intercept = ' + string(coeff(0), format='(f10.3)') + ' DN^2'
    plottimes, .1, .74, 'X-Intercept = ' + string(-coeff(0)/coeff(1), format='(f10.3)') + ' DN'
    ;
    if (keyword_set(tit2)) then plottimes, .1, .9, tit2
    plottime, 0, 0, progver
    plottime
end
;
if (barr(1)) then begin
    !p.multi = [0,1,2]
    plot, out2.avg, gain2, psym=2, tit=infil_ref, /ynoz, $
	xtit = 'Signal (DN)', ytit = 'Gain (e/DN)'
    plottimes, .1, .85, 'Average Gain = ' + string(total(gain)/n_elements(gain), format='(f8.3)') + ' e/DN'

    plot_hist, gain2, bin=.1, $
	xtit = 'Gain (e/DN)'

    if (keyword_set(tit2)) then plottimes, .1, .9, tit2
    plottime, 0, 0, progver
    plottime
end
;
if (barr(2)) then begin		;multiple fix signal vs noise^2 plot
    !p.multi = 0
    ;
    nnn = n_elements(out(0,*))
    if (nnn ne 1) then begin
	tmp = out(*,(nnn-3)>0:*,*).noise
;;	yr = [0, median(tmp)*1.15]
;;	ss33 = avoid_spikes(out.noise)
	yr = [0, max(out(ss33).noise)*1.5]
    end
    ;
    ;plot, out.avg, out.noise, psym=2, tit = infil_ref, $
    ;		xtit='Signal (DN)', ytit = 'Noise^2 (DN^2)', yr=yr
    plot, out.avg, out.noise, psym=3, tit = infil_ref, $
		xtit='Signal (DN)', ytit = 'Noise^2 (DN^2)', yr=yr
    oplot, out(ss33).avg, out(ss33).noise, psym=2
    ;
    x0 = .1
    y0 = .9
    nsteps = 5
    if (keyword_set(info)) then plottimes, .03, 0.98, info.tit
    plottimes, x0+nsteps/2.*.08, y0+.03*2,         ' Fit Lower Threshold
    plottimes, x0-.05,            y0-nsteps*.03, ' Fit Upper Threshold', 1, 90
    nsteps = 5
    fmax2 = max(out2.avg)<fmax
    iimin = min(out2.avg)
    range = fmax2 - iimin
    for ix=0,nsteps-1 do begin
	for iy=0,nsteps-1 do begin
	    imin = ix*range/nsteps + iimin
	    imax = (iy+1)*range/nsteps + iimin
	    if (iy eq 0) then plottimes, x0+ix*.1, y0+.03,    string(long(imin), format='(i10)')
	    if (ix eq 0) then plottimes, x0-   .1, y0-iy*.03, string(long(imax), format='(i10)')
	    ss = where((out2.avg ge imin) and (out2.avg le imax), nss)
	    if (nss ge 6) then begin
		out4 = out2(ss)
		coeff=poly_fit(out4.avg, out4.noise, 1, yfit)
		oplot, out4.avg, yfit
		plottimes, x0+ix*.1, y0-iy*.03, string(1/coeff(1), format='(f10.2)')
		if (n_elements(mgain) eq 0) then mgain = 1/coeff(1) $
					else mgain = [mgain, 1/coeff(1)]
	    end
	end
    end
    if (keyword_set(mgain)) then begin
	idev = stdev(mgain, iavg)
	strout = string('Avg=', iavg, ' Dev=', idev, ' (', idev/iavg*100, ' %)  Range=', $
			max(mgain)-min(mgain), format='(4(a,f10.3))')
	plottimes, x0, y0-nsteps*.03, strout
	outstr.avg_gain	= iavg
	outstr.dev_gain = idev
	outstr.range_gain = max(mgain)-min(mgain)
    end
    if (keyword_set(tit2)) then plottimes, .1, .9, tit2
    plottime, 0, 0, progver
    plottime
end
;
end
;-------------------------------------------------------
pro sag_ltc_1pair, infil1, infil2, sel_struct, out, ipair, npair
;
;
;
print, 'Now processing: ', infil1, ' and ', infil2
;
img1 = rfits(infil1, h=h1) & img1 = unsign(img1)
img2 = rfits(infil2, h=h2) & img2 = unsign(img2)
;
if (data_type(sel_struct) ne 8) then begin
    sel_struct = {nx_block: 30, ny_block: 30, $
		x0: 10, x1: data_chk(img1, /xsize)-10, $
		y0: 10, y1: data_chk(img1, /ysize)-10, $
		b0: 3, b1: 6, $
		offset: 'NONE'}
end

nx_block 	= sel_struct.nx_block
ny_block 	= sel_struct.ny_block
x0		= sel_struct.x0
x1		= sel_struct.x1
y0		= sel_struct.y0
y1		= sel_struct.y1
b0		= sel_struct.b0
b1		= sel_struct.b1    
;
tmp = img1(b0:b1, y0:y1)
bls_avg0 = total(tmp)/n_elements(tmp)
tmp = img1(x0:x1, y0:y1)
img1_avg = total(tmp)/n_elements(tmp)
tmp = img2(x0:x1, y0:y1)
img2_avg = total(tmp)/n_elements(tmp)
print, 'BLS_AVG = ', bls_avg0, ' IMG1_AVG = ', img1_avg, ' IMG2_AVG = ', img2_avg
bls_avg = 0
case sel_struct.offset of
    'USE_BLS': bls_avg = bls_avg0
    'NONE':
endcase

nx_steps = fix((x1-x0)/nx_block)
ny_steps = fix((y1-y0)/ny_block)
;
if (n_elements(out) eq 0) then begin
    out0 = {avg: 0., dev: 0., noise: 0., min: 0., max: 0., $
	avg1: 0., avg2: 0., dev1: 0., dev2: 0., $
	x00: 0, y00: 0, stat: 0}
    out = replicate(out0, nx_steps, ny_steps, npair)
end
;
for ix=0,nx_steps-1 do begin
    for iy=0,ny_steps-1 do begin
	x00 = x0 + ix*nx_block
	y00 = y0 + iy*ny_block
	img1b = img1(x00:x00+nx_block-1, y00:y00+ny_block-1)
	img2b = img2(x00:x00+nx_block-1, y00:y00+ny_block-1)

	out(ix, iy, ipair).dev1 = stdev(img1b, avg_temp1)
	out(ix, iy, ipair).avg1 = avg_temp1
	out(ix, iy, ipair).dev2 = stdev(img2b, avg_temp2)
	out(ix, iy, ipair).avg2 = avg_temp2

	out(ix, iy, ipair).x00 = x00
	out(ix, iy, ipair).y00 = y00

	diff = img1b - img2b
	dev = stdev(diff, avg_temp)
	out(ix, iy, ipair).avg = (avg_temp1 + avg_temp2) / 2. - bls_avg
	out(ix, iy, ipair).dev = dev
	out(ix, iy, ipair).noise = dev^2 / 2
	out(ix, iy, ipair).min = min(diff)
	out(ix, iy, ipair).max = max(diff)
	;
	if (min(img1b) eq 0) or (min(img2b) eq 0) then out(ix, iy, ipair).stat=1
    end
end
;
if (n_elements(out2) ne 0) then out2 = [[out], [out2]] else out2 = out
end
;--------------------------------------------------------
pro sag_ltc, infil_arr, nexp, npe, fmax, sel_struct, qstop=qstop, $
		vmask=vmask, out=out, info=info, $
		logfil=logfil, hc=hc, gif=gif
;
;
delvarx, out
if (n_elements(npe) eq 0) then npe = 3		;number of images per exposure
if (n_elements(nexp) eq 0) then nexp = 8	;number of different exposure durations
if (n_elements(fmax) eq 0) then fmax = 45000.

;
npair = nexp	;later to 3*nexp

for ie=0,nexp-1 do begin
    infil1 = infil_arr(ie*npe)
    infil2 = infil_arr(ie*npe+1)
    sag_ltc_1pair, infil1, infil2, sel_struct, out, ie, npair
    if (ie eq 0) then infil_ref = infil1
end
gain = out.avg / out.noise

v = out.noise
v = reform(v, n_elements(v))
junk = filter_1d(v,ix_good=ix_good, smth_wid = 100<(.2*n_elements(v)))
bad = intarr(n_elements(v))+1
bad(ix_good) = 0
ss = where(bad, nss)
if (nss ne 0) then out(ss) = out((ss-1)>0)
;
sag_ltc_plots, out, infil_ref, vmask, tit2, fmax=fmax, $
	outstr=outstr, progver=progver, info=info
;
;---- Optionally write results to log file
;
if (keyword_set(logfil)) then begin
    break_file, infil_ref, dsk_log, dir00, filnam00
    outlog = string(filnam00, outstr.gain, outstr.avg_gain, outstr.dev_gain, outstr.range_gain, $
	progver, format='(a, 4(1x, f10.3), 3x, a)')
    print, outlog
    file_append, logfil, outlog, /uniq
end
;
;---- Optionally write results to GIF file
;
if (keyword_set(gif)) then begin
    save_dev = !d.name
    set_plot,'z'
    sag_ltc_plots, out, infil_ref, vmask, tit2, fmax=fmax, $
	outstr=outstr, progver=progver, info=info
    zbuff2file, gif
    set_plot, save_dev
end
;
;---- Optionally write results to GIF file
;
if (keyword_set(hc)) then begin
    save_dev = !d.name
    setps,/port
    if (data_type(hc) eq 7) then psfil = hc else psfil = 'idl.ps'
    device, file=psfil
    sag_ltc_plots, out, infil_ref, vmask, tit2, fmax=fmax, $
	outstr=outstr, progver=progver, info=info
    device, /close
    sprint, psfil
    set_plot, save_dev
end
;

if (keyword_set(qstop)) then stop
end