function rowmed, data, mask

rowm = make_array(/float, (dim(data))[1], value=0)

for i=0, (dim(data))[1]-1 do begin
	row = data[*, i]
	rmask = mask[*, i]
	good = where(rmask)
	if dim(good) gt 0 then rowm[i] = median(row[good])
endfor

return, rowm

end


pro gain

!except = 0
!quiet = 1
dw
device, retain=2

n1 = 20
i1 = 955
n2 = 11
i2 = 975
n3 = 5
i3 = 986

data = make_array(/float, 512, 512, n1 + n2 + n3, value=0)
for i = 0, n1 - 1 do begin
	fname = 'images/data0' + string(i + i1, format='(i3)') + '.a.fits'
	j = i
	print, j, ' Reading ', fname
	data[*, *, j] = float(readfits(fname)) > 0
endfor
for i = 0, n2 - 1 do begin
	fname = 'images/data0' + string(i + i2, format='(i3)') + '.a.fits'
	j = i + n1
	print, j, ' Reading ', fname
	data[*, *, j] = float(readfits(fname)) > 0
endfor
for i = 0, n3 - 1 do begin
	fname = 'images/data0' + string(i + i3, format='(i3)') + '.a.fits'
	j = i + n1 + n2
	print, j, ' Reading ', fname
	data[*, *, j] = float(readfits(fname)) > 0
endfor

x = make_array(/float, n1 + n2 + n3, value=0)
x[0:n1 - 1] = 0.25 * 4
x[n1:n1 + n2 - 1] = 0.25 * 40
x[n1 + n2:n1 + n2 + n3 - 1] = 0.25 * 80
print, 'Exposure = ', x

a = make_array(/float, 512, 512, value=0)
b = make_array(/float, 512, 512, value=0)
c2 = make_array(/float, 512, 512, value=0)

print, 'Fitting...'
for i = 0, 511 do begin
	for j = 0, 511 do begin
		y = data[i, j, *]
		z = linfit(x, y, chisq=c)
		a[i, j] = z[0]
		b[i, j] = z[1]
		c2[i, j] = c
	endfor
endfor

window, 0, xs=512, ys=512, title='Intercept'
tvscl, a
print, 'max(a), min(a), mean(a), median(a), stddev(a)'
print, max(a), min(a), mean(a), median(a), stddev(a)
window, 1, xs=512, ys=512, title='Slope'
tvscl, b
print, 'max(b), min(b), mean(b), median(b), stddev(b)'
print, max(b), min(b), mean(b), median(b), stddev(b)
window, 2, xs=512, ys=512, title='Chi^2'
tvscl, c2
print, 'max(c2), min(c2), mean(c2), median(c2), stddev(c2)'
print, max(c2), min(c2), mean(c2), median(c2), stddev(c2)

stop, 'a, b, & chi^2 displayed. .c to continue.'

mask1 = abs(a - mean(a)) lt 2.5 * stddev(a)
mask2 = abs(b - mean(b)) lt 2.5 * stddev(b)
mask3 = abs(c2 - mean(c2)) lt 2.5 * stddev(c2)
mask = mask1 * mask2 * mask3

a = a * mask
b = b * mask
c2 = c2 * mask

window, 3, xs=512, ys=512, title='Mask'
tvscl, mask
wset, 0
tvscl, a
print, 'max(a), min(a), mean(a), median(a), stddev(a)'
print, max(a), min(a), mean(a), median(a), stddev(a)
wset, 1
tvscl, b
print, 'max(b), min(b), mean(b), median(b), stddev(b)'
print, max(b), min(b), mean(b), median(b), stddev(b)
wset, 2
tvscl, c2
print, 'max(c2), min(c2), mean(c2), median(c2), stddev(c2)'
print, max(c2), min(c2), mean(c2), median(c2), stddev(c2)

writefits, 'images/a.fits', a
writefits, 'images/b.fits', b
writefits, 'images/c2.fits', c2
writefits, 'images/mask2.fits', mask
stop, 'Masked a, b, & chi^2 displayed.  .c to continue.'

print, 'mask[200,100:109]'
print, mask[200,100:109]
window, 4, title='Counts vs. exposure'
plot, x, data[200,100,*], psym=7, yrange=[0,3e4]
plot, x, a[200, 100] + b[200, 100] * x, yrange=[0,3e4], /noerase
for i = 100,109 do begin
	plot, x, data[200,i,*], yrange=[0,3e4], psym=7, /noerase
	plot, x, a[200, i] + b[200, i] * x, yrange=[0,3e4], /noerase
endfor

stop, 'Counts and lines plotted.  .c to continue.'

flat = a + b
flat = flat / median(flat)
print, 'max(flat), min(flat), mean(flat), median(flat), stddev(flat)'
print, max(flat), min(flat), mean(flat), median(flat), stddev(flat)
window, 5, xs=512, ys=512, title='Flat'
tvscl, flat
writefits, 'images/flat.fits', flat

stop, 'Flat displayed.  .c to continue.'

d = float(readfits('images/data1136.a.fits')) > 0
good = where(mask)
d2 = make_array(/float, 512, 512, value=0)
d2[good] = d[good] / flat[good]

window, 6, xs=512, ys=512, title='Raw data'
tvscl, d
window, 7, xs=512, ys=512, title='Normalized data'
tvscl, d2
writefits, 'images/data2.fits', d2
print, 'max(d), min(d), mean(d), median(d), stddev(d)'
print, max(d), min(d), mean(d), median(d), stddev(d)
print, 'max(d2), min(d2), mean(d2), median(d2), stddev(d2)'
print, max(d2), min(d2), mean(d2), median(d2), stddev(d2)

stop, 'Raw and normalized data displayed.  .c to continue.'

wset, 6
tv, 255 * 10 * d / max(d)
wset, 7                  
tv, 255 * 10 * d2 / max(d2)

stop, 'Stretched data displayed.  .c to continue.'

rm = rowmed(d2, mask)
window, 8, title='Rowmeds'
plot, rm

stop, 'Row medians displayed.  .c to continue.'

mask = mask * float(readfits('images/mask1.fits'))
d2 = d2 * mask
rm = rowmed(d2, mask)
plot, rm
wset, 7
tvscl, d2

stop, 'Masked data and row medians displayed.  --All done--.  .c to return.'

!except = 1
!quiet = 0

end