pro welldepth_reduce,sec=sec,offset=offset,tv=tv


; get the device type and part number from the pathname
ff=file_expand_path('ramp1.fits')
dev=strmid(ff,23)
dev=strmid(dev,0,strpos(dev,'/'))
dettemp=strcompress(string(fix(float(getkeyval('ramp1.fits','DETTEMP'))+0.5)),/rem)
wheel1=strcompress(getkeyval('ramp1.fits','WHEEL1'),/rem)
wheel2=strcompress(getkeyval('ramp1.fits','WHEEL2'),/rem)

fits_read,'ramp1.fits',r1
sz=size(r1)
fits_read,'ramp2.fits',r2
sz=size(r2)

; image section for stats
if not keyword_set(sec) then sec=[4,sz(1)-5,4,sz(2)-5]

; get the exposure times for each read of the ramp
n=sz(3)
exptime=(float(getkeyval('ramp2.fits','EXP_TIME'))/n)*findgen(n)

; find where images are off both rails. Be conservative
m1=dblarr(n)
for i=0,n-1 do m1(i)=median(r1(sec(0):sec(1),sec(2):sec(3),i))
u1=where(m1 gt -31000.0 and m1 lt 31000.0)
m2=dblarr(n)
for i=0,n-1 do m2(i)=median(r2(sec(0):sec(1),sec(2):sec(3),i))
u2=where(m2 gt -31000.0 and m2 lt 31000.0)
uboth=where(m2 gt -31000.0 and m2 lt 31000.0 and m1 gt -31000.0 and m1 lt 31000.0)

xx=findgen(n)

; get the mean offset between the two ramps in DN

; if user supplied offset, use this
;offset=62200.0
;offset=61800.0
if keyword_set(offset) then begin
	offset=offset
	; assume no ratedif
	ratedif=1.0
	uboth(0)=(-1)
endif

; if there is overlap between the ramps, use this
if uboth(0) ne -1 then begin
	ratedif=mean((m1(uboth(1:*))-m1(uboth(0))) / (m2(uboth(1:*))-m2(uboth(0))))
	delta=(((m1(uboth)-m1(0))/ratedif)+m1(0))  -m2(uboth)
	offset=mean(delta)
endif

; if no overlap between the ramps, use this
if uboth(0) eq -1 then begin
	; use a linear fit to the last 3 good points of the first ramp and the first 3 
	; points of the second ramp.
	p1=poly_fit(exptime(u1(n_elements(u1)-3:*)),m1(u1(n_elements(u1)-3:*)),1,/double)
	p2=poly_fit(exptime(u2(0:2)),m2(u2(0:2)),1,/double)
	; match to the midpoint between the ends of the two ramps
	midtime=mean([exptime(u1(n_elements(u1)-1)),exptime(u2(0))])
	offset=poly(midtime,p1)-poly(midtime,p2)
	; assume no ratedif?
	ratedif=1.0
endif

	
print,'Offset between ramps: '+strcompress(string(offset),/rem)+' ADU'

; apply the offset to ramp1
for i=0,n-1 do m1(i)=m1(i)-offset


ndif1=[9999.9,m1(1:*)-m1(0:*)]
udif1=where(ndif1(u1) lt 200.)
ndif2=[9999.9,m2(1:*)-m2(0:*)]
udif2=where(ndif2(u2) lt 200.)

count=0
topend=dblarr(sz(1),sz(2))
if udif2(0) ne -1 then begin
	topend=double(r2(*,*,u2(udif2(n_elements(udif2)-1))))
	ind=udif2(n_elements(udif2)-1)
	usat=where(abs(double(r2(*,*,u2(ind))) - double(r2(*,*,u2(ind-1)))) lt 200)
	count=1
endif


fname=dev+'_linsat.ps'
if keyword_set(tv) then window,0,xs=700,ys=900
if not keyword_set(tv) then begin
	set_plot,'ps
	device,filename=fname,/inches,/port,xsize=8.0,ysize=10.5,xoffset=0.25,yoffset=0.25,/color,bits=8
	!p.font=0       ; default font is Times
	;psinit,/full
	;device,filename=fname
endif

!p.multi=[0,1,2]

if count eq 0 then begin
	plot,findgen(100),xst=5,yst=5,/nodata
	med=60.
	xyouts,10,45,'Saturation not reached in these ramps',/data
endif

; bottom end is just the first good read of the first ramp, as long at that first read is in
; range
if u1(0) eq 0 then bottomend=double(r1(*,*,u1(0)))-offset
if u1(0) ne 0 then begin
	!p.multi=[0,1,2]
	plot,finden(100),xst=5,yst=5,/nodata
	xyouts,10,35,'0th read level at bottom rail in first ramp',/data
	med=60.
endif

if count gt 0 then begin
topend=topend/count
dif=topend-bottomend

satmask=dif-dif
satmask(usat)=1
dif=dif(sec(0):sec(1),sec(2):sec(3))
difmask=satmask(sec(0):sec(1),sec(2):sec(3))
uu=where(difmask eq 1)
dif=dif(uu)


binsize=100
npix=double(n_elements(uu))

