pro ipsim2
; 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,rotate(p,2),0,'Original Frame 0',0,800
answer=""
print,'Enter the area to be processed (hri, mri, sli).'
ymin=0
ymax=255
read,answer
case answer of
  'hri': begin
    xmin=364
    xmax=523
    end
  'mri': begin
    xmin=54
    xmax=229
    end
  'sli': begin
    xmin=234
    xmax=361
    end
else: print,'Invalid area'
endcase
xsize=xmax-xmin+1
ysize=ymax-ymin+1
p1=intarr(xsize,ysize)
p1=p(xmin:xmax,ymin:ymax)
p1=rotate(p1,2)
wdelete,0
tvfit,rebin(p1,2*xsize,2*ysize,/sample),1,answer+' 1',0,1100-2*270
p_histo,p1,x,h,2,answer+' histogram',0,70
new_sqrt,p1,f,b
p1s=byte(f(p1))
tvfit,rebin(p1s,2*xsize,2*ysize,/sample),3,'Square-rooted 3',2*xsize+10,1100-2*270
p_histo,p1s,x,h,4,'Square-rooted Image Histogram',430,70
openw,1,'compress.inp'
writeu,1,p1s
close,1
cmd='/users/ldoose/TUB/darks_test/disrcomp compress.inp compress.tmp 256 '+ $
  string(xsize)+' 4'
spawn,cmd
cmd='/users/ldoose/TUB/software/dcs compress.tmp compress.out '+ $
  string(ysize)+string(xsize)
spawn,cmd
wdelete,0
openr,1,'compress.out'
pc=bytarr(xsize,ysize)
readu,1,pc
close,1
spawn,'rm compress.inp compress.tmp compress.out tmp.dec'
tvfit,rebin(pc,2*xsize,2*ysize,/sample),0,'Decompressed 0',4*xsize+20,1100-2*270
pcb=intarr(xsize,ysize)
pcb=b(pc)
tvfit,rebin(pcb,2*xsize,2*ysize,/sample),5,'Unsquare-rooted 5',6*xsize+30,1100-2*270
dark=fltarr(256)
dark=p(523-291,ymin:ymax)
pcbd=pcb
for i=0,xsize-1 do pcbd(i,*)=pcb(i,*)-dark
tvfit,rebin(pcbd,2*xsize,2*ysize,/sample),6,'Dark Correction 6',20,1044-2*270-20
print,'Select a flat-field image'
d_wread,hf,pf
pf1=intarr(xsize,ysize)
pf1=pf(xmin:xmax,ymin:ymax)
pf1=rotate(pf1,2)
tvfit,rebin(pf1,2*xsize,2*ysize,/sample),7,'Flat Field 7',2*xsize+40,1044-2*270-20
p_histo,pf1,x,h,8,'Flat Field Histogram',860,70
new_sqrt,pf1,f,b
pf1min=min(pf1)
pf1max=max(pf1)
pfs=f(pf1)
pfs=byte(pfs)
openw,1,'compress.inp'
writeu,1,pfs
close,1
cmd='/users/ldoose/TUB/darks_test/disrcomp compress.inp compress.tmp 256 '+ $
  string(xsize)+' 4'
spawn,cmd
cmd='/users/ldoose/TUB/software/dcs compress.tmp compress.out '+ $
  string(ysize)+string(xsize)
spawn,cmd
openr,1,'compress.out'
pfsc=bytarr(xsize,ysize)
readu,1,pfsc
close,1
spawn,'rm compress.inp compress.tmp compress.out tmp.dec'
pf1w=fix(b(pfsc))
tvfit,rebin(pf1w,2*xsize,2*ysize,/sample),9,'Decompressed Flat Field 9', $
  4*xsize+60,1044-2*270-20
p_histo,fix(pf1w),x,h,10,'Decompressed Flat Field Histogram',20,0
darkf=pf(10,ymin:ymax)
pfd=pf1w
for i=0,xsize-1 do pfd(i,*)=pf1w(i,*)-darkf
tvfit,rebin(pfd,2*xsize,2*ysize,/sample),11,'Dark Corrected Flat Field 11', $
  6*xsize+80,1044-2*270-20
p_histo,fix(pfd),x,h,12,'Dark-corrected Flat Field Histogram',450,0
normimg=pfd(10:xsize-10,10:ysize-10)
std=stdev(normimg,flat)
pfd=float(pfd)/flat
pcbdf=pcbd
pcbdf=0>fix(float(pcbd)/(pfd>0.1)+0.5)<4095
tvfit,rebin(pcbdf,2*xsize,2*ysize,/sample),13,'Final Image 13',40,1044-2*270-50
p_histo,pcbdf,x,h,14,'Final Image Histogram',880,0
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
normimg=pf1dc(10:xsize-10,10:ysize-10)
std=stdev(normimg,oflat)
pf1dc=float(pf1dc)/oflat
p1dcf=0>fix(float(p1dc)/(pf1dc>0.1)+0.5)<4095
tvfit,rebin(p1dcf,2*xsize,2*ysize,/sample),15,'Original Image Flat-fielded 15', $
  2*xsize+60,1044-2*270-50
bestflat=fltarr(xsize,ysize)
bestflat=0>float(pcbd)*oflat/float(p1dcf>1)<4095
tvfit,rebin(bestflat,2*xsize,2*ysize,/sample),16,'Best Flat Field 16', $
  4*xsize+80,1044-2*270-50
stop
return
end