pro ipsim8,ratio,min,max,bad_pix_map
; 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
d_read,'landsat',h,p
answer=""
ymin=0
ymax=255
answer="hri"
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)
image_scale_tmp,min,max,p1,pscaled
p2=pscaled
p1=rotate(p2,2)
minp1=min(pscaled) & maxp1=max(pscaled)
new_sqrt,p1,f,b
; Evaluate sqrt damage
p1s=f(p1)
p1new=b(p1s)
s=snr(p1,p1new)
print,'From square rooting alone snr = ',s(0),', mean/rmse = ',s(1), $
  ',and variance = ',s(2)
sqrt_dmg=s(1)
p1s=byte(f(p1))
openw,1,'compress.inp'
writeu,1,p1s
close,1
cmd='/users/ldoose/TUB/darks_test/disrcomp compress.inp compress.tmp 256 '+ $
  string(xsize)+' '+string(ratio)
spawn,cmd
cmd='/users/ldoose/TUB/software/dcs compress.tmp compress.out '+ $
  string(ysize)+string(xsize)
spawn,cmd
openr,1,'compress.out'
pc=bytarr(xsize,ysize)
readu,1,pc
close,1
spawn,'rm compress.inp compress.tmp compress.out tmp.dec'
; Evaluate compression damage
s=snr(p1s,pc)
print,'From compression alone snr = ',s(0),', mean/rmse = ',s(1), $
  ', and variance = ', s(2)
pcb=intarr(xsize,ysize)
pcb=b(pc)
dark=fltarr(256)
dark=p(523-291,ymin:ymax)
pcbd=pcb
for i=0,xsize-1 do pcbd(i,*)=pcb(i,*)-dark
s=snr(p1,pcbd)
print,'Without blemishes snr = ',s(0),', mean/rmse = ',s(1), $
  ', and variance = ', s(2)
stop
return
end