;; use gauss fit to histogram
;;med=median(dif)
;a=gaussfithist(dif/100.,bin=1)
;med=a(1)*100.

; use max of histogram
y=((histogram(dif,bin=binsize,min=0)/npix)*100.)/binsize
;maxx=y(a(1))
maxx=max(y)
xx=findgen(n_elements(y))*binsize
u=where(y eq max(y))
med=xx(u(0))

cumulative=dblarr(n_elements(y))
cumulative(0)=y(0)
for i=1,n_elements(y)-1 do cumulative(i)=cumulative(i-1)+y(i)

plot,xx,y,xra=[med-10000,med+10000],psym=10,xtitle='Well Depth (ADU)',ytitle='% Number of Pixels',title=dev+' Well Depth Histogram at '+dettemp+' K ('+wheel1+' + '+wheel2+')',thick=5,xthick=6,ythick=6,charsize=1.2,yra=[0,maxx*1.1],yst=9,xmar=[8,6],xst=1

plots,[med,med],[0,maxx],lines=2,thick=5
xyouts,med+600.,maxx*0.1,strcompress(string(med),/rem),orient=90
axis,yaxis=1,yrange=[0,100],/save,ytitle='Cumulative %',charsize=1.2,ytickv=indgen(11)*10.,yticks=10,yminor=5,color=black,ythick=6
oplot,xx,cumulative*100.,lines=3,thick=5

endif



; get voltages from thea header
if strmid(dev,0,1) eq 'S' then ss=['VDETCOM','VDDUC','VDDCL','VROWOFF','VSSUC','ISLEWREF','VDDOUT','VSUB','VSSOUT1','VP','VNROW','VGGCL','VNCOL']
if strmid(dev,0,1) eq 'H' then ss=['DSUB','VRESET','VBIASPWR','VBIASGTE','DRAIN','VPULLA7','VPULLA15','OUTWINA','VREF','VPULLA23','VPULLA31','VDDA','VDD']


st=strarr(n_elements(ss))
for i=0,n_elements(ss)-1 do st(i)=strcompress(getkeyval('ramp2.fits',ss(i)),/rem)
for i=0,12 do xyouts,med+4000.,85.-(i*5),ss(i)+' = '+st(i)+' V',charsize=0.9


; fit for linearity

vec=[(m1(u1(1:*))-m1(u1(0)))/ratedif,m2(u2(1:*))-m1(u1(0))]
vec1=(vec/[exptime(u1(1:*)),exptime(u2(1:*))])

ufit=where(vec lt 0.8*max(vec))
s=sort(vec)
vec=vec(s)
vec1=vec1(s)

p=poly_fit(vec(ufit),vec1(ufit),2,/double)
p1=p/p(0)

ugood=where(vec gt 0.0)
plot,vec(ugood),vec1(ugood)/p(0),yra=[0.6,1.05],xra=[0,max(vec)*1.1],thick=5,xthick=6,ythick=6,xst=1,yst=1,title=dev+' Linearity',xtitle='Total DN in Well',ytitle='Fractional Countrate',charsize=1.2,xmar=[8,6]
circsym,1,1
oplot,vec(ugood),vec1(ugood)/p(0),psym=8,thick=5
circsym,0.6,1
oplot,vec(ugood),vec1(ugood)/p(0),psym=8,thick=5,color=255
circsym,1,1
oplot,vec(ufit),vec1(ufit)/p(0),psym=8,thick=5
oplot,findgen(max(vec)*1.1),poly(findgen(max(vec)*1.1),p1),lines=3,thick=4
plots,[0,max(vec)*1.1],[1.0,1.0],lines=2,thick=3

a1=strcompress(string(p1(1),format='$(e20.10)'),/rem)
a1=strmid(a1,0,5)+strmid(a1,strpos(a1,'e'))
a2=strcompress(string(p1(2),format='$(e20.10)'),/rem)
a2=strmid(a2,0,5)+strmid(a2,strpos(a2,'e'))

xyouts,max(vec)*0.1,0.8,'fit = 1 + ('+a1+' * DN) + ('+a2+' * DN^2)'
xyouts,max(vec)*0.2,0.77,'= included in fit',/data
xyouts,max(vec)*0.2,0.74,'= not included in fit',/data
circsym,1,1
oplot,[max(vec)*0.17],[0.775],psym=8,thick=5
circsym,1,1
oplot,[max(vec)*0.17],[0.745],psym=8,thick=5
circsym,0.6,1
oplot,[max(vec)*0.17],[0.745],psym=8,thick=5,color=255


; for some single pixels
;plot,r2(1000,1100,*)-float(r2(1000,1100,0)),((r2(1000,1100,*)-float(r2(1000,1100,0)))/(findgen(10)*30.))/280.,yra=[0.6,1.1]
;for i=1000,1100 do oplot,r2(i,1100,*)-float(r2(i,1100,0)),((r2(i,1100,*)-float(r2(i,1100,0)))/(findgen(10)*30.))/280.


!p.multi=0

if not keyword_set(tv) then begin
	device,/close
	set_plot,'x
	;psterm,/noplot,/nostamp
	a=execute("spawn,'ps2gif -antialias -density 150 '+fname")
endif

	!p.font=(-1)

end