pro ipsim
; variables:
; p = the original ccd image
; p1 = the subarray of p to be processed
; xmin,ymin,xmax,ymax = the coordinates of the subarray corners
; f = the forward square root lookup table
; b = the inverse square root lookup table
; p1s = the square-rooted image
; pc = the decompressed, square-rooted image
; pcb = the unsquare-rooted, decompressed image
; dark = the dark current vector from columns 6 - 9
; pcbd = the dark-corrected, unsquare-rooted, decompressed image
; pf = the raw flat-field image pixels
; pf1 = the subarray of the flat-field to be processed
; pf1c = the decompressed flat-field image
; pf1w = the word-sized decompressed flat-field image
; darkf = the dark vector of the flat field
; pfd = the dark-corrected, decompressed flat-field image
; pcbdf = the final flat-fielded, dark-corrected, unsquare-rooted,
;         decompressed image
; p1dc = the original subarray with dark correction
; pf1dc = the original flat field with dark correction
; p1dcf = the original subarray with dark correction and flat field correction
; bestflat = pcbd/p1dcf, the ideal flat field
print,'Select the image input file.'
d_wread,h,p
tvfit,p,0,'Original Frame'
print,'Enter colmin,rowmin,colmax,rowmax for the area to be processed.'
read,xmin,ymin,xmax,ymax
xsize=xmax-xmin+1
ysize=ymax-ymin+1
p1=intarr(xsize,ysize)
p1=p(xmin:xmax,ymin:ymax)
tvfit,rebin(p1,2*xsize,2*ysize,/sample),1,'Subframe'
p_histo,p1,x,h,2,'Subframe Histogram'
read_sqrt,f,b
p1s=byte(f(p1))
tvfit,rebin(p1s,2*xsize,2*ysize,/sample),3,'Square-rooted'
p_histo,p1s,x,h,4,'Square-rooted Image Histogram'
openw,1,'compress.inp'
writeu,1,p1s
close,1
print,'Please compress file compress.inp.'
print,'Save the decompressed image as compress.out.'
print,'Type a return to proceed'
s=' '
read,s
wdelete,0
openr,1,'compress.out'
pc=bytarr(xsize,ysize)
readu,1,pc
close,1
tvfit,rebin(pc,2*xsize,2*ysize,/sample),0,'Decompressed'
pcb=intarr(xsize,ysize)
pcb=b(pc)
tvfit,rebin(pcb,2*xsize,2*ysize,/sample),5,'Unsquare-rooted'
dark=fltarr(256)
dark=(p(6,*)+p(7,*)+p(8,*)+p(9,*))/4.
dark=fix(dark(ymin:ymax)+0.5)
pcbd=pcb
for i=0,xsize-1 do pcbd(i,*)=pcb(i,*)-dark
tvfit,rebin(pcbd,2*xsize,2*ysize,/sample),6,'Dark Correction'
print,'Select a flat-field image'
d_wread,hf,pf
pf1=intarr(xsize,ysize)
pf1=pf(xmin:xmax,ymin:ymax)
tvfit,rebin(pf1,2*xsize,2*ysize,/sample),7,'Flat Field'
p_histo,pf1,x,h,8,'Flat Field Histogram'
openw,1,'flat.inp'
writeu,1,bytscl(pf1)
close,1
pf1min=min(pf1)
pf1max=max(pf1)
print,'Please compress file flat.inp.'
print,'Save the decompressed image as flat.out.'
print,'Type a return to proceed'
read,s
openr,1,'flat.out'
pf1c=bytarr(xsize,ysize)
readu,1,pf1c
close,1
pf1w=fltarr(xsize,ysize)
pf1w=float(pf1max-pf1min)*pf1c/256.+pf1min
pf1w=fix(pf1w)
tvfit,rebin(pf1w,2*xsize,2*ysize,/sample),9,'Decompressed Flat Field'
p_histo,fix(pf1w),x,h,10,'Decompressed Flat Field Histogram'
darkf=fltarr(256)
darkf=(pf(6,*)+pf(7,*)+pf(8,*)+pf(9,*))/4.
darkf=fix(darkf(ymin:ymax)+0.5)
pfd=pf1w
for i=0,xsize-1 do pfd(i,*)=pf1w(i,*)-darkf
tvfit,rebin(pfd,2*xsize,2*ysize,/sample),10,'Dark Corrected Flat Field'
p_histo,fix(pfd),x,h,11,'Dark-corrected Flat Field Histogram'
flat=total(pfd)/(float(xsize)*float(ysize))
pfd=float(pfd)/flat
pcbdf=pcbd
pcbdf=fix(float(pcbd)/pfd+0.5)
tvfit,rebin(pcbdf,2*xsize,2*ysize,/sample),12,'Final Image'
p_histo,pcbdf,x,h,13,'Final Image Histogram'
p1dc=p1
for i=0,xsize-1 do p1dc(i,*)=p1(i,*)-dark
pf1dc=pf1
for i=0,xsize-1 do pf1dc(i,*)=pf1(i,*)-darkf
oflat=pf1dc/(total(pf1dc)/(xsize*ysize))
p1dcf=pf1
p1dcf=fix(p1dc/oflat+0.5)
tvfit,rebin(p1dcf,2*xsize,2*ysize,/sample),13,'Original Image Flat-fielded'
bestflat=fltarr(xsize,ysize)
bestflat=float(pcbd)*oflat/float(p1dc>1)
window,14,xsize=2*xsize,ysize=2*ysize,title='Best Flat Field'
tv,rebin(bytscl(bestflat,min=0.2,max=1.5),2*xsize,2*ysize,/sample)
stop
return
